Hi everyone,
I’m running a shear test simulation on granular materials and encountering an unexpected drop in shear stress right after continuing from a restart file. I’m wondering if I’ve made a mistake in my restart script.
First, I run an initial simulation to set up and start the shearing process. Here is the input script that generates the initial state and saves a restart file called restart.200s.
variable name string large_200s # Name of the simulation
atom_style sphere # Use spherical particles
units si # Use Lennard-Jones (reduced) units
###############################################
# Geometry and box parameters
###############################################
variable scale equal 0.001
# Define the dimensions of the simulation box
variable boxx equal 61.64*${scale}
variable boxy equal 61.64*${scale}
variable boxz equal 64.4*${scale}
# Variables for applying shear deformation
variable xy_shear equal 0 # Shear in the xy plane
variable xz_shear equal 0 # Shear in the xz plane
variable yz_shear equal 0 # Shear in the yz plane
###############################################
# Simulation box and boundary conditions
###############################################
# Set periodic boundaries in all three dimensions
boundary p p p
# Define a triclinic box with shear capability
region boxreg prism 0 ${boxx} 0 ${boxy} 0 ${boxz} ${xy_shear} ${xz_shear} ${yz_shear}
# Create a simulation box within the defined region
create_box 3 boxreg
# Set a random seed for reproducibility of random processes
variable rand_seed equal 1234
###############################################
# Particle properties
###############################################
# Define particle size range and related properties
variable rlo equal 0.45*${scale} # Minimum particle radius
variable rhi equal 0.5*${scale} # Maximum particle radius
variable dlo equal 2.0*${rlo} # Minimum particle diameter
variable dhi equal 2.0*${rhi} # Maximum particle diameter
variable wall_d equal 3*${scale} # wall particle diameter
variable skin equal 0.4*${rhi} # Skin distance for neighbor list
###############################################
# Define regions for applying shear
###############################################
variable Lywall equal 4*${scale} # Thickness of the walls
variable walltopystart equal ${boxz}-${Lywall} # Starting position of the top wall
variable wallskin equal 0.01*${Lywall}
# Define top and bottom wall regions
region top_region block 0 ${boxx} 0 ${boxy} ${walltopystart} ${boxz} units box
region bottom_region block 0 ${boxx} 0 ${boxy} ${wallskin} ${Lywall} units box
region mid_region block 0 ${boxx} 0 ${boxy} ${Lywall} ${walltopystart} units box
###############################################
# initialization
###############################################
variable lattice_length equal 3.2*${scale}
# Define a custom lattice with basis points
lattice custom ${lattice_length} &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.5 0.5 0.5 &
basis 0.0 0.0 0.0
# Randomly create particles within the box
create_atoms 1 random 280328 ${rand_seed} mid_region #360000*0.92^3
create_atoms 2 region top_region basis 2 2
create_atoms 2 region bottom_region basis 2 2
# Assign random diameters to particles within the defined range
variable rand_diam atom random(${dlo},${dhi},${rand_seed})
set type 1 diameter v_rand_diam
set type 2 diameter ${wall_d}
set group all density 500000
###############################################
# Group particles and assign wall properties
###############################################
# Group particles based on regions
group topwallparticles region top_region
group bottomwallparticles region bottom_region
group wallparticles union topwallparticles bottomwallparticles
group midparticles region mid_region
group notwalls subtract all wallparticles
###############################################
# Interparticle interactions
###############################################
# Set the granular pair style for particle interactions
pair_style granular
pair_coeff * * hertz/material 1e7 0.4 0.3 tangential mindlin NULL 1.0 0.3 damping tsuji
# Apply the NVE integrator for particles
fix 1 notwalls nve/sphere
###############################################
# Apply shearing forces and rigid body motion
###############################################
# Fix rigid motion of top wall with force applied along z-axis
fix Fontopwalls topwallparticles rigid/nve single force * off off on torque * off off off
# Fix rigid motion of bottom wall with no external forces
#fix Fonbotwalls bottomwallparticles rigid/nve single force * off off off torque * off off off
###############################################
# Neighbor list and communication settings
###############################################
# Enable communication of velocities between processes
comm_modify vel yes
# Define neighbor list settings with appropriate skin distance
neighbor ${skin} bin
neigh_modify delay 0 every 1 check yes
###############################################
# Output settings
###############################################
#Tensor componenet are listed as: xx, yy, zz, xy, xz, yz
compute stressTensor all stress/atom thermo_temp virial
compute per_atom_vol all voronoi/atom
compute mid_volume midparticles reduce sum c_per_atom_vol[1]
compute stressXZ_vol midparticles reduce sum c_stressTensor[5]
# Dump particle information to a file
dump 1 all custom 10000 out.${name}.dump id type diameter x y z vx vy vz fx fy fz c_stressTensor[5] c_per_atom_vol[1]
###############################################
# Main simulation steps
###############################################
# Set the timestep for integration
timestep 1e-5
# Apply a shear force to the top wall and run the simulation
variable desire_pressure equal 5 / 1000000 * step
variable desire_force equal -1*v_desire_pressure*${boxx}*${boxy}
fix stresstop topwallparticles addforce 0.0 0.0 v_desire_force
###############################################
# Thermodynamic output
###############################################
# Define thermo output format for monitoring pressure and energy
thermo_style custom step temp ke pe press c_stressXZ_vol c_mid_volume v_desire_pressure
thermo_modify lost warn
thermo 100
run 1000000
write_restart restart.stay10s
# Apply shear by setting velocities for top and bottom walls
variable desire_pressure equal 5
variable desire_force equal -1*${desire_pressure}*${boxx}*${boxy}
fix stresstop topwallparticles addforce 0.0 0.0 ${desire_force}
variable vtop equal 1e-2*${scale} # Velocity of top wall
variable vbottom equal -1*${vtop} # Velocity of bottom wall
velocity topwallparticles set ${vtop} 0.0 0.0 units box
fix move_bot bottomwallparticles move linear ${vbottom} 0.0 0.0 units box
# Run the simulation for an additional period under shear
run 19000000
write_restart restart.200s
Next, I use the following script to continue the simulation from restart.200s:
# Simulate a shear test on granular materials using a triclinic box.
# Spherical particles are randomly distributed inside the box and interact via granular mechanics.
# Shear deformation is applied in the xz direction, with controlled interactions between walls and particles.
variable name string large_400s # Name of the simulation
atom_style sphere # Use spherical particles
units si # Use Lennard-Jones (reduced) units
###############################################
# Read Restart File
###############################################
read_restart restart.200s
###############################################
# Restore Necessary Settings
###############################################
# Define the same scale variable as the original script (keeping for potential use)
variable scale equal 0.001
# Define the dimensions of the simulation box
variable boxx equal 61.64*${scale}
variable boxy equal 61.64*${scale}
variable boxz equal 64.4*${scale}
# Set the granular pair style for particle interactions (ensure consistency with the restart)
pair_style granular
pair_coeff * * hertz/material 1e7 0.4 0.3 tangential mindlin NULL 1.0 0.3 damping tsuji # Use the original coefficients
# Fixes are NOT stored in restart, so we need to reapply them
fix 1 notwalls nve/sphere
fix Fontopwalls topwallparticles rigid/nve single force * off off on torque * off off off
# Restore neighbor list settings (adjust skin distance if needed based on particle size)
variable rhi equal 0.5*${scale}
variable skin equal 0.4*${rhi}
neighbor ${skin} bin
neigh_modify delay 0 every 1 check yes
# Restore computes (computes are NOT stored in restart)
compute stressTensor all stress/atom thermo_temp virial
compute per_atom_vol all voronoi/atom
compute mid_volume midparticles reduce sum c_per_atom_vol[1]
compute stressXZ_vol midparticles reduce sum c_stressTensor[5]
# Restore output settings (dumps are NOT stored in restart)
dump 1 all custom 10000 out.${name}.dump id type diameter x y z vx vy vz fx fy fz c_stressTensor[5] c_per_atom_vol[1]
# Restore thermo settings (also not stored in restart)
thermo_style custom step temp ke pe press c_stressXZ_vol c_mid_volume
thermo_modify lost warn
thermo 100
###############################################
# Continue Simulation 1
###############################################
# Set the timestep for integration
timestep 1e-5
# Apply a shear force to the top wall and run the simulation
variable desire_pressure equal 5
variable desire_force equal -1*${desire_pressure}*${boxx}*${boxy}
fix stresstop topwallparticles addforce 0.0 0.0 ${desire_force}
#run 0
#write_restart restart.130s
###############################################
# Continue Simulation 2
###############################################
# Apply shear by setting velocities for top and bottom walls (using variables)
variable vtop equal 1e-2*${scale} # Velocity of top wall
variable vbottom equal -1*${vtop} # Velocity of bottom wall
velocity topwallparticles set ${vtop} 0.0 0.0 units box
fix move_bot bottomwallparticles move linear ${vbottom} 0.0 0.0 units box
run 20000000
write_restart restart.400s
However, when I plot the shear stress (c_stressXZ_vol) over time, I see a sharp, artificial drop at the point of the restart.
I’ve read the LAMMPS documentation which mentions that granular pair styles might not restart exactly the same due to their dependence on half-step velocities. However, the drop seems too immediate to be just a statistical fluctuation.
Is there a mistake in my restart script that could be causing this? Any help would be greatly appreciated! Thanks!