[gmx-users] How to implement dihedral restraints

chris.neale at utoronto.ca chris.neale at utoronto.ca
Thu Dec 14 18:57:07 CET 2006

I am not sure what you are asking since I believe that I have already  
given you this answer on the gromacs mailing list.  

I certainly do not mind personal emails and I am happy to help, but  
there are three very good reasons to keep everything on the gromacs  
mailing list:
1. Answers may be read by people who disagree with portions of the  
answer -- such as Mark's suggestion that my system's stability using a  
dihedral restraint around 180deg was fortuitous and that you should  
redefine your dihedral to avoid that.  
2. I may not know the answer to your question but somebody else might.
3. The questions and answers are maintained for somebody else to  
search through years later.

In fact, I don't know how to do exactly what you are suggesting. Do  
you mean that at every integration step the dihedral feels a pull  
toward it's 180deg opposite position? Or perhaps you want to do this  
(using phiAtTimeZero instead of phi on the right side of the equation):

{ X [ A (1+cos (nT-phi)) ] } - { (1-X)  [ A (1+cos (nT-( phiAtTimeZero  
+ 180 ))) ] }

in place of the equation that you have written (below). In any event,  
if what you are asking is how to slowly change the umbrella then you  
could look into the afm pulling code and see if that is possible for  
dihedrals. Otherwise, you could do what I have suggested previously  
except that you should periodically stop and change the phi value of  
your dihedral restraint in the .top file and then grompp to get the  
updated tpr and then mdrun again. If you are asking how to use a cos  
function and thus be able to use 180deg... then you have me stumped. I  
also tried to do this but I think that it is not currently possible.
If you do use my previous suggestion then please note that there is no  
guarantees that anything you do around a dihedral of about 180 using  
the dihedral restraint will work (since it is not a cos function).  
Therefore you should redefine the value. For example, I am working  
similarly with a chi1 dihedral. For this I have the options of  
N-CA-CB-CG or C'-CA-CB-CG. If I want N-CA-CB-CG to be at 180deg then  
one could instead require that C'-CA-CB-CG be at 60 deg (check this  
actual value for yourself though to ensure that I haven't slipped a  
sign or rotated the wrong direction... this is all from my memory of  
what I set up a while ago).

Please don't hesitate to send another email to the list if you have  
some more questions on this topic.

Quoting Prasad Gajula <prasad.gajula at uos.de>:

> Dear Chris,
> Thanks a lot for your answer. I dont know if iam disturbing you. Iam  
>  very sorry for that.
> Hope you can help me to solve my prolem.
> i need to rotate a bond simultaneously from its present orientation   
> to opposite orientation.
> my idea is like this...
> the dihedral angle formula is
> [ A (1+cos (nT-phi)) ]  right.
> now, i want to do like this...
> { X [ A (1+cos (nT-phi)) ] } - { (1-X)  [ A (1+cos (nT-( phi + 180 ))) ] }
> I change the value of X from 1 down to 0, say 0.9, 0.8,.....
> that means i completely change the orientation of the bond.
> Now, my problem is how can i implement this special term in gromacs   
> code only for my particular dihedral atoms say B-C (bonded atoms   
> A-B-C-D) with simulation continues as usual normally for all other   
> bonds.
> Please do not mind for sending mail to you personally.
> looking forward for your answer.
> Thank you very much!
> Best regards
> Prasad
> Message: 9
> Date: Wed, 13 Dec 2006 15:54:28 -0500
> From: chris.neale at utoronto.ca
> Subject: [gmx-users] How to implement dihedral restraints
> To: gmx-users at gromacs.org
> Message-ID: <20061213155428.gcel5p932o2scs04 at webmail.utoronto.ca>
> Content-Type: text/plain; charset=ISO-8859-1; DelSp="Yes";
> format="flowed"
>> 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.

More information about the gromacs.org_gmx-users mailing list