[EXTERNAL] Re: Lammps MD simulation help

I would take a step back and check that both LAMMPS and your other code
are running the same force field. ReaxFF has quite a few variants. You
should be able to match forces and energies from an interesting
configuration of atoms to high precision. For best results, use the same
cell and position input data for both codes. ***Don't run MD***, as the
trajectories will generally be different for many reasons that have
nothing to do with the force field.

Please direct future correspondence on this to the LAMMPS mailing list.


Dear Aidan,

Thank you very much for the advise of making sure the ReaxFF stand-alone code and LAMMPS can reproduce the same energy.
After I checked, I did use the same force field for the two calculations and my test calculations for the same input structure do result in the same energy for the two methods:
From Lammps:
Step Temp E_pair E_mol TotEng Press
       0 0 -662210.87 0 -662210.87 1240.9704

From ReaxFF:
Iteration Nmol Time(fs) Epot(kcal) Vol(A^3) T(K) Pres(MPa) Dens(kg/dm3)
        0 6000 0.00 -662210.88 6623488.15 4000.00 0.00 0.02

So it seems I am using the correct force field, but as far as I carry on MD simulations, the two methods will give completely different pictures.
I have been trying to adjust the parameters in the reax/c control files but haven't got any luck yet. I would greatly appreciate it if you can give me some suggestions on running the MD correctly.



Hi Qing,

It is encouraging that you get the same energy, however the pressure
is different. From LAMMPS you got 1240 atm or 124 MPa, while with
ReaxFF the pressure was 0.0 MPa; this indicates some initial
setting(s) is/are still different. Otherwise first step
energy/pressure should be the same.

The default settings in the control file, given in this page:
http://lammps.sandia.gov/doc/pair_reax_c.html, are the same as the
stand-alone code, so you don't have to play with the settings in the
control file. Instead you can use "pair_style reax/c NULL".

Try examine your input script and try dump the forces on each atom and
examine if they are exceptional large.


That is good agreement, however you should still worry about how well this
snapshot is exercising the different parts of the potential. Look at the
examples/reax examples to see how to look at the itemized contributions to
the energy:

compute reax all pair reax/c

variable eb equal c_reax[1]
variable ea equal c_reax[2]

Also, for LAMMPS, the control file only affects the potential, not the MD
simulation. If you think the potential is correct, then you should examine
whether you are performing the LAMMPS MD simulation correctly. Check
timestep, energy conservation, etc.

I assumed the value of zero was due to the pressure not being calculated.


Dear Aidan,

I have extracted the itemized contributions and they are all very consistent with each other:

Step eb ea elp emol ev epen ecoa ehb et eco ew ep efi eqeq
0 -1145437.6 -29722.567 -2.3441655e-06 0 0 0 0 0 0 0 519229.29 2414.7532 0 -8694.7081

Iter. Ebond Eatom Elp Emol Eval Ecoa Ehbo Etors Econj Evdw Ecoul Echarge Efield
0 -1145437.64 -29722.57 0.00 0.00 0.00 0.00 0.00 0.00 0.00 519229.29 2414.76 -8694.72 0.00

So I would think the problem is not from the potential part.

I did, however, notice that the pressure term is being funky (also pointed by Ray)

Following is the pressure output from Lammps for a NVT simulation:

Step Temp Press
0 0 1240.9704
100 543.04078 -3271.3535
200 1881.6707 6332.5587
300 4064.0868 -807.75083
400 3347.3628 -625.46872
500 4046.5729 883.19246
600 4329.614 523.96399
700 3902.9782 124.86305
800 3480.9631 1130.7066
900 3980.2941 230.469
1000 4434.3128 680.8794
1100 3854.9525 1294.6761
1200 3889.6867 -64.187006
1300 3986.0691 523.91264
1400 3924.022 663.16824
1500 3987.5709 197.08486
1600 3941.0818 1373.6556
1700 3985.9305 282.61069
1800 4124.7859 608.76089

ADF/Reaxff always outputs pressure as 0 for NVT/NVE simulations, which I assumed as they are not calculated.
But it is quite puzzle here how Lammps has outputted such unstable pressure and even negative pressures sometimes.
My simulation stepsize is at 0.1 fs, and I used "fix 1 all nvt temp 4000.0 4000.0 1.0". Is there anything else I need to be aware of to set the MD correctly?



The pressure from LAMMPS looks perfectly normal. Nothing funky there. Your
thermostat damping constant is 1.0 fs. This is very small. A better choice
would be 10.0 fs or 100.0 fs. There may be other problems with your
simulation set-up too. I suggest that you spend some time familiarizing
yourself with MD and LAMMPS. Try running some of the examples in the
examples folder.