Inconsistency of Etot caused by restart

Dear Lammps users,

I am using LAMMPS version 29 Oct 2020.

I am simulating the evaporation of a nickel nitrate nanodroplet (consists of Ni2+, NO3-, H2O). The droplet is surrounded by bathgas N2. I used fix nvt command to control the temperature of the droplet and the surrounding N2. “fix evaporate” command is used to deleted evaporated water molecules.
I tried to restart my simulation using read_restart command at 600ps and 1200ps. However, I noticed that E_tot is discontinuous at the restart timestep.
Following is my input file for the 1st run:

variable            DP      equal 20 #nm, droplet diameter
variable            totH2O  equal 420426 #total number of atoms in water
variable            T0      equal 300  #K
variable            T1      equal 450  #K
variable            P0      equal 0.0  #atm
variable            P1      equal 1.0  #atm
variable            cutlj   equal   10.0  # LJ cutoff
variable            cutcoul equal   10.0  # Coulomb cutoff
variable            RELAXTIME equal 20000
variable            PRODTIME  equal 2000000
variable            dt1 equal 0.1
variable            dt2 equal 1.0
variable            dmp1 equal 100000
variable            dmp2 equal 10000
variable            evapRadius equal $(v_DP*10/2.0*3.0) #A, region to delete H2O
variable            RDN  equal 4975

boundary            p p p
atom_style          full
units               real
bond_style          harmonic
angle_style         harmonic
improper_style      harmonic
pair_style          lj/cut/coul/long ${cutlj} ${cutcoul}
special_bonds       lj/coul 0.0 0.0 0.5
kspace_style        pppm 0.0001

read_data           Ni_d20_cmax_450K_N2_cfg1.data

pair_coeff 1 1 0.0118 2.4460
pair_coeff 1 2 0.0486 3.1730
pair_coeff 1 3 0.0428 2.8000
pair_coeff 1 4 0.0428 2.8060
pair_coeff 1 5 0.0000 1.2230
pair_coeff 1 6 0.0292 2.8830
pair_coeff 2 2 0.2001 3.9000
pair_coeff 2 3 0.1762 3.5270
pair_coeff 2 4 0.1763 3.5330
pair_coeff 2 5 0.0000 1.9500
pair_coeff 2 6 0.1204 3.6100
pair_coeff 3 3 0.1551 3.1540
pair_coeff 3 4 0.1552 3.1600
pair_coeff 3 5 0.0000 1.5770
pair_coeff 3 6 0.1060 3.2370
pair_coeff 4 4 0.1553 3.1660
pair_coeff 4 5 0.0000 1.5830
pair_coeff 4 6 0.1060 3.2430
pair_coeff 5 5 0.0000 0.0000
pair_coeff 5 6 0.0000 1.6600
pair_coeff 6 6 0.0724 3.3200

bond_coeff          1  525.0   1.27  # N-O bond in NO3-
angle_coeff         1  105.0   120.0  # O-N-O angle in NO3-
improper_coeff      1  60.0    0.0  # NO3-
bond_coeff          2  600.0   1.0  #SPC/E water, O-H bond
angle_coeff         2  75.0    109.47  #SPC/E water, H-O-H angle
bond_coeff          3  600.0   1.106  #N2 bond length fixed

group               metalion type 1
group               ntr  type 2 3
group               spce type 4 5
group               rgd  type 2 3 4 5 6
group               drop type 1 2 3 4 5
group               bathgas type 6 
group               ON type 3
group               OW type 4
region              OUTSPHERE sphere 0.0 0.0 0.0 ${evapRadius} side out units box
region              INSPHERE  sphere 0.0 0.0 0.0 ${evapRadius} side in  units box
fix                 blc all balance 1000 1.1 shift xyz 10 1.05
neigh_modify        every 1 delay 0 check yes

variable            a equal xcm(drop,x)*(-1.0)
variable            b equal xcm(drop,y)*(-1.0)
variable            c equal xcm(drop,z)*(-1.0)
displace_atoms      all move $(v_a) $(v_b) $(v_c) units box

compute             1 drop com
compute             2 metalion group/group ntr kspace no molecule inter
compute             3 metalion group/group spce kspace no molecule inter
compute             4 metalion coord/atom cutoff 2.5 group ON
compute             5 metalion coord/atom cutoff 2.5 group OW
compute             6 metalion reduce ave c_4  # CN of Ni-ON
compute             7 metalion reduce ave c_5  # CN of Ni-OW
compute             8 drop chunk/atom bin/sphere 0.0 0.0 0.0 1.0 $(v_DP*5.0+11.0) $(v_DP*5+10.0) discard yes
compute             TDROP drop temp
compute             TBATH bathgas temp

thermo_style        custom step time atoms temp c_TDROP c_TBATH pe ke etotal evdwl ecoul epair &
ebond eangle eimp elong c_1[1] c_1[2] c_1[3] c_2 c_3
thermo              1000

fix                 freeze drop setforce 0.0 0.0 0.0
minimize            1.0e-12 1.0e-12 10000 10000
unfix               freeze
reset_timestep      0

velocity            bathgas create ${T1} ${RDN} dist gaussian
velocity            drop create ${T0} ${RDN} dist gaussian
timestep            ${dt2}
fix                 1 bathgas nvt temp ${T1} ${T1} $(100.0*v_dt2)

dump                1 all custom ${dmp2} bathgas_relax.lammpstrj &
id mol type q x y z vx vy vz
dump_modify         1 sort id
run                 ${RELAXTIME}
#run                 1000

reset_timestep      0
variable            numH2O equal count(spce,INSPHERE)
undump              1
dump                1 all custom ${dmp2} Ni_d${DP}_cmax_evap${T1}K.lammpstrj &
id mol type element q x y z vx vy vz
dump_modify         1 sort id element Ni N O O H N
compute             NiOW drop rdf 200 1 4 cutoff 10.0
compute             NiON drop rdf 200 1 3 cutoff 10.0

fix                 massDens drop ave/chunk 2 50 1000 8 density/mass file massDensity.txt
fix                 numON ON ave/chunk 2 50 1000 8 density/number file ON.num
fix                 numOW OW ave/chunk 2 50 1000 8 density/number file OW.num
fix                 numMetal metalion ave/chunk 2 50 1000 8 density/number file metal.num

fix                 grNiOW drop ave/time 1 1000 10000 c_NiOW[1] c_NiOW[2] c_NiOW[3] &
file NiOW_rdf.txt mode vector
fix                 grNiON drop ave/time 1 1000 10000 c_NiON[1] c_NiON[2] c_NiON[3] &
file NiON_rdf.txt mode vector
fix                 fCN metalion ave/time 1 100 1000 c_6 c_7 file Ni_O.coord

fix                 2 drop nvt temp ${T1} ${T1} $(100.0*v_dt2)
fix                 3 spce evaporate 1000 ${totH2O} OUTSPHERE 56987 molecule yes
fix                 fRattle rgd rattle 0.0001 20 0 b 1 2 3 a 1 2
fix                 CENTER drop recenter 0.0 0.0 0.0 shift all
fix                 4 all print 1000 "${numH2O}" file numH2O_inDrop screen no

restart             100000  Ni_d${DP}_cmax_evap${T1}K.restart
run                 ${PRODTIME}
#run                 1000
write_restart       Ni_d${DP}_cmax_evap${T1}K_2ns.restart

The 1st run have written a restart file at 600ps – Ni_d20_cmax_evap450K.restart.600000.

Then I used read_restart command to continue the 2nd run from 600ps. The input is as follows:

variable            DP      equal 20 #nm, droplet diameter
variable            totH2O  equal 420426 #total number of atoms in water
variable            T0      equal 300  #K
variable            T1      equal 450  #K
variable            P0      equal 0.0  #atm
variable            P1      equal 1.0  #atm
variable            cutlj   equal   10.0  # LJ cutoff
variable            cutcoul equal   10.0  # Coulomb cutoff
variable            RELAXTIME equal 20000
variable            PRODTIME  equal 2000000
variable            dt1 equal 0.1
variable            dt2 equal 1.0
variable            dmp1 equal 100000
variable            dmp2 equal 10000
variable            evapRadius equal $(v_DP*10/2.0*3.0) #A, region to delete H2O
variable            RDN  equal 4975

boundary            p p p
atom_style          full
units               real
bond_style          harmonic
angle_style         harmonic
improper_style      harmonic
pair_style          lj/cut/coul/long ${cutlj} ${cutcoul}
special_bonds       lj/coul 0.0 0.0 0.5
kspace_style        pppm 0.0001

read_restart        Ni_d20_cmax_evap450K.restart.600000

pair_coeff 1 1 0.0118 2.4460
pair_coeff 1 2 0.0486 3.1730
pair_coeff 1 3 0.0428 2.8000
pair_coeff 1 4 0.0428 2.8060
pair_coeff 1 5 0.0000 1.2230
pair_coeff 1 6 0.0292 2.8830
pair_coeff 2 2 0.2001 3.9000
pair_coeff 2 3 0.1762 3.5270
pair_coeff 2 4 0.1763 3.5330
pair_coeff 2 5 0.0000 1.9500
pair_coeff 2 6 0.1204 3.6100
pair_coeff 3 3 0.1551 3.1540
pair_coeff 3 4 0.1552 3.1600
pair_coeff 3 5 0.0000 1.5770
pair_coeff 3 6 0.1060 3.2370
pair_coeff 4 4 0.1553 3.1660
pair_coeff 4 5 0.0000 1.5830
pair_coeff 4 6 0.1060 3.2430
pair_coeff 5 5 0.0000 0.0000
pair_coeff 5 6 0.0000 1.6600
pair_coeff 6 6 0.0724 3.3200

bond_coeff          1  525.0   1.27  # N-O bond in NO3-
angle_coeff         1  105.0   120.0  # O-N-O angle in NO3-
improper_coeff      1  60.0    0.0  # NO3-
bond_coeff          2  600.0   1.0  #SPC/E water, O-H bond
angle_coeff         2  75.0    109.47  #SPC/E water, H-O-H angle
bond_coeff          3  600.0   1.106  #N2 bond length fixed

group               metalion type 1
group               ntr  type 2 3
group               spce type 4 5
group               rgd  type 2 3 4 5 6
group               drop type 1 2 3 4 5
group               bathgas type 6 
group               ON type 3
group               OW type 4
region              OUTSPHERE sphere 0.0 0.0 0.0 ${evapRadius} side out units box
region              INSPHERE  sphere 0.0 0.0 0.0 ${evapRadius} side in  units box
fix                 blc all balance 1000 1.1 shift xyz 10 1.05
neigh_modify        every 1 delay 0 check yes

compute             1 drop com
compute             2 metalion group/group ntr kspace no molecule inter
compute             3 metalion group/group spce kspace no molecule inter
compute             4 metalion coord/atom cutoff 2.5 group ON
compute             5 metalion coord/atom cutoff 2.5 group OW
compute             6 metalion reduce ave c_4  # CN of Ni-ON
compute             7 metalion reduce ave c_5  # CN of Ni-OW
compute             8 drop chunk/atom bin/sphere 0.0 0.0 0.0 1.0 $(v_DP*5.0+11.0) $(v_DP*5+10.0) discard yes
compute             TDROP drop temp
compute             TBATH bathgas temp

thermo_style        custom step time atoms temp c_TDROP c_TBATH pe ke etotal evdwl ecoul epair &
ebond eangle eimp elong c_1[1] c_1[2] c_1[3] c_2 c_3
thermo              1000
dump                1 all custom ${dmp2} Ni_d${DP}_cmax_evap${T1}K_600ps_cont.lammpstrj &
id mol type element q x y z vx vy vz
dump_modify         1 sort id element Ni N O O H N

timestep            ${dt2}
variable            numH2O equal count(spce,INSPHERE)
compute             NiOW drop rdf 200 1 4 cutoff 10.0
compute             NiON drop rdf 200 1 3 cutoff 10.0

fix                 massDens drop ave/chunk 2 50 1000 8 density/mass file massDensity.txt
fix                 numON ON ave/chunk 2 50 1000 8 density/number file ON.num
fix                 numOW OW ave/chunk 2 50 1000 8 density/number file OW.num
fix                 numMetal metalion ave/chunk 2 50 1000 8 density/number file metal.num

fix                 grNiOW drop ave/time 1 1000 10000 c_NiOW[1] c_NiOW[2] c_NiOW[3] &
file NiOW_rdf.txt mode vector
fix                 grNiON drop ave/time 1 1000 10000 c_NiON[1] c_NiON[2] c_NiON[3] &
file NiON_rdf.txt mode vector
fix                 fCN metalion ave/time 1 100 1000 c_6 c_7 file Ni_O.coord

fix                 1 bathgas nvt temp ${T1} ${T1} $(100.0*v_dt2)
fix                 2 drop nvt temp ${T1} ${T1} $(100.0*v_dt2)
fix                 3 spce evaporate 1000 ${totH2O} OUTSPHERE 56987 molecule yes
fix                 fRattle rgd rattle 0.0001 20 0 b 1 2 3 a 1 2
fix                 CENTER drop recenter 0.0 0.0 0.0 shift all
fix                 4 all print 1000 "${numH2O}" file numH2O_inDrop screen no

restart             100000  Ni_d${DP}_cmax_evap${T1}K.restart
run                 ${PRODTIME}
write_restart       Ni_d${DP}_cmax_evap${T1}K_final.restart

Then I found that Etot(600ps) = -9173432.6, Etot(601ps) = -9172760.1 in the 1st run. In the 2nd run, Etot(600ps, 2nd run) = -9173432.6 = Etot(600ps, 1st run). However, Etot(601ps, 2nd run) = -9230152.1 which is significantly different from Etot(601ps, 1st run) = -9172760.1.

Is this energy discontinuity a bug of LAMMPS or if I have done something wrong?
Any suggestion will be of great help.
Thank you very much!

Best wishes,
Dingyu

Hi Dingyu,

Please next time submit your inputs as upload or external link, it is difficult to read as it is.

Indeed I don’t think that the energy should be different after a restart… Can you try to use write_restart to generate plain text file instead of binary ? (according to the doc binary restart are more exact, so it should not improve your result, but it’s worth the try)

Simon

This is difficult to say given the complexity of your inputs. There can be some divergence after a restart, but it should not grow as quickly as you are seeing. The first thing to check would be to create dump files with rather frequent output after the restart point from both the original run and the restarted run and compare those. Then I would do the same after removing fix balance and fix evaporate. Both of those are not guaranteed to restart perfectly and thus may result in small differences and those may become bigger quickly.