src/primitives/Bend.cpp

#include "primitives/Bend.hpp"

namespace oopse {

/**@todo still a lot left to improve*/
void Bend::calcForce() {
    Vector3d pos1 = atom1_->getPos();
    Vector3d pos2 = atom2_->getPos();
    Vector3d pos3 = atom3_->getPos();

    Vector3d r12 = pos1 - pos2;
    double d12 = r12.length();

    double d12inv = 1. 0 / d12;

    Vector3d r32 = pos3 - pos2;
    double d32 = r32.length();

    double d32inv = 1. 0 / d32;

    double cosTheta = (r12 * r32) / (d12 * d32);

    //check roundoff     
    if (cosTheta > 1.0) {
        cosTheta = 1.0;
    } else if (cosTheta < -1.0) {
        cosTheta = -1.0;
    }

    double theta = acos(cosTheta);

    double firstDerivative;

    bendType_->calcForce(theta, firstDerivative, potential_);

    double sinTheta = sqrt(1.0 - cosTheta * cosTheta);

    if (fabs(sinTheta) < 1.0E - 12) {
        sinTheta = 1.0E - 12;
    }

    double commonFactor1 = -firstDerivative / sinTheta * d12inv;
    double commonFactor2 = -firstDerivative / sinTheta * d32inv;

    Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);

    Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);

    //total force in current bend is zero
    Vector3d force2 = force1 + force3;
    force2 *= -1.0;

    atom1_->addFrc(force1);
    atom2_->addFrc(force2);
    atom3_->addFrc(force3);
}

    double k = value->k * scale;
    double theta0 = value->theta0;

    double diff = theta - theta0;

    double energy = k * diff * diff;

    double firstDerivative = 2.0 * k * diff;

} //end namespace oopse
Revision:	1746
Committed:	Wed Nov 17 06:37:56 2004 UTC (20 years, 11 months ago) by tim
File size:	1557 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/Bend.hpp"
2
3			namespace oopse {
4
5			/*@todo still a lot left to improve/
6			void Bend::calcForce() {
7			Vector3d pos1 = atom1_->getPos();
8			Vector3d pos2 = atom2_->getPos();
9			Vector3d pos3 = atom3_->getPos();
10
11			Vector3d r12 = pos1 - pos2;
12			double d12 = r12.length();
13
14			double d12inv = 1. 0 / d12;
15
16			Vector3d r32 = pos3 - pos2;
17			double d32 = r32.length();
18
19			double d32inv = 1. 0 / d32;
20
21			double cosTheta = (r12 * r32) / (d12 * d32);
22
23			//check roundoff
24			if (cosTheta > 1.0) {
25			cosTheta = 1.0;
26			} else if (cosTheta < -1.0) {
27			cosTheta = -1.0;
28			}
29
30			double theta = acos(cosTheta);
31
32			double firstDerivative;
33
34			bendType_->calcForce(theta, firstDerivative, potential_);
35
36			double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
37
38			if (fabs(sinTheta) < 1.0E - 12) {
39			sinTheta = 1.0E - 12;
40			}
41
42			double commonFactor1 = -firstDerivative / sinTheta * d12inv;
43			double commonFactor2 = -firstDerivative / sinTheta * d32inv;
44
45			Vector3d force1 = commonFactor1(r12(d12invcosTheta) - r32d32inv);
46
47			Vector3d force3 = commonFactor2(r32(d32invcosTheta) - r12d12inv);
48
49			//total force in current bend is zero
50			Vector3d force2 = force1 + force3;
51			force2 *= -1.0;
52
53			atom1_->addFrc(force1);
54			atom2_->addFrc(force2);
55			atom3_->addFrc(force3);
56			}
57
58			double k = value->k * scale;
59			double theta0 = value->theta0;
60
61			double diff = theta - theta0;
62
63			double energy = k * diff * diff;
64
65			double firstDerivative = 2.0 * k * diff;
66
67			} //end namespace oopse