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