Direct continuation of pulling after changing mdp parameters

User discussions
,
1

GROMACS version: 2024.4
GROMACS modification: No

Dear community,

I am attempting a direct continuation of a pulling experiment in which I change the pulling rate midway. As oriented in the user guide, I am calling grompp to generate a new .tpr file with the updated pulling rate, and then calling mdrun with the -cpi flag for continuation. Since I am changing the pulling rate, I have also included init-step in the .mdp file.

For context, this is a coarse-grained protein-protein simulation. My goal is to accelerate the initial force loading by applying a pulling rate of 5e-07 nm/ps (which I refer to as the pre-stretch phase) and then switching to 5e-08 nm/ps to save time during this initial phase.

What I’m observing is an immediate loss of tension following the continuation, after which the force increases as expected until rupture. Is it possible to achieve a continuation without this loss of tension, or is it an inherent behavior by the physics and the code?

Below are the .mdp for both stages.

5e-07 mdp

define                   = -DPOSRES_ANCHOR
integrator               = md
dt                       = 0.02
nsteps                   = 575000000
nstcomm                  = 100
nstxout                  = 1000000000
nstvout                  = 1000000000
nstfout                  = 1000000000
nstlog                   = 1000000000
nstenergy                = 92000
nstxout-compressed       = 9200
compressed-x-precision   = 100
cutoff-scheme            = Verlet
nstlist                  = 40
ns_type                  = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005
coulombtype              = reaction-field
rcoulomb                 = 1.1
epsilon_r                = 15
epsilon_rf               = 0
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 1.1
tcoupl                   = v-rescale
tc-grps                  = system
tau_t                    = 1.0
ref_t                    = 310
Pcoupl                   = no
Pcoupltype               = isotropic
tau_p                    = 12.0 
compressibility          = 3e-4 
ref_p                    = 1.0  
gen_vel                  = no
gen_temp                 = 310
gen_seed                 = 473529
refcoord-scaling         = COM
constraints              = none
constraint_algorithm     = Lincs
pull                     = yes
pull_ncoords             = 1          
pull_ngroups             = 2          
pull_group1_name         = anchor     
pull_group2_name         = pull       
pull_coord1_type         = umbrella
pull_coord1_geometry     = direction-periodic
pull-coord1-vec          = 0 0 1 
pull_coord1_groups       = 1 2
pull_coord1_start        = yes
pull_coord1_rate         = 5e-07
pull_coord1_k            = 418
pull_nstxout             = 9200
pull_nstfout             = 9200

5e-08 mdp

define                   = -DPOSRES_ANCHOR
integrator               = md
dt                       = 0.02
nsteps                   = 2500000000
nstcomm                  = 100
init-step                = 5750000000
nstxout                  = 10000000000
nstvout                  = 10000000000
nstfout                  = 10000000000
nstlog                   = 10000000000
nstenergy                = 40000
nstxout-compressed       = 4000
compressed-x-precision   = 100
cutoff-scheme            = Verlet
nstlist                  = 40
ns_type                  = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005
coulombtype              = reaction-field
rcoulomb                 = 1.1
epsilon_r                = 15
epsilon_rf               = 0
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 1.1
tcoupl                   = v-rescale
tc-grps                  = system
tau_t                    = 1.0
ref_t                    = 310
Pcoupl                   = no
Pcoupltype               = isotropic
tau_p                    = 12.0 
compressibility          = 3e-4 
ref_p                    = 1.0  
gen_vel                  = no
gen_temp                 = 310
gen_seed                 = 473529
refcoord-scaling         = COM
constraints              = none
constraint_algorithm     = Lincs
pull                     = yes
pull_ncoords             = 1          
pull_ngroups             = 2          
pull_group1_name         = anchor     
pull_group2_name         = pull       
pull_coord1_type         = umbrella
pull_coord1_geometry     = direction-periodic
pull-coord1-vec          = 0 0 1 
pull_coord1_groups       = 1 2
pull_coord1_start        = yes
pull_coord1_rate         = 5e-08
pull_coord1_k            = 418
pull_nstxout             = 4000
pull_nstfout             = 4000

Thank you for the assistance!

2

Honestly I think this might be the correct behavior.

You are pulling with a spring. After a while (depending on the spring constant), the reference of the pulling will be far away ‘enough’ from the group you are pulling to actually exert a force and start deforming/dragging/etc. At a certain point you stop the simulation and compile a new tpr. This will start with the spring again at ‘position 0’ and will have to start loading again before exerting force on the groups you are pulling.

I do not know how you can circumvent this, but I am pretty sure that specifying init-step won’t change this behavior. Probably you can use some mdp options like pull-coord1-init to set the starting distance in a way that gives instantly rise to the force you want!