Change in Electrostatics/Pressure Upon Restart of NPT Simulation

LAMMPS LAMMPS General Discussion
1

Hi all,

Relatively new LAMMPS user. Doing some tests, I noticed that when I restart from an NPT (fix NPT) simulation (either by reading a restart file, or continuing a run in the same input script with a second run command), the electrostatics energy (e_coul+e_long) and subsequently the instantaneous “pressure” of the system, both slightly change (the electrostatics energy differs in the tenths digit, the pressure in the ones digit) for the same system state. This seems to be due to a change in the portioning of the electrostatics energy calculation between the short and long range portions, as e_coul and e_long both change by ~10x more than what their sum changes by. All other energy terms, the instantaneous “temperature”, and the system volume/dimensions are all identical before/after restart.

It only matters what type of simulation I restart from, not what type I restart to, e.g. restarting from an NVT simulation and starting an NPT simulation, this behavior does not occur, but restarting from an NPT simulation and starting an NVT simulation, it does.

Not running shake or rattle during these tests

Also, when I force the neighbor list to update every time step (neigh_modify check no delay 0 every 1), the behavior still occurs

I understand that if the partitioning of the electrostatics interactions between the short-range/real space and long-range/k space contributions changes, this will change the total electrostatics energy, due to various approximations/interpolations/floating point math numerics performed during the calculation.

My questions are:

  • Why does this only occur upon restarting from an NPT simulation, and not from an NVT (fix NVT) or NVE (fix NVE) simulation as well?

  • What is the root cause of this change in apportioning of interactions between short- and long-range electrostatics upon restart from an NPT simulation?

Simulation Details below. I’m running LAMMPS version 29 Oct 2020. My test system is 433 benzene molecules using the Amber Gaff2 force field

Any input would be helpful, thanks

Initial Interaction Energy Setup:

pair_style lj/cut/coul/long 10.0
pair_modify mix arithmetic
dielectric 1.0
special_bonds amber
bond_style harmonic
angle_style harmonic
dihedral_style charmm
improper_style cvff

Example Input Script for a Restart

units real
atom_style full
boundary p p p

read_restart restart.rst # restart from previous NPT run

reset_timestep 0 # compare energy at time step 0 vs end of run restarting from
kspace_style pppm 1E-4
group mobile union all
timestep 1
thermo_style custom etotal ke pe emol ebond eangle edihed eimp evdwl ecoul elong temp press vol lx ly lz
thermo_modify line multi
thermo 250
thermo_modify flush yes

neigh_modify check no delay 0 every 1 # issue occurs with or without this command

fix 1 mobile npt temp 293 293 $(100.0*dt) iso 1.000000 1.000000 $(1000.0*dt)
dump 2 all dcd 500 test.dcd
run 10000
write_restart restart.A.rst

run 10000 # compare energy at start of this run step vs end of previous run step
write_restart restart.B.rst

2

There are three factors that need to be considered here (with decreasing impact):

  • when using kspace style pppm you have multiple factors that affect accuracy and the choice of the alpha parameter which determines the damping. With a variable cell (fix npt) those parameters are chosen at the beginning of a run based on the coulomb cutoff and the convergence parameter and default settings (that can be changed with kspace_modify). With different cell volumes you can get slightly different settings and with PPPM also, this can affect the FFT grid. The latter cannot be changed continuously, but in discrete steps. Please note that you are using a rather relaxed convergence of 1.0e-4. That is leaving a rather large error margin. It will result in reasonably well converged forces (still rather aggressive), but stress converges less well and thus for NPT you should using a tighter convergence, say 1.0e-6. Based on that, if the volume of a NPT run changes significantly, it is strongly recommended to restart and thus re-initialize the settings so that they will remain within the desired error margin.
  • when restarting you enforce a neighbor list rebuild, that may change the order of atoms on individual MPI ranks and that in turn affects the summing of forces and energies, which - due to the use of floating point math and its non-associative behavior - can result in subtle changes in energies and forces.
  • when restarting fixes, the (thermostat state) data stored in the restart is only restored when restarting the same fix with the same fix ID.
3

Thanks for the quick response!

More details below around mostly your first bullet point. But before that, for practical purposes, is the solution:

  • Increase PPPM relative accuracy to 1E-6 or 1E-8?
  • Treat the small changes in electrostatics energy/pressure upon restart as negligible/“within the noise”, and to be ignored (similar to effects of shake/rattle on system pressure upon restart)?**

Regarding PPPM

I have tried tests with tighter PPPM accuracy (1E-6 and 1E-8), and the same behavior occurs. As an example, below is output from two sequential run 10000 commands, using 1E-6 kspace error tolerance, with initial restart file coming from an NVT simulation

The PPPM output before and after restart are identical (i.e. when writing restart file at end of first run 10000 command and at beginning of second run 10000 command), but the e_coul+e_long sum (and the pressure) are different ever so slightly for the two step 10000 thermo outputs:

Pre-Restart:

  • E_coul: 5390.4638
  • E_long: -4662.9674
  • E_coul+E_long: 727.4964
  • Pressure: 191.4731

Post-Restart:

  • E_coul: 5371.9158
  • E_long: -4644.4216
  • E_coul+E_long: 727.4942
  • Pressure: 191.4579

(full output at end of post)

So the electrostatics energy is changing for the same system configuration for what I think (?) is the same PPPM parameter set (see output below) between just before and right after the reset. The short range and long range terms change considerably, but their sum (the important number) only changes in thousandths place. And I’m assuming the pressure changes due to the small change in the electrostatics energy/forces. The box size itself is constant during the restart (as it should be). The box side length changes by ~2% during the first 10000 step run, however the PPPM grid size is still 30 30 30 at beginning and end (and constant throughout?)

I’m curious as to why this only occurs when restarting from an NPT simulation?

For your second two points, yes, agree completely. I always keep my fix numbers the same between scripts when doing restarts for the reason you specify.

Simulation Results

First run 10000 command

PPPM Setup

PPPM initialization …
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.29382554
grid = 30 30 30
stencil order = 5
estimated absolute RMS force accuracy = 0.0002539759
estimated relative force accuracy = 7.6484088e-07
using double precision FFTW3
3d grid and FFT values/proc = 6776 1800

Step 0 Thermo Output

---------------- Step 0 ----- CPU = 0.0000 (sec) ----------------
TotEng = 7361.4504 KinEng = 4534.2625 PotEng = 2827.1879
E_mol = 3707.4750 E_bond = 1484.8443 E_angle = 1114.8403
E_dihed = 1023.3749 E_impro = 84.4155 E_vdwl = -1601.4767
E_coul = 5387.2808 E_long = -4666.0912 Temp = 292.8106
Press = 425.4471 Volume = 64000.0000 Lx = 40.0000
Ly = 40.0000 Lz = 40.0000

Step 10000 Thermo Output

---------------- Step 10000 ----- CPU = 43.8061 (sec) ----------------
TotEng = 7474.9852 KinEng = 4439.8798 PotEng = 3035.1053
E_mol = 3777.9161 E_bond = 1488.6982 E_angle = 1110.9299
E_dihed = 1087.3544 E_impro = 90.9336 E_vdwl = -1470.3072
E_coul = 5390.4638 E_long = -4662.9674 Temp = 286.7156
Press = 191.4731 Volume = 67619.0440 Lx = 40.7402
Ly = 40.7402 Lz = 40.7402

PPPM Info for Restart

System init for write_restart …
PPPM initialization …
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.29243697
grid = 30 30 30
stencil order = 5
estimated absolute RMS force accuracy = 0.00026828132
estimated relative force accuracy = 8.0792124e-07
using double precision FFTW3
3d grid and FFT values/proc = 6776 1800

Second run 10000 command

PPPM Info

PPPM initialization …
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.29243697
grid = 30 30 30
stencil order = 5
estimated absolute RMS force accuracy = 0.00026828132
estimated relative force accuracy = 8.0792124e-07
using double precision FFTW3
3d grid and FFT values/proc = 6776 1800

Step 10000 Thermo Output After Restarting

---------------- Step 10000 ----- CPU = 0.0000 (sec) ----------------
TotEng = 7474.9829 KinEng = 4439.8798 PotEng = 3035.1031
E_mol = 3777.9161 E_bond = 1488.6982 E_angle = 1110.9299
E_dihed = 1087.3544 E_impro = 90.9336 E_vdwl = -1470.3072
E_coul = 5371.9158 E_long = -4644.4216 Temp = 286.7156
Press = 191.4579 Volume = 67619.0440 Lx = 40.7402
Ly = 40.7402 Lz = 40.7402

4

It is not clear to me from your description what exactly you are comparing, but please note the following difference between:

and:

Because every time PPPM is being re-initialized without explicitly setting all related parameters, you will get slight differences if the volume is not identical. That won’t happen with fix nve or fix nvt (unless you use some other fix deforming the cell).

You could try adding kspace_modify gewald 0.2935 to your input and see if that will make your restarting results more consistent.

5

Thanks, that makes sense.

These differences seem suitably small that with appropriate kspace error tolerance (1E-6 or lower) that they can be ignored?

6

What matters is not so much the energies but the forces and stresses. With forces there is quite a bit of error cancellation relative to the change in energy.

The only way to be certain is to make some tests to determine the error margins for your specific settings and system.

Also, keep in mind that the KSpace style is not fully re-initialized during a run, so the actual errors will vary as the volume changes. This is particularly problematic when a system is expanding since that lowers the accuracy.