Title: Question on fix smd equilibrium position persistence across run and restart

LAMMPS LAMMPS General Discussion
,
1

Hello everyone,

I’m running a simple SMD test case in LAMMPS, where a spring pulls a single particle. The relevant commands are:

create_atoms 1 single 0.03 0.03 0
fix pull tlattice smd cvel 100 0.1 tether 0.02 0.03 0.0 0

According to the documentation, for cvel mode, the spring’s equilibrium position is incremented linearly over time based on the specified velocity.

However, I’m encountering an issue when splitting the simulation into multiple runs:

run 3000
run 3000

or when using a restart:

run 3000
restart 3000 restart/restart*.bishear

Then in a new input script:

read_restart restart\\restart3000.bishear
fix pull tlattice smd cvel 100 0.1 tether 0.02 0.03 0.0 0
run 3000

In both cases, it seems the spring equilibrium position is re-initialized after each new fix smd command, and computed from the current position of the particle and the tether point, instead of continuing from the previous state. This causes the spring force to reset to zero, breaking the continuity of the force evolution.

In the output (see attached force plot), you can clearly see that after the second run or restart, the spring force drops back to zero instead of continuing smoothly.

My question:

  • Is there a way to persist the internal state of fix smd (specifically the current equilibrium position) across multiple run commands or after a restart?
  • Or is there a recommended way to perform continuous SMD pulling across runs or simulations without force discontinuity?

Any suggestions or workarounds would be greatly appreciated.

Thanks in advance!
Here are the output curves of the spring force and equilibrium position.
force
l0

The following script can reproduce my question：

#! bin/bash
# file header
atom_style  sphere
atom_modify	map array
dimension   3
boundary	p p f
newton		off
comm_modify	vel yes
units		si

region		 reg block 0 0.06 0 0.06 -0.04 0.046 units box
create_box	 1 reg
neighbor	 0.001 bin
neigh_modify every 1 delay 0

# 
pair_style granular
pair_coeff * * hertz 36630036630.03663 0.9 tangential mindlin 45248868778.28054 1.0 0.5 damping coeff_restitution 
                
pair_style      gran/hertz/history 37e9 45e9 1000 1000 0.5 1
pair_coeff	* *

create_atoms 1 single 0.03 0.03 0 

group tlattice id 1:1
set         group tlattice diameter 0.005 density 2499.9999999999995
velocity    tlattice zero linear
velocity    tlattice zero angular
fix		    trlc tlattice aveforce NULL NULL NULL

fix pull tlattice smd cvel 100 0.1 tether 0.02 0.03 0.0 0

timestep   0.00001

compute    1 tlattice com

thermo_style    custom step atoms f_trlc[1] c_1[1]
thermo    10
thermo_modify    lost ignore norm no

shell    mkdir post
shell    mkdir restart

variable m_time equal time

fix integr all nve/sphere

fix ave_data all ave/time 1 1 10 v_m_time f_trlc[1] c_1[1] f_pull[1] f_pull[5] f_pull[6] file plate.txt title1 "" title2 "step time sf lx pullf l0 dx"

dump dmp all cfg 1000 post/dump*.cfg mass type xs ys zs id type x y z vx vy vz fx fy fz radius

restart 3000 restart/restart*.bishear

run 3000

run 3000
2

Hi @jing-lee,

If the fix command is new, then it is different and this is expected behavior. However looking into the code it looks likes the r0 is indeed reset at the beginning of each run through its init method. So what you see is normal but not explicitely documented as far as I can tell.

However, from the read_restart doc page I can read:

As indicated in the above list, the fixes used for a simulation are not stored in the restart file. This means the new input script should specify all fixes it will use. However, note that some fixes store an internal “state” which is written to the restart file. This allows the fix to continue on with its calculations in a restarted simulation. To re-enable such a fix, the fix command in the new input script must be of the same style and use the same fix-ID as was used in the input script that wrote the restart file.

If a match is found, LAMMPS prints a message indicating that the fix is being re-enabled. If no match is found before the first run or minimization is performed by the new script, the “state” information for the saved fix is discarded. At the time the discard occurs, LAMMPS will also print a list of fixes for which the information is being discarded. See the doc pages for individual fixes for info on which ones can be restarted in this manner. Note that fixes which are created internally by other LAMMPS commands (computes, fixes, etc) will have style names which are all-capitalized, and IDs which are generated internally.

This means that the fix definition must be done before reading the restart file. Have you also tried that? This is not what you show in your script extracts.

2 Likes
3

You’re absolutely right. However, in my script, I actually have to use the fix command after read_restart, because I need to read back the simulation box and atom configurations first before applying any fixes.

Interestingly, when I run the script this way, LAMMPS correctly recognizes the fix pull and outputs:

Resetting global fix info from restart file:
fix style: smd, fix ID: pull
So it appears that the fix is successfully matched and restored.

That said, based on the documentation, it seems that the fix smd does not store the spring equilibrium position (r0) in the restart file. It only saves the pulling direction, current target distance, and the running PMF:

“The fix stores the direction of the spring, current pulling target distance and the running PMF to binary restart files.”
So even though the fix is restored, it looks like r0 is re-initialized each time, which explains the behavior I’m seeing.

In any case, I’ve found a workaround by manually computing and applying the spring force. It’s a bit rudimentary, but works for now — and I hope it might be helpful for others facing similar issues. Here’s a snippet:

variable k_s equal 100.0         # spring stiffness
variable c_vel equal 0.1         # pulling velocity
variable l0 equal 5.0            # initial pulling position

variable x_target equal time*v_c_vel + v_l0
variable dx equal v_x_target - c_1[1]
variable pull_force equal -v_k_s * v_dx

fix pull tlattice addforce v_pull_force 0.0 0.0
4

R0 is not the equilibrium position but an input parameter and should be set to 0.0 for almost all applications. Fix smd has inherited that from fix spring that it is based on.
The distance to the tether or couple point and the pulling direction are stored.
See the source code.