1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2008 The Regents of the University of California 4 // 5 // This file is part of Qbox 6 // 7 // Qbox is distributed under the terms of the GNU General Public License 8 // as published by the Free Software Foundation, either version 2 of 9 // the License, or (at your option) any later version. 10 // See the file COPYING in the root directory of this distribution 11 // or <http://www.gnu.org/licenses/>. 12 // 13 //////////////////////////////////////////////////////////////////////////////// 14 // 15 // TorsionConstraint.h 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 19 #ifndef TORSIONCONSTRAINT_H 20 #define TORSIONCONSTRAINT_H 21 22 #include "Constraint.h" 23 #include "D3vector.h" 24 #include <cassert> 25 class AtomSet; 26 27 class TorsionConstraint : public Constraint 28 { 29 std::string name1_, name2_, name3_, name4_; 30 int ia1_, ia2_, ia3_, ia4_, is1_, is2_, is3_, is4_; 31 double m1_, m2_, m3_, m4_, m1_inv_, m2_inv_, m3_inv_, m4_inv_; 32 double angle_, velocity_, force_, weight_, tol_, sin_angle_, cos_angle_; 33 double sigma(D3vector a, D3vector b, 34 D3vector c, D3vector d) const; 35 void grad_sigma(const D3vector &r1, const D3vector &r2, 36 const D3vector &r3, const D3vector &r4, 37 D3vector &g1, D3vector &g2,D3vector &g3,D3vector &g4) const; 38 double torsion_angle(D3vector a, D3vector b, 39 D3vector c, D3vector d) const; 40 41 public: 42 TorsionConstraint(std::string name,std::string name1,std::string name2,std::string name3,std::string name4,double angle,double velocity,double tolerance)43 TorsionConstraint(std::string name, std::string name1, std::string name2, 44 std::string name3, std::string name4, 45 double angle, double velocity, double tolerance): 46 name1_(name1), name2_(name2), name3_(name3), name4_(name4), 47 velocity_(velocity), 48 tol_(tolerance), m1_(0.0), m2_(0.0), m3_(0.0), m4_(0.0) 49 { 50 set_value(angle); 51 name_ = name; 52 names_.resize(4); 53 names_[0] = name1_; 54 names_[1] = name2_; 55 names_[2] = name3_; 56 names_[3] = name4_; 57 force_ = 0.0; 58 weight_ = 1.0; 59 } ~TorsionConstraint(void)60 ~TorsionConstraint(void) {} 61 type(void)62 std::string type(void) const { return "torsion"; } value(void)63 double value(void) const { return angle_; } velocity(void)64 double velocity(void) const { return velocity_; } force(void)65 double force(void) const { return force_; } weight(void)66 double weight(void) const { return weight_; } tolerance(void)67 double tolerance(void) const { return tol_; } set_value(double value)68 void set_value(double value) 69 { 70 angle_ = value; 71 if ( angle_ < -180.0 ) angle_ = 180.0; 72 if ( angle_ > 180.0 ) angle_ = 180.0; 73 sin_angle_ = sin((M_PI/180.0)*angle_); 74 cos_angle_ = cos((M_PI/180.0)*angle_); 75 } set_velocity(double velocity)76 void set_velocity(double velocity) 77 { 78 velocity_ = velocity; 79 } 80 81 void setup(const AtomSet& atoms); dofs(void)82 int dofs(void) const { return 1; } 83 void update(double dt); 84 bool enforce_r(const std::vector<std::vector<double> > &r0, 85 std::vector<std::vector<double> > &rp) const; 86 bool enforce_v(const std::vector<std::vector<double> > &r0, 87 std::vector<std::vector<double> > &v0) const; 88 void compute_force(const std::vector<std::vector<double> > &r0, 89 const std::vector<std::vector<double> > &f); 90 std::ostream& print( std::ostream& os ); 91 92 }; 93 #endif 94