Bug #32

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

Added by David Mobley almost 15 years ago. Updated almost 15 years ago.

Target version:
Affected version - extra info:
Affected version:


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;
ddphi = -kfac*ddp;

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.


#1 Updated by David van der Spoel almost 15 years ago

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

#2 Updated by David Mobley almost 15 years ago

(In reply to comment #1)

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

It seems to work fine in the sense that I don't see the problem I was seeing
anymore. But of course my system is hardly a test suite. I will let you know if
I notice any more problems but it looks good right now. (And I agree with your
changes; I implemented my own while I was waiting to hear back from you and it
was essentially the same).

Also available in: Atom PDF