Pulling seven reaction coordinates

User discussions
1

GROMACS version: GROMACS/2021.5-foss-2021b-CUDA-11.4.1
GROMACS modification: No
Hi,
My structure contains a flat particle standing on an infinite surface. I am trying to pull this flat particle against the surface along the X-axis to test friction. I divided this flat particle into 6 parts, set up 6 reaction coordinates, and pulled these 6 parts simultaneously with ‘pull_coord1_type = constraint’ ‘pull_coord1_geometry = direction-periodic’. Furthermore, I added a 7th reaction coordinate that put loading on the flat particle from the z-axis. After running MD simulation, I got pullf.xvg and pullx.xvg.

  1. I am wondering what pullf.xvg means. It is easy to understand if there is only one reaction coordinate, but I currently have seven reaction coordinates.

  2. I mainly want to know the force constant changes of the spring that pulls this flat particle during the whole pulling, and then know the friction of this flat particle against the surface. It would be appreciated a lot if someone can help. Thanks!

2

Your pullf.xvg file should contain 8 columns. Maybe your plotting program only shows column 2 versus column 1.

I don’t see how the pulling force constant is an interesting parameter. Maybe the pulling speed is of interest?

3

Thank you very much for your reply.

  1. Yes, there are 8 columns when I used xmgrace -nxy pullf.xvg, but now I tried to use only one spring to pull the particle as a whole to calculate the friction of the particle on the surface.
  2. The reason why I am concerned about pullf.xvg is that the change in f can reflect the friction of this particle on the surface. I tested different pull_coord1_k and pull_coord1_rate with umbrella and restraint respectively. When visualizing pullf.xvg (f <here I am a little confused the Y axis f is force constant or force?> vs time) and pullx.xvg (distance vs time), both f vs time and distance vs time are linear. What I am expecting is that when the particle starts to move and then reaches a constant velocity, would f in pullf.xvg increase first and then reach a plateau?
    I am looking forward to your reply!
4

I think I misread your question. Pulling one molecule with one that one reaction coordinate leads to results that are very difficult to interpret.

The pull force can show different behavior, depending on the system, pull parameters and initial conditions. I wouldn’t focus to much on the initial behavior. You need to simulate long enough that the system reaches a steady state. For this you need to pull long enough that your (force) signal is significantly higher than the noise.

5

Thanks for the explanation. I am trying to simulate friction as people always do in Lammps. I tested different pulling rates and spring force constants. For example, I used 0.001 nm/ps, 602 kJ mol^-1 nm^-2 for the constraint to pull the particle in the X direction on a surface. In total, my simulation is 12ns with 0.002 ps as the timestep, so I pulled the particle for 12nm. I attached the figures of pullf.xvg and pullx.xvg from this simulation. The force and distance increase linearly with respect to the simulation time. The pullx.xvg shows that the particle moves at a steady speed, but I am confused about the pullf.xvg curve, why the force is increasing throughout the simulation?

pullf
pullf792×612 5.07 KB

pullx
pullx760×588 3.65 KB

6

Are these curves for one dimensional pulling? Then there should only be one curve in pullf.

Can you post your pull mdp parameters?

7

In addition to the pulling reaction coordinate, I had two other reaction coordinates that press the particle in the z-direction and prevent the surface from moving. I also used enforced rotation to prevent my flat particle from falling down, so I can investigate the friction between my particle edge and surface. Below is my MDP file.

;title = OPLS Lysozyme MD simulation
;define = -DPOSRES ; position restrain the CNT
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 6000000 ; 2 * 50000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
; No output for .trr trjectory.
nstxout = 10000 ; save coordinates every 200.0 ps
nstvout = 10000 ; save velocities every 200.0 ps
nstfout = 10000
nstenergy = 10000 ; save energies every 200.0 ps
nstlog = 10000 ; update log file every 200.0 ps
nstxout-compressed = 10000 ; save compressed coordinates every 200.0 ps
; nstxout-compressed replaces nstxtcout
compressed-x-grps = System ; replaces xtc-grps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tc-grps = System ; Treat the whole system as one
bd_fric = 0 ; Brownian dynamics friction coefficient. When bd-fric is 0, the friction coefficient for each particle is calculated as mass/ tau-t.
tau_t = 1 ; time constant, in ps
ref_t = 0 ; reference temperature, one for each group, in K

;tcoupl = V-rescale ; modified Berendsen thermostat
;tc-grps = System ; Treat the whole system as one group.
;tau_t = 0.1 ; time constant, in ps
;ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
;pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
;pcoupltype = isotropic ; uniform scaling of box vectors
;tau_p = 2.0 ; time constant, in ps
;ref_p = 1.0 ; reference pressure, in bar
;compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
;refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; Pull code
pull = yes
pull_ncoords = 3 ; only one reaction coordinate
pull_ngroups = 3 ; two groups defining one reaction coordinate
pull_nstxout = 10000
pull_nstfout = 10000
pull_xout_average = no
pull_fout_average = no

pull_group1_name = ILL
;pull_group1_pbcatom = 24294
pull_coord1_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord1_geometry = direction ; simple distance increase
pull_coord1_dim = Y N N ; pull along x
pull_coord1_groups = 0 1 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull-coord1-origin = 1.033 0.0 0.0
pull_coord1_init = 0.0
pull_coord1_vec = 1.0 0.0 0.0
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.001 ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k = 602 ; kJ mol^-1 nm^-2

pull_group2_name = ILL
;pull_group2_pbcatom = 23864
;pull_group2_name = MON
;pull_group2_pbcatom = 10359
pull_coord2_type = constant-force ; a linear potential and therefore a constant force. For this option there is no reference ;position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord2_geometry = direction ; simple distance increase
pull_coord2_groups = 0 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord2_dim = N N Y ; pull along z
pull_coord2_origin = 0.0 0.0 7.874
pull_coord2_init = 0.0
pull_coord2_vec = 0.0 0.0 1.0
pull_coord2_start = yes ; define initial COM distance > 0
;pull_coord2_rate = -0.00001 ; 0.001 nm per ps = 1 nm per ns
pull_coord2_k = 602 ; kJ mol^-1 nm^-1

pull_group3_name = MON
pull_group3_pbcatom = 41021
;pull_group3_name = MON
;pull_group3_pbcatom = 10359
pull_coord3_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord3_geometry = direction ; simple distance increase
pull_coord3_groups = 0 3 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord3_dim = Y N N ; pull along z
pull_coord3_origin = -0.222 0.0 0.0
pull_coord3_init = 0.0
pull_coord3_vec = 1.0 0.0 0.0
pull_coord3_start = yes ; define initial COM distance > 0
pull_coord3_rate = 0 ; 0.001 nm per ps = 1 nm per ns
pull_coord3_k = 100000 ; kJ mol^-1 nm^-2

; Enforced rotation
rotation = yes
rot_nstrout = 10000
rot_nstsout = 10000
rot_ngroups = 1
rot_group0 = ILL
rot_type0 = iso
rot_massw0 = no
rot_vec0 = 0 1 0
;rot_pivot0 = 0 0 0
rot_rate0 = 0
rot_k0 = 10
rot_fit_method0 = rmsd

freezegrps = MON
freezedim = N N Y
;pull_coord1_init = 0.0
pull_pbc_ref_prev_step_com = yes

8

There are several issues here, but I am not sure what is causing the increase in pull force. I suspect that the rotation code might reset the group to the original location, try pulling without the rotation code active to check if that is the cause.

Adding a constant force between the molecule and substrate along z increases the friction, is that really what you want?

You should not use an SD integrator when measuring friction, unless you are sure that the friction that SD adds is much smaller than the friction you want to measure.

9

If I do not use a rotation code, the standing flat particle will fall down on the surface at the very beginning when I pull the COM.

I want to calculate friction coefficients (friction coefficient=friction force/loading) under different loadings, so I add a force from z directions.

The reason why I use SD is that I want to use Langevin thermostat. My simulation box includes the particle, the surface, and a few water molecules, so most of the box is empty. I got to know that in this case, the Langevin thermostat is better than other thermostats such as v-rescale that is more suitable for homogeneous systems.

I can try to add restraints on the z coordinates of the particle to prevent falling rather than using enforced rotation, but in this way, the z restraints that prevent falling will compete with the loading I add on the particle, further affecting the friction simulation. It is worth a try to see if the pull force keeps increasing without using enforced rotation.

10

I meant “enforced rotation”, that’s on in your mdp settings.

That’s a valid reason for using SD. But then you should choose the friction for SD small enough so it doesn’t affect the friction you want to measure.

11

I tried to do the pulling without using enforced rotation, and the pullf.xvg and pullx.xvg look like below. I froze the surface in xyz directions to avoid its moving and the particle in z-direction to prevent falling. The MDP file is also attached.
I will try to change the friction for SD to minimize the noise. I am also wondering if there are other ways to make the flat particle stand in addition to ‘enforced rotation’ or ‘add restraint on z coordinates of all atoms in particle’?

pullf-freeze
pullf-freeze792×612 21.8 KB

pullx-freeze
pullx-freeze792×612 15.8 KB

MDP:
;title = OPLS Lysozyme MD simulation
;define = -DPOSRES ; position restrain the CNT
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 6000000 ; 2 * 50000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
; No output for .trr trjectory.
nstxout = 10000 ; save coordinates every 200.0 ps
nstvout = 10000 ; save velocities every 200.0 ps
nstfout = 10000
nstenergy = 10000 ; save energies every 200.0 ps
nstlog = 10000 ; update log file every 200.0 ps
nstxout-compressed = 10000 ; save compressed coordinates every 200.0 ps
; nstxout-compressed replaces nstxtcout
compressed-x-grps = System ; replaces xtc-grps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tc-grps = System ; Treat the whole system as one
bd_fric = 0 ; Brownian dynamics friction coefficient. When bd-fric is 0, the friction coefficient for each particle is calculated as mass/ tau-t.
tau_t = 1 ; time constant, in ps
ref_t = 0 ; reference temperature, one for each group, in K

;tcoupl = V-rescale ; modified Berendsen thermostat
;tc-grps = System ; Treat the whole system as one group.
;tau_t = 0.1 ; time constant, in ps
;ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
;pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
;pcoupltype = isotropic ; uniform scaling of box vectors
;tau_p = 2.0 ; time constant, in ps
;ref_p = 1.0 ; reference pressure, in bar
;compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
;refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; Pull code
pull = yes
pull_ncoords = 2 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_nstxout = 10000
pull_nstfout = 10000
pull_xout_average = no
pull_fout_average = no

pull_group1_name = ILL
;pull_group1_pbcatom = 24294
pull_coord1_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord1_geometry = direction ; simple distance increase
pull_coord1_dim = Y N N ; pull along x
pull_coord1_groups = 0 1 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull-coord1-origin = 1.033 0.0 0.0
pull_coord1_init = 0.0
pull_coord1_vec = 1.0 0.0 0.0
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.001 ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k = 602 ; kJ mol^-1 nm^-2

pull_group2_name = ILL
;pull_group2_pbcatom = 23864
;pull_group2_name = MON
;pull_group2_pbcatom = 10359
pull_coord2_type = constant-force ; a linear potential and therefore a constant force. For this option there is no reference ;position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord2_geometry = direction ; simple distance increase
pull_coord2_groups = 0 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord2_dim = N N Y ; pull along z
pull_coord2_origin = 0.0 0.0 7.874
pull_coord2_init = 0.0
pull_coord2_vec = 0.0 0.0 1.0
pull_coord2_start = yes ; define initial COM distance > 0
;pull_coord2_rate = -0.00001 ; 0.001 nm per ps = 1 nm per ns
pull_coord2_k = 602 ; kJ mol^-1 nm^-1

;pull_group3_name = MON
;pull_group3_pbcatom = 41021
;pull_group3_name = MON
;pull_group3_pbcatom = 10359
;pull_coord3_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
;pull_coord3_geometry = direction ; simple distance increase
;pull_coord3_groups = 0 3 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
;pull_coord3_dim = Y N N ; pull along z
;pull_coord3_origin = -0.222 0.0 0.0
;pull_coord3_init = 0.0
;pull_coord3_vec = 1.0 0.0 0.0
;pull_coord3_start = yes ; define initial COM distance > 0
;pull_coord3_rate = 0 ; 0.001 nm per ps = 1 nm per ns
;pull_coord3_k = 100000 ; kJ mol^-1 nm^-2

; Enforced rotation
;rotation = yes
;rot_nstrout = 10000
;rot_nstsout = 10000
;rot_ngroups = 1
;rot_group0 = ILL
;rot_type0 = iso
;rot_massw0 = no
;rot_vec0 = 0 1 0
;rot_pivot0 = 0 0 0
;rot_rate0 = 0
;rot_k0 = 10
;rot_fit_method0 = rmsd

freezegrps = MON ILL
freezedim = Y Y Y N N Y
;freezegrps = ILL
;freezedim = N N Y
;pull_coord1_init = 0.0
pull_pbc_ref_prev_step_com = yes

12

Those jumps in the force look suspicious. I don’t know if they can be explained by stick-slip effects or rotational jumps.

You might be able to use pull angle coordinates to restrain the orientation of your molecule.

13

Hello again,
I tried to add one pull angle coordinate to restrain the orientation of the flat particle, which can keep the particle perpendicular to the surface, but the particle drifted in the xy plane. I tried to set the reference atom ‘pull_group1_pbcatom’ for the pulling as centered as possible, but the particle still drifted around 45 degrees in the xy plane. Is there a way to make the particle move along the x-axis perfectly while standing vertically along the z-axis?

One good thing is that the pulling force did not increase linearly when I use the pull angle coordinate but fluctuate around zero. I used four reaction coordinates in total. 1. Pulling with constraint in the x-direction. 2. Pressing with constant force in the z-direction. 3. Large constraint in the x-direction to prevent the surface from moving when pulling the particle. 4. Pulling angle coordinate with umbrella to make the particle parallel to the z-axis.

You mentioned that ‘you should choose the friction for SD small enough so it doesn’t affect the friction you want to measure’. Does this mean that I need to increase tau_t since the friction coefficient in the Langevin’s equations is defined as ‘m/tau_t’?

pullf-angle
pullf-angle792×612 5.34 KB

MDP parameters:
;title = OPLS Lysozyme MD simulation
;define = -DPOSRES ; position restrain the CNT
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 5000000 ; 2 * 50000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
; No output for .trr trjectory.
nstxout = 100000 ; save coordinates every 200.0 ps
nstvout = 100000 ; save velocities every 200.0 ps
nstfout = 100000
nstenergy = 100000 ; save energies every 200.0 ps
nstlog = 100000 ; update log file every 200.0 ps
nstxout-compressed = 100000 ; save compressed coordinates every 200.0 ps
; nstxout-compressed replaces nstxtcout
compressed-x-grps = System ; replaces xtc-grps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tc-grps = System ; Treat the whole system as one
bd_fric = 0 ; Brownian dynamics friction coefficient. When bd-fric is 0, the friction coefficient for each particle is calculated as mass/ tau-t.
tau_t = 2 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K

;tcoupl = V-rescale ; modified Berendsen thermostat
;tc-grps = System ; Treat the whole system as one group.
;tau_t = 0.1 ; time constant, in ps
;ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
;pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
;pcoupltype = isotropic ; uniform scaling of box vectors
;tau_p = 2.0 ; time constant, in ps
;ref_p = 1.0 ; reference pressure, in bar
;compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
;refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; Pull code
pull = yes
pull_ncoords = 4 ; only one reaction coordinate
pull_ngroups = 5 ; two groups defining one reaction coordinate
pull_nstxout = 100000
pull_nstfout = 100000
pull_xout_average = no
pull_fout_average = no

pull_group1_name = ILL
pull_group1_pbcatom = 24344
pull_coord1_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord1_geometry = direction-periodic ; simple distance increase
pull_coord1_dim = Y N N ; pull along x
pull_coord1_groups = 0 1 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull-coord1-origin = 7.632 0.0 0.0
pull_coord1_init = 0.0
pull_coord1_vec = 1.0 0.0 0.0
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.001 ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k = 602 ; kJ mol^-1 nm^-2

pull_group2_name = ILL
pull_group2_pbcatom = 24344
;pull_group2_name = MON
;pull_group2_pbcatom = 10359
pull_coord2_type = constant-force ; a linear potential and therefore a constant force. For this option there is no reference ;position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord2_geometry = direction ; simple distance increase
pull_coord2_groups = 0 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord2_dim = N N Y ; pull along z
pull_coord2_origin = 0.0 0.0 6.681
pull_coord2_init = 0.0
pull_coord2_vec = 0.0 0.0 1.0
pull_coord2_start = yes ; define initial COM distance > 0
;pull_coord2_rate = -0.00001 ; 0.001 nm per ps = 1 nm per ns
pull_coord2_k = 602 ; kJ mol^-1 nm^-1

pull_group3_name = MON
pull_group3_pbcatom = 10392
;pull_group3_name = MON
;pull_group3_pbcatom = 10359
pull_coord3_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord3_geometry = direction ; simple distance increase
pull_coord3_groups = 0 3 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord3_dim = Y N N ; pull along z
pull_coord3_origin = 7.479 0.0 0.0
pull_coord3_init = 0.0
pull_coord3_vec = 1.0 0.0 0.0
pull_coord3_start = yes ; define initial COM distance > 0
pull_coord3_rate = 0 ; 0.001 nm per ps = 1 nm per ns
pull_coord3_k = 100000 ; kJ mol^-1 nm^-2

pull_group4_name = part1
pull_group4_pbcatom = 26575
pull_group5_name = part3
pull_group5_pbcatom = 22078
pull_coord4_type = umbrella ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord4_geometry = angle-axis ; simple distance increase
pull_coord4_groups = 4 5 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
;pull_coord4_dim = N Y Y ; pull along z
pull_coord4_origin = 0.0 0.0 0.0
pull_coord4_init = 0.0
pull_coord4_vec = 0.0 0.0 1.0
pull_coord4_start = yes ; define initial COM distance > 0
pull_coord4_rate = 0 ; 0.001 nm per ps = 1 nm per ns
pull_coord4_k = 602000 ; kJ mol^-1 nm^-2

; Enforced rotation
;rotation = yes
;rot_nstrout = 10000
;rot_nstsout = 10000
;rot_ngroups = 1
;rot_group0 = ILL
;rot_type0 = iso
;rot_massw0 = no
;rot_vec0 = 0 1 0
;rot_pivot0 = 0 0 0
;rot_rate0 = 0
;rot_k0 = 10
;rot_fit_method0 = rmsd

;freezegrps = MON ILL
;freezedim = Y Y Y N N Y
;freezegrps = ILL
;freezedim = N N Y
;pull_coord1_init = 0.0
pull_pbc_ref_prev_step_com = yes

14

I can’t follow your orientation issues with just text. You can add a second angle constraint, if that would help.

I would not use a pull coordinate for translating the molecule and one for keeping the surface in place. Instead you can pull to molecule with respect to the surface using a single pull coordinate.

The friction of the SD integrator will affect the measured friction when made small enough. I don’t know if the current SD friction is so small that it does. You can increase or decrease tau_t and see if this affects the measured friction.

15

I made a drawing that might look more clear. Yes, I can try to add a second angle constraint in the y direction to prevent the drift.

The reason why I tether the surface with a strong spring is that the bottom of the particle has a strong interaction with the surface. During pulling, the surface also moves in the x direction. Yes, I can try to use a single pull coordinate to pull the molecule with respect to the surface.

IMG_0372
IMG_03721668×1437 170 KB

16

Hi, I hope you had a nice summer vacation. As you suggested, I tried to use two ‘pull angle’ restraints to keep the particle standing and not drifting in the xy plane. In total, I have four reaction coordinates. 1. Pulling in the x direction. 2. Pressing in the z direction. 3. Pulling angle in the z direction. 4. Pulling angle in the y direction. I tried to simulate in vacuum (no water in the simulation box) and 0k temperature to eliminate other effects and to see if I can get a clean pulling force curve. The pullf.xvg is shown below. I am still confused why the pulling force goes to negative. Another question is that is the pulling force related to the two ‘pull angle’ restraints (reaction coordinates 3 and 4)? Thanks!

image
image797×851 9.62 KB

MDP file:
;title = OPLS Lysozyme MD simulation
;define = -DPOSRES ; position restrain the CNT
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 5000000 ; 2 * 50000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
; No output for .trr trjectory.
nstxout = 50000 ; save coordinates every 200.0 ps
nstvout = 50000 ; save velocities every 200.0 ps
nstfout = 50000
nstenergy = 50000 ; save energies every 200.0 ps
nstlog = 50000 ; update log file every 200.0 ps
nstxout-compressed = 50000 ; save compressed coordinates every 200.0 ps
; nstxout-compressed replaces nstxtcout
compressed-x-grps = System ; replaces xtc-grps
; Bond parameters
continuation = no ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tc-grps = System ; Treat the whole system as one
bd_fric = 0 ; Brownian dynamics friction coefficient. When bd-fric is 0, the friction coefficient for each particle is calculated as mass/ tau-t.
tau_t = 1 ; time constant, in ps
ref_t = 0 ; reference temperature, one for each group, in K
ld_seed = -1

; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is off
; Pull code
pull = yes
pull_ncoords = 4 ; only one reaction coordinate
pull_ngroups = 7 ; two groups defining one reaction coordinate
pull_nstxout = 50000
pull_nstfout = 50000
pull_xout_average = no
pull_fout_average = no

pull_group1_name = ILL
pull_group1_pbcatom = 45007
pull_group2_name = MON
pull_group2_pbcatom = 31791
pull_coord1_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord1_geometry = direction-periodic ; simple distance increase
pull_coord1_dim = Y N N ; pull along x
pull_coord1_groups = 2 1 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull-coord1-origin = 0.0 0.0 0.0
pull_coord1_init = 0
pull_coord1_vec = 1.0 0.0 0.0
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.001 ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k = 6020 ; kJ mol^-1 nm^-2

pull_group3_name = ILL
pull_group3_pbcatom = 45007
;pull_group2_name = MON
;pull_group2_pbcatom = 10359
pull_coord2_type = constant-force ; a linear potential and therefore a constant force. For this option there is no reference ;position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord2_geometry = direction ; simple distance increase
pull_coord2_groups = 0 3 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord2_dim = N N Y ; pull along z
pull_coord2_origin = 0.0 0.0 0.0
pull_coord2_init = 0
pull_coord2_vec = 0.0 0.0 1.0
pull_coord2_start = yes ; define initial COM distance > 0
;pull_coord2_rate = -0.00001 ; 0.001 nm per ps = 1 nm per ns
pull_coord2_k = 6020 ; kJ mol^-1 nm^-1

;pull_group3_name = MON
;pull_group3_pbcatom = 10392
;pull_group3_name = MON
;pull_group3_pbcatom = 10359
;pull_coord3_type = constraint ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
;pull_coord3_geometry = direction ; simple distance increase
;pull_coord3_groups = 0 3 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
;pull_coord3_dim = Y N N ; pull along z
;pull_coord3_origin = 7.479 0.0 0.0
;pull_coord3_init = 0.0
;pull_coord3_vec = 1.0 0.0 0.0
;pull_coord3_start = yes ; define initial COM distance > 0
;pull_coord3_rate = 0 ; 0.001 nm per ps = 1 nm per ns
;pull_coord3_k = 100000 ; kJ mol^-1 nm^-2

pull_group4_name = part1
pull_group4_pbcatom = 47281
pull_group5_name = part3
pull_group5_pbcatom = 42784
pull_coord3_type = umbrella ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord3_geometry = angle-axis ; simple distance increase
pull_coord3_groups = 4 5 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
;pull_coord3_dim = N Y Y ; pull along z
pull_coord3_origin = 0.0 0.0 0.0
pull_coord3_init = 0.0
pull_coord3_vec = 0.0 0.0 1.0
pull_coord3_start = yes ; define initial COM distance > 0
pull_coord3_rate = 0 ; 0.001 nm per ps = 1 nm per ns
pull_coord3_k = 602000 ; kJ mol^-1 nm^-2

pull_group6_name = part4
pull_group6_pbcatom = 43496
pull_group7_name = part6
pull_group7_pbcatom = 46568
pull_coord4_type = umbrella ; a linear potential and therefore a constant force. For this option there is no reference position and therefore the parameters pull-coord1-init and pull-coord1-rate are not used.
pull_coord4_geometry = angle-axis ; simple distance increase
pull_coord4_groups = 6 7 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
;pull_coord4_dim = N Y Y ; pull along z
pull_coord4_origin = 0.0 0.0 0.0
pull_coord4_init = 0.0
pull_coord4_vec = 0.0 1.0 0.0
pull_coord4_start = yes ; define initial COM distance > 0
pull_coord4_rate = 0 ; 0.001 nm per ps = 1 nm per ns
pull_coord4_k = 602000 ; kJ mol^-1 nm^-2

; Enforced rotation
;rotation = yes
;rot_nstrout = 100000
;rot_nstsout = 100000
;rot_ngroups = 1
;rot_group0 = ILL
;rot_type0 = iso
;rot_massw0 = no
;rot_vec0 = 0 1 0
;rot_pivot0 = 0 0 0
;rot_rate0 = 0
;rot_k0 = 10
;rot_fit_method0 = rmsd

;freezegrps = MON ILL
;freezedim = Y Y Y N N Y
freezegrps = mon_middle_half
freezedim = Y Y Y
;pull_coord1_init = 0.0
pull_pbc_ref_prev_step_com = yes

17

The pull force fluctuates because of fluctuations in the system. I can’t say if the average is significantly positive or negative.

I don’t understand your second question.

18

Thanks. The second question is that if we look at the green curve (the pull angle restraint in the z direction), most of the time it was negative, which means the particle tends to lean forward in the pulling, so the pull angle restraint tries to prevent this process by adding a negative force. Will this negative force counteract some pulling force?