[gmx-developers] much faster dihedrals

[Alexei Podtelezhnikov]

Hi,

Please consider this much faster diheral angle algorithm, it uses atan2 wich is much safer and faster than acos. In fact, acos is slow because it is computed via atan2. It also gets rid of if-condition and a division.

real dih_angle(matrix box,
               rvec xi,rvec xj,rvec xk,rvec xl,
               rvec r_ji,rvec r_kj,rvec r_lk,rvec m,rvec n,
               int *t1,int *t2,int *t3)
{
  real ipr,phi;

  *t1 = pbc_rvec_sub(xj,xi,r_ji);                       /*  3           */
  *t2 = pbc_rvec_sub(xk,xj,r_kj);                       /*  3           */
  *t3 = pbc_rvec_sub(xl,xk,r_lk);                       /*  3           */

  oprod(r_ji,r_kj,m);                   /*  9           */
  oprod(r_kj,r_lk,n);                   /*  9           */
  phi=atan2(iprod(m,r_lk)*sqrt(iprod(r_kj,r_kj)),iprod(m,n));

  return phi;
}

Also, even a simple bond angle is a lot faster with atan2. Here is on way of doing it:

real bond_angle(matrix box,
                rvec xi,rvec xj,rvec xk,        /* in  */
                rvec r_ij,rvec r_kj,rvec m
                int *t1,int *t2)                /* out */
/* Return value is the angle between the bonds i-j and j-k */
{
  real th;

  *t1 = pbc_rvec_sub(xi,xj,r_ij);                       /*  3           */
  *t2 = pbc_rvec_sub(xk,xj,r_kj);                       /*  3           */

  oprod(r_ij,r_kj,m)

  th = atan2(sqrt(iprod(m,m)),iprod(r_ij,r_kj));

  return th;
}

Note that there is no slow divisions. We are just cutting though directly to atan2 in both cases. My point is that all acos'es should die if you want speed.

Thanks,
Alexei Podtelezhnikov, PhD