Energies from previous run do not match while using restart files

LAMMPS LAMMPS General Discussion
, ,
1

Hi all,
I am running a npt simulation of a system consisting a nanoparticle (NP) immeresed in a solvent mixture (contains salts and organic solvents). I noticed that when I continue the simulation using restart files from the previous run, the potential energy and columbic pairwise energy do not start from the same value from the previous run. I have attached a short input script and the log files. I also noticed that when I exclude the interactions between the NP using the neigh_modify command, the differnces between the energies from the restart file and the previous run is larger. I believe the problem lies in the kspace solver but not sure what it is exactly. Hope someone can help. Thanks.

Vaishali

include "../system.in.init"

read_data "../system.data"
# read_restart test.restart.500

include "../system.in.settings"


group NP type 1 2 3 4 5 6
group LP30 type 7 8 9 10 11 12 13 14 15 16 17 18 19 20
group OH type 1 4

neigh_modify exclude group NP NP

print "--------- beginning simulation (using fix npt) ---------"

velocity        LP30 create 300 49245 rot yes mom yes dist gaussian
dump            1 all custom 1000 traj_npt2.lammpstrj id mol type x y z ix iy iz


thermo_style    custom step pe etotal evdwl ecoul temp press vol 
thermo          500  
restart 500 test.restart

fix fxMobile LP30 npt temp 300 300 100 iso 1.0 1.0 1000 drag 1 dilate LP30

compute tempMobile LP30 temp
fix_modify fxMobile temp tempMobile


timestep 1.0   
run		500

the log file:

PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
  G vector (1/distance) = 0.27818548
  grid = 90 90 90
  stencil order = 5
  estimated absolute RMS force accuracy = 0.00034430813
  estimated relative force accuracy = 1.0368737e-06
  using double precision FFTW3
  3d grid and FFT values/proc = 140608 97200
Generated 190 of 190 mixed pair_coeff terms from arithmetic mixing rule
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (src/neighbor.cpp:654)
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 14
  ghost atom cutoff = 14
  binsize = 7, bins = 12 12 12
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 58.98 | 59.77 | 61.3 Mbytes
   Step         PotEng         TotEng         E_vdwl         E_coul          Temp          Press          Volume    
         0  -124691.18     -85059.223      195945.96      316626.21      280.73691      35184.126      555412.25    
       500  -299504.45     -186230.92      5556.2468      286587.09      802.38428     -53097.14       449896.5

with restart file generated from the above run:

include "../system.in.init"

# read_data "../system.data"
read_restart test.restart.500

include "../system.in.settings"


group NP type 1 2 3 4 5 6
group LP30 type 7 8 9 10 11 12 13 14 15 16 17 18 19 20
group OH type 1 4

neigh_modify exclude group NP NP

print "--------- beginning simulation (using fix npt) ---------"

# velocity        LP30 create 300 49245 rot yes mom yes dist gaussian
dump            1 all custom 1000 traj_npt2.lammpstrj id mol type x y z ix iy iz


thermo_style    custom step pe etotal evdwl ecoul temp press vol 
thermo          500  
restart 500 test.restart

fix fxMobile LP30 npt temp 300 300 100 iso 1.0 1.0 1000 drag 1 dilate LP30

compute tempMobile LP30 temp
fix_modify fxMobile temp tempMobile


timestep 1.0   
run		0

the log file from run using restart file:

PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
  G vector (1/distance) = 0.28200344
  grid = 90 90 90
  stencil order = 5
  estimated absolute RMS force accuracy = 0.00028092041
  estimated relative force accuracy = 8.4598347e-07
  using double precision FFTW3
  3d grid and FFT values/proc = 140608 97200
Generated 190 of 190 mixed pair_coeff terms from arithmetic mixing rule
All restart file global fix info was re-assigned
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (src/neighbor.cpp:654)
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 14
  ghost atom cutoff = 14
  binsize = 7, bins = 11 11 11
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 61.87 | 62.35 | 63.2 Mbytes
   Step         PotEng         TotEng         E_vdwl         E_coul          Temp          Press          Volume    
       500  -304443.06     -191169.54      5556.2468      291640.05      802.38428     -53100.458      449896.5
2

The vdW energy tells that the configuration is indeed the same:

But without looking at the system.in and system.in.settings files is impossible to say. This is a hint of the difference observed in the electrostatic energy:

WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (src/neighbor.cpp:654)
3

Your volume has changed significantly from the beginning to the end of your first run (555412.25 vs 449896.5) that results in a change of the (automatically) chosen k-space cutoff (see G-vector) and thus the pre-computed structure factors when PPPM is re-initialized. For efficiency reasons, LAMMPS does not recompute those values during a run, so the change in energy is expected since with PPPM there are only certain choices possible due to the use of grids that must be decomposible for the 3d-Fourier transform. This change in energy should be smaller with a tighter PPPM convergence. For fix npt this should be tighter than the commonly used 1.0e-4 since that only leads to reasonably well converged forces, but not to a similarly well converged pressure tensor. Try with 1.0e-6.

1 Like
4

For completeness, if you force your restart input to reuse the exact same long-ranged electrostatics settings with

kspace_modify gewald  0.27818548 grid 90 90 90

you should get much closer results.

2 Likes
5

Thanks for your reply. I was using a pppm convergence of 1.0e-4 in the actual runs which is when I noticed this difference. The volume did equilibrate eventually by 10 ns, however, even after a run of 80 ns my system energy didn’t seem to converge. This is the plot of TotEnergy vs time.

image
image888×525 21.1 KB

And this the plot of the energy gradient with time, wherein I noticed these unusual energy dips wherever I started a run with the restart files.

image
image866×530 14.3 KB

Then I ran some short test runs with different pppm tolerence and the log files that I shared above are from using a convergence of 1.0e-6. It did seem to improve the energy differences between subsequent runs and the least difference was observed when I didn’t include the neigh_exclude command. However for my system I do not require the interactions within the NP as I treat it as a rigid body. So should I include those interactions as it also warns me about inconsistent Coulombic energies. Since the energy differences does get shorter with simuation time, and probably will be less by the end of the equlibration run, should I be bothered about it?

This is the sys.in.init:

# Intialization
units          real
dimension      3
boundary       p p p
atom_style     full
neighbor 2 bin
neigh_modify delay 0 every 1 check yes
kspace_style    pppm 0.000001

# Atom Definition
pair_style lj/cut/coul/long 12.0 12.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_modify mix arithmetic
special_bonds amber

sys.in.settings:

pair_coeff 1 1 0 0
pair_coeff 2 2 0.0958 3.42
pair_coeff 3 3 0.0977 3.51
pair_coeff 4 4 0.0929 3.29
pair_coeff 5 5 3.295 1.9
pair_coeff 6 6 3.295 1.9
bond_coeff 1 780.83 1.0
bond_coeff 2 1912.046 1.9
bond_coeff 3 1912.046 1.9
bond_coeff 4 1912.046 1.9
bond_coeff 5 1912.046 1.9
bond_coeff 6 1912.046 1.9
bond_coeff 7 1912.046 1.9
angle_coeff 1 129.78 114.85
angle_coeff 2 129.78 114.85
pair_coeff 7 7 0.2000000 3.740
pair_coeff 8 8 0.0610000 3.118
pair_coeff 9 9 0.0005000 2.870
bond_coeff 8 370.4589 1.606
angle_coeff 3 139.2209 90.00
angle_coeff 4 34.77535 180.00
pair_coeff 10 10 0.0005000 2.870
pair_coeff 11 11 0.21 2.96
pair_coeff 12 12 0.105 3.75
pair_coeff 13 13 0.17 3.0
pair_coeff 14 14 0.066 3.5
pair_coeff 15 15 0.015 2.42
bond_coeff 9 570.0 1.229
bond_coeff 10 214.0 1.327
bond_coeff 11 320.0 1.41
bond_coeff 12 340.0 1.09
angle_coeff 5 83.0 123.4
angle_coeff 6 33.0 107.8
angle_coeff 7 35.0 109.5
angle_coeff 8 83.0 116.9
angle_coeff 9 92.6 111.55
dihedral_coeff 1 0.0 5.124 0.0 0.0
dihedral_coeff 2 0.0 0.0 0.198 0.0
dihedral_coeff 3 -0.375 1.358 0.004 0.0
improper_coeff 1 10.5 180.0
pair_coeff 16 16 0.21 2.96
pair_coeff 17 17 0.105 3.75
pair_coeff 18 18 0.17 3.0
pair_coeff 19 19 0.066 3.5
pair_coeff 20 20 0.015 2.42
bond_coeff 13 570.0 1.229
bond_coeff 14 214.0 1.327
bond_coeff 15 268.0 1.529
bond_coeff 16 320.0 1.41
bond_coeff 17 340.0 1.09
angle_coeff 10 83.0 123.4
angle_coeff 11 50.0 109.5
angle_coeff 12 33.0 107.8
angle_coeff 13 35.0 109.5
angle_coeff 14 37.5 110.7
angle_coeff 15 83.0 116.9
angle_coeff 16 92.6 111.55
dihedral_coeff 4 0.0 5.124 0.0 0.0
dihedral_coeff 5 -0.55 0.0 0.0 0.0
dihedral_coeff 6 0.0 0.0 0.468 0.0
dihedral_coeff 7 0.0 0.0 0.3 0.0
dihedral_coeff 8 -1.22 -0.126 0.422 0.0
dihedral_coeff 9 0.0 0.0 0.198 0.0
dihedral_coeff 10 -0.375 1.358 0.004 0.0
improper_coeff 2 10.5 180.0
set type 1 charge 0.417
set type 2 charge -1.035
set type 3 charge -1.124
set type 4 charge -0.913
set type 5 charge 2.159
set type 6 charge 2.248
set type 7 charge 1.34
set type 8 charge -0.39
set type 9 charge 1.0
set type 10 charge 1.00076
set type 11 charge -0.632
set type 12 charge 0.918
set type 13 charge -0.374
set type 14 charge -0.072
set type 15 charge 0.101
set type 16 charge -0.611
set type 17 charge 0.891
set type 18 charge -0.425
set type 19 charge 0.149
set type 20 charge 0.068