Open jchodera opened 8 years ago
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.
CustomGBForce.SingleParticle
) is scaled by (lambda_electrostatics*alchemical+(1-alchemical))
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
.CustomGBForce.SingleParticle
) remains unmodifiedCustomGBForce.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 =
https://github.com/pandegroup/openmm/pull/1641 has been merged, so I can move ahead with implementing this.
Looks like https://github.com/pandegroup/openmm/pull/1311 may have our GBSA support? Looking into this...