--- trunk/src/primitives/Bend.cpp 2005/01/12 22:41:40 246 +++ trunk/src/primitives/Bend.cpp 2005/11/16 23:10:02 749 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -43,8 +43,8 @@ namespace oopse { namespace oopse { -/**@todo still a lot left to improve*/ -void Bend::calcForce() { + /**@todo still a lot left to improve*/ + void Bend::calcForce(double& angle) { Vector3d pos1 = atom1_->getPos(); Vector3d pos2 = atom2_->getPos(); Vector3d pos3 = atom3_->getPos(); @@ -63,9 +63,9 @@ void Bend::calcForce() { //check roundoff if (cosTheta > 1.0) { - cosTheta = 1.0; + cosTheta = 1.0; } else if (cosTheta < -1.0) { - cosTheta = -1.0; + cosTheta = -1.0; } double theta = acos(cosTheta); @@ -73,11 +73,13 @@ void Bend::calcForce() { double dVdTheta; bendType_->calcForce(theta, potential_, dVdTheta); + //std::cout << atom1_->getType() << "\t" << atom2_->getType() << "\t" << atom3_->getType() << "\t"; + //std::cout << "theta = " << theta/M_PI * 180.0 <<", potential = " << potential_ << std::endl; double sinTheta = sqrt(1.0 - cosTheta * cosTheta); if (fabs(sinTheta) < 1.0E-6) { - sinTheta = 1.0E-6; + sinTheta = 1.0E-6; } double commonFactor1 = dVdTheta / sinTheta * d21inv; @@ -93,6 +95,8 @@ void Bend::calcForce() { atom1_->addFrc(force1); atom2_->addFrc(force2); atom3_->addFrc(force3); -} + angle = theta /M_PI * 180.0; + } + } //end namespace oopse