Enforced rotation of domain

User discussions
1

GROMACS version: 2018.5

Hi all,

I have to make an enforced rotation simulation, where I have to rotate a domain in my protein. However, I have no experience with such simulations. I have made an mdp file where I am guessing a bit.

; Enforced rotation
rotation = yes ; apply rotation potential (yes/no)
rot-nstrout = 1 ; (100) output frequency (in steps) for the angle of the rotation group, as well as for the torque and the rotation potential energy
rot-nstout = 10 ; (1000) output frequency for per-slab data of the flexible axis potentials, i.e. angles, torques and slab centers
rot-nqroups = 1 ; number of rotation groups
rot-group0 = ; nome of rotation group 0 in the index file
rot-type0 =iso ; (iso) type of rotation protential applied to rotation group 0 (iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t)
rot-mass0 = no ; (no) use mass weighted rotation group positions
rot-vec0 = 1 0 0 ; rotation vector, will get normalised
rot-pivot0 = 0 0 0 ; pivot point (nm) for potentials iso, pm, rm and rm2
rot-rate0 = 10.0 ; reference rotation rate (degree/ps) for group 0
rot-k0 = 500.0 ; force constance (kJ/(mol*nm^2)) for group 0
rot-fit-method0 = norm ; fetting method to determine angle of rotation group (rmsd, norm, or potential)

I have come to the conclusion after reading about enforced rotation in Gromacs manual that I should use iso as my rotation type as I don’t want a flexible rotation.

How do I decide on the rot-vec0 and rot-pivot0?

Additionally, what else do I need in my mdp file to perform my enforced rotation? And how should I submit it?

Best regards,
Irene

2

You have to specify a rot-vec because that is the axis about which GROMACS will rotate the group for you. The length of the vector does not matter, but the direction will. You could use pymold or VMD to determine a rotation vector that suits your needs, or calculate one from the difference of two atomic positions, e.g. Same is true for the pivot, but you could also use the pivot-free variant of the potential. In that case, no pivot is needed and the center of mass of the rotation group will be used as pivot.

Hope that helps,
Carsten

3

Thank you very much Carsten.

I have tried using https://pymolwiki.org/index.php/RotationAxis to determine my axis, where I use draw_axis between my the domain I wish to rotate and the rest of the protein. I get the following:

/Applications/PyMOL.app/Contents/lib/python2.7/site-packages/pymol/__init__.py:165: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
   _pymol.lock_api_status = threading.RLock() # mutex for PyMOL status info
Transformation (TTT) matrix
-0.97,     0.25,     0.05,    23.97
 0.25,     0.87,     0.43,   -19.72
 0.07,     0.43,    -0.90,    71.74
 0.00,     0.00,     0.00,     1.00
The direction cosines of the rotation axis is: 0.13, 0.97, 0.22
The angle of rotation is 179.43 degrees
The lenght of the translation vector along the rotation axis is 0.03 Angstroms
The distance between mass centers is 15.25 Angstroms
Lines to be used in a pml script to generate the axis
CYLINDER, 20.53, 55.38, 50.98, 5.02, -60.56, 24.23, 0.60, 1.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.0
cmd.load_cgo(obj, 179.43)

I have illustrated my protein after the draw_axis

rot_axis
rot_axis1387×1052 521 KB

Am I doing this right? From this my rot-vec would be 0.13, 0.97, 0.22.

If this is right, which gmx command should I use to submit this?

Best regards,
Irene

4

You simply enter the x, y, and z values of your rotation vector in rot-vec in the .mdp file, e.g. in your case rot-vec0 = 0.13 0.97 0.22 and then use grompp to produce your .tpr file.

5

Thank you very much!
Should I use some specific flags for mdrun after using grompp, or should it just be a normal mdrun submission? Also, how long do you suggest the simulation to be?

gmx grompp -f $mdp -c $gro -p $top -r $gro -n $ndx -t $cpt -o $tprout -maxwarn 5
gmx mdrun -deffnm enforced_rot -s $tprrot -mp topol.top

With

title = Enforced Rot
; Run parameters
integrator = md
nsteps = 2000000
dt = 0.002

I have thinking about something like this?

Also, should this be present in the enforced_rotation mdp?

define = -DPOSRES_B
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 250000 ; 500 ps
nstcomm = 10
; Output parameters
nstxout = 5000 ; every 10 ps
nstvout = 5000
nstfout = 500
nstxtcout = 500 ; every 1 ps
nstenergy = 500
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds
continuation = yes ; continuing from NPT
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; 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
; Berendsen temperature coupling is on in two groups
Tcoupl = Nose-Hoover
tc_grps = Protein Non-Protein
tau_t = 1.0 1.0
ref_t = 310 310
; Pressure coupling is on
Pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 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

6

No, you don’t need any extra flags on the mdrun command line for enforced rotation. Also for the rest of the .mdp input file apart from the enforced rotation options just have it idential to your ‘normal’ simulation without enforced rotation.

Best,
Carsten

7

I have tried to run it a few times now. It seems there are no trouble with grompp, but when I perform my production run I get this error:

Fatal error:
2 particles communicated to PME rank 4 are more than 2/3 times the cut-off out
of the domain decomposition cell of their charge group in dimension x.
This usually means that your system is not well equilibrated.
I have tried to make my simulation box bigger, and smaller. And I have tried to run my NVT and NPT steps for twice as long as I usually do.

Do you have any ideas what I can do differently?

8

Hi

maybe you need to pull slower.

Best,
Carsten

9

Hi Carsten,

That sounds plausible, I will try some different values.

Is the pull rate is decided by rot-k0 or should I also change rot-rate0?

; reference rotation rate (degree/ps) for group 0
rot-rate0 = 10.0
; force constance (kJ/(mol*nm^2)) for group 0
rot-k0 = 100.0 ;I have tried with 500 as well

10

Hi

k is the spring constant and rate is the rotation rate.

Carsten

11

If k is the spring constant, then which flag do I need to change to pull slower?

12

To pull slower, set a smaller value for rot-rate (given in degrees per picosecond).

Carsten