--- trunk/src/primitives/GhostTorsion.cpp 2005/04/15 22:04:00 507 +++ trunk/src/primitives/GhostTorsion.cpp 2006/05/17 21:51:42 963 @@ -46,7 +46,7 @@ namespace oopse { GhostTorsion::GhostTorsion(Atom *atom1, Atom *atom2, DirectionalAtom* ghostAtom, TorsionType *tt) : Torsion(atom1, atom2, ghostAtom, ghostAtom, tt) {} - void GhostTorsion::calcForce() { + void GhostTorsion::calcForce(RealType& angle) { DirectionalAtom* ghostAtom = static_cast(atom3_); Vector3d pos1 = atom1_->getPos(); @@ -59,28 +59,25 @@ namespace oopse { // Calculate the cross products and distances Vector3d A = cross(r21, r32); - double rA = A.length(); + RealType rA = A.length(); Vector3d B = cross(r32, r43); - double rB = B.length(); + RealType rB = B.length(); Vector3d C = cross(r32, A); - double rC = C.length(); + RealType 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); + RealType cos_phi = dot(A, B) ; - double dVdPhi; - torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi); + RealType dVdcosPhi; + torsionType_->calcForce(cos_phi, potential_, dVdcosPhi); Vector3d dcosdA = (cos_phi * A - B) /rA; Vector3d dcosdB = (cos_phi * B - A) /rB; - double dVdcosPhi = -dVdPhi / sin_phi; - Vector3d f1 = dVdcosPhi * cross(r32, dcosdA); Vector3d f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA)); Vector3d f3 = dVdcosPhi * cross(dcosdB, r32); @@ -92,6 +89,8 @@ namespace oopse { f3.negate(); ghostAtom->addTrq(cross(r43, f3)); + + angle = acos(cos_phi) /M_PI * 180.0; } }