Enable support for all OpenMM GBSA models

[choderalab]

Enable support for all OpenMM GBSA models

[jchodera]

In order to support the additional GB models implemented via CustomGBForce, we need to rewrite the computed value and energy expressions using these rules:

• lambda_electrostatics is added as a global parameter:

custom_force.addGlobalParameter("lambda_electrostatics", 1.0)

• alchemical is added as a per-particle parameter:

custom_force.addPerParticleParameter("alchemical");

All atoms in the alchemical group have this parameter set to 1; otherwise 0.

• Any single-particle energy term (CustomGBForce.SingleParticle) is scaled by (lambda_electrostatics*alchemical+(1-alchemical))
• Any two-particle energy term (CustomGBForce.ParticlePairNoExclusions) has charge 1 (q1) replaced by (lambda_electrostatics*alchemical1+(1-alchemical1))*q1 and charge 2 (q2) replaced by (lambda_electrostatics*alchemical2+(1-alchemical2))*q2.
• Any single-particle computed value (CustomGBForce.SingleParticle) remains unmodified
• Any two-particle computed value (CustomGBForce.ParticlePairNoExclusions) is scaled by (lambda_electrostatics*alchemical2 + (1-alchemical2))

Scaling of a term should always prepend and capture the value with an intermediate variable. For example, this complex expression

        self.addComputedValue("I", "Ivdw+neckScale*Ineck;"
                                   "Ineck=step(radius1+radius2+neckCut-r)*getm0(radindex1,radindex2)/(1+100*(r-getd0(radindex1,radindex2))^2+"
                                   "0.3*1000000*(r-getd0(radindex1,radindex2))^6);"
                                   "Ivdw=select(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), 0);"
                                   "U=r+sr2;"
                                   "L=max(or1, D);"
                                   "D=abs(r-sr2);"
                                   "radius1=or1+offset; radius2=or2+offset;"
                                   "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)

can be scaled by prependingscaling * unscaled; unscaled =

[jchodera]

openmm/openmm#1641 has been merged, so I can move ahead with implementing this.