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

[Justin A. Lemkul]

Daniel L. Ensign wrote:
> Justin,
> 
> I appreciate your helpful reply. To summarize the additional points I 
> make below:
> (1) the code is right as far as calculating the dihedral restraint 
> potentials, but the manual has a bug, and

Yes, perhaps some additional information in the manual would be useful.

> (2) the settings for dihedral restraints should be documented, 
> preferably in the wiki, which I'm not able to edit.
> 

I've updated it:

http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints

-Justin

> Justin wrote:
>> 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;
> 
> Thanks for pointing that out. As I look closer, I also see
> 
>       if (dp > dphi)
>     ddp = dp-dphi;
>       else if (dp < -dphi)
>     ddp = dp+dphi;
>       else
>     ddp = 0;
> 
> followed by
> 
>     vtot += 0.5*kfac*ddp*ddp;
> 
> Your code snippet indicates that the code is correct (absolute value of 
> dp when calculating the angle modulus) and mine also indicates that the 
> code is correct (absolute value of dp when calculating the restraint 
> energy).
> 
> The manual, however, should have absolute values in each place that 
> phi-phi0 appears. How can the manual be corrected expending the least 
> effort?
> 
>>> 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.
> 
> snip
> 
>>> 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.
> 
> So then it seems there are two ways to set force constants -- by setting 
> dihre_fc in the mdp file and by setting kfac in the dihedral restraints 
> itp file. Good to know.
> 
>>> 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.
> 
> For now I'll just cross my fingers and put a 2 ... but it would be good 
> for all of these doohickeys to be documented.
> 
>>> 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.
> 
> sniip
> 
>> I might get around to updating this page with the above information, 
>> unless
>> there is more information to be considered, or I'm wrong :)
> 
> I don't appear to have the appropriate privileges for editing, but 
> that's okay, as it seems like a lot of responsibility.
> 
> Have fun,
> 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

========================================