{
if (ir->efep != efepNO)
{
+ real bondA,bondB;
dt_2 = 1/sqr(ir->delta_t);
dvdl = 0;
for (i = 0; i < ncons; i++)
{
type = idef->il[F_CONSTR].iatoms[3*i];
- dvdl += lagr[i]*dt_2*
- (idef->iparams[type].constr.dB-idef->iparams[type].constr.dA);
+
+ /* dh/dl contribution from constraint force is dh/dr (constraint force) dot dr/dl */
+ /* constraint force is -\sum_i lagr_i* d(constraint)/dr, with constrant = r^2-d^2 */
+ /* constraint force is -\sum_i lagr_i* 2 r */
+ /* so dh/dl = -\sum_i lagr_i* 2 r * dr/dl */
+ /* However, by comparison with lincs and with
+ comparison with a full thermodynamics cycle (see
+ redmine issue #1255), this is off by a factor of
+ two -- the 2r should apparently just be r. Further
+ investigation should be done at some point to
+ understand why and see if there is something deeper
+ we are missing */
+
+ bondA = idef->iparams[type].constr.dA;
+ bondB = idef->iparams[type].constr.dB;
+ dvdl += lagr[i] * dt_2 * ((1.0-lambda)*bondA + lambda*bondB) * (bondB-bondA);
}
*dvdlambda += dvdl;
}