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.
-The simulation cell corresponds to a triclinic 4x4x6 supercell of crystalline TATB (192 molecules, 4608 atoms) where
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).
-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.
University of Missouri - Columbia
Department of Chemistry
energy_comparison.pdf (392 KB)
in.thermalize (3.54 KB)