src/primitives/GhostBend.cpp

#include "primitives/GhostBend.hpp"
#include "primitives/DirectionalAtom.hpp"
namespace oopse {

/**@todo still a lot left to improve*/
void GhostBend::calcForce() {
    DirectionalAtom* ghostAtom = static_cast<DirectionalAtom*>(atom2_);
    
    Vector3d pos1 = atom1_->getPos();
    Vector3d pos2 = ghostAtom->getPos();

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

    double d12inv = 1.0 / d12;

    Vector3d r32 = ghostAtom->getElectroFrame().getColumn(2);
    double d32 = r32.length();

    double d32inv = 1.0 / d32;

    double cosTheta = dot(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);
    atom1_->addFrc(force1);
    ghostAtom->addFrc(-force1);
    /**@todo test correctness */
    ghostAtom->addTrq(cross(r32, force3) );

}

} //end namespace oopse

Revision:	1832
Committed:	Thu Dec 2 16:04:19 2004 UTC (19 years, 7 months ago) by tim
File size:	1439 byte(s)
Log Message:	still have two linking problem
#	User	Rev	Content
1	tim	1832	#include "primitives/GhostBend.hpp"
2	tim	1782	#include "primitives/DirectionalAtom.hpp"
3	tim	1748	namespace oopse {
4
5			/*@todo still a lot left to improve/
6	tim	1832	void GhostBend::calcForce() {
7	tim	1782	DirectionalAtom* ghostAtom = static_cast<DirectionalAtom*>(atom2_);
8
9	tim	1748	Vector3d pos1 = atom1_->getPos();
10	tim	1782	Vector3d pos2 = ghostAtom->getPos();
11	tim	1748
12			Vector3d r12 = pos1 - pos2;
13			double d12 = r12.length();
14
15	tim	1782	double d12inv = 1.0 / d12;
16	tim	1748
17	tim	1816	Vector3d r32 = ghostAtom->getElectroFrame().getColumn(2);
18	tim	1748	double d32 = r32.length();
19
20	tim	1782	double d32inv = 1.0 / d32;
21	tim	1748
22	tim	1782	double cosTheta = dot(r12, r32) / (d12 * d32);
23	tim	1748
24			//check roundoff
25			if (cosTheta > 1.0) {
26			cosTheta = 1.0;
27			} else if (cosTheta < -1.0) {
28			cosTheta = -1.0;
29			}
30
31			double theta = acos(cosTheta);
32
33			double firstDerivative;
34
35			bendType_->calcForce(theta, firstDerivative, potential_);
36
37			double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
38
39	tim	1782	if (fabs(sinTheta) < 1.0E-12) {
40			sinTheta = 1.0E-12;
41	tim	1748	}
42
43			double commonFactor1 = -firstDerivative / sinTheta * d12inv;
44			double commonFactor2 = -firstDerivative / sinTheta * d32inv;
45
46			Vector3d force1 = commonFactor1(r12(d12invcosTheta) - r32d32inv);
47	tim	1782	Vector3d force3 = commonFactor2(r32(d32invcosTheta) - r12d12inv);
48	tim	1748	atom1_->addFrc(force1);
49	tim	1782	ghostAtom->addFrc(-force1);
50			/*@todo test correctness /
51			ghostAtom->addTrq(cross(r32, force3) );
52
53			}
54
55	tim	1748	} //end namespace oopse
56	tim	1782