AngleConstraint.h 2.33 KB
Newer Older
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
//  AngleConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: AngleConstraint.h,v 1.3 2007-03-17 01:14:00 fgygi Exp $
7 8 9 10 11 12 13

#ifndef ANGLECONSTRAINT_H
#define ANGLECONSTRAINT_H

#include "Constraint.h"
#include "D3vector.h"
#include <cassert>
14

15 16 17 18
class AtomSet;

class AngleConstraint : public Constraint
{
19
  std::string name1_, name2_, name3_;
20 21
  int    ia1_, ia2_, ia3_, is1_, is2_, is3_;
  double m1_, m2_, m3_, m1_inv_, m2_inv_, m3_inv_;
Francois Gygi committed
22
  double angle_, velocity_, force_, weight_, tol_;
23 24 25 26 27 28 29
  double sigma(D3vector a, D3vector b, D3vector c) const;
  void grad_sigma(const D3vector &r1, const D3vector &r2, const D3vector &r3,
                  D3vector &g1, D3vector &g2,D3vector &g3) const;
  double bond_angle(D3vector a, D3vector b, D3vector c) const;
  
  public:
  
30 31 32
  AngleConstraint(std::string name, std::string name1, 
                  std::string name2, std::string name3,
                  double angle, double velocity, double tolerance):
33 34 35 36 37
  name1_(name1), name2_(name2), name3_(name3), 
  velocity_(velocity),
  tol_(tolerance), m1_(0.0), m2_(0.0), m3_(0.0)
  { 
    set_value(angle);
Francois Gygi committed
38
    name_ = name;
39 40 41 42
    names_.resize(3);
    names_[0] = name1_;
    names_[1] = name2_;
    names_[2] = name3_;
Francois Gygi committed
43 44
    force_ = 0.0;
    weight_ = 1.0;
45 46
  }
  
47
  std::string type(void) const { return "angle"; }
48 49
  double value(void) const { return angle_; }
  double velocity(void) const { return velocity_; }
Francois Gygi committed
50 51
  double force(void) const { return force_; }
  double weight(void) const { return weight_; }
52 53 54 55 56 57 58 59 60 61 62 63 64 65
  double tolerance(void) const { return tol_; }
  void set_value(double value)
  {
    angle_ = value;
    if ( angle_ <   0.0 ) angle_ =   0.0;
    if ( angle_ > 180.0 ) angle_ = 180.0;
  }
  void set_velocity(double velocity)
  {
    velocity_ = velocity;
  }
  
  void setup(const AtomSet& atoms);
  void update(double dt);
66 67 68 69 70 71 72
  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 );
73 74 75
  
};
#endif