# Enforced rotation of domain

**URL:** <https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516>\
**Category:** User discussions\
**Created:** [July 16, 2020, 7:05pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516 "2020-07-16T19:05:23Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 16, 2020, 7:05pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/1 "2020-07-16T19:05:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [July 17, 2020, 8:32am UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/2 "2020-07-17T08:32:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 17, 2020, 11:30am UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/3 "2020-07-17T11:30:39Z")

</div>

Thank you very much Carsten.

I have tried using [https://pymolwiki.org/index.php/RotationAxis](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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/1X/5d488b0a0313fcda107fb09e56e4285a2be07b26.jpeg)

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

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [July 17, 2020, 12:25pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/4 "2020-07-17T12:25:05Z")

</div>

> [@Irene](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 17, 2020, 12:32pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/5 "2020-07-17T12:32:40Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [July 17, 2020, 2:27pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/6 "2020-07-17T14:27:56Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 19, 2020, 4:50pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/7 "2020-07-19T16:50:39Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [July 19, 2020, 6:19pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/8 "2020-07-19T18:19:29Z")

</div>

Hi

maybe you need to pull slower.

Best,  
Carsten

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 19, 2020, 6:26pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/9 "2020-07-19T18:26:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [July 20, 2020, 5:43am UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/10 "2020-07-20T05:43:23Z")

</div>

Hi

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

Carsten

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [July 20, 2020, 10:41am UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/11 "2020-07-20T10:41:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\
**Post date:** [August 31, 2020, 1:45pm UTC](https://gromacs.bioexcel.eu/t/enforced-rotation-of-domain/516/12 "2020-08-31T13:45:02Z")

</div>

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

Carsten