Questions about LAMMPS and Materials Studio Calculations

Case · October 29, 2015, 2:24am

Hi LAMMPS users and developers:

I am learning to use LAMMPS in my research now and encountered some problems which I searched for the solution but found none online.
The main problem is the discrepancy of energy between LAMMPS and Materials Studio/Forcite (referred MS in the following text). The system studied is a PVDF chain with 10 units)

LAMMPS structural optimization (log file is shown below and the calculation is using data file converted from MS with msi2lmp code):
LAMMPS (10 Aug 2015)
#----------------------------- Initialization ------------------------------
units real
boundary p p p
atom_style full
#--------------------------- pcff potential information ------------------------------
pair_style lj/class2/coul/long 12.0 10.0
kspace_style ewald 1.0e-6
bond_style class2
angle_style class2
improper_style class2
dihedral_style class2
#neighbor 0.4 bin
#neigh_modify every 10 one 10000
read_data pvdf.data
orthogonal box = (-2.17853 -2.15272 -42.2431) to (37.8215 37.8473 -2.24307)
1 by 1 by 1 MPI processor grid
reading atoms …
62 atoms
scanning bonds …
4 = max bonds/atom
scanning angles …
6 = max angles/atom
scanning dihedrals …
18 = max dihedrals/atom
scanning impropers …
4 = max impropers/atom
reading bonds …
61 bonds
reading angles …
120 angles
reading dihedrals …
171 dihedrals
reading impropers …
80 impropers
4 = max # of 1-2 neighbors
6 = max # of 1-3 neighbors
12 = max # of 1-4 neighbors
16 = max # of special neighbors
#velocity all create 300.0 1231
#fix 1 all nve/limit 0.05
#fix 2 all langevin 300.0 300.0 10.0 904297
#thermo_style custom step temp pe
#thermo 10000
#timestep 1
#run 1000000
#unfix 1
#unfix 2
#write_restart restart.${simname}.dreiding1
#minimize 1.0e-25 1.0e-25 50000 100000
#thermo 1
thermo_style custom step press etotal pe evdwl ecoul elong epair ebond eangle edihed
dump 1 all cfg 6 dump.comp_*.cfg mass type xs ys zs
#dump 1 all xyz 1 file.xyz
#min_style cg
minimize 1.0e-25 1.0e-25 50000 100000
WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:168)
Ewald initialization …
G vector (1/distance) = 0.266803
estimated absolute RMS force accuracy = 0.000338812
estimated relative force accuracy = 1.02032e-06
KSpace vectors: actual max1d max3d =2084 10 4630
kxmax kymax kzmax = 10 10 10
Neighbor list info …
1 neighbor list requests
update every 1 steps, delay 0 steps, check yes
master list distance cutoff = 14
ghost atom cutoff = 14
Memory usage per processor = 17.6165 Mbytes
Step Press TotEng PotEng E_vdwl E_coul E_long E_pair E_bond E_angle E_dihed
0 162.18374 325.43216 325.43216 -1.9378116 525.56826 -178.78222 344.84823 0.24822394 1.7470875 -21.506128
1454 1.8019873 309.60019 309.60019 -4.2304778 501.14871 -179.52404 317.39419 3.4584588 15.982549 -27.897851
Loop time of 4.21872 on 1 procs for 1454 steps with 62 atoms
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
325.43215768 309.600187407 309.600187407
Force two-norm initial, final = 24.486 0.155976
Force max component initial, final = 6.42681 0.0519857
Final line search alpha, max atom move = 4.07521e-07 2.11853e-08
Iterations, force evaluations = 1454 2945
Pair time () = 0.188569 (4.46981) Bond time () = 0.567602 (13.4544)
Kspce time () = 3.41002 (80.8306) Neigh time () = 0.00050211 (0.0119019)
Comm time () = 0.00234032 (0.0554746) Outpt time () = 0.0364492 (0.863987)
Other time (%) = 0.0132387 (0.313808)
Nlocal: 62 ave 62 max 62 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 35 ave 35 max 35 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1531 ave 1531 max 1531 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1531
Ave neighs/atom = 24.6935
Ave special neighs/atom = 11.3548
Neighbor list builds = 7
Dangerous builds = 0
MS structural optimization:
---- Initial structure ----
Total enthalpy : -294.754158 kcal/mol
External pressure term : 0.000000 kcal/mol
Total energy : -294.754158 kcal/mol
Contributions to total energy (kcal/mol):
Valence energy (diag. terms) : 70.990
Bond : 87.073
Angle : 5.026
Torsion : -21.109
Inversion : 0.000
Valence energy (cross terms) : 2.622
Stretch-Stretch : -0.010
Stretch-Bend-Stretch : -0.400
Stretch-Torsion-Stretch : 0.183
Separated-Stretch-Stretch : 0.000
Torsion-Stretch : 0.036
Bend-Bend : -0.084
Torsion-Bend-Bend : 1.176
Bend-Torsion-Bend : 1.721
Non-bond energy : -368.366
van der Waals : 21.532
Long range correction : -0.008
Electrostatic : -389.890
rms force : 4.041E+001 kcal/mol/A
max force : 1.646E+002 kcal/mol/A
Cell parameters: a: 40.000000 A b: 40.000000 A c: 40.000000 A
alpha: 90.000 deg beta: 90.000 deg gamma: 90.000 deg
---- Final structure ----
Total enthalpy : -386.499890 kcal/mol
External pressure term : 0.000000 kcal/mol
Total energy : -386.499890 kcal/mol
Contributions to total energy (kcal/mol):
Valence energy (diag. terms) : -18.996
Bond : 0.251
Angle : 1.871
Torsion : -21.118
Inversion : 0.000
Valence energy (cross terms) : -0.498
Stretch-Stretch : 0.001
Stretch-Bend-Stretch : -0.124
Stretch-Torsion-Stretch : 0.022
Separated-Stretch-Stretch : 0.000
Torsion-Stretch : -0.119
Bend-Bend : 0.014
Torsion-Bend-Bend : 0.015
Bend-Torsion-Bend : -0.308
Non-bond energy : -367.005
van der Waals : 11.404
Long range correction : -0.008
Electrostatic : -378.402
rms force : 8.172E-002 kcal/mol/A
max force : 3.807E-001 kcal/mol/A
Cell parameters: a: 40.000000 A b: 40.000000 A c: 40.000000 A
alpha: 90.000 deg beta: 90.000 deg gamma: 90.000 deg
SOME SETTINGS for MS calculation: PCFF(force field), Force Field Assigned charge, quality is medium, van der Walls is atom-based, electrostatic is using Ewald summation.
Things I’d like to point out:
The energies have different sign for MS and LAMMPS claculations, while MS has negative energies(325.43216 and 309.60019 for initial and final structures) , LAMPPS has postivie energies (-294.754158 and -386.499890 for initial and final structures).
The energy term Ebond is different from each other for LAMMPS and MS initial structural calculations. MS has a Ebond of 87.073 while LAMMPS gives 0.24822394.
The final optimized structure of MS and LAMMPS calculations are quite different shown as follow: The curved one is from LAMMPS and linear one is from MS.
With the difference between LAMMPS and MS results presented above, I am wondering if the difference is made by mistake or it truly exists. Suppose the calculation process is correct (can’t guarantee this since I am a new user of LAMMPS and new to CLASS2 force field), what is the origin of the difference? Shouldn’t MS and
LAMMPS produce the same results? The positive energy from LAMMPS seems not reasonable to me? Maybe there is something wrong with my data file. By the way, when I modify the charge in data file to zero (no charge for all atoms), the result gives a negative total energy and the optimized structure is more linear-like, close to that optimized from MS.
Anyway, above-mentioned difference has been confusing me for quite a while. Hope you guys can give comments and suggestions and thanks a lot.Attached is the data file for LAMMPS.

pvdfchain.data (30.9 KB)

Ray_Shan · October 29, 2015, 4:09am

If you did not get the same energy/force at the zeroth step, then there could be something wrong with the converted data file. From my understanding the msi2lmp tool has been unmaintained for a while. You need to verify that you have a correct data file for LAMMPS.

Ray

akohlmey · October 29, 2015, 1:28pm

Hi LAMMPS users and developers:

I am learning to use LAMMPS in my research now and encountered some
problems which I searched for the solution but found none online.
The main problem is the discrepancy of energy between LAMMPS and Materials
Studio/Forcite (referred MS in the following text).

one major problem in this context is, that you will need to consult with
an expert in both, materials studio and LAMMPS. however, it is my
impression that there are rather few people around here that have that. i,
for example, have never used MS and have no access to it.

the second issue is, that absolute energies are as such arbitrary with
classical potentials. different MD codes have different conventions and to
some degree, they can be shifted (check out "pair_modify shift yes" ).
using long-range electrostatics makes things even more complicated, as also
the division between real space and kspace is dependent on several
settings, that are usually determined through some heuristics in the
respective codes and those may not be the same.

thus to reliable identify differences, you need to map out for each code
each individual energy/force term, e.g. set up two atoms with only a bond
interaction and compute the energy and forces for several distances.
similarly, it is far easier to compare coulomb with cutoff and compute
forces/energies for VdW and Coulomb separately by either setting epsilon or
the charge to zero for a pair of atoms. in each case you *have* to have
multiple points, so you can determine the functional form used with
sufficient accuracy and compare the shapes of the energy/force curves.

[...]

1. The energies have different sign for MS and LAMMPS claculations, while

MS has negative energies(325.43216 and 309.60019 for initial and final
structures) , LAMPPS has postivie energies (-294.754158 and -386.499890 for
initial and final structures).

as stated above, you should not pay much attention to total energies.

2. The energy term Ebond is different from each other for LAMMPS and MS
initial structural calculations. MS has a Ebond of 87.073 while LAMMPS
gives 0.24822394.

please make sure you are looking numbers with the right units and check
what kind of normalization is used on output. it is more important to
compare forces. thus it is best to compare just a single interaction (build
a minimal structure and then turn contributions selectively off. keep in
mind that non-bonded terms are usually not computed between pairs of atoms
that have a bond between then).

3. The final optimized structure of MS and LAMMPS calculations are quite
different shown as follow: The curved one is from LAMMPS and linear one is
from MS.

given the asymmetric structure of your compound, the result from the
LAMMPS minimization looks more likely to me. but for a compound like yours,
both of these structures are obviously local minima and thus have little
relevance as to whether your force field parameters are translated/used
correctly in LAMMPS.

With the difference between LAMMPS and MS results presented above, I am
wondering if the difference is made by mistake or it truly exists. Suppose
the calculation process is correct (can’t guarantee this since I am a new
user of LAMMPS and new to CLASS2 force field), what is the origin of the
difference? Shouldn’t MS and
LAMMPS produce the same results? The positive energy from LAMMPS seems not
reasonable to me? Maybe there is something wrong with my data file. By the
way, when I modify

so far, nothing you have reported is convincing evidence toward either
correctness or incorrectness of the LAMMPS simulation. as mentioned before,
short of having around somebody with extensive experience in MS, you have
to drill it down to the smallest possible pieces and check out each piece
separately.

there are some known issues, though. the msi2lmp tool has only seen minor
updates for the last >10 years, mostly to improve portability and
compatibility with current versions of LAMMPS. it was written for some
specific purposes and does not support *all* features of the .frc file
format, as became evident from a recent posting to this mailing list. in
particular, the .frc file format, as read by MS, allows to have multiple
functional forms for the same interaction, e.g. use a morse or a harmonic
potential for bond interactions, or different ways to represent LJ
parameters. the msi2lmp tool as bundled with LAMMPS does not check for any
of that and will have some fixed assumptions what the functional form is
based on the class of force field. thus, if you are using settings in MS
that result in a selection of different parameters and/or functional forms
than what msi2lmp does, you will get differences. from the information you
provide and the test you made, it is inconclusive. setting up a simulation
results into something that looks mostly reasonable, but the differences
you note seem a bit drastic. but then again, in a class2 force field the
balance between LJ and Coulomb is different and thus it is not automatic to
have a total negative energy.

the charge in data file to zero (no charge for all atoms), the result
gives a negative total energy and the optimized structure is more
linear-like, close to that optimized from MS.
Anyway, above-mentioned difference has been confusing me for quite a
while. Hope you guys can give comments and suggestions and thanks a
lot.Attached is the data file for LAMMPS.

yes, it *is* confusing. however, the general attitude towards such
observations is, that the person best suited to work it out, is the person
that has a vested interest to make this work (well), i.e. you. if you can
provide more conclusive proof on which particular details are not working
correctly and make proper test cases available, people way be able to look
into it and - at least - add code to detect inconsistencies when using
msi2lmp. there currently is not much desire among the LAMMPS developers to
enable msi2lmp to be fully compatible with all options that MS offers. like
most codes in the "tools" tree, they come with no claim to completeness or
guarantee to be fully tested or compatible with the current version of
LAMMPS.

axel.

Andrew_Jewett · October 29, 2015, 11:34pm

as stated above, you should not pay much attention to total energies.

2. The energy term Ebond is different from each other for LAMMPS and MS
initial structural calculations. MS has a Ebond of 87.073 while LAMMPS gives
0.24822394.

it is best to compare just a single interaction (build a
minimal structure and then turn contributions selectively off. keep in mind
that non-bonded terms are usually not computed between pairs of atoms that
have a bond between then).

Agreed. Instead of a PVDF chain with 62 atoms, start with something
even simpler: a single molecule of either O2, CO2, methane, ethylene,
ethanol, or ethane.

the general attitude towards such
observations is, that the person best suited to work it out, is the person
that has a vested interest to make this work (well), i.e. you. if you can
provide more conclusive proof on which particular details are not working
correctly and make proper test cases available, people way be able to look
into it and - at least - add code to detect inconsistencies when using
msi2lmp.

If you want to help fix this, then dig deeper into the source of
the discrepancy. Ray mentioned electrostatic effects, and Axel
mentioned different pair cutoff schemes.

In addition to those issues, I would check the 3-body and 4-body
bonded interactions: count number of (3-body) angles, (4-body)
dihedral interactions, and (4-body) improper interactions for your
molecule. Do this in both LAMMPS and in MS. And do it for a simple
molecule (like ethylene). (In LAMMPS, this is easy. Just open the
DATA file in a text-editor and look at the first few lines of the
file.) Make sure that, for your small molecule, MS and LAMMPS are
calculating the same number of these interactions. (As a reference,
ethylene has 4 diheral interactions and 6 improper interactions, when
using the OPLS force-field. This might be different for you since you
are using a class2 force-field.)

If the number of interactions is the same, then compare their
energies. Calculate the (2-body) bond-energy, (3-body) angle-energy,
(4body) dihedral-energy, and (4body) improper-ergy for your molecule,
using the same starting coordinates.

To do that, add "eimp" to your thermo_style:

thermo 1
thermo_style custom step pe evdwl ecoul elong epair ebond eangle edihed eimp
run 1 # run a single timestep and print out the energy

...As Axel mentioned, absolute energy numbers are physically
meaningless. (Different software may add different offsets to each
term in the energy.) However the energy differences (ie. between two
different conformations of the same molecule), are not meaningless,
and should be the same in both MS and LAMMPS. So if the energies for
a given structure are different between LAMMPS and MS, then move the
coordinates of one of the atoms in the molecule slightly, so that you
have two different structures, calculate the energy differences (in
epair, ebond, eangle, edihed, eimp), using LAMMPS and MS and compare.
If there are still differences, then this means there is a bug in
msi2lmp or lammps (but more likely msi2lmp). Knowing which kind of
energy has the bug will help you (or someone else) fix that bug. (But
keep in mind, it will probably fall on your shoulders.)

Finally, to rule out the effect of electrostatics, before you
calculate the energies, set the charge of all the atoms to zero (for
both MS and LAMMPS).

Cheers

Andrew

Case · November 5, 2015, 5:39am

Thank Ray，Axel，Andrew，and others for your kind reply and good suggestions. I did some tests as suggested.
PART I: First, for the small molecule test. CH4 was chosen to clarify the difference between MS and LAMMPS.
The setting is the same as before. The following are the energies for the optimization of CH4 by LAMMPS and MS:
LAMMPS:

akohlmey · November 5, 2015, 9:22am

Thank Ray，Axel，Andrew，and others for your kind reply and good suggestions. I did some tests as suggested.

however, those tests are not looking at the individual contributions
and are still only comparing energies and not forces. thus they are
essentially useless.

[...]

As seen above, the energies of CH4(initial and final) for LAMMPS are 0.41983146 and 0.222921 and they are0.408397 and 0.211576 for MS. The energies are close to each other for MS and LAMMPS, but some decomposed terms of energy are different, e.g. eangle though the ebond’s for MS and LAMMPS are very close. At this point, I am wondering if the eangle in LAMMPS is corresponding to the angle term in MS. In general, how the terms in MS such as Angle, Torsion, Inversion, Stretch-Stretch relate to the energy terms in LAMMPS?

the various "complex" force field terms, that are one of the
differences between class 1 and class 2 force fields are in LAMMPS
internally mapped to angle, bond, dihedral and improper styles,
depending on the number of atoms involved.

[...]

Termination status : Normal
From above energies, we can see that single-PVDF is more different for LAMMPS and MS. In particular , for MS results, the Electrostatic energy is non-zero. The negative value of electrostatic seems not reasonable to me since I
Built a cell(20x20x20) and the electrostatic energy should be very small unless the energy includes not only the interaction between two-PVDF molecule from another periodic cell.

you don't seem to be aware of the implications of using ewald
summation for electrostatic interactions. ...and - again - energies
can have arbitrary offsets. what matters for MD are the forces. if you
want to reduce unknowns and complexity, use cutoff coulomb instead of
long-range electrostatics.

[...]

Anyway, these are what I learned about the issue so far. Since I am new to LAMMPS and force field calculations, I most likely missed some critical points in my analyses, pls enlighten me. Thank you very much for your time and attention.

you have critical comments above. most importantly, what you are
asking the mailing list to do here, is beyond what can be done and is
something that you need to discuss with your adviser and/or more
experienced colleagues, or failing sufficient experience in either,
your adviser needs to find you are suitable tutor. this is not some
minor details that you are getting wrong, but you need to study and
learn some fundamentals of the methods you want to use. it is probably
best to dedicate some time to that, before you revisit the particular
problem at hand. what you are currently trying to do is backwards and
extremely inefficient.

axel.

Case · November 6, 2015, 6:16am

Axel，thanks for reply。I'll learn more basics.

In the meantime,if anyone encounters this problem or has experience in using class2 force field in LAMMPS, it would be very nice of you to take a look at the following questions:
1. Do class2 force field parameters ,(i.e. BondBond, bondangle, middlebondtorsion, endbondtorsion, angle torsion, angleangletorsion,bondbond13,angleangle coeffs ) affect the electrostatic energy? Seems to me, it does not affect much.
2. How to calculate by approximate formula the electrostatic energy given the charges of all atoms in a molecule?
3. Any recommendations of study materials for learning to use
Class2 force field in LAMMPS?

Thanks a lot.

akohlmey · November 6, 2015, 11:51am

Axel，thanks for reply。I'll learn more basics.

In the meantime,if anyone encounters this problem or has experience in using class2 force field in LAMMPS, it would be very nice of you to take a look at the following questions:

we are going in circles here. if you *do* know my basics, you would
not need to ask any of the questions below.

1. Do class2 force field parameters ,(i.e. BondBond, bondangle, middlebondtorsion, endbondtorsion, angle torsion, angleangletorsion,bondbond13,angleangle coeffs ) affect the electrostatic energy? Seems to me, it does not affect much.

no and yes.. all contributions to the potential function in a
conventional force field are additive, i.e. computed independently. so
no, when you change any of those parameters for a given component, it
will not change what you compute for another component. however, if
you change the total potential function, you also change the
equilibrium location, thus after some MD steps the overall geometry
will change, and you will get differences. thus yes, if you change
*any* parameter, it will affect the overall structure.

2. How to calculate by approximate formula the electrostatic energy given the charges of all atoms in a molecule?

no need to approximate. coulomb's law is simple and can be easily
computed exactly for coulomb with cutoff. the expressions for doing
coulomb with ewald summation is also explained in any good text book
on MD (and many papers).

but i have to repeat. absolute energies don't matter, if you want to
compare different MD programs, you need to compare the *forces*.

3. Any recommendations of study materials for learning to use
Class2 force field in LAMMPS?

none of your problems are really class2 specific, but due to lack of
understanding of force fields and MD codes in general (and lack of
following advice).

axel.

Case · November 9, 2015, 6:56am

Axel，thanks for reply。I'll learn more basics.

In the meantime,if anyone encounters this problem or has experience in using class2 force field in LAMMPS, it would be very nice of you to take a look at the following questions:

we are going in circles here. if you *do* know my basics, you would
not need to ask any of the questions below.

1. Do class2 force field parameters ,(i.e. BondBond, bondangle, middlebondtorsion, endbondtorsion, angle torsion, angleangletorsion,bondbond13,angleangle coeffs ) affect the electrostatic energy? Seems to me, it does not affect much.

no and yes.. all contributions to the potential function in a
conventional force field are additive, i.e. computed independently. so
no, when you change any of those parameters for a given component, it
will not change what you compute for another component. however, if
you change the total potential function, you also change the
equilibrium location, thus after some MD steps the overall geometry
will change, and you will get differences. thus yes, if you change
*any* parameter, it will affect the overall structure.

Thanks, "components are independent" is good enough for my question.

2. How to calculate by approximate formula the electrostatic energy given the charges of all atoms in a molecule?

no need to approximate. coulomb's law is simple and can be easily
computed exactly for coulomb with cutoff. the expressions for doing
coulomb with ewald summation is also explained in any good text book
on MD (and many papers).

but i have to repeat. absolute energies don't matter, if you want to
compare different MD programs, you need to compare the *forces*.

3. Any recommendations of study materials for learning to use
Class2 force field in LAMMPS?

none of your problems are really class2 specific, but due to lack of
understanding of force fields and MD codes in general (and lack of
following advice).

axel.

........ .............