[gmx-users] How to implement dihedral restraints

[chris.neale at utoronto.ca]

> For example, my peptide is rotating 50 degrees about a bond B-C, out of four
> bonded atoms A-B-C-D.  I want to apply additional force on this
> dihedral(B-C) during my simulation(only for this particular atoms to make it
> rotate further till 180 degrees). All the remaining atoms in the peptide
> will have normal force as usual.
> May I ask you, how can I apply this in gromacs code.
> Any help will be appreciated!

The current manual is good at describing what the options do for  
dihedral restraints, but (to my knowlegde) doesn't explain at all how  
to get the implementation up and running. Here is how I have done it  
based on searching the mailing list for answers. I suggest that  
something like this is added by way of example to the manual for  
future releases.

Also, the manual is a bit unclear about whether this type of dihedral  
restraint is stable for use near 180deg. I have found that for my  
system everything appears to behave normally and as expected over the  
entire range of dihedral angles including 180deg.

In your .top file:

; Include forcefield parameters
#include "ffoplsaa.itp"
; Include topologies
#include "myprotein.itp"
[ dihedral_restraints ]
; ai   aj    ak    al  type  label  phi  dphi  kfac  power
    A    B     C     D     1      1  180     0     1      2
#include "tip4p.itp"
etc...
Make sure that the dihedral_restraints section comes immediately after  
the inclusion of the protein topology. If you have the protein  
topology directly in your .top file then just include  
dihedral_restraints after the protein listing but before any mention  
of things that are not that protein molecule.

Add this to your .mdp
;dihedral restraints
dihre               =  simple
dihre_fc            =  100     ; or whatever value you desire
dihre_tau           =  0.0
nstdihreout         =  50

For more information about what these options mean, please refer to  
the manual.