LAMMPS version: 2 Aug 2023, built from source using CMake.
I have encountered an unexpected issue when writing restart files during a simulation.
I start from the same initial configuration (State A) and run an NVT simulation for 10000 timesteps. Without writing restart files, the system evolves to a final configuration (State B). However, when I enable periodic writing of restart files during the run, the simulation no longer reaches State B after 10000 timesteps, but instead ends in a different configuration (State C).
To investigate this, I prepared the minimal input script shown below and ran it on 8 MPI processes.
When use_restart = 0, the simulation starts with a potential energy of approximately 1.25 (after minimization) and ends with a potential energy of approximately 1.75 after 10000 steps.
When use_restart = 1, the only difference is that restart files are written periodically during the run. The simulation still starts from the same initial state with a potential energy of approximately 1.25, but after 10000 steps the final potential energy is approximately 1.70 instead of 1.75.
My understanding is that writing restart files should not alter the trajectory of the simulation. Has anyone encountered similar behavior, or can suggest a possible explanation for this discrepancy?
Minimal input script:
variable use_restart equal 0
units lj
dimension 2
boundary p p p
atom_style atomic
atom_modify map yes
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes
newton on
variable rho equal 1.2
variable N equal 500
variable Temp equal 0.5
variable L equal ({N}/{rho})^(1.0/2)
region box block 0 {L} 0 {L} -0.5 0.5
create_box 1 box
create_atoms 1 random ${N} 12345 box
mass 1 1.0
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
fix enforce2d all enforce2d
velocity all create ${Temp} 987654 dist gaussian
timestep 0.005
variable elapsed_time equal elapsed
minimize 1e-6 1e-8 1000 10000
reset_timestep 0
compute pressure all pressure NULL virial
variable SXY equal c_pressure[4]
variable PE equal pe
variable KE equal ke
thermo 100
thermo_style custom step pe vol ke temp press atoms v_SXY
fix nvt1 all nvt temp {Temp} {Temp} 0.5
fix extra all print 1000 ‘${elapsed_time} {PE} {KE} {SXY}’ file Dynamics\_{use_restart}.dat
if “${use_restart} == 0” then “restart 1000 restart_${use_restart}_*.equil”
run 10000
print “Final PE = {PE}, Final KE={KE}, Final SXY = ${SXY}”
unfix nvt1
unfix enforce2d
uncompute pressure