[gmx-users] g_angle calculation

[Bogdan Costescu]

On Wed, Dec 21, 2011 at 18:44, SebastianWaltz
<sebastian.waltz at physik.uni-freiburg.de> wrote:
> Hi all,
>
> I try to reproduce the results from g_angle using my own c++ code. Using
> the formula
>
> θ=acos(r_ij * r_kj/|r_ij||r_kj|)
>
> with the acos from the cmath library, I obtain an angle which is always
> ~5% larger when the angel calculated by g_angle.
> Does g_angle use the same formula or does it calculate θ in a different way?

The angles are calculated through src/gmxlib/bondfree.c::bond_angle()
which does (after taking care of PBC):

  *costh=cos_angle(r_ij,r_kj);
  th=acos(*costh);

cos_angle() is defined in include/vec.h to do  r_ij *
r_kj/|r_ij||r_kj|; acos() comes from the math lib, it's not a GROMACS
function. So it's essentially what you seem to do. I guess
approximations and precision loss might play a role, f.e. cos_angle()
uses GROMACS' own gmx_invsqrt() - do you use math lib sqrt() to
calculate the vector norm ? You could try cut-and-pasting the few
pieces of code making up the GROMACS implementation of these functions
- they are well contained in include/vec.h.

Cheers,
Bogdan