Optimize and simplify the calculation of dihedral gradients

This commit uses the equations (27i) ~ (27l) from
https://onlinelibrary.wiley.com/doi/10.1002/(SICI)1096-987X(19960715)17:9%3C1132::AID-JCC5%3E3.0.CO;2-T
to simplify the calculation of dihedral gradients.

src/colvarcomp_angles.cpp

@@ -267,74 +267,22 @@ void colvar::dihedral::calc_value()
 
 void colvar::dihedral::calc_gradients()
 {
-  cvm::rvector A = cvm::rvector::outer(r12, r23);
-  cvm::real   rA = A.norm();
-  cvm::rvector B = cvm::rvector::outer(r23, r34);
-  cvm::real   rB = B.norm();
-  cvm::rvector C = cvm::rvector::outer(r23, A);
-  cvm::real   rC = C.norm();
-
-  cvm::real const cos_phi = (A*B)/(rA*rB);
-  cvm::real const sin_phi = (C*B)/(rC*rB);
-
-  cvm::rvector f1, f2, f3;
-
-  rB = 1.0/rB;
-  B *= rB;
-
-  if (cvm::fabs(sin_phi) > 0.1) {
-    rA = 1.0/rA;
-    A *= rA;
-    cvm::rvector const dcosdA = rA*(cos_phi*A-B);
-    cvm::rvector const dcosdB = rB*(cos_phi*B-A);
-    // rA = 1.0;
-
-    cvm::real const K = (1.0/sin_phi) * (180.0/PI);
-
-        f1 = K * cvm::rvector::outer(r23, dcosdA);
-        f3 = K * cvm::rvector::outer(dcosdB, r23);
-        f2 = K * (cvm::rvector::outer(dcosdA, r12)
-                   +  cvm::rvector::outer(r34, dcosdB));
-  }
-  else {
-    rC = 1.0/rC;
-    C *= rC;
-    cvm::rvector const dsindC = rC*(sin_phi*C-B);
-    cvm::rvector const dsindB = rB*(sin_phi*B-C);
-    // rC = 1.0;
-
-    cvm::real    const K = (-1.0/cos_phi) * (180.0/PI);
-
-    f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
-              - r23.x*r23.y*dsindC.y
-              - r23.x*r23.z*dsindC.z);
-    f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
-              - r23.y*r23.z*dsindC.z
-              - r23.y*r23.x*dsindC.x);
-    f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
-              - r23.z*r23.x*dsindC.x
-              - r23.z*r23.y*dsindC.y);
-
-    f3 = cvm::rvector::outer(dsindB, r23);
-    f3 *= K;
-
-    f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
-              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
-              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
-              +dsindB.z*r34.y - dsindB.y*r34.z);
-    f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
-              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
-              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
-              +dsindB.x*r34.z - dsindB.z*r34.x);
-    f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
-              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
-              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
-              +dsindB.y*r34.x - dsindB.x*r34.y);
-  }
+  // Eqs. (27i) ~ (27l) from https://doi.org/10.1002/(SICI)1096-987X(19960715)17:9<1132::AID-JCC5>3.0.CO;2-T.
+
+  const cvm::rvector A = cvm::rvector::outer(r12, r23);
+  const cvm::rvector B = cvm::rvector::outer(r23, r34);
+  const cvm::real   nG = r23.norm();
+  const cvm::real   A2 = A.norm2();
+  const cvm::real   B2 = B.norm2();
+
+  const cvm::real     K = 180.0/PI;
+  const cvm::rvector f1 = K * nG / A2 * A;
+  const cvm::rvector f2 = K * ((r12 * r23 / (A2 * nG)) * A + (r34 * r23 / (B2 * nG)) * B);
+  const cvm::rvector f3 = K * nG / B2 * B;
 
   group1->set_weighted_gradient(-f1);
-  group2->set_weighted_gradient(-f2 + f1);
-  group3->set_weighted_gradient(-f3 + f2);
+  group2->set_weighted_gradient( f2 + f1);
+  group3->set_weighted_gradient(-f3 - f2);
   group4->set_weighted_gradient(f3);
 }