Commit 8836760a by Francois Gygi

New constraint-related functions


git-svn-id: http://qboxcode.org/svn/qb/trunk@414 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 9b68a022
////////////////////////////////////////////////////////////////////////////////
//
// AngleConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AngleConstraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef ANGLECONSTRAINT_H
#define ANGLECONSTRAINT_H
#include "Constraint.h"
#include "D3vector.h"
#include <cassert>
class AtomSet;
class AngleConstraint : public Constraint
{
string name1_, name2_, name3_;
int ia1_, ia2_, ia3_, is1_, is2_, is3_;
double m1_, m2_, m3_, m1_inv_, m2_inv_, m3_inv_;
double angle_, velocity_, tol_;
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:
AngleConstraint(string name1, string name2, string name3,
double angle, double velocity, double tolerance):
name1_(name1), name2_(name2), name3_(name3),
velocity_(velocity),
tol_(tolerance), m1_(0.0), m2_(0.0), m3_(0.0)
{
set_value(angle);
names_.resize(3);
names_[0] = name1_;
names_[1] = name2_;
names_[2] = name3_;
}
string type(void) const { return "angle"; }
double value(void) const { return angle_; }
double velocity(void) const { return velocity_; }
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);
bool enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const;
bool enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const;
double projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const;
ostream& print( ostream& os );
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// Constraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Constraint.C,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#include "Constraint.h"
#include <iostream>
using namespace std;
ostream& operator << ( ostream &os, Constraint &c )
{
return c.print(os);
}
////////////////////////////////////////////////////////////////////////////////
//
// Constraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Constraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef CONSTRAINT_H
#define CONSTRAINT_H
#include <string>
#include <vector>
#include <cassert>
using namespace std;
class AtomSet;
class Constraint
{
protected:
vector<string> names_; // names of atoms involved in the constraint
public:
virtual string type(void) const = 0;
virtual double value(void) const = 0;
virtual double velocity(void) const = 0;
virtual double tolerance(void) const = 0;
virtual void set_value(double value) = 0;
virtual void set_velocity(double velocity) = 0;
virtual bool enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const = 0;
virtual bool enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const = 0;
virtual double projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const = 0;
virtual void update(double dt) = 0;
virtual void setup(const AtomSet& atoms) = 0;
virtual ostream& print(ostream &os) = 0;
string names(int i)
{
assert( i >= 0 && i < names_.size() );
return names_[i];
}
};
ostream& operator << ( ostream &os, Constraint &c );
#endif
////////////////////////////////////////////////////////////////////////////////
//
// ConstraintCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintCmd.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef CONSTRAINTCMD_H
#define CONSTRAINTCMD_H
#include "UserInterface.h"
#include "Sample.h"
class ConstraintCmd : public Cmd
{
public:
Sample *s;
ConstraintCmd(Sample *sample) : s(sample) {};
char *name(void) const { return "constraint"; }
char *help_msg(void) const
{
return
"\n constraint\n\n"
" syntax:\n\n"
" constraint distance name1 name2 distance [velocity]\n"
" The constraint command defines a new constraint and adds it to\n"
" the constraint list. Constraints are enforced at each MD step\n"
" if ions are allowed to move. The optional velocity parameter\n"
" must be given in atomic units.\n\n";
//" constraint multidistance weight name1 name2 \n"
//" [weight name1 name2] distance [velocity]\n"
}
int action(int argc, char **argv)
{
string subcmd(argv[1]);
if ( subcmd == "distance" || subcmd == "angle" || subcmd == "torsion" )
return s->constraints.set_constraint(s->atoms,argc,argv);
else if ( subcmd == "enforce" )
{
s->constraints.enforce(s->atoms);
// Note: should catch exception here if constraints could not be enforced
}
return 0;
}
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// ConstraintSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef CONSTRAINTSET_H
#define CONSTRAINTSET_H
#include <vector>
#include <string>
using namespace std;
class Atom;
class AtomSet;
class Constraint;
class Context;
class ConstraintSet
{
private:
const Context& ctxt_;
vector<Constraint *> constraint_list;
public:
ConstraintSet(const Context& ctxt) : ctxt_(ctxt) {}
bool set_constraint(AtomSet &atoms, int argc, char **argv);
bool del_constraint(int argc, char **argv);
void list_constraints(ostream &os);
int size(void) const { return constraint_list.size(); }
void enforce(AtomSet& atoms);
void enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp);
void enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0);
double projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x);
void update_constraints(double dt);
void setup(AtomSet& atoms);
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// DistanceConstraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceConstraint.C,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#include "DistanceConstraint.h"
#include "AtomSet.h"
#include "Atom.h"
#include "Species.h"
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
void DistanceConstraint::setup(const AtomSet& atoms)
{
// find position in tau array corresponding to atom name1
is1_ = atoms.is(name1_);
ia1_ = atoms.ia(name1_);
assert(is1_>=0);
assert(ia1_>=0);
m1_ = atoms.species_list[is1_]->mass() * 1822.89;
assert(m1_>0.0);
m1_inv_ = 1.0/m1_;
is2_ = atoms.is(name2_);
ia2_ = atoms.ia(name2_);
assert(is2_>=0);
assert(ia2_>=0);
m2_ = atoms.species_list[is2_]->mass() * 1822.89;
assert(m2_>0.0);
m2_inv_ = 1.0/m2_;
}
////////////////////////////////////////////////////////////////////////////////
void DistanceConstraint::update(double dt)
{
// cout << " DistanceConstraint::update" << endl;
if ( distance_ + velocity_ * dt > 0.0 )
distance_ += velocity_ * dt;
}
////////////////////////////////////////////////////////////////////////////////
bool DistanceConstraint::enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const
{
const double* pr1 = &r0[is1_][3*ia1_];
const double* pr2 = &r0[is2_][3*ia2_];
D3vector r1(pr1);
D3vector r2(pr2);
double* pr1p = &rp[is1_][3*ia1_];
double* pr2p = &rp[is2_][3*ia2_];
D3vector r1p(pr1p);
D3vector r2p(pr2p);
// compute gradient at r
D3vector r12(r1-r2);
D3vector g1,g2;
g1 = 2.0 * r12;
g2 = -g1;
double ng = g1*g1 + g2*g2;
assert(ng>=0.0);
// compute gradient at rp
D3vector r12p(r1p-r2p);
D3vector g1p,g2p;
g1p = 2.0 * r12p;
g2p = -g1p;
const double ngp = g1p*g1p + g2p*g2p;
assert(ngp>=0.0);
// test alignment of the gradient at r and at rp
// compute scalar product of normalized gradients
const double sp = ( g1*g1p + g2*g2p ) / sqrt( ng * ngp );
if ( fabs(sp) < 0.5*sqrt(3.0) )
{
g1 = g1p;
g2 = g2p;
#if DEBUG_CONSTRAINTS
cout << " <!-- g and gp nearly orthogonal, use gp only -->" << endl;
#endif
}
const double sigma = r12p*r12p - distance_*distance_;
#if DEBUG_CONSTRAINTS
cout << " <!-- ";
cout << " DistanceConstraint::enforce_r: " << name1_ << " " << name2_ << endl;
cout << " DistanceConstraint::enforce_r: r1 = " << r1 << endl;
cout << " DistanceConstraint::enforce_r: r2 = " << r2 << endl;
cout << " DistanceConstraint::enforce_r: tol = " << tol_ << endl;
cout << " DistanceConstraint::enforce_r: err = " << sigma << endl;
cout << " DistanceConstraint::enforce_r: g1 = " << g1 << endl;
cout << " DistanceConstraint::enforce_r: g2 = " << g2 << endl;
cout << " DistanceConstraint::enforce_r: g1p = " << g1p << endl;
cout << " DistanceConstraint::enforce_r: g2p = " << g2p << endl;
cout << " -->" << endl;
#endif
if ( fabs(sigma) < tol_ ) return true;
// make one shake iteration
const double den = m1_inv_ * g1 * g1p + m2_inv_ * g2 * g2p;
const double lambda = -sigma / den;
pr1p[0] += m1_inv_ * lambda * g1.x;
pr1p[1] += m1_inv_ * lambda * g1.y;
pr1p[2] += m1_inv_ * lambda * g1.z;
pr2p[0] += m2_inv_ * lambda * g2.x;
pr2p[1] += m2_inv_ * lambda * g2.y;
pr2p[2] += m2_inv_ * lambda * g2.z;
return false;
}
////////////////////////////////////////////////////////////////////////////////
bool DistanceConstraint::enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const
{
const double* pr1 = &r0[is1_][3*ia1_];
const double* pr2 = &r0[is2_][3*ia2_];
D3vector r1(pr1);
D3vector r2(pr2);
double* pv1 = &v0[is1_][3*ia1_];
double* pv2 = &v0[is2_][3*ia2_];
D3vector v1(pv1);
D3vector v2(pv2);
// compute gradient at r
D3vector r12(r1-r2);
D3vector g1,g2;
g1 = 2.0 * r12;
g2 = -g1;
const double norm2 = g1*g1 + g2*g2;
assert(norm2>=0.0);
// if the gradient is too small, do not attempt correction
if ( norm2 < 1.e-6 ) return true;
const double proj = v1 * g1 + v2 * g2;
const double err = fabs(proj)/sqrt(norm2);
#if DEBUG_CONSTRAINTS
cout << " <!-- ";
cout << " DistanceConstraint::enforce_v: " << name1_ << " " << name2_ << endl;
cout << " DistanceConstraint::enforce_v: r1 = " << r1 << endl;
cout << " DistanceConstraint::enforce_v: r2 = " << r2 << endl;
cout << " DistanceConstraint::enforce_v: v1 = "
<< v1[0] << " " << v1[1] << " " << v1[2] << endl;
cout << " DistanceConstraint::enforce_v: v2 = "
<< v2[0] << " " << v2[1] << " " << v2[2] << endl;
cout << " DistanceConstraint::enforce_v: tol = " << tol_ << endl;
cout << " DistanceConstraint::enforce_v: err = " << err
<< endl;
cout << " -->" << endl;
#endif
if ( err < tol_ ) return true;
// make one shake iteration
const double den = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2;
assert(den>0.0);
const double eta = -proj / den;
pv1[0] += m1_inv_ * eta * g1.x;
pv1[1] += m1_inv_ * eta * g1.y;
pv1[2] += m1_inv_ * eta * g1.z;
pv2[0] += m2_inv_ * eta * g2.x;
pv2[1] += m2_inv_ * eta * g2.y;
pv2[2] += m2_inv_ * eta * g2.z;
return false;
}
////////////////////////////////////////////////////////////////////////////////
double DistanceConstraint::projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const
{
const double* pr1 = &r0[is1_][3*ia1_];
const double* pr2 = &r0[is2_][3*ia2_];
D3vector r1(pr1);
D3vector r2(pr2);
const double* px1 = &x[is1_][3*ia1_];
const double* px2 = &x[is2_][3*ia2_];
D3vector x1(px1);
D3vector x2(px2);
// compute gradient at r
D3vector r12(r1-r2);
D3vector g1,g2;
g1 = 2.0 * r12;
g2 = -g1;
const double norm2 = g1*g1;
assert(norm2>=0.0);
if ( norm2 == 0.0 )
return 0.0;
const double norm = sqrt(norm2);
D3vector e1(g1/norm);
D3vector e2(-e1);
const double proj1 = x1*e1;
const double proj2 = x2*e2;
return 0.5*(proj1+proj2);
}
////////////////////////////////////////////////////////////////////////////////
ostream& DistanceConstraint::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <!-- constraint ";
os << type() << " ";
os << setw(4) << name1_ << " ";
os << setw(4) << name2_ << " ";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << setw(10) << setprecision(6) << distance_ << " "
<< setw(10) << setprecision(6) << velocity_ << " -->";
return os;
}
////////////////////////////////////////////////////////////////////////////////
//
// DistanceConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceConstraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef DISTANCECONSTRAINT_H
#define DISTANCECONSTRAINT_H
#include "Constraint.h"
#include <cassert>
#include <cmath> // fabs
class AtomSet;
class DistanceConstraint : public Constraint
{
string name1_, name2_;
int ia1_, ia2_, is1_, is2_;
double m1_, m2_, m1_inv_, m2_inv_;
double distance_, velocity_, tol_;
public:
DistanceConstraint(string name1, string name2,
double distance, double velocity, double tolerance):
name1_(name1), name2_(name2), distance_(distance), velocity_(velocity),
tol_(tolerance), m1_(0.0), m2_(0.0)
{
names_.resize(2);
names_[0] = name1_;
names_[1] = name2_;
}
string type(void) const { return "distance"; }
double value(void) const { return distance_; }
double velocity(void) const { return velocity_; }
double tolerance(void) const { return tol_; }
void set_value(double value)
{
distance_ = fabs(value);
}
void set_velocity(double velocity)
{
velocity_ = velocity;
}
void setup(const AtomSet& atoms);
void update(double dt);
bool enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const;
bool enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const;
double projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const;
ostream& print( ostream& os );
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// ListConstraintsCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ListConstraintsCmd.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef LISTCONSTRAINTSCMD_H
#define LISTCONSTRAINTSCMD_H
#include <iostream>
using namespace std;
#include "UserInterface.h"
#include "Sample.h"
class ListConstraintsCmd : public Cmd
{
public:
Sample *s;
ListConstraintsCmd(Sample *sample) : s(sample) {};
char *name(void) const { return "list_constraints"; }
char *help_msg(void) const
{
return
"\n list_constraints\n\n"
" syntax: list_constraints\n\n"
" The list_constraints command prints the list"
" of active constraints.\n\n";
}
int action(int argc, char **argv)
{
if ( s->ctxt_.onpe0() ) s->constraints.list_constraints(cout);
return 0;
}
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// TorsionConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: TorsionConstraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
#ifndef TORSIONCONSTRAINT_H
#define TORSIONCONSTRAINT_H
#include "Constraint.h"
#include "D3vector.h"
#include <cassert>
class AtomSet;
class TorsionConstraint : public Constraint
{
string name1_, name2_, name3_, name4_;
int ia1_, ia2_, ia3_, ia4_, is1_, is2_, is3_, is4_;
double m1_, m2_, m3_, m4_, m1_inv_, m2_inv_, m3_inv_, m4_inv_;
double angle_, velocity_, tol_, sin_angle_, cos_angle_;
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;
public:
TorsionConstraint(string name1, string name2, string name3, string name4,
double angle, double velocity, double tolerance):
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);
names_.resize(4);
names_[0] = name1_;
names_[1] = name2_;
names_[2] = name3_;
names_[3] = name4_;
}
string type(void) const { return "torsion"; }
double value(void) const { return angle_; }
double velocity(void) const { return velocity_; }
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;
}
void setup(const AtomSet& atoms);
void update(double dt);
bool enforce_r(const vector<vector<double> > &r0,
vector<vector<double> > &rp) const;
bool enforce_v(const vector<vector<double> > &r0,
vector<vector<double> > &v0) const;
double projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const;
ostream& print( ostream& os );
};
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment