move part of dihedral angle calculation from mol to geom

modules/geom/pymod/export_vecmat3_op.cc

@@ -49,6 +49,7 @@ void export_VecMat3_op()
   def("CompDivide",Vec3CompDivide);
   def("Distance",Vec3Distance);  
   def("Equal",Mat3Equal);
+  def("DihedralAngle", &DihedralAngle);
   def("Dot",Mat3Dot);
   def("Det",Mat3Det);
   def("Cross",Vec3Cross);

modules/geom/src/vecmat3_op.cc

@@ -171,4 +171,17 @@ Mat3 AxisRotation(const Vec3& axis, Real angle)
           xz-ca*xz-sa*y, yz-ca*yz+sa*x,zz+ca-ca*zz);
 }
 
+
+Real DihedralAngle(const Vec3& p1, const Vec3& p2, const Vec3& p3, 
+                   const Vec3&p4)
+{
+  Vec3 r1=p2-p1;
+  Vec3 r2=p3-p2;
+  Vec3 r3=p4-p3;
+  Vec3 r12cross = Cross(r1, r2);
+  Vec3 r23cross = Cross(r2, r3);
+  return atan2(Dot(r1*Length(r2), r23cross),
+               Dot(r12cross, r23cross));
+}
+
 } // ns

modules/geom/src/vecmat3_op.hh

@@ -163,6 +163,10 @@ Mat3 DLLEXPORT_OST_GEOM AxisRotation(const Vec3& axis, Real angle);
 /// The returned vector is of length 1
 Vec3 DLLEXPORT_OST_GEOM OrthogonalVector(const Vec3& axis);
 
+
+Real DLLEXPORT_OST_GEOM DihedralAngle(const Vec3& p1, const Vec3& p2, 
+                                      const Vec3& p3, const Vec3&p4);
+
 //! returns std::min of each component
 inline Vec3 Min(const Vec3& v1, const Vec3& v2)
 {

modules/mol/base/src/impl/dihedral.cc

@@ -46,13 +46,8 @@ Dihedral::Dihedral(const AtomImplList& atoms)
 
 
 Real Dihedral::GetAngleXCS() const {
-  Vec3 r1=atoms_[1]->GetPos()-atoms_[0]->GetPos();
-  Vec3 r2=atoms_[2]->GetPos()-atoms_[1]->GetPos();
-  Vec3 r3=atoms_[3]->GetPos()-atoms_[2]->GetPos();
-  Vec3 r12cross = Cross(r1, r2);
-  Vec3 r23cross = Cross(r2, r3);
-  return atan2(Dot(r1*Length(r2), r23cross),
-               Dot(r12cross, r23cross));
+  return geom::DihedralAngle(atoms_[0]->GetPos(), atoms_[1]->GetPos(), 
+                             atoms_[2]->GetPos(), atoms_[3]->GetPos());
 }
 
 Real Dihedral::GetAngleICS() const {