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