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)