[gmx-users] Pulling ligand - Different Profiles (Force vs time)

[Steven Neumann]

Dear Gmx Users,

I obtained a protein-ligand complex from 100ns simulation. Now I am
pulling my ligand away from the protein after the energy minimzation
in water and equilibration of 100ps (two coupling baths: Protein,
LIG_Water_and_ions).
Then I proceed my pulling :

grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr

mdrun -s pull.tpr -deffnm pull


title       = Umbrella pulling simulation
define      = -DPOSRES
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 500000    ; 1 ns
nstcomm     = 10
; Output parameters
nstxout     = 0
nstvout     = 0
nstfout     = 500
nstxtcout   = 1000       ; every 1 ps
nstenergy   = 500
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes       ; continuing from NPT
; Single-range cutoff scheme
nstlist     = 5
ns_type     = grid
rlist       = 1.4
rcoulomb    = 1.4
rvdw        = 1.2
vdwtype     = Switch
rvdw-switch = 1.0
; PME electrostatics parameters
coulombtype     = PME
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Temperature coupling is on
tcoupl      = V-rescale                                  ; modified
Berendsen thermostat
tc_grps     = Protein LIG_Water_and_ions   ; two coupling groups - more accurate
tau_t       = 0.1   0.1                                 ; time constant, in ps
ref_t       = 298   298                                  ; reference
temperature, one for each group, in K
; Pressure coupling is on
Pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 2.0
compressibility = 4.5e-5
ref_p           = 1.0
; Generate velocities is off
gen_vel     = no
; Periodic boundary conditions are on in all directions
pbc     = xyz
; Long-range dispersion correction
DispCorr    = EnerPres
; Pull code
pull            = umbrella
pull_geometry   = distance  ; simple distance increase
pull_dim        = N N Y
pull_start      = yes       ; define initial COM distance > 0
pull_ngroups    = 1
pull_group0     = Protein
pull_group1     = LIG
pull_rate1      = 0.004      ; 0.004 nm per ps = 4 nm per ns
pull_k1         = 500      ; kJ mol^-1 nm^-2

I run 3 pulling simulations with the same mdp  and I obtain 3
different profiles (Force vs time). Then I used 2xlonger pulling with
the same pulling distance and I run 3 simulations again. Each time I
obtain different profile. Can anyone explain me this? I am using
velocities from npt simulation as above (gen_vel = no and continuation
= yes) so I presume the output should be similar. Please, advice.

Steven