TorsionConstraint.h 3.21 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17
//  TorsionConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
18
// $Id: TorsionConstraint.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
19 20 21 22 23 24 25 26 27 28 29

#ifndef TORSIONCONSTRAINT_H
#define TORSIONCONSTRAINT_H

#include "Constraint.h"
#include "D3vector.h"
#include <cassert>
class AtomSet;

class TorsionConstraint : public Constraint
{
30
  std::string name1_, name2_, name3_, name4_;
31 32
  int    ia1_, ia2_, ia3_, ia4_, is1_, is2_, is3_, is4_;
  double m1_, m2_, m3_, m4_, m1_inv_, m2_inv_, m3_inv_, m4_inv_;
Francois Gygi committed
33
  double angle_, velocity_, force_, weight_, tol_, sin_angle_, cos_angle_;
34 35 36 37 38 39 40
  double sigma(D3vector a, D3vector b,
               D3vector c, D3vector d) const;
  void grad_sigma(const D3vector &r1, const D3vector &r2,
                    const D3vector &r3, const D3vector &r4,
                    D3vector &g1, D3vector &g2,D3vector &g3,D3vector &g4) const;
  double torsion_angle(D3vector a, D3vector b,
                       D3vector c, D3vector d) const;
41

42
  public:
43 44

  TorsionConstraint(std::string name, std::string name1, std::string name2,
45
                    std::string name3, std::string name4,
Francois Gygi committed
46
                    double angle, double velocity, double tolerance):
47 48 49 50 51
  name1_(name1), name2_(name2), name3_(name3), name4_(name4),
  velocity_(velocity),
  tol_(tolerance), m1_(0.0), m2_(0.0), m3_(0.0), m4_(0.0)
  {
    set_value(angle);
Francois Gygi committed
52
    name_ = name;
53 54 55 56 57
    names_.resize(4);
    names_[0] = name1_;
    names_[1] = name2_;
    names_[2] = name3_;
    names_[3] = name4_;
Francois Gygi committed
58 59
    force_ = 0.0;
    weight_ = 1.0;
60
  }
61
  ~TorsionConstraint(void) {}
62

63
  std::string type(void) const { return "torsion"; }
64 65
  double value(void) const { return angle_; }
  double velocity(void) const { return velocity_; }
Francois Gygi committed
66 67
  double force(void) const { return force_; }
  double weight(void) const { return weight_; }
68 69 70 71 72 73 74 75 76 77 78 79 80
  double tolerance(void) const { return tol_; }
  void set_value(double value)
  {
    angle_ = value;
    if ( angle_ < -180.0 ) angle_ = 180.0;
    if ( angle_ >  180.0 ) angle_ = 180.0;
    sin_angle_ = sin((M_PI/180.0)*angle_);
    cos_angle_ = cos((M_PI/180.0)*angle_);
  }
  void set_velocity(double velocity)
  {
    velocity_ = velocity;
  }
81

82
  void setup(const AtomSet& atoms);
83
  int dofs(void) const { return 1; }
84
  void update(double dt);
85 86 87 88 89 90 91
  bool enforce_r(const std::vector<std::vector<double> > &r0,
                 std::vector<std::vector<double> > &rp) const;
  bool enforce_v(const std::vector<std::vector<double> > &r0,
                 std::vector<std::vector<double> > &v0) const;
  void compute_force(const std::vector<std::vector<double> > &r0,
                     const std::vector<std::vector<double> > &f);
  std::ostream& print( std::ostream& os );
92

93 94
};
#endif