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 r21 = pos1 - pos2;
    Vector3d r32 = pos2 - pos3;
    Vector3d r43 = pos3 - pos4;

    //  Calculate the cross products and distances
    Vector3d A = cross(r21, r32);
    double rA = A.length();
    Vector3d B = cross(r32, r43);
    double rB = B.length();
    Vector3d C = cross(r32, 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(r32, dcosdA);
        f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
        f3 = dVdcosPhi * cross(r32, 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*((r32.y()*r32.y() + r32.z()*r32.z())*dsindC.x() - r32.x()*r32.y()*dsindC.y() - r32.x()*r32.z()*dsindC.z());

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

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

        f2.x() = dVdsinPhi*(-(r32.y()*r21.y() + r32.z()*r21.z())*dsindC.x() + (2.0*r32.x()*r21.y() - r21.x()*r32.y())*dsindC.y()
        + (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());

        f2.y() = dVdsinPhi*(-(r32.z()*r21.z() + r32.x()*r21.x())*dsindC.y() + (2.0*r32.y()*r21.z() - r21.y()*r32.z())*dsindC.z()
        + (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());

        f2.z() = dVdsinPhi*(-(r32.x()*r21.x() + r32.y()*r21.y())*dsindC.z() + (2.0*r32.z()*r21.x() - r21.z()*r32.x())*dsindC.x()
        +(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());

        f3 = dVdsinPhi * cross(dsindB, r32);

    }

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

}
Revision:	1857
Committed:	Mon Dec 6 20:15:02 2004 UTC (19 years, 9 months ago) by tim
File size:	3215 byte(s)
Log Message:	short range interaction for butane is correct.Still something wrong with long range one
#	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	tim	1857	Vector3d r21 = pos1 - pos2;
16			Vector3d r32 = pos2 - pos3;
17			Vector3d r43 = pos3 - pos4;
18	tim	1746
19			// Calculate the cross products and distances
20	tim	1857	Vector3d A = cross(r21, r32);
21	tim	1746	double rA = A.length();
22	tim	1857	Vector3d B = cross(r32, r43);
23	tim	1746	double rB = B.length();
24	tim	1857	Vector3d C = cross(r32, A);
25	tim	1746	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	tim	1857	double dVdcosPhi = -dVdPhi / sin_phi;
53	tim	1746
54	tim	1857	f1 = dVdcosPhi * cross(r32, dcosdA);
55			f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
56			f3 = dVdcosPhi * cross(r32, dcosdB);
57	tim	1746
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	tim	1857	double dVdsinPhi = dVdPhi /cos_phi;
63	tim	1746	Vector3d dsindB = (sin_phi * B - C) /rB;
64			Vector3d dsindC = (sin_phi * C - B) /rC;
65
66	tim	1857	f1.x() = dVdsinPhi((r32.y()r32.y() + r32.z()r32.z())dsindC.x() - r32.x()r32.y()dsindC.y() - r32.x()r32.z()dsindC.z());
67	tim	1746
68	tim	1857	f1.y() = dVdsinPhi((r32.z()r32.z() + r32.x()r32.x())dsindC.y() - r32.y()r32.z()dsindC.z() - r32.y()r32.x()dsindC.x());
69	tim	1746
70	tim	1857	f1.z() = dVdsinPhi((r32.x()r32.x() + r32.y()r32.y())dsindC.z() - r32.z()r32.x()dsindC.x() - r32.z()r32.y()dsindC.y());
71	tim	1746
72	tim	1857	f2.x() = dVdsinPhi(-(r32.y()r21.y() + r32.z()r21.z())dsindC.x() + (2.0r32.x()r21.y() - r21.x()r32.y())dsindC.y()
73			+ (2.0r32.x()r21.z() - r21.x()r32.z())dsindC.z() + dsindB.z()r43.y() - dsindB.y()r43.z());
74	tim	1746
75	tim	1857	f2.y() = dVdsinPhi(-(r32.z()r21.z() + r32.x()r21.x())dsindC.y() + (2.0r32.y()r21.z() - r21.y()r32.z())dsindC.z()
76			+ (2.0r32.y()r21.x() - r21.y()r32.x())dsindC.x() + dsindB.x()r43.z() - dsindB.z()r43.x());
77	tim	1746
78	tim	1857	f2.z() = dVdsinPhi(-(r32.x()r21.x() + r32.y()r21.y())dsindC.z() + (2.0r32.z()r21.x() - r21.z()r32.x())dsindC.x()
79			+(2.0r32.z()r21.y() - r21.z()r32.y())dsindC.y() + dsindB.y()r43.x() - dsindB.x()r43.y());
80	tim	1746
81	tim	1857	f3 = dVdsinPhi * cross(dsindB, r32);
82	tim	1746
83			}
84
85			atom1_->addFrc(f1);
86			atom2_->addFrc(f2 - f1);
87			atom3_->addFrc(f3 - f2);
88			atom4_->addFrc(-f3);
89			}
90
91			}