65 RealType d21 = r21.length();
67 RealType d21inv = 1.0 / d21;
70 Vector3d r23 = ghostAtom->getA().transpose().getColumn(2);
71 RealType d23 = r23.
length();
73 RealType d23inv = 1.0 / d23;
75 RealType cosTheta =
dot(r21, r23) / (d21 * d23);
80 }
else if (cosTheta < -1.0) {
84 RealType theta = acos(cosTheta);
88 bendType_->calcForce(theta, potential_, dVdTheta);
90 RealType sinTheta = sqrt(1.0 - cosTheta * cosTheta);
92 if (fabs(sinTheta) < 1.0E-6) { sinTheta = 1.0E-6; }
94 RealType commonFactor1 = dVdTheta / sinTheta * d21inv;
95 RealType commonFactor2 = dVdTheta / sinTheta * d23inv;
97 Vector3d force1 = commonFactor1 * (r23 * d23inv - r21 * d21inv * cosTheta);
98 Vector3d force3 = commonFactor2 * (r21 * d21inv - r23 * d23inv * cosTheta);
102 atoms_[0]->addFrc(force1);
103 ghostAtom->addFrc(-force1);
105 ghostAtom->addTrq(
cross(r23, force3));
107 atoms_[0]->addParticlePot(potential_);
108 ghostAtom->addParticlePot(potential_);
111 angle = theta / Constants::PI * 180.0;