--- trunk/src/restraints/ObjectRestraint.cpp 2009/11/25 20:02:06 1390 +++ branches/development/src/restraints/ObjectRestraint.cpp 2013/02/20 15:39:39 1850 @@ -35,8 +35,9 @@ * * [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). - * [4] Vardeman & Gezelter, in progress (2009). + * [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). */ #include "restraints/ObjectRestraint.hpp" @@ -46,13 +47,12 @@ namespace OpenMD { void ObjectRestraint::calcForce(Vector3d struc) { pot_ = 0.0; - if (restType_ & rtDisplacement) { Vector3d del = struc - refPos_; - RealType r = del.length(); + RealType r = del.length(); Vector3d frc = -kDisp_ * del; RealType p = 0.5 * kDisp_ * del.lengthSquare(); - pot_ += p; + pot_ = p; force_ = frc * scaleFactor_; restInfo_[rtDisplacement] = std::make_pair(r,p); } @@ -71,13 +71,12 @@ namespace OpenMD { Quat4d quat = temp.toQuaternion(); - RealType twistAngle, swingAngle; + RealType twistAngle; Vector3d swingAxis; - RealType tw, swingX, swingY; + RealType swingX, swingY; - quat.getTwistSwingAxisAngle(twistAngle, swingAngle, swingAxis); - quat.toSwingTwist(tw, swingX, swingY); - + quat.toSwingTwist(swingX, swingY, twistAngle); + RealType dVdtwist, dVdswingX, dVdswingY; RealType dTwist, dSwingX, dSwingY; RealType p;