diff --git a/modules/geom/pymod/export_vecmat3_op.cc b/modules/geom/pymod/export_vecmat3_op.cc index cdafc821ce62797a555c1c2b98e29f661208e3ce..5f0939b1a2f2441e646ba4dd0e1593cbba1ddfea 100644 --- a/modules/geom/pymod/export_vecmat3_op.cc +++ b/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); diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc index 4c6a2732a70ecb9f02b0325300bd534aaa64e7fb..4194cac02880ea0d4fb7d3cd1fd528c247dacec9 100644 --- a/modules/geom/src/vecmat3_op.cc +++ b/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 diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh index e4ca5a65b5eb0b2a9d11717fb68fa7a2c2c14a9d..76bd719d282ffc4261fdcc15d3104c9f2e4ff079 100644 --- a/modules/geom/src/vecmat3_op.hh +++ b/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) { diff --git a/modules/mol/base/src/impl/dihedral.cc b/modules/mol/base/src/impl/dihedral.cc index ff484d4e694d70256e60525d5000897dfd75dcba..dab6b84d8542d4c8a34330c8f52e3788cb125028 100644 --- a/modules/mol/base/src/impl/dihedral.cc +++ b/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 {