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)