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 phi = -atan2(sin_phi, cos_phi);

    double dVdPhi;
    torsionType_->calcForce(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:	1746
Committed:	Wed Nov 17 06:37:56 2004 UTC (19 years, 9 months ago) by tim
File size:	3048 byte(s)
Log Message:	add PolynomialBondType, PolynomialBendType, PolynomialTorsionType, HarmonicBendType and CharmmTorsionType. Need to refine the design and add document for them
#	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 phi = -atan2(sin_phi, cos_phi);
39
40			double dVdPhi;
41			torsionType_->calcForce(phi, potential_, dVdPhi);
42
43			Vector3d f1;
44			Vector3d f2;
45			Vector3d f3;
46
47			// Next, we want to calculate the forces. In order
48			// to do that, we first need to figure out whether the
49			// sin or cos form will be more stable. For this,
50			// just look at the value of phi
51			if (fabs(sin_phi) > 0.1) {
52			// use the sin version to avoid 1/cos terms
53
54			Vector3d dcosdA = (cos_phi * A - B) /rA;
55			Vector3d dcosdB = (cos_phi * B - A) /rB;
56
57			double dVdcosPhi = dVdPhi / sin_phi;
58
59			f1 = dVdcosPhi * cross(r23, dcosdA);
60			f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
61			f3 = dVdcosPhi * cross(r23, dcosdB);
62
63			} else {
64			// This angle is closer to 0 or 180 than it is to
65			// 90, so use the cos version to avoid 1/sin terms
66
67			double dVdsinPhi = -dVdPhi /cos_phi;
68			Vector3d dsindB = (sin_phi * B - C) /rB;
69			Vector3d dsindC = (sin_phi * C - B) /rC;
70
71			f1.x = dVdsinPhi((r23.yr23.y + r23.zr23.z)dsindC.x - r23.xr23.ydsindC.y - r23.xr23.zdsindC.z);
72
73			f1.y = dVdsinPhi((r23.zr23.z + r23.xr23.x)dsindC.y - r23.yr23.zdsindC.z - r23.yr23.xdsindC.x);
74
75			f1.z = dVdsinPhi((r23.xr23.x + r23.yr23.y)dsindC.z - r23.zr23.xdsindC.x - r23.zr23.ydsindC.y);
76
77			f2.x = dVdsinPhi(-(r23.yr12.y + r23.zr12.z)dsindC.x + (2.0r23.xr12.y - r12.xr23.y)dsindC.y
78			+ (2.0r23.xr12.z - r12.xr23.z)dsindC.z + dsindB.zr34.y - dsindB.yr34.z);
79
80			f2.y = dVdsinPhi(-(r23.zr12.z + r23.xr12.x)dsindC.y + (2.0r23.yr12.z - r12.yr23.z)dsindC.z
81			+ (2.0r23.yr12.x - r12.yr23.x)dsindC.x + dsindB.xr34.z - dsindB.zr34.x);
82
83			f2.z = dVdsinPhi(-(r23.xr12.x + r23.yr12.y)dsindC.z + (2.0r23.zr12.x - r12.zr23.x)dsindC.x
84			+(2.0r23.zr12.y - r12.zr23.y)dsindC.y + dsindB.yr34.x - dsindB.xr34.y);
85
86			f3 = dVdsinPhi * cross(dsindB, r23);
87
88			}
89
90			atom1_->addFrc(f1);
91			atom2_->addFrc(f2 - f1);
92			atom3_->addFrc(f3 - f2);
93			atom4_->addFrc(-f3);
94			}
95
96			}