Commit b305d27b by Francois Gygi

Redesign of constraints.


git-svn-id: http://qboxcode.org/svn/qb/trunk@422 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 3dcc6e96
......@@ -3,7 +3,7 @@
// AngleConstraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AngleConstraint.C,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: AngleConstraint.C,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#include "AngleConstraint.h"
#include "AtomSet.h"
......@@ -249,15 +249,15 @@ vector<vector<double> > &v0) const
}
////////////////////////////////////////////////////////////////////////////////
double AngleConstraint::projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const
void AngleConstraint::compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f)
{
const double* pr1 = &r0[is1_][3*ia1_];
const double* pr2 = &r0[is2_][3*ia2_];
const double* pr3 = &r0[is3_][3*ia3_];
const double* px1 = &x[is1_][3*ia1_];
const double* px2 = &x[is2_][3*ia2_];
const double* px3 = &x[is3_][3*ia3_];
const double* pf1 = &f[is1_][3*ia1_];
const double* pf2 = &f[is2_][3*ia2_];
const double* pf3 = &f[is3_][3*ia3_];
D3vector r1(pr1);
D3vector r2(pr2);
......@@ -266,16 +266,40 @@ double AngleConstraint::projection(const vector<vector<double> > &r0,
grad_sigma(r1,r2,r3,g1,g2,g3);
D3vector x1(px1);
D3vector x2(px2);
D3vector x3(px3);
D3vector f1(pf1);
D3vector f2(pf2);
D3vector f3(pf3);
const double norm2 = g1*g1 + g2*g2 + g3*g3;
if ( norm2 == 0.0 )
return 0.0;
{
force_ = 0.0;
return;
}
const double proj = x1*g1 + x2*g2 + x3*g3;
return proj/sqrt(norm2);
const double proj = f1*g1 + f2*g2 + f3*g3;
force_ = -proj/norm2;
// compute weight
const double z = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2 +
m3_inv_ * g3 * g3;
assert(z > 0.0);
weight_ = 1.0 / sqrt(z);
#if DEBUG_CONSTRAINTS
// check value of z
const double r12s = norm(r1-r2);
const double r32s = norm(r3-r2);
const double fac = 180.0/M_PI;
const double cos_theta = normalized(r1-r2)*normalized(r3-r2);
const double zcheck = fac*fac *
( m1_inv_ / r12s +
m2_inv_ * ( (r12s+r32s-2*sqrt(r12s*r32s)*cos_theta) /(r12s*r32s) ) +
m3_inv_ / r32s
);
cout << " <!-- AngleConstraint: z=" << z << " zcheck=" << zcheck << " -->"
<< endl;
#endif
}
////////////////////////////////////////////////////////////////////////////////
......@@ -371,14 +395,15 @@ double AngleConstraint::bond_angle(D3vector a, D3vector b, D3vector c) const
ostream& AngleConstraint::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <!-- constraint ";
os << type() << " ";
os << setw(4) << name1_ << " ";
os << setw(4) << name2_ << " ";
os << setw(4) << name3_ << " ";
os << " <constraint name=\"" << name();
os << "\" type=\"" << type();
os << "\" atoms=\"" << name1_ << " ";
os << name2_ << " " << name3_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << setw(10) << setprecision(6) << angle_ << " "
<< setw(10) << setprecision(6) << velocity_ << " -->";
os << " value=\"" << setprecision(6) << angle_;
os << "\" velocity=\"" << setprecision(6) << velocity_ << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
return os;
}
......@@ -3,7 +3,7 @@
// AngleConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AngleConstraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: AngleConstraint.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef ANGLECONSTRAINT_H
#define ANGLECONSTRAINT_H
......@@ -18,7 +18,7 @@ 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 angle_, velocity_, force_, weight_, 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;
......@@ -26,22 +26,27 @@ class AngleConstraint : public Constraint
public:
AngleConstraint(string name1, string name2, string name3,
AngleConstraint(string name, 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);
name_ = name;
names_.resize(3);
names_[0] = name1_;
names_[1] = name2_;
names_[2] = name3_;
force_ = 0.0;
weight_ = 1.0;
}
string type(void) const { return "angle"; }
double value(void) const { return angle_; }
double velocity(void) const { return velocity_; }
double force(void) const { return force_; }
double weight(void) const { return weight_; }
double tolerance(void) const { return tol_; }
void set_value(double value)
{
......@@ -60,8 +65,8 @@ class AngleConstraint : public Constraint
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;
void compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f);
ostream& print( ostream& os );
};
......
......@@ -3,7 +3,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.24 2005-06-27 22:31:39 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.25 2005-09-16 23:08:11 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -255,20 +255,14 @@ void BOSampleStepper::step(int niter)
if ( compute_forces )
{
#if 0
if ( s_.constraints.size() > 0 )
{
const double projected_force =
s_.constraints.projection(ionic_stepper->r0(), fion);
s_.constraints.compute_forces(ionic_stepper->r0(), fion);
if ( onpe0 )
{
cout << " <constraint_projected_force> "
<< projected_force
<< " </constraint_projected_force>"
<< endl;
s_.constraints.list_constraints(cout);
}
}
#endif
// move atoms to new position: r0 <- r0 + v0*dt + dt2/m * fion
ionic_stepper->compute_r(energy,fion);
ef_.atoms_moved();
......
......@@ -3,7 +3,7 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.10 2005-06-27 22:30:28 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.11 2005-09-16 23:08:11 fgygi Exp $
#include "CPSampleStepper.h"
#include "SlaterDet.h"
......@@ -190,17 +190,13 @@ void CPSampleStepper::step(int niter)
}
}
}
#if 0
#if 1
if ( s_.constraints.size() > 0 )
{
const double projected_force =
s_.constraints.projection(mdionic_stepper->r0(), fion);
s_.constraints.compute_forces(mdionic_stepper->r0(), fion);
if ( onpe0 )
{
cout << " <constraint_projected_force> "
<< projected_force
<< " </constraint_projected_force>"
<< endl;
s_.constraints.list_constraints(cout);
}
}
#endif
......
......@@ -3,7 +3,7 @@
// Constraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Constraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: Constraint.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef CONSTRAINT_H
#define CONSTRAINT_H
......@@ -19,6 +19,7 @@ class Constraint
{
protected:
string name_; // constraint name
vector<string> names_; // names of atoms involved in the constraint
public:
......@@ -26,6 +27,8 @@ class Constraint
virtual string type(void) const = 0;
virtual double value(void) const = 0;
virtual double velocity(void) const = 0;
virtual double force(void) const = 0;
virtual double weight(void) const = 0;
virtual double tolerance(void) const = 0;
virtual void set_value(double value) = 0;
virtual void set_velocity(double velocity) = 0;
......@@ -33,12 +36,13 @@ class Constraint
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 compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f) = 0;
virtual void update(double dt) = 0;
virtual void setup(const AtomSet& atoms) = 0;
virtual ostream& print(ostream &os) = 0;
string names(int i)
string name(void) const { return name_; }
string names(int i) const
{
assert( i >= 0 && i < names_.size() );
return names_[i];
......
......@@ -3,7 +3,7 @@
// ConstraintCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintCmd.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: ConstraintCmd.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef CONSTRAINTCMD_H
#define CONSTRAINTCMD_H
......@@ -25,25 +25,57 @@ class ConstraintCmd : public Cmd
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 define distance name atom1 atom2 distance [velocity]\n"
" constraint define angle name atom1 atom2 atom3 angle [velocity]\n"
" constraint define torsion name atom1 atom2 atom3 atom4 angle [velocity]\n"
" constraint set name value [velocity]\n"
" constraint delete name\n"
" constraint list\n"
" constraint enforce\n\n"
" Constraints are enforced at each MD step if ions are allowed to move.\n"
" Velocity parameters are optional.\n\n";
//" constraint multidistance weight name1 name2 \n"
//" [weight name1 name2] distance [velocity]\n"
}
int action(int argc, char **argv)
{
const bool onpe0 = s->ctxt_.onpe0();
if ( argc < 2 )
{
if ( onpe0 )
cout << help_msg();
return 1;
}
string subcmd(argv[1]);
if ( subcmd == "distance" || subcmd == "angle" || subcmd == "torsion" )
return s->constraints.set_constraint(s->atoms,argc,argv);
if ( subcmd == "define" )
{
return s->constraints.define_constraint(s->atoms,argc,argv);
}
else if ( subcmd == "set" )
{
return s->constraints.set_constraint(argc,argv);
}
else if ( subcmd == "delete" )
{
return s->constraints.delete_constraint(argc,argv);
}
else if ( subcmd == "enforce" )
{
s->constraints.enforce(s->atoms);
// Note: should catch exception here if constraints could not be enforced
}
else if ( subcmd == "list" )
{
if ( onpe0 )
s->constraints.list_constraints(cout);
}
else
{
if ( onpe0 )
cout << help_msg();
}
return 0;
}
};
......
......@@ -3,7 +3,7 @@
// ConstraintSet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: ConstraintSet.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef CONSTRAINTSET_H
#define CONSTRAINTSET_H
......@@ -27,8 +27,9 @@ class ConstraintSet
public:
ConstraintSet(const Context& ctxt) : ctxt_(ctxt) {}
bool set_constraint(AtomSet &atoms, int argc, char **argv);
bool del_constraint(int argc, char **argv);
bool define_constraint(AtomSet &atoms, int argc, char **argv);
bool set_constraint(int argc, char **argv);
bool delete_constraint(int argc, char **argv);
void list_constraints(ostream &os);
int size(void) const { return constraint_list.size(); }
void enforce(AtomSet& atoms);
......@@ -36,8 +37,8 @@ class ConstraintSet
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 compute_forces(const vector<vector<double> > &r0,
const vector<vector<double> > &f);
void update_constraints(double dt);
void setup(AtomSet& atoms);
};
......
......@@ -3,7 +3,7 @@
// DistanceCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $
// $Id: DistanceCmd.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef DISTANCECMD_H
#define DISTANCECMD_H
......@@ -66,7 +66,7 @@ class DistanceCmd : public Cmd
cout.setf(ios::fixed,ios::floatfield);
cout << " <!-- distance " << name1 << "-" << name2 << ": "
<< setprecision(3)
<< d << " (a.u.) / " << 0.529177*d << " (Ang)" << endl;
<< d << " (a.u.) / " << 0.529177*d << " (Ang) -->" << endl;
}
return 0;
}
......
......@@ -3,7 +3,7 @@
// DistanceConstraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceConstraint.C,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: DistanceConstraint.C,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#include "DistanceConstraint.h"
#include "AtomSet.h"
......@@ -177,17 +177,17 @@ vector<vector<double> > &v0) const
}
////////////////////////////////////////////////////////////////////////////////
double DistanceConstraint::projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const
void DistanceConstraint::compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f)
{
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);
const double* pf1 = &f[is1_][3*ia1_];
const double* pf2 = &f[is2_][3*ia2_];
D3vector f1(pf1);
D3vector f2(pf2);
// compute gradient at r
D3vector r12(r1-r2);
......@@ -197,30 +197,34 @@ double DistanceConstraint::projection(const vector<vector<double> > &r0,
const double norm2 = g1*g1;
assert(norm2>=0.0);
if ( norm2 == 0.0 )
return 0.0;
{
force_ = 0.0;
return;
}
const double norm = sqrt(norm2);
D3vector e1(g1/norm);
D3vector e2(-e1);
const double proj1 = x1*e1;
const double proj2 = x2*e2;
const double proj1 = f1*e1;
const double proj2 = f2*e2;
return 0.5*(proj1+proj2);
force_ = -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 << " <constraint name=\"" << name();
os << "\" type=\"" << type();
os << "\" atoms=\"" << name1_ << " " << name2_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << setw(10) << setprecision(6) << distance_ << " "
<< setw(10) << setprecision(6) << velocity_ << " -->";
os << " value=\"" << setprecision(6) << distance_;
os << "\" velocity=\"" << setprecision(6) << velocity_ << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
return os;
}
......
......@@ -3,7 +3,7 @@
// DistanceConstraint.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceConstraint.h,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: DistanceConstraint.h,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#ifndef DISTANCECONSTRAINT_H
#define DISTANCECONSTRAINT_H
......@@ -18,23 +18,28 @@ class DistanceConstraint : public Constraint
string name1_, name2_;
int ia1_, ia2_, is1_, is2_;
double m1_, m2_, m1_inv_, m2_inv_;
double distance_, velocity_, tol_;
double distance_, velocity_, force_, weight_, tol_;
public:
DistanceConstraint(string name1, string name2,
DistanceConstraint(string name, 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)
{
name1_(name1), name2_(name2), distance_(distance),
velocity_(velocity), tol_(tolerance), m1_(0.0), m2_(0.0)
{
name_ = name;
names_.resize(2);
names_[0] = name1_;
names_[1] = name2_;
force_ = 0.0;
weight_ = 1.0;
}
string type(void) const { return "distance"; }
double value(void) const { return distance_; }
double velocity(void) const { return velocity_; }
double force(void) const { return force_; }
double weight(void) const { return weight_; }
double tolerance(void) const { return tol_; }
void set_value(double value)
{
......@@ -51,8 +56,8 @@ class DistanceConstraint : public Constraint
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;
void compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f);
ostream& print( ostream& os );
};
......
#-------------------------------------------------------------------------------
# $Id: Makefile,v 1.35 2005-06-27 22:37:32 fgygi Exp $
# $Id: Makefile,v 1.36 2005-09-16 23:08:11 fgygi Exp $
#------------------------------------------------------------------------------
#
include $(TARGET).mk
......@@ -258,9 +258,6 @@ LDAFunctional.o: XCFunctional.h
ListAtomsCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h
ListAtomsCmd.o: D3vector.h ConstraintSet.h Wavefunction.h UnitCell.h
ListAtomsCmd.o: Control.h
ListConstraintsCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h
ListConstraintsCmd.o: D3vector.h ConstraintSet.h Wavefunction.h UnitCell.h
ListConstraintsCmd.o: Control.h
ListSpeciesCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h
ListSpeciesCmd.o: D3vector.h ConstraintSet.h Wavefunction.h UnitCell.h
ListSpeciesCmd.o: Control.h
......@@ -327,14 +324,13 @@ PSDWavefunctionStepper.o: UnitCell.h Control.h Timer.h
qb.o: isodate.h release.h qbox_xmlns.h Context.h UserInterface.h Sample.h
qb.o: AtomSet.h Atom.h D3vector.h ConstraintSet.h Wavefunction.h UnitCell.h
qb.o: Control.h Timer.h AngleCmd.h AtomCmd.h ConstraintCmd.h DistanceCmd.h
qb.o: HelpCmd.h ListAtomsCmd.h ListConstraintsCmd.h ListSpeciesCmd.h
qb.o: LoadCmd.h MoveCmd.h PrintCmd.h QuitCmd.h RandomizeWfCmd.h RunCmd.h
qb.o: SaveCmd.h SetCmd.h SpeciesCmd.h StatusCmd.h TorsionCmd.h AtomsDyn.h
qb.o: Cell.h CellDyn.h SlaterDet.h Basis.h Matrix.h CellLock.h CellMass.h
qb.o: ChargeMixCoeff.h ChargeMixRcut.h Debug.h Ecut.h Ecutprec.h Ecuts.h
qb.o: Emass.h ExtStress.h FermiTemp.h Dt.h Nempty.h NetCharge.h Nrowmax.h
qb.o: RefCell.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h WfDiag.h
qb.o: WfDyn.h Xc.h
qb.o: HelpCmd.h ListAtomsCmd.h ListSpeciesCmd.h LoadCmd.h MoveCmd.h
qb.o: PrintCmd.h QuitCmd.h RandomizeWfCmd.h RunCmd.h SaveCmd.h SetCmd.h
qb.o: SpeciesCmd.h StatusCmd.h TorsionCmd.h AtomsDyn.h Cell.h CellDyn.h
qb.o: SlaterDet.h Basis.h Matrix.h CellLock.h CellMass.h ChargeMixCoeff.h
qb.o: ChargeMixRcut.h Debug.h Ecut.h Ecutprec.h Ecuts.h Emass.h ExtStress.h
qb.o: FermiTemp.h Dt.h Nempty.h NetCharge.h Nrowmax.h RefCell.h Stress.h
qb.o: Thermostat.h ThTemp.h ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
QuitCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
QuitCmd.o: ConstraintSet.h Wavefunction.h UnitCell.h Control.h
......
......@@ -3,7 +3,7 @@
// SampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleStepper.C,v 1.19 2005-06-27 22:19:12 fgygi Exp $
// $Id: SampleStepper.C,v 1.20 2005-09-16 23:08:11 fgygi Exp $
#include "SampleStepper.h"
#include "Species.h"
......@@ -13,7 +13,7 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
SampleStepper::SampleStepper(Sample& s) : s_(s), atoms_(s_.atoms)
SampleStepper::SampleStepper(Sample& s) : s_(s)
{
fion.resize(s_.atoms.nsp());
for ( int is = 0; is < fion.size(); is++ )
......@@ -121,13 +121,13 @@ void SampleStepper::compute_sigma(void)
{
sigma_kin = 0.0;
// compute kinetic contribution to stress using velocities at time t0
for ( int is = 0; is < atoms_.atom_list.size(); is++ )
for ( int is = 0; is < s_.atoms.atom_list.size(); is++ )
{
int i = 0;
double mass = atoms_.species_list[is]->mass() * 1822.89;
for ( int ia = 0; ia < atoms_.atom_list[is].size(); ia++ )
double mass = s_.atoms.species_list[is]->mass() * 1822.89;
for ( int ia = 0; ia < s_.atoms.atom_list[is].size(); ia++ )
{
Atom* pa = atoms_.atom_list[is][ia];
Atom* pa = s_.atoms.atom_list[is][ia];
D3vector v = pa->velocity();
const double vx = v.x;
const double vy = v.y;
......
......@@ -3,7 +3,7 @@
// SampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleStepper.h,v 1.8 2004-03-11 21:52:31 fgygi Exp $
// $Id: SampleStepper.h,v 1.9 2005-09-16 23:08:11 fgygi Exp $
#ifndef SAMPLESTEPPER_H
#define SAMPLESTEPPER_H
......@@ -22,11 +22,6 @@ class SampleStepper
protected:
Sample& s_;
AtomSet& atoms_;
int nsp_;
int ndofs_;
vector<int> na_; // number of atoms per species na_[nsp_]
vector<double> pmass_; // pmass_[nsp_]
vector<vector<double> > fion;
valarray<double> sigma_eks, sigma_kin, sigma_ext, sigma;
......
......@@ -3,7 +3,7 @@
// SpeciesCmd.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SpeciesCmd.C,v 1.6 2005-06-27 22:18:35 fgygi Exp $
// $Id: SpeciesCmd.C,v 1.7 2005-09-16 23:08:11 fgygi Exp $
#include "SpeciesCmd.h"
#include "SpeciesReader.h"
......@@ -37,9 +37,9 @@ int SpeciesCmd::action(int argc, char **argv)
s->atoms.addSpecies(sp,argv[1]);
if ( s->ctxt_.onpe0() )
{
cout << endl << " species " << sp->name() << ":" << endl;
cout << endl << " <!-- species " << sp->name() << ":" << endl;
sp->info(cout);
cout << endl;
cout << " -->" << endl;
}
}
......
......@@ -3,7 +3,7 @@
// TorsionConstraint.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: TorsionConstraint.C,v 1.1 2005-06-27 22:34:46 fgygi Exp $
// $Id: TorsionConstraint.C,v 1.2 2005-09-16 23:08:11 fgygi Exp $
#include "TorsionConstraint.h"
#include "AtomSet.h"
......@@ -269,27 +269,27 @@ bool TorsionConstraint::enforce_v(const vector<vector<double> > &r0,
}
////////////////////////////////////////////////////////////////////////////////
double TorsionConstraint::projection(const vector<vector<double> > &r0,
const vector<vector<double> > &x) const
void TorsionConstraint::compute_force(const vector<vector<double> > &r0,
const vector<vector<double> > &f)
{
const double* pr1 = &r0[is1_][3*ia1_];
const double* pr2 = &r0[is2_][3*ia2_];
const double* pr3 = &r0[is3_][3*ia3_];
const double* pr4 = &r0[is4_][3*ia4_];
const double* px1 = &x[is1_][3*ia1_];
const double* px2 = &x[is2_][3*ia2_];
const double* px3 = &x[is3_][3*ia3_];
const double* px4 = &x[is4_][3*ia4_];
const double* pf1 = &f[is1_][3*ia1_];
const double* pf2 = &f[is2_][3*ia2_];
const double* pf3 = &f[is3_][3*ia3_];
const double* pf4 = &f[is4_][3*ia4_];
D3vector r1(pr1);
D3vector r2(pr2);
D3vector r3(pr3);
D3vector r4(pr4);
D3vector x1(px1);
D3vector x2(px2);
D3vector x3(px3);
D3vector x4(px4);
D3vector f1(pf1);
D3vector f2(pf2);
D3vector f3(pf3);
D3vector f4(pf4);
const double h = 0.001;
const double fac = 0.5 / h;
......@@ -303,26 +303,36 @@ double TorsionConstraint::projection(const vector<vector<double> > &r0,
const double norm2 = g1*g1 + g2*g2 + g3*g3 +g4*g4;
assert(norm2>0.0);
const double proj = x1*g1 + x2*g2 + x3*g3 + x4*g4;
const double proj = f1*g1 + f2*g2 + f3*g3 + f4*g4;
if ( norm2 == 0.0 )
return 0.0;
return proj/sqrt(norm2);
{
force_ = 0.0;
return;
}
force_ = -proj/norm2;
// compute weight
const double z = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2 +
m3_inv_ * g3 * g3 +
m4_inv_ * g4 * g4;
assert(z > 0.0);
weight_ = 1.0 / sqrt(z);
}
////////////////////////////////////////////////////////////////////////////////
ostream& TorsionConstraint::print( ostream &os )
{
os.setf(ios::left,ios::adjustfield);
os << " <!-- constraint ";
os << type() << " ";
os << setw(4) << name1_ << " ";