Merge branch 'develop'

src/AngleConstraint.C

@@ -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 +

src/BOSampleStepper.C

@@ -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 )

src/DistanceConstraint.C

@@ -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);
 }
 
 ////////////////////////////////////////////////////////////////////////////////

src/PairExtForce.C

@@ -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);
 }
 
 ////////////////////////////////////////////////////////////////////////////////

src/TorsionConstraint.C

@@ -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 +

src/release.C

@@ -19,5 +19,5 @@
 #include "release.h"
 std::string release(void)
 {
-  return std::string("rel1_70_0");
+  return std::string("rel1_70_0dev");
 }