...
 
Commits (5)
......@@ -285,7 +285,8 @@ void AngleConstraint::compute_force(const vector<vector<double> > &r0,
}
const double proj = f1*g1 + f2*g2 + f3*g3;
force_ = -proj/norm2;
// A positive force on the angle constraint tends to increase the angle
force_ = proj/norm2;
// compute weight
const double z = m1_inv_ * g1 * g1 +
......
......@@ -506,6 +506,7 @@ void BOSampleStepper::step(int niter)
{
if ( s_.constraints.size() > 0 )
{
s_.constraints.update_constraints(dt);
s_.constraints.compute_forces(ionic_stepper->r0(), fion);
if ( onpe0 )
{
......@@ -1241,9 +1242,6 @@ void BOSampleStepper::step(int niter)
}
}
if ( atoms_move )
s_.constraints.update_constraints(dt);
// if using force_tol or stress_tol, check if maxforce and maxstress
// within tolerance
if ( force_tol > 0.0 )
......
......@@ -216,7 +216,8 @@ void DistanceConstraint::compute_force(const vector<vector<double> > &r0,
const double proj1 = f1*e1;
const double proj2 = f2*e2;
force_ = -0.5*(proj1+proj2);
// A positive force on the constraint tends to increase the distance
force_ = 0.5*(proj1+proj2);
}
////////////////////////////////////////////////////////////////////////////////
......
......@@ -50,16 +50,18 @@ double PairExtForce::energy(const vector<vector<double> > &r,
D3vector r1(pr1);
D3vector r2(pr2);
// A positive value of force_ corresponds to a repulsive force
// i.e. a positive force tends to increase the distance
D3vector e12 = normalized(r2-r1);
f[is1_][3*ia1_+0] += force_ * e12.x;
f[is1_][3*ia1_+1] += force_ * e12.y;
f[is1_][3*ia1_+2] += force_ * e12.z;
f[is1_][3*ia1_+0] -= force_ * e12.x;
f[is1_][3*ia1_+1] -= force_ * e12.y;
f[is1_][3*ia1_+2] -= force_ * e12.z;
f[is2_][3*ia2_+0] -= force_ * e12.x;
f[is2_][3*ia2_+1] -= force_ * e12.y;
f[is2_][3*ia2_+2] -= force_ * e12.z;
f[is2_][3*ia2_+0] += force_ * e12.x;
f[is2_][3*ia2_+1] += force_ * e12.y;
f[is2_][3*ia2_+2] += force_ * e12.z;
return 2 * force_ * length(r2-r1);
return -2 * force_ * length(r2-r1);
}
////////////////////////////////////////////////////////////////////////////////
......
......@@ -314,7 +314,8 @@ void TorsionConstraint::compute_force(const vector<vector<double> > &r0,
force_ = 0.0;
return;
}
force_ = -proj/norm2;
// A positive force on the torsion constraint tends to increase the angle
force_ = proj/norm2;
// compute weight
const double z = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2 +
......
--------------------------------------------------------------------------------
rel1_71_0
--------------------------------------------------------------------------------
a26cbfa Update constraints before ionic step
88181a8 Change sign of force on constraints and pair force
--------------------------------------------------------------------------------
rel1_70_0
--------------------------------------------------------------------------------
c790914 Fix xsd links in upf2qso dir
......
......@@ -19,5 +19,5 @@
#include "release.h"
std::string release(void)
{
return std::string("rel1_70_0");
return std::string("rel1_71_0");
}