## Bug #33

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

**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 over 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 over 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.