Hi, just upgraded to the 29 Sep 2021 release of Lammps, and am playing around with the new Thermo outputs Econserve and Ecouple. More details & background below, but does anyone know if the thermo outputs Ecouple and Econserve account for the effects of Rattle or Shake on the system?
Background
System is a collection of 433 benzene molecules (box size ~40 Angstroms) with the Amber gaff2 forcefield. Seeing the following behavior:
- With Rattle on, constraining all the C-H bond lengths, under NPT conditions, the Econserve value rapidly decreases with time in an essentially linear fashion (simulation time of 15 ns w/ a 2 fs timestep). This is opposite what the documentation states:
The econserve keyword is the sum of the potential and kinetic energy of the system as well as the energy that has been transferred by thermostatting or barostatting to their coupling reservoirs. I.e. it is pe + ke + econserve . Ideally, for a simulation in the NVE, NPH, or NPT ensembles, the econserve quantity should remain constant over time.
-
I ran the same system under NVE, NVT, and NPH conditions to double check everything. I see similar “runaway” behavior of the Econserve and Ecouple values for the NPH simulation. Conversely, I see good energy conservation for the NVE simulation, with a small drift of ~1-2 kcal/mol over 15 ns, and local energy fluctuations of ~3-4 kcal/mol. The Econserve value for the NVT simulation behaves similarly to that from the NVE simulation (drift of ~2-3 kcal/mol and local energy fluctuations of ~5 kcal/mol).
-
In the above simulations, the thermodynamic values of the benzene system are at equilibrium and not trending anywhere, e.g. if I plot the different energies of the system, or the temperature, pressure, volume, or cell side length vs simulation time, they are not trending in any direction, just fluctuating around average values. The simulations all ran for 15 ns, which should be enough time to see long-time behavior given the approximate temperature (~20 C). Similarly, when visualizing the trajectories, nothing looks off.
-
Changing the Rattle tolerance from the original value I used of 1e-6 to values between 1e-2 to 1e-14 had no effect on the simulation results.
-
I reran the above cases with Shake instead of Rattle. Seeing similar behavior, except the Econserve value rapidly increase with time for the NPT & NPH simulations instead of rapidly decreasing like it did with Rattle. Behavior outside that was qualitatively the same
-
I also reran without either Rattle or Shake turned on, with an appropriately reduced time step (0.5 fs, down from 2 fs). Here, the Econserve value for the NPT and NPH simulations look more like a random walk (albeit with rather large excursions). Everything else is qualitatively similar to what I was seeing for the Rattle or Shake simulations.
-
Last, I looked at changing the value of the barostat damping parameter. I ran the above simulations using a value of 1e3. In addition, I tried values of 1e2, 1e4, and 1e5 for the Rattle, Shake, and no-Rattle-or-Shake NPT and NPH simulations. Redcuing the value to 1e2 intensifies the problem of “runaway” of Econserve values. However, increasing the value to 1e4 or 1e5, eliminated the issue, but at the expense of introducing long-time sinusoidal correlations (“ringing”) in the time autocorrelation functions of the volume and non-bonded potential energy terms. So this is not a viable solution due to it introducing these unphysical long-time autocorrelations in the thermo properties.
-
System volumes are close to experimental values, and the RDF of the system is as well, so I don’t believe I implemented the FF incorrectly. I’ve also been able to reproduce Amber energy values using this FF in Lammps for the same single- and two-benzene systems.
-
As a result of the above, I’m thinking that the Econserve (and Ecouple) values may not be accounting for the effects of Rattle or Shake on the system. Is this true?
Example NPT Code
units real
atom_style full
boundary p p p
variable n file npt.nstep
variable q file npt.lstep
variable t file heat.temp # temperature to heat to/run dynamics at
read_restart npt.$q.rst
if “$n==1” then “reset_timestep 0”
kspace_style pppm 1E-6
group mobile union all
timestep 0.5
#thermo_style custom etotal ke pe emol ebond eangle edihed eimp evdwl ecoul elong temp press vol lx ly lz
thermo_style custom elapsed elaplong etotal ke pe emol ebond eangle edihed eimp evdwl ecoul etail ecouple econserve elong temp press vol lx ly lz
thermo_modify line multi
thermo 1000
thermo_modify flush yes
fix 1 mobile npt temp $t $t $(100.0*dt) iso 1.000000 1.000000 $(1000.0*dt)
#fix rattle all rattle 1E-2 500 50 t 2
dump 2 all dcd 1000 npt.$n.dcd
run 2000000
write_restart npt.$n.rst