[gmx-users] manual eq. 4.74-4.75 (dihedral restraints) head scratcher

[Berk Hess]

> Date: Sat, 17 Apr 2010 14:58:22 -0400
> From: jalemkul at vt.edu
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] manual eq. 4.74-4.75 (dihedral restraints)	head	scratcher
> 
> 
> 
> Daniel L. Ensign wrote:
> > Hello gmx-users, you rock and rollers,
> > 
> > Equations 4.74 and 4.75 in my copy of the manual have (please pardon my 
> > pseudo-LaTeX):
> > 
> > (4.74)
> > \phi' = (\phi - \phi_0) MOD 2\pi
> > 
> > (4.75)
> > V(\phi') = \frac{1}{2}k ( \phi' - \phi_0 - \Delta \phi )^2 if \phi' > 
> > \Delta \phi
> > or
> > V(\phi') = 0 if \phi' \leq \Delta \phi
> > 
> > but should there be absolute values around all of the \phi-\phi_0? Which 
> > way is it in the code -- with the absolute distance between \phi and 
> > \phi_0 or the directed distance?
> > 
> 
> It looks like absolute values are considered.  From src/gmxlib/dihres.c:
> 
>      /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
>       * force changes if we just apply a normal harmonic.
>       * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
>       * This means we will never have the periodicity problem, unless
>       * the dihedral is Pi away from phiO, which is very unlikely due to
>       * the potential.
>       */
>      dp = phi-phi0;
>      if (fabs(dp) > dphi) {
>        /* dp cannot be outside (-2*pi,2*pi) */
>        if (dp >= M_PI)
>      dp -= 2*M_PI;
>        else if(dp < -M_PI)
>      dp += 2*M_PI;
> 
> 
> > Also, as far as I can tell (and some mornings I definitely don't read 
> > too good) neither the manual nor 
> > http://wiki.gromacs.org/index.php/Dihedral_Restraints define the fields 
> > in [ dihedral_restraints ], although the latter does name them. There, I 
> > see
> > 
> > [ dihedral_restraints ]
> > ; ai   aj    ak    al  type  label  phi  dphi  kfac  power
> >     5    7     9    15     1      1  180     0     1      2
> > 
> > ai, aj, ak, al = atom numbers, obviously
> > type = ?, but I'm guessing there's only one type anyway
> 
> Probably so.
> 
> > label = what is this one?
> 
> Looks to be bookkeeping.  The code doesn't seem to use it other than to print 
> debug information, but I could be wrong since I haven't surfed around it very long.
> 

power and label have been copied from distance restraints, but not used at all.
I have already removed them for the next release.

Berk

> > phi means phi0 ?
> 
> Yes.
> 
> > dphi means Delta phi, I guess
> 
> Yes.
> 
> > kfac is the force constant, probably
> 
> Indirectly.  This term is equivalent to the fac value in distance restraints. 
> Since the force constant is specified in the .mdp file, different restraints 
> would otherwise have to be restrained with equivalent force constants.  The 
> value of kfac is multiplied by the value of dihre_fc in the .mdp file, so that 
> different restraints could have different force constants.
> 
> > power = what is this one? Does 2 give me harmonic constraints?
> > 
> 
> Not a clue on this one.  Also doesn't seem to be used in the code, but maybe 
> it's somewhere outside of dihres.c.
> 
> > I'd be happy to translate any answers given to the wiki, assuming I get 
> > answered and that I'm allowed to edit the wiki.
> > 
> 
> The wiki link you posted above doesn't actually exist any more, and actually 
> points to an empty page, even though the content is still online:
> 
> http://www.gromacs.org/index.php?title=Documentation/How-tos/Dihedral_Restraints
> 
> The old wiki doesn't actually exist, but has been merged into the Gromacs site, 
> so if you register as a user you can make contributions.  If you have troubles, 
> I might get around to updating this page with the above information, unless 
> there is more information to be considered, or I'm wrong :)
> 
> -Justin
> 
> > May the Force be with you (as long as it's calculated in hand-tuned 
> > assembly loops),
> > Dan
> > 
> 
> -- 
> ========================================
> 
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> ========================================
> -- 
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
 		 	   		  
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100418/bc646ad4/attachment.html>