Commit 88181a82 by Francois Gygi

Change sign of force on constraints and pair force

The change of sign makes the value of the force consistent with
the definition of the force ( -grad(E) )
A positive value of the constraint force tends to increase the
constraint value.
The change of sign of PairForce makes it such that a positive
PairForce value tends to increase the bond length.
The value of the PairForce needed to compensate for a distance
constraint is opposite to the force acting on the constraint.
parent 5e226066
......@@ -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 +
......
......@@ -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 +
......
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