Lennard-Jones model and Energy Conservation

[onizuka]

Hi,
i've been doing a simulation using LJ model, but I'm having a troublesome time on figuring out, what's happening with the energy :uhh:

This is the LJ potencial:
$$E = 4\epsilon\left[\left(\frac{\sigma}{R}\right)^{12} - \left(\frac{\sigma}{R}\right)^6\right]$$

but I've been using this simplified form given on Fosdick's book:
$$E = \left[\left(\frac{1}{R}\right)^{12} - \left(\frac{2}{R}\right)^6\right]$$

And for time integration I'm using the velocity Verlet method.

I am using noble gases only, and in this example, i have 3 atoms forming an equilaterum triangle. The system starts with KE = 0, and PE < 0.

What happens is that while the atoms approach each other, the KE increases (as one would expect), and the PE decreases (it is always a negative value, so it increases in absolute value).
When they are too close, they repel themselfs and the opposite happens.

So i have a kind-of "harmonic" motion, where the $$E$$ increases and decreases. :yuck:
To give you an idea of the values i have, i initially have an $$E$$ = -0.00051 and can go up to $$E$$ = -0.3619


My suspection of why this doesn't work migth be because of the values I'm using in the simulation. (and yes... i really don't know what is the magnitude of the SI units i am using here... much less the ones i should be using :shy: )
time step = 0.01
particle mass = 1
radius = 0.6
distance between each particle = 5


I've read that using the conservation of energy, one can do some corrections on the values obtained... but uhh, i think what i have here i something completely different.

Thanks in advance for any help.

 

[CharlesConcord]

From what you described, I cannot tell whether the energy in your simulation is conserved or not. The total energy (i.e. p.e. plus k.e.) should conserve
nearly perfectly if you use a smoothing function to avoid the error caused by the potential cutoff. A classical book for learning MD simulation is Allen and Tildesley's "Computer Simulation of Liquids".

I myself tend not to use the reduced unit notation for the LJ potential because it would then create headache when you add new types of potentials such as chemical bonds, electricstatic force fields and so on. I don't think one would have a pleasant time in converting the energy of these potentials into the reduced unit, which is based solely on the LJ potential.

The Molecular Workbench software I am developing might be a useful reference for you to learn MD simulations. Check it out at: http://mw.concord.org/modeler/index.html

 

[quetzalcoatl9]

when i started my PhD i wrote MD code to do exactly what you describe.

you should see energy conservation to at least 4 significant figures. check to make sure that you are calculating the energy correctly first of all. also check your periodic boundaries and your verlet integrator. visualize the output to make sure that the motions of the LJ particles looks correct.

long-range corrections should make only a very small difference to the energy conservation.

also remember that you should equilibrate your system, but this still should not affect the energy conservation.

once you have it working, it is interesting to apply some stat.mech. to calculate things like g(r), velocity autocorrelation, MSD, etc. and calculate the pressure, diffusion constant, etc.

as the other poster mentioned, it is benefitial to work in reduced units. refer to the seminal paper by Tildesley and Allen (I forget which journal it was published in) where they give the equation of state and a table of parameters they used.