[lammps-users] Miscalculation of Dihedrals in n-Alkane system

Hello all,

I am using the 16Mar18 version.

I have been attempting to run an n-Dodecane simulation using the United Atom model along with the OPLS model to restrain the dihedrals in the molecules. Recently, I was testing each individual intra-molecular interaction (bond, angles, and dihedrals) to test if there was any problem with the parameters I was inputting into lammps since I kept getting discrepancies in the vapor density of my system compared to results found in literature. I ran a simulation containing a single dihedral, 2 bending angles, and 3 stretching bonds to compare the energy computed by lammps for each interaction to hand calculations. The details of each interaction are the following:

Stretching : Harmonic Style, (Energy = 95.88199 kcal/molA^2), (Ro = 1.54 Angsts.)
Bending : Harmonic Style, (Energy = 62.00997 kcal/mol
rad^2), (Theta_o = 114 degrees)
Dihedral: OPLS Style, (K1 = 1.411, K2 = -0.271, K3 = 3.145, K4 = 0.0, all in kcal/mol)

My method of comparing these calculations was to output the individual energies using thermo and the ebond, eangle, and edih keywords. Also, I have outputted the bond distances, bending angles, and dihedral angles using the local computes to use in my hand calculations.

The results obtained from lammps area the following:

E_Bond = 8.9944e-5 kcal/mol

E_Angle = 0.00079077 kcal/mol

E_Dih = 1.3299e-8 kcal/mol

My hand results are the following:

E_Bond = 8.09435e-5 kcal/mol

E_Angle = 0.00079005 kcal/mol

E_Dih = 8.721836e-9 kcal/mol

As seen by the results, the bond and angle energies are pretty close with a 1% and 0.09% errors respectively. The issue arises with the dihedral which deviated by a 30% from what my hand calculations suggest. These energies were obtained from the initial structure at timestep 0, in order to avoid any integration of the motion equations that could affect my calculations.

So my question is, is there something that I am missing in terms of how lammps calculates these values? I looked around the source code for the OPLS dihedral and did not find any mayor difference in the way I calculated it.

Any insight, critic, or suggestion in this manner is appreciated.

Thank you.

Jesus Gutierrez
[email protected]

Hello all,

I am using the 16Mar18 version.
I have been attempting to run an n-Dodecane simulation using the United Atom model along with the OPLS model to restrain the dihedrals in the molecules. Recently, I was testing each individual intra-molecular interaction (bond, angles, and dihedrals) to test if there was any problem with the parameters I was inputting into lammps since I kept getting discrepancies in the vapor density of my system compared to results found in literature. I ran a simulation containing a single dihedral, 2 bending angles, and 3 stretching bonds to compare the energy computed by lammps for each interaction to hand calculations. The details of each interaction are the following:

Stretching : Harmonic Style, (Energy = 95.88199 kcal/molA^2), (Ro = 1.54 Angsts.)
Bending : Harmonic Style, (Energy = 62.00997 kcal/mol
rad^2), (Theta_o = 114 degrees)
Dihedral: OPLS Style, (K1 = 1.411, K2 = -0.271, K3 = 3.145, K4 = 0.0, all in kcal/mol)

My method of comparing these calculations was to output the individual energies using thermo and the ebond, eangle, and edih keywords. Also, I have outputted the bond distances, bending angles, and dihedral angles using the local computes to use in my hand calculations.

The results obtained from lammps area the following:

E_Bond = 8.9944e-5 kcal/mol

E_Angle = 0.00079077 kcal/mol

E_Dih = 1.3299e-8 kcal/mol

My hand results are the following:

E_Bond = 8.09435e-5 kcal/mol

E_Angle = 0.00079005 kcal/mol

E_Dih = 8.721836e-9 kcal/mol

As seen by the results, the bond and angle energies are pretty close with a 1% and 0.09% errors respectively. The issue arises with the dihedral which deviated by a 30% from what my hand calculations suggest. These energies were obtained from the initial structure at timestep 0, in order to avoid any integration of the motion equations that could affect my calculations.

So my question is, is there something that I am missing in terms of how lammps calculates these values? I looked around the source code for the OPLS dihedral and did not find any mayor difference in the way I calculated it.

Assuming you are using “dihedral_style opls”, the code that calculates the energy is located near line 196 of the “dihedral_opls.cpp” file (in the DihedralOPLS::compute() function), which you can find here:
https://github.com/lammps/lammps/blob/5ea9d9702424a2d1086fc6a483c11d25b463e7bf/src/MOLECULE/dihedral_opls.cpp#L196

For testing, I suggest running LAMMPS using a simple DATA file that only contains only one dihedral interaction with known angles and 4 atoms (which are not coplanar or collinear). I have attached a test file that I use for testing other dihedral styles in LAMMPS. (“4atoms.data”. Bond lengths are 1, bond angles are 120 degrees, and the dihedral angle is 60 degrees.) Download the “4atoms.data” and “test_dihedral.in” files, and run

lmp_mpi -i test_dihedral.in

(I am using the same OPLS parameters you are using.)

The number on the far right is the dihedral energy. (I ran the simulation for 0 timesteps, so this is the initial dihedral energy.)

If you are familiar with gdb, you can halt the code at this point and query “phi” and other variables to see where the discrepancy might have occurred. Otherwise, you can edit this file and use printf() to print out this information.

Hope this helps.

4atoms.data (540 Bytes)

test_dihedral.in (1.32 KB)

I just read your email more carefully.
Is a difference between 8.721836e-9 and 1.3299e-8 kcal/mol significant for your application? What kind of simulation are you running where you need to have energy accuracy less than 5e-9 kcal/mol. (This is more than 100 million times smaller than the thermal energy at 300K).

Andrew