LAMMPS energy conservation error for triclinic systems using pppm k-space solver

Hello All,

I have been modeling molecular, triclinic, crystalline systems using LAMMPS and I have been having problems with energy conservation using the new pppm k-space solver that was recently implemented for triclinic systems. The specific version that I have been having problems with is the June 17, 2013 version of LAMMPS (lammps-17Jun13).

The attached PDF shows plots of total system energy for six systems that have identical initial conditions, but with different options for k-space solvers.

The plots show NVT equilibration to T=300K for the first 405,000 time steps during which the system total energy fluctuates, followed by 2,000,000 time steps of NVE integration (using a step size of 0.1 fs giving 200 ps of NVE trajectory). Using NVE integration, the ewald and ewald/disp methods give consistent energy conservation, whereas the pppm solver gives a steady, monotonically increasing total system energy. The increase in total energy found using the pppm solver is on the order of the fluctuation in the system kinetic energy. A visualization of the trajectory does not show any problems with pppm.

I have verified (shown in PDF) that the ewald/disp k-space solver conserves energy whether one computes the r-6 term to infinity or simply truncates it at a cutoff of Rc=11 angstroms, indicating that the issue with pppm is not due to neglecting the dispersion beyond the cutoff. Moreover, varying the cutoff distance used with the pppm solver does not give a physically expected result–one would expect the total system energy to decrease with increasing cutoff distance as more of the negative-valued dispersion tail is considered in the calculation, but instead one sees an increase.

System Details
-The simulation cell corresponds to a triclinic 4x4x6 supercell of crystalline TATB (192 molecules, 4608 atoms) where
lx=36.010
ly=31.333
lz=38.427
xy=-18.93
xz=-1.608
yz=-14.333
Thus the cutoff distances are appropriate for the simulation cell.
-All bonds, angles, and impropers are modeled as harmonic oscillators and dihedrals are modeled with single cosine terms.
-Intramolecular nonbonded interactions are excluded.
-Intramolecular potential energies are identical at t=0 for the six systems considered, but nonbonded terms differ.
-The accuracy threshold for all k-space calculations was set to 10-6 and I allowed LAMMPS to automatically set the pppm mesh parameters.
-The pppm jobs were run using a single node (8 CPUs) and the ewald jobs were run using three nodes (24 CPUs).

Attached Files
-PDF of plots of total system energy for six systems
-Sample input script corresponding to the top-middle figure in the PDF.

Thank you for any assistance you can provide in resolving these issues with energy conservation.

Matt Kroonblawd
University of Missouri - Columbia
Department of Chemistry

energy_comparison.pdf (392 KB)

in.thermalize (3.54 KB)

Hi,

the pppm with the ik differentiation, the default in LAMMPS, can not conserve
energy. You could try to use the

kspace_modify diff ad

command to use the ad differentiation. You can find more on that in

@article{doi:10.1021/ct2001792,
author = {Ballenegger, V. and Cerdà, J. J. and Holm, C.},
title = {How to Convert SPME to P3M: Influence Functions and Error Estimates},
journal = {Journal of Chemical Theory and Computation},
volume = {8},
number = {3},
pages = {936-947},
year = {2012},
doi = {10.1021/ct2001792},

URL = {http://pubs.acs.org/doi/abs/10.1021/ct2001792},
eprint = {http://pubs.acs.org/doi/pdf/10.1021/ct2001792}
}

Best,

Rolf

Stan can comment on this.

Steve