Improve PME overlap by computing reciprocal-space exceptions

[jchodera]

We can improve the PME overlap by subtracting out the exceptions computed in reciprocal space.

The relevant OpenMM code that implements this is located here.

To implement this, we would add a CustomBondForce that computes all exceptions and exclusions as bond terms using the direct-space Ewald energy:

-(lambda_electrostatics^softcore_d)*ONE_4PI_EPS0*chargeprod*erf(alpha_ewald*reff_electrostatics)/reff_electrostatics

(note the minus sign, since we are subtracting these).

chargeprod would be computed as charge1*charge2, and would NOT match the chargeprod used in the exception/exclusion. This is because we are subtracting off the PME reciprocal space contribution, which uses the product of charges.

We would iterate over all exceptions and exclusions in the original NonbondedForce.

I think we can just use r instead of bothering to compute softcore reff_electrostatics.

[andrrizzi]

Just a typo above, the formula actually include erf, not erfc.

Are we sure we are not introducing artifacts by including the correction for the reciprocal space exceptions without including the reciprocal space energy? My first guess is the our potential will be "more correct" if we don't. We can certainly try them both and see how the overlap changes.

[jchodera]

Are we sure we are not introducing artifacts by including the correction for the reciprocal space exceptions without including the reciprocal space energy? My first guess is the our potential will be "more correct" if we don't. We can certainly try them both and see how the overlap changes.

We are aiming to include any terms that are very sensitive to ligand conformation or position. If the terms are relatively insensitive, then we don't need to include them.

My physical intuition is that many of these terms are exceptions due to (1,2) directly bonded atoms. Those will be pretty constant, but the (1,3) or (1,4) exclusions/exceptions will likely be very sensitive to ligand conformation. So I'm guessing we will get a significant improvement in overlap.

[andrrizzi]

My doubts arise from the fact that this correction is not a physical force, but it's added only to cancel the incorrect contribution from the reciprocal space, so the reference system should not feel this.

I can try to make a quick experiment and test this locally.

[Lnaden]

I'm in agreement with @jchodera on this, since the force we are comparing to IS the full PME force, we should try to minimize the number of things we have to exclude so we are not changing the energy in the wrong direction relative to the reference state. You saw in my script how bad the direct space comparison is.

[andrrizzi]

so we are not changing the energy in the wrong direction relative to the reference state

I think that we'll do exactly that if we do include the correction. I'll try to make this clearer so that you can tell me if I'm misunderstanding something. If we decompose reciprocal space and direct space energy that OpenMM gives us in

reciprocal_energy = true_reciprocal_energy + exclusion_reciprocal_energy
direct_energy = true_direct_energy - exclusion_reciprocal_energy

then the reference system energy will be

reference_energy = true_direct_energy + true_reciprocal_energy + exclusion_reciprocal_energy - exclusion_reciprocal_energy
                 = true_direct_energy + true_reciprocal_energy

currently our alchemical system energy is true_direct_energy, if we added the exclusion correction we would get

alchemical_energy = true_direct_energy - exclusion_reciprocal_energy

but the force arising from the exclusion_reciprocal_energy term is not in the reference system at all, so I don't understand why this should increase our overlap with the reference system.

We can also talk about this in person on Monday, if it's easier to explain.

[jchodera]

It's totally not a physical force, I agree. And @andrrizzi may actually be right now that I think about it: if we are adding some f(x) but the full PME includes -f(x) in reciprocal space and +f(x) in direct space, we aren't helping anything.

[jchodera]

Put another way: We may be making what we compute closer to the direct space part of real PME but farther from the total PME energy because we are adding a term that cancels out in full PME.

[Lnaden]

Good point. If its easy to test, we can try, but if not, go for @andrrizzi suggestion