Add the missing term from Alch OBC2
Adds the missing `C` term to the alchemical OBC2 GB Force
Partial resolution to
@@ -1776,11 +1776,11 @@ class AbsoluteAlchemicalFactory(object):
custom_force = openmm.CustomGBForce()
# Add per-particle parameters.
- custom_force.addGlobalParameter("lambda_electrostatics", 1.0);
- custom_force.addPerParticleParameter("charge");
- custom_force.addPerParticleParameter("radius");
- custom_force.addPerParticleParameter("scale");
- custom_force.addPerParticleParameter("alchemical");
+ custom_force.addGlobalParameter("lambda_electrostatics", 1.0)
+ custom_force.addPerParticleParameter("charge")
+ custom_force.addPerParticleParameter("radius")
+ custom_force.addPerParticleParameter("scale")
+ custom_force.addPerParticleParameter("alchemical")
# Set nonbonded method.
custom_force.setNonbondedMethod(reference_force.getNonbondedMethod())
@@ -1791,8 +1791,9 @@ class AbsoluteAlchemicalFactory(object):
custom_force.addGlobalParameter("soluteDielectric", reference_force.getSoluteDielectric())
custom_force.addGlobalParameter("offset", 0.009)
- custom_force.addComputedValue("I", "(lambda_electrostatics*alchemical2 + (1-alchemical2))*step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+ custom_force.addComputedValue("I", "(lambda_electrostatics*alchemical2 + (1-alchemical2))*step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r+C);"
"U=r+sr2;"
+ "C=2*(1/or1-1/L)*step(sr2-r-or1)"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
