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