| 35 | 
  | 
 *                                                                       | 
| 36 | 
  | 
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).              | 
| 37 | 
  | 
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).           | 
| 38 | 
< | 
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).           | 
| 39 | 
< | 
 * [4]  Vardeman & Gezelter, in progress (2009).                         | 
| 38 | 
> | 
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).           | 
| 39 | 
> | 
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | 
> | 
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 | 
  | 
 */ | 
| 42 | 
  | 
 | 
| 43 | 
  | 
#include "restraints/ObjectRestraint.hpp" | 
| 47 | 
  | 
  void ObjectRestraint::calcForce(Vector3d struc) { | 
| 48 | 
  | 
 | 
| 49 | 
  | 
    pot_ = 0.0; | 
| 49 | 
– | 
     | 
| 50 | 
  | 
    if (restType_ & rtDisplacement) { | 
| 51 | 
  | 
      Vector3d del = struc - refPos_; | 
| 52 | 
< | 
      RealType r = del.length();   | 
| 52 | 
> | 
      RealType r = del.length(); | 
| 53 | 
  | 
      Vector3d frc = -kDisp_ * del; | 
| 54 | 
  | 
      RealType p = 0.5 * kDisp_ * del.lengthSquare(); | 
| 55 | 
< | 
      pot_ += p; | 
| 55 | 
> | 
 | 
| 56 | 
> | 
      pot_ = p; | 
| 57 | 
  | 
      force_ = frc * scaleFactor_; | 
| 58 | 
  | 
      restInfo_[rtDisplacement] = std::make_pair(r,p); | 
| 59 | 
  | 
    } | 
| 72 | 
  | 
 | 
| 73 | 
  | 
      Quat4d quat = temp.toQuaternion(); | 
| 74 | 
  | 
 | 
| 75 | 
< | 
      RealType twistAngle, swingAngle; | 
| 75 | 
> | 
      RealType twistAngle; | 
| 76 | 
  | 
      Vector3d swingAxis; | 
| 77 | 
< | 
      RealType tw, swingX, swingY; | 
| 77 | 
> | 
      RealType swingX, swingY; | 
| 78 | 
  | 
       | 
| 79 | 
< | 
      quat.getTwistSwingAxisAngle(twistAngle, swingAngle, swingAxis); | 
| 80 | 
< | 
      quat.toSwingTwist(swingX, swingY, tw); | 
| 80 | 
< | 
       | 
| 81 | 
< | 
      RealType dVdtwist, dVdswingX, dVdswingY; | 
| 82 | 
< | 
      RealType dTwist, dSwingX, dSwingY; | 
| 79 | 
> | 
      quat.toSwingTwist(swingX, swingY, twistAngle); | 
| 80 | 
> | 
 | 
| 81 | 
  | 
      RealType p; | 
| 82 | 
  | 
      Vector3d tTwist, tSwing; | 
| 83 | 
  | 
 | 
| 84 | 
  | 
      if (restType_ & rtTwist){ | 
| 85 | 
< | 
        dTwist = twistAngle - twist0_; | 
| 86 | 
< | 
        dVdtwist = kTwist_ * sin(dTwist); | 
| 85 | 
> | 
        RealType dTwist = twistAngle - twist0_; | 
| 86 | 
> | 
        RealType dVdtwist = kTwist_ * sin(dTwist); | 
| 87 | 
  | 
        p = kTwist_ * (1.0 - cos(dTwist) ); | 
| 88 | 
  | 
        pot_ += p; | 
| 89 | 
  | 
        tBody -= dVdtwist * V3Z; | 
| 91 | 
  | 
      } | 
| 92 | 
  | 
 | 
| 93 | 
  | 
      if (restType_ & rtSwingX){ | 
| 94 | 
< | 
        dSwingX = swingX - swingX0_; | 
| 95 | 
< | 
        dVdswingX = kSwingX_ * 0.5 * sin(2.0 * dSwingX); | 
| 94 | 
> | 
        RealType dSwingX = swingX - swingX0_; | 
| 95 | 
> | 
        RealType dVdswingX = kSwingX_ * 0.5 * sin(2.0 * dSwingX); | 
| 96 | 
  | 
        p = 0.25 * kSwingX_ * (1.0 - cos(2.0 * dSwingX)); | 
| 97 | 
  | 
        pot_ += p; | 
| 98 | 
  | 
        tBody -= dVdswingX * V3X; | 
| 100 | 
  | 
      } | 
| 101 | 
  | 
 | 
| 102 | 
  | 
      if (restType_ & rtSwingY){ | 
| 103 | 
< | 
        dSwingY = swingY - swingY0_; | 
| 104 | 
< | 
        dVdswingY = kSwingY_ * 0.5 * sin(2.0 * dSwingY); | 
| 103 | 
> | 
        RealType dSwingY = swingY - swingY0_; | 
| 104 | 
> | 
        RealType dVdswingY = kSwingY_ * 0.5 * sin(2.0 * dSwingY); | 
| 105 | 
  | 
        p = 0.25 * kSwingY_ * (1.0 - cos(2.0 * dSwingY)); | 
| 106 | 
  | 
        pot_ += p; | 
| 107 | 
  | 
        tBody -= dVdswingY * V3Y; |