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

[Daniel L. Ensign]

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
(2) the settings for dihedral restraints should be documented,  
preferably in the wiki, which I'm not able to edit.

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