# Pulling seven reaction coordinates **URL:** \ **Category:** User discussions\ **Tags:** pull-code\ **Created:** [July 12, 2023, 3:08pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818 "2023-07-12T15:08:12Z")\ **Posts on this page:** 18\ **Page:** 1
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 12, 2023, 3:08pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/1 "2023-07-12T15:08:12Z")
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! ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 24, 2023, 2:28pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/2 "2023-07-24T14:28:49Z")
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? ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 24, 2023, 4:47pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/3 "2023-07-24T16:47:21Z")
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 \ 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! ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 25, 2023, 10:24am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/4 "2023-07-25T10:24:48Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 25, 2023, 2:29pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/5 "2023-07-25T14:29:24Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/902473a10fae74de9d9e20cd0b1c0ae0e6cc8a89.png) ![pullx](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/e/e6465f59c65af9c9dee5ae3c236179c6cc17702a.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 25, 2023, 2:55pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/6 "2023-07-25T14:55:28Z")
Are these curves for one dimensional pulling? Then there should only be one curve in pullf. Can you post your pull mdp parameters? ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 25, 2023, 3:17pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/7 "2023-07-25T15:17:56Z")
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 ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 25, 2023, 3:37pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/8 "2023-07-25T15:37:15Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 25, 2023, 4:06pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/9 "2023-07-25T16:06:40Z")
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. ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 26, 2023, 10:11am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/10 "2023-07-26T10:11:15Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [July 26, 2023, 3:04pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/11 "2023-07-26T15:04:31Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/4/4926a056dcd4708511de12217ac7fabcde81c642.jpeg) ![pullx-freeze](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/99b2644bfc5d069e834c529f71a0cb59ceab0823.jpeg) 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 ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [July 26, 2023, 3:24pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/12 "2023-07-26T15:24:13Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [August 3, 2023, 10:45am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/13 "2023-08-03T10:45:08Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/8/86ceb3e9e082eb0808858f77caf7dbec82009671.png) 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 ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 3, 2023, 11:18am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/14 "2023-08-03T11:18:48Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [August 3, 2023, 12:43pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/15 "2023-08-03T12:43:38Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/b/b780dcfded28d301b82d593535cc58a0f97a7450.jpeg) ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [August 21, 2023, 7:58pm UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/16 "2023-08-21T19:58:32Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/6/6104deaddd2c944f26ec810c6d8a8437c1cc3cff.png) 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 ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 22, 2023, 6:52am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/17 "2023-08-22T06:52:18Z")
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. ---
**Author:** ![lige1018](https://avatars.discourse-cdn.com/v4/letter/l/cab0a1/32.png) [@lige1018](https://gromacs.bioexcel.eu/u/lige1018)\ **Post date:** [August 22, 2023, 7:49am UTC](https://gromacs.bioexcel.eu/t/pulling-seven-reaction-coordinates/6818/18 "2023-08-22T07:49:16Z")
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?