Flat-bottom pull

User discussions
1

GROMACS version: 2024
GROMACS modification: No
Dear GROMACS users,

I have two molecules to be pulled toward each other. I need to the pulling force to stop so I set the a flat-bottom pull, thus the force can become zero once the two molecules are close enough. I set the pull-coord1-init as the target distance that considered to be close, for example 1.5 nm.

However, when running simulation, I noticed that the force drops to zero even far beyond the pull-coord1-init distance 1.5 nm. The force become zero when the two molecules are still 1.8 nm away. I used a K value 500 and rate 0.001. How is the pull force adjusted according to the pull-coord1-init? Should I set other parameters accordingly?

Thanks in advance!

2

One thing springs to mind, are you using pull-coord1-start = yes? Would you be able to copy-paste the pulling settings of your mdp file?

3

Hi MagnusL, thanks for your reply. I have definitely set the pull-coord1-start=no. Below please see my copy pasted .mdp file. This is just one example. I have removed my previous problematic trials. There are some variable parameters that are changed for different conditions, and those are the lines above the “end of variable” comment.

Also I later noticed that if using simple umbrella pull, that is, setting pull-coord1-type=umbrella, if pull rate is positive, the force direction also reverses at a position not equal to reference position. If using zero pull rate, the force direction seems fine. I am not sure why pull rate affects the force this way.

title = Umbrella pulling simulation
define = -DPOSRES_DNA ;-DPOSRES_B

nsteps = 1000000
pull_coord1_k = 300 ; kJ mol^-1 nm^-2
pull_coord1_init = 1.5095794205794206 ; set init position for window center
pull_coord1_rate = 0
pull-group1-pbcatom =335
pull-group2-pbcatom =335

;END_OF_VARIABLES
;-------------------------------
; NOTE: calculate energygrps slow down simulation significantly

; Run parameters
integrator = md
dt = 0.002
tinit = 0
;nstcomm (100) number of steps that elapse between calculating the energies
nstcomm = 100 ; 10 nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider setting nstcomm equal to nstcalcenergy for less overhead

; Output parameters
nstxout-compressed = 5000 ;
nstenergy = 5000

; Bond parameters
constraint_algorithm = lincs
constraints = h-bonds ; all-bonds, “the forcefield has been parametrized only with constraints involving hydrogen atoms”
continuation = yes

; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist = 10 ; 20
ns_type = grid
rlist = 1.1 ;1.4
rcoulomb = 1.0 ;1.4
rvdw = 1.0 ;1.4

; PME electrostatics parameters
coulombtype = PME
fourierspacing = 0.125 ;0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes

; Berendsen temperature coupling is on in two groups
Tcoupl = V-rescale
tc_grps = system
tau_t = 1.0
ref_t = 310

; Pressure coupling is on
Pcoupl = C-rescale ;Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0 ;1.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = com

; 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 = yes
pull_ncoords = 1
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = DNA
pull_group2_name = DBD

pull_coord1_groups = 1 2
pull_coord1_type = flat-bottom ;umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase.

;should confine the pull direction to z only, otherwise the box size may contradict the pull distance on y or x direction
; pull_coord1_dim = Y Y Y
pull_coord1_dim = N N Y

; set ref point as restoring location
pull_coord1_start = no ; define initial COM distance > 0
;pull_coord1_rate = 0

pull-pbc-ref-prev-step-com = yes

4

Strange. I don’t have any ideas right now. If you set pull-print-com = yes and pull-print-ref-value = yes and then look at pullx.xvg and pullf.xvg at the positions where the pull force is not as you expect, does that help you at all? Posting a few such lines here might give some hints.

5

Hi MagnusL,

Thanks for your suggestion. I tried printing out the ref value. Now what I have in mdp is:
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_init = 1.5095794205794206 ; set init position for window center
pull_coord1_rate = 0.001

pull_coord1_groups = 1 2
pull_coord1_type = umbrella
pull_coord1_geometry = distance
pull_coord1_dim = N N Y
pull_coord1_start = no ; define initial COM distance > 0

Here is what I get in pullx.xvg, the reference position moves (the third column).
|0.0000|2.99189|1.50958|1.9934|4.98529|
|0.1000|2.98964|1.50968|1.99483|4.98447|
|0.2000|2.98467|1.50978|1.99692|4.9816|
|0.3000|2.981|1.50988|1.99563|4.97663|
|0.4000|2.9755|1.50998|1.99427|4.96977|
|0.5000|2.97129|1.51008|1.99298|4.96427|
|0.6000|2.96862|1.51018|1.99216|4.96077|
|0.7000|2.96705|1.51028|1.99257|4.95962|
…
|99.3000|1.92075|1.60888|1.99266|3.9134|
|99.4000|1.92049|1.60898|1.99231|3.91279|
|99.5000|1.924|1.60908|1.99106|3.91506|
|99.6000|1.9274|1.60918|1.9901|3.9175|
|99.7000|1.93027|1.60928|1.99036|3.92063|
|99.8000|1.93062|1.60938|1.99108|3.9217|
|99.9000|1.92915|1.60948|1.99355|3.9227|
|100.0000|1.92872|1.60958|1.9945|3.92322|

It seems the problem can be fixed by using zero pull rate. But I don’t understand why this happens. Thanks.

6

The reference position is the position you are pulling towards (i.e., what the first column - the actual coordinate - is trying to reach). The reference coordinate should change at exactly the same rate as the pull-rate.

7

And the fourth and fifth columns show the positions of the two pull groups. Since your pull coordinate is a distance their difference is the same as the coordinate in the second column.

8

I see. Thanks a lot for your explaination!