OFF-EP 0007 — Add long-range dispersion attribute in vdW section

[mattwthompson]

Implements a proposal as discussed recently; follow links in the EP for more information.

Supersedes #22

cc: @mikemhenry

[j-wags]

Thanks for the excellent writeup. This looks good to me!

[davidlmobley]

Looks good, thanks for getting this in place/removing the ambiguity.

[jchodera]

I've just taken a pass through this. I'm definitely in favor of removing ambiguity, but this appears to be too simplified.

For example, it is very likely we will want to use LJPME instead of a simple isotropic correction.

In addition, isotropic is ambiguous. Does it mean that the entire vdW energy function is integrated out to infinity against the 4*pi*r^2 surface term, as \int_0^\infty dr 4*pi*r^2 U(r)? Or just the dispersive (attractive) part? If a switching function S(r) is in use, is a factor of 1-S(r) included?

[mattwthompson]

Thanks for the feedback!

For example, it is very likely we will want to use LJPME instead of a simple isotropic correction.

This seems to be where we'd all like things to end up, but until then we need to provide guidance on the non-LJPME approach. As I understand it, these are mutually exclusive choices; one could either

                    non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.PME)
                    non_bonded_force.setUseDispersionCorrection(True)

or

                    non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.LJPME)

If it's a matter of clarity, and my understanding is accurate, I'd be happy to add in a comment along the lines of "If LJ-PME is used, long_range_dispersion should be ignored."

In addition, isotropic is ambiguous.

I completely agree and I admit I don't know what it's supposed to mean. It's been in the spec for years and as far as I understand it, the downstream effect is simply to flip non_bonded_force.setUseDispersionCorrection(True). Answering your questions would involve digging into the OpenMM source code and/or your paper from 15 years ago and figuring out how exactly it's implemented and/or how it's suggested to be implemented. The OpenMM docs provide good user-facing guidance but (necessarily) not this level of detail:

Set whether to add a contribution to the energy that approximately represents the effect of Lennard-Jones interactions beyond the cutoff distance. The energy depends on the volume of the periodic box, and is only applicable when periodic boundary conditions are used. When running simulations at constant pressure, adding this contribution can improve the quality of results.

Perhaps it would be simplest to ask if you recall what the intended meaning of "isotropic" was when it was originally added into the spec?

[davidlmobley]

I agree with @jchodera 's questions, but I think these have to be dealt with separately from this EP. In particular, right now the spec is ambiguous or makes an unstated assumption -- that we are using an isotropic long-range dispersion correction for LJ. The immediate "stop the bleeding" response to this should be to state this assumption, which is what this EP proposes to do, and I think it's the right path forward.

John's questions here are good:

In addition, isotropic is ambiguous. Does it mean that the entire vdW energy function is integrated out to infinity against the 4pir^2 surface term, as \int_0^\infty dr 4pir^2 U(r)? Or just the dispersive (attractive) part? If a switching function S(r) is in use, is a factor of 1-S(r) included?

But... we can't address this by changing the spec because, as Matt notes, this is not even addressed in the OpenMM docs, and even if it were, our spec couldn't control the form of the isotropic correction. I propose:
a) We add a note here saying how we achieve this when using OpenMM (non_bonded_force.setUseDispersionCorrection(True)) and
b) We open an issue on the OpenMM issue tracker to get John's questions answered/addressed in the OpenMM docs

And then relating to LJPME:

For example, it is very likely we will want to use LJPME instead of a simple isotropic correction.

Yes, but that's a separate issue from the change here, which is to remove the unstated assumption our implementation of the spec is currently making (that it should be using an isotropic LJ dispersion correction). So any changes in that regard should come in a separate EP.

I agree with @mattwthompson 's proposed solution:

If it's a matter of clarity, and my understanding is accurate, I'd be happy to add in a comment along the lines of "If LJ-PME is used, long_range_dispersion should be ignored."

[jchodera]

As I understand it, these are mutually exclusive choices; one could either

                    non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.PME)
                    non_bonded_force.setUseDispersionCorrection(True)

or

                    non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.LJPME)

Let's not make the mistake of confusing the actual logical choices here with how OpenMM chooses to implement things via its API (especially the convoluted NonbondedForce).

I propose the following:

1. We should allow the following very popular options or dispersion correction:
   a. Do not add a long-range dispersion correction
   b. Use an isotropic long-range dispersion correction that corresponds to the most popular and sensible choice
   c. Use LJPME for vdW, which should either be specified here or in the nonbonded method for vdW (in which case the dispersion correction is ignored). This can only apply to vdW potentials that include sums of inverse even powers of r
2. Be as explicit as possible about what is meant. For (b), this would be "To compute the dispersion correction, the vdW potential U_{ij}(r) between pairs of LJ sites i and j is integrated as

U_dispersion = [N (N-1)]^{-1} \sum_{i < j} \int_{0}^{\infty} dr [4 \pi r^2] [1 - S(r)] U_{ij}(r)

where 4 \pi r^2 is the surface area of a sphere of radius r and S(r) is the switching function or function representing a cutoff, where S(r) = 0 outside the cutoff and S(r) = 1 inside the cutoff before the switching function has turned on.

[mattwthompson]

Thanks for all of the feedback. I have implemented the two changes suggested by @davidlmobley above and also fixed a few typos.

At this point I consider task I set out to do - contribute to the specification a detail that has been implicit for quite some period of time but in a way that enables future extensions - as accomplished. I advocate that further improvements should be handled in future EPs brought forth by members of the committee or the community more broadly. I understand that means this EP may not be accepted.

[davidlmobley]

@jchodera I propose we approve/merge this EP and migrate your suggestions to a separate EP. I don't yet understand (nor does @mattwthompson , I think) how to translate them to a specific change to the spec, which to me is a strong argument to making that a SEPARATE change to the spec.

This EP fixes the problem it was intended to fix, in my view; the fact that there are other, related problems doesn't mean they MUST be dealt with in the same EP.

[jchodera]

I don't yet understand (nor does @mattwthompson , I think) how to translate them to a specific change to the spec, which to me is a strong argument to making that a SEPARATE change to the spec.

@davidlmobley : I've clarified my points in an alternative EP: #44

#44 should provide a cleaner path forward that also allows us to start doing force field science in cases where LJPME is important for accurately modeling heterogeneous systems like membrane proteins.