Project

General

Profile

Bug #33

Inconsistent sign convention for calculation of angles between g_angle and angle restraints

Added by David Mobley about 14 years ago. Updated about 14 years ago.

Status:
Closed
Priority:
Low
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

This is a bug either in the documentation or code for g_angle and angle restraints; I've corresponded
about this on the mailing list but I want to make sure that it gets fixed. I discovered the problem
(described below) by calculating an angle between three atoms (A, B, C) using g_angle, then imposing a
restraint on this angle by setting up angle restraints in my toplogy file restraining the angle between
ABBC. This angle restraint forced the restrained angle AWAY from the preferred value.

I've now identified the exact problem. It seems that the routine angres computes angles using
essentially the cosine of the angle cos_angle(r_ij,r_kl), while g_angle (which is what I use to compute the
angles to restrain to) computes cos_angle(r_ij, r_kj). For the angres case, when j=k, this is cos_angle
(r_ij,r_jl). Renumber l as k and this is cos_angle(r_ij,r_jk).

So g_angle computes the cosine of the angle between r_ij and r_kj while the angle restraints (for the
case where the two middle atoms are the same) uses the cosine of the angle between r_ij and r_jk,
which turns out in this case to give different answers. It must be that cos_angle computes the cosine by
taking the dot product of the two vectors without taking the magnitude. That is, one can write cos
(theta)=|r_ij dot r_jk|/|r_ij||r_jk|, and if cos_angle were using this to compute the angle, there would be
no problem. But I think (based also on the documentation of the angle restraints, around p. 64) it's
using cos(theta)=r_ij dot r_jk/|r_ij||r_jk|, so the result is dependent on the sign of r_ij dot r_jk. And r_ij
dot r_jk has a different sign from r_ij dot r_kj.

In other words, g_angle uses r_ij dot r_kj, and angres uses r_ij dot r_jk, the two of which have different
signs.

It's not clear to me that this is exactly a bug, since it's documented how the angres routine computes
the angle. It's not, however, documented how g_angle does it, so at the very least that should be
documented. And it seems like a rather poor choice to use a different sign convention for the two when
what I'm doing is an obvious use of angle restraints. So my preferred fix would be to fix g_angle or
angres so they use the same sign convention. Alternatively, you could fix the documentation.

Once I knew this was the problem I simply changed the atom ordering in my topology file (ABBC ->
ABCB) to switch the sign of the vector between the last two atoms and now angres is consistent with
g_angle. But this, of course, is not documented anywhere.

Severity: Minor once you know what the issue is, but if people don't know that this is going on, angle
restraints will do absolutely the wrong thing for the case where one restrains an angle computed using
g_angle and involving three atoms.

Thanks,
David

History

#1 Updated by Berk Hess about 14 years ago

This is not a bug.

I would say that the definition of an angle ABC between particles
A, B and C is a 'universal' definition that does not need to be
stated in the manual.

The angle restraints are, as indeed stated in the manual,
restraints for the relative orientation of two vectors AB and CD.
I see no reason why one should expect that the angle for the restraint
between AB and BC should match the angle ABC.
To get this to correspond, one should use AB and CB (or BA and BC),
which, according to me, is not more and not less logical than AB and BC.

Berk.

#2 Updated by David Mobley about 14 years ago

(In reply to comment #1)

This is not a bug.

I would say that the definition of an angle ABC between particles
A, B and C is a 'universal' definition that does not need to be
stated in the manual.

The angle restraints are, as indeed stated in the manual,
restraints for the relative orientation of two vectors AB and CD.
I see no reason why one should expect that the angle for the restraint
between AB and BC should match the angle ABC.
To get this to correspond, one should use AB and CB (or BA and BC),
which, according to me, is not more and not less logical than AB and BC.

Berk.

It's obvious now, and I agree it's not a bug per se. But you might just make a
note in the documentation that points out that the form used for the angle
restraints means that AB-BC is not equivalent to the angle ABC. There have been
other people on the list who mention restraining an angle AB-BC, and it seems
like a logical mistake to make to assume that restraining an angle AB-BC is the
same as restraining the angle ABC.

In other words, I know now it's not a bug... But just for the sake of helping
people avoid making errors, you might want to point this out in the
documentation. A simple statement like "If you want to restrain the angle ABC,
note that this functional form means that you should use angle restraints on
AB-CB or BA-BC, not AB-BC," would do, if put in the angle restraints section.

It's inevitable that people will want to restrain the angle between three atoms
ABC, I think. It seems like it would be really helpful if you could make it
quite clear in the documentation how to do this.

Also available in: Atom PDF