src/primitives/Torsion.cpp

#include "primitives/Torsion.hpp"

namespace oopse {

Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
                 TorsionType *tt) :
    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }

void Torsion::calcForce() {
    Vector3d pos1 = atom1_->getPos();
    Vector3d pos2 = atom2_->getPos();
    Vector3d pos3 = atom3_->getPos();
    Vector3d pos4 = atom4_->getPos();

    Vector3d r12 = pos1 - pos2;
    Vector3d r23 = pos2 - pos3;
    Vector3d r34 = pos3 - pos4;

    //  Calculate the cross products and distances
    Vector3d A = cross(r12, r23);
    double rA = A.length();
    Vector3d B = cross(r23, r34);
    double rB = B.length();
    Vector3d C = cross(r23, A);
    double rC = C.length();

    A.normalize();
    B.normalize();
    C.normalize();
    
    //  Calculate the sin and cos
    double cos_phi = dot(A, B) ;
    double sin_phi = dot(C, B);

    double dVdPhi;
    torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);

    Vector3d f1;
    Vector3d f2;
    Vector3d f3;

    //  Next, we want to calculate the forces.  In order
    //  to do that, we first need to figure out whether the
    //  sin or cos form will be more stable.  For this,
    //  just look at the value of phi
    if (fabs(sin_phi) > 0.1) {
        //  use the sin version to avoid 1/cos terms

        Vector3d dcosdA = (cos_phi * A - B) /rA;
        Vector3d dcosdB = (cos_phi * B - A) /rB;

        double dVdcosPhi = dVdPhi / sin_phi;

        f1 = dVdcosPhi * cross(r23, dcosdA);
        f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
        f3 = dVdcosPhi * cross(r23, dcosdB);

    } else {
        //  This angle is closer to 0 or 180 than it is to
        //  90, so use the cos version to avoid 1/sin terms

        double dVdsinPhi = -dVdPhi /cos_phi;
        Vector3d dsindB = (sin_phi * B - C) /rB;
        Vector3d dsindC = (sin_phi * C - B) /rC;

        f1.x() = dVdsinPhi*((r23.y()*r23.y() + r23.z()*r23.z())*dsindC.x() - r23.x()*r23.y()*dsindC.y() - r23.x()*r23.z()*dsindC.z());

        f1.y() = dVdsinPhi*((r23.z()*r23.z() + r23.x()*r23.x())*dsindC.y() - r23.y()*r23.z()*dsindC.z() - r23.y()*r23.x()*dsindC.x());

        f1.z() = dVdsinPhi*((r23.x()*r23.x() + r23.y()*r23.y())*dsindC.z() - r23.z()*r23.x()*dsindC.x() - r23.z()*r23.y()*dsindC.y());

        f2.x() = dVdsinPhi*(-(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() = dVdsinPhi*(-(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() = dVdsinPhi*(-(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());

        f3 = dVdsinPhi * cross(dsindB, r23);

    }

    atom1_->addFrc(f1);
    atom2_->addFrc(f2 - f1);
    atom3_->addFrc(f3 - f2);
    atom4_->addFrc(-f3);
}

}
Revision:	1850
Committed:	Sat Dec 4 19:29:54 2004 UTC (19 years, 9 months ago) by tim
File size:	3215 byte(s)
Log Message:	rewind to old copy
#	User	Rev	Content
1	tim	1746	#include "primitives/Torsion.hpp"
2
3			namespace oopse {
4
5			Torsion::Torsion(Atom atom1, Atom atom2, Atom atom3, Atom atom4,
6			TorsionType *tt) :
7	tim	1847	atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
8	tim	1746
9			void Torsion::calcForce() {
10			Vector3d pos1 = atom1_->getPos();
11			Vector3d pos2 = atom2_->getPos();
12			Vector3d pos3 = atom3_->getPos();
13			Vector3d pos4 = atom4_->getPos();
14
15			Vector3d r12 = pos1 - pos2;
16			Vector3d r23 = pos2 - pos3;
17			Vector3d r34 = pos3 - pos4;
18
19			// Calculate the cross products and distances
20			Vector3d A = cross(r12, r23);
21			double rA = A.length();
22			Vector3d B = cross(r23, r34);
23			double rB = B.length();
24			Vector3d C = cross(r23, A);
25			double rC = C.length();
26
27			A.normalize();
28			B.normalize();
29			C.normalize();
30
31			// Calculate the sin and cos
32			double cos_phi = dot(A, B) ;
33			double sin_phi = dot(C, B);
34
35			double dVdPhi;
36	tim	1782	torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);
37	tim	1746
38			Vector3d f1;
39			Vector3d f2;
40			Vector3d f3;
41
42			// Next, we want to calculate the forces. In order
43			// to do that, we first need to figure out whether the
44			// sin or cos form will be more stable. For this,
45			// just look at the value of phi
46			if (fabs(sin_phi) > 0.1) {
47			// use the sin version to avoid 1/cos terms
48
49			Vector3d dcosdA = (cos_phi * A - B) /rA;
50			Vector3d dcosdB = (cos_phi * B - A) /rB;
51
52			double dVdcosPhi = dVdPhi / sin_phi;
53
54			f1 = dVdcosPhi * cross(r23, dcosdA);
55			f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
56			f3 = dVdcosPhi * cross(r23, dcosdB);
57
58			} else {
59			// This angle is closer to 0 or 180 than it is to
60			// 90, so use the cos version to avoid 1/sin terms
61
62			double dVdsinPhi = -dVdPhi /cos_phi;
63			Vector3d dsindB = (sin_phi * B - C) /rB;
64			Vector3d dsindC = (sin_phi * C - B) /rC;
65
66	tim	1782	f1.x() = dVdsinPhi((r23.y()r23.y() + r23.z()r23.z())dsindC.x() - r23.x()r23.y()dsindC.y() - r23.x()r23.z()dsindC.z());
67	tim	1746
68	tim	1782	f1.y() = dVdsinPhi((r23.z()r23.z() + r23.x()r23.x())dsindC.y() - r23.y()r23.z()dsindC.z() - r23.y()r23.x()dsindC.x());
69	tim	1746
70	tim	1782	f1.z() = dVdsinPhi((r23.x()r23.x() + r23.y()r23.y())dsindC.z() - r23.z()r23.x()dsindC.x() - r23.z()r23.y()dsindC.y());
71	tim	1746
72	tim	1782	f2.x() = dVdsinPhi(-(r23.y()r12.y() + r23.z()r12.z())dsindC.x() + (2.0r23.x()r12.y() - r12.x()r23.y())dsindC.y()
73			+ (2.0r23.x()r12.z() - r12.x()r23.z())dsindC.z() + dsindB.z()r34.y() - dsindB.y()r34.z());
74	tim	1746
75	tim	1782	f2.y() = dVdsinPhi(-(r23.z()r12.z() + r23.x()r12.x())dsindC.y() + (2.0r23.y()r12.z() - r12.y()r23.z())dsindC.z()
76			+ (2.0r23.y()r12.x() - r12.y()r23.x())dsindC.x() + dsindB.x()r34.z() - dsindB.z()r34.x());
77	tim	1746
78	tim	1782	f2.z() = dVdsinPhi(-(r23.x()r12.x() + r23.y()r12.y())dsindC.z() + (2.0r23.z()r12.x() - r12.z()r23.x())dsindC.x()
79			+(2.0r23.z()r12.y() - r12.z()r23.y())dsindC.y() + dsindB.y()r34.x() - dsindB.x()r34.y());
80	tim	1746
81			f3 = dVdsinPhi * cross(dsindB, r23);
82
83			}
84
85			atom1_->addFrc(f1);
86			atom2_->addFrc(f2 - f1);
87			atom3_->addFrc(f3 - f2);
88			atom4_->addFrc(-f3);
89			}
90
91			}