--- trunk/src/primitives/Inversion.cpp 2012/08/22 02:28:28 1782 +++ trunk/src/primitives/Inversion.cpp 2013/06/18 21:06:25 1890 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -50,7 +50,9 @@ namespace OpenMD { Inversion::Inversion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4, InversionType *it) : atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), - inversionType_(it) { } + inversionType_(it) { + inversionKey_ = inversionType_->getKey(); + } void Inversion::calcForce(RealType& angle, bool doParticlePot) { @@ -86,7 +88,22 @@ namespace OpenMD { if (cos_phi < -1.0) cos_phi = -1.0; RealType dVdcosPhi; - inversionType_->calcForce(cos_phi, potential_, dVdcosPhi); + switch (inversionKey_) { + case itCosAngle: + inversionType_->calcForce(cos_phi, potential_, dVdcosPhi); + break; + case itAngle: + RealType phi = acos(cos_phi); + RealType dVdPhi; + inversionType_->calcForce(phi, potential_, dVdPhi); + RealType sin_phi = sqrt(1.0 - cos_phi * cos_phi); + if (fabs(sin_phi) < 1.0E-6) { + sin_phi = 1.0E-6; + } + dVdcosPhi = dVdPhi / sin_phi; + break; + } + Vector3d f1 ; Vector3d f2 ; Vector3d f3 ;