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) { }

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