## Bug #32

### bug in dihedral restraints code causes restraint force often to be in the wrong direction

Status:
Closed
Priority:
High
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized

Description

I ran a set of simulations with dihedral restraints of varying strengths, and find that as I make the
dihedral restraints stronger, the restraints force the restrained dihedrals away from the preferred
orientation. It looks like the mistake is in dihres.c in an if/else statement, such that the computed
ddphi value always has the same sign except in the case where dp is wrapped. The net result is that phi
is forced away from phi0 until it reaches a difference of pi from phi0, where the rapid change in the
force at dp=pi traps it.

Consider the following code snippet from dihres.c:

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;

``ddp = dp-dphi;      vtot += 0.5*kfac*ddp*ddp;      if (phi < phi0)        ddphi = kfac*ddp;      else        ddphi = -kfac*ddp;``
``do_dih_fup(ai,aj,ak,al,(real)(-ddphi),r_ij,r_kj,r_kl,m,n,                 f,fshift,pbc,g,x,t1,t2,t3);``

Consider the following cases (in degrees; here, dphi=0, so ddp=-dp):
phi phi0 dp ddp ddphi (related to force)
160 170 -10 10 k(10)
180 170 10 -10 -k(-10)=k(10)
-170 170 20 -20 k(-20)
170 -170 -20 20 -k(20)=k(-20)
0 180 180->-180 180 k(180)
0 179 179 -179 k(-179)

Notice that the last two cases look more or less OK, but the first four are bad: regardless of whether phi
or phi0 is smaller, the force ends up being in the same direction. I think the problem is the if/else
statement; it looks to me like it might work fine just to use ddphi=kfac*ddp regardless (unless I have
an overall sign error), since ddp already has a sign dependent on phi-phi0. Or maybe it would work if
you check whether ddp is negative rather than checking if phi<phi0.

This should be an easy fix, but it's a critical mistake, as it means that dihedral restraints ordinarily do
exactly the opposite of what is intended, except where phi and phi0 are different by at least 180
degrees (pi).

Please fix ASAP; it should be an easy fix and I need working dihedral restraints for my project.

### History

#### #1 Updated by David van der Spoelalmost 15 years ago

I've put a new version in CVS. Could yuo please test it?