Atomic simulation and lots of questions

[Zach]

First, hi. This is my first post, and after careful perusing of the forum I *think* this post is in the correct place. Feel free to move it around if not, apologies.

The meat of my questions. I'm a High School student that would like to simulate atoms on the computer. I've been getting help from a lot of various sources, and therefore have a real hodgepodge of information. I'm not quite sure how it all fits together. And to top it off, I think there may be some math needed that I don't know yet (calculus, which I have next trimester). But I'm willing to learn :)

I have been corresponding with a few chemistry professors. They have told me to start off small using the Lennard Jones potential. I understand the theory behind the potential quite well now (at medium distance, the two atoms in question are mildly attractive, and at short distances they are extremely repulsive. As the distance increases, the attraction rapidly drops off to about nothing). The LJ potential is as follows:

$$E = 4\epsilon\left[\left(\frac{\sigma}{R}\right)^{12} - \left(\frac{\sigma}{R}\right)^6\right]$$

Where E is the energy in joules/mol, $$\epsilon$$ is the depth of the well in joules/mol, $$\sigma$$ is the distance in angstroms for the interaction to occur. R is distance in angstroms between the two atoms. Both $$\epsilon$$ and $$\sigma$$ are constants

With that in mind, I was told this:

The basic way to use the L-J potential, any potential really, is to take the derivative of it with respect to the distance between the atoms to get the force between the atoms. Since you know the mass of the atoms, you can then determine the acceleration acting on them. From the acceleration and current velocity you can get their new velocity and ultimately the new position of the atoms some amount of time later.

And from here my questions start. Do I need to differentiate that equation? I don't know calculus, so I'm not sure if that is needed. Or does that merely mean plug in the distance to get the energy.

From there, I know the mass, as the quote says. Finding the acceleration is trivial (F=ma). The next part of my questions has to deal with the "ultimately finding the new positions" part. I have been told by a few people I will need to "integrate with respect to time". I've also been told by someone else that I don't need classical calculus, but instead need to use Euler's Integration, as that is the only way to do such a thing on the computer (I believe I understand this. Calculus is continuous, while the Euler integration assumes a constant value in between timesteps. If I were to use classical calc, I would have to reduce my timestep infinitely to achive the results). But again, I have no idea, and could be talking out of my rear.

Thanks, any help would be greatly appreciated. I am swimming in a sea of knowledge, but no real way to piece it together.

BTW, this LaTex math function stuff is extremely cool :)


EDIT: I forgot to mention, I've also heard of using the Verlet Algorithm for integratin instead of the others. Again, no idea if this is what I should aim for or not. Thanks again.

 

[pervect]

Well, I think you'll have to learn to differentiate at least a*x^n to proceed much further. Fortunately, it's pretty easy.

$$\frac{d}{dt} a x^n = a n x^{n-1}$$

So to get the forces out of your potential, you'll get

$$F = \frac{dE}{dt} = \frac{4\epsilon}{R}\left[12\left(\frac{\sigma}{R}\right)^{11} - 6\left(\frac{\sigma}{R}\right)^5\right]$$

hmmm, I realized you'd need the chain rule to understand where the 1/R in the above came from. Well, you need the 1/R, just take my word for it :-)

Anyway, the Verlett algorithm appears to be commonly used in molecular dynamics, a reasonably clear implmenetation is at

http://www.ph.ed.ac.uk/~graeme/compmeth/verlet.html

I'm *guessing* that the reason the Verlett agorithim is popular is because it's symplectic, but i don't really know for sure. At the pre-calculus level, I don't think you should worry too much about the algorithm selection. If you can get Verlett to work, great.

 

[Zach]

Great! Thanks for the help. Couple of simple questions. Why does the equation need to be differentiated (just out of curiosity)? What is it in the undifferentiated form if it is not usable to get the force?

I'll look into the Verlett algorithm more. I believe I have a book around here that also describes it as well.

 

[pervect]

Work (energy) is force * distance.

So to get the force, you need to take the derivative of the energy with respect to distance.

Consider the force of gravity, for instance

the potential energy is E = -GmM/r, where r is the radial distance.

The force is dE/dt = GmM/r^2, directed in a radial direction.

You have a similar situation - your prof gave you the energy, you need to find the force.

You know that like the case of gravity, the force is radial, you also know that it gets repulsive at close range and attractive at long range, so you should be able to get the direction and sign worked out correctly.

 

[kdkdkd]

first of all I am amazed a high school student has taken up a project on simulation of atoms. that's shows you are an inspired guy.
but simulation of atoms is much beyond just simulation of the gravitational-electrostatic force equilibrium right?
you'll have to take the characteristics of atoms into consideration, the temperature at which you want to get them closer( because the atoms might actually form a bonding to form a molecule if they have sufficient energy). bonding distances come into picture. no discouragement intended!
accomodating all these will strengthen your model. you can also check out the Runge-kutta numerical algorithm for integration, its used in many applications.

 

[Zach]

Pervect: Ahh, I understand now. Thanks for working that out for me. I didn't realize (should have, but guess I wasnt thinking) that energy and force were two separate things.

kdkdkd: Yeah, that seems to be the general response (although usually more negative, such as "Go and work on something easier, you don't stand a chance"). I plan on modeling, at first, only noble gasses. The LJ potential works with all atoms, like you said, but works exceptionally well with noble gasses because they are inert. I won't have to worry about any reactions taking place.

But I do plan, later, to add in additional complexity. But first, the noble gasses. Baby steps. :)

Thanks again everyone for the help!

 

[Zach]

First off, sorry to revive this ancient thread. Basically, I would like to look more into this subject, and now that I have had a good (well, better) handle on basica calculus, I think I might have a shot.

So, a few clarifications if I want to model this on the computer. I need the first derivative to the LJ potential $$\frac{dE}{dt}$$. Given starting position and mass of all the atoms, I use the $$f=ma$$ (force coming from the derived LJ potential, mass given) to find acceleration. I then use acceleration (or perhaps determine velocity from the acceleration) over a predetermined discrete timestep to find displacement. This is done for all the atoms in the system. Time is incremented, all atoms moved. Process repeats.

Is this process correct?