CutoffNonPeriodic should not use reaction-field electrostatics

[jchodera]

When a cutoff is used in a non-periodic system, reaction-field electrostatics should not be used.

I may be reading the code wrong, but this seems to suggest that reaction-field electrostatics are used whenever a cutoff is employed:
https://github.com/SimTk/openmm/blob/master/platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp#L590-L591
Also, the tests (specifically, testCutoff() in TestReferenceNonbondedForce) test the CutoffNonPeriodic condition against reaction-field electrostatics:
https://github.com/SimTk/openmm/blob/master/platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp#L590-L591

[peastman]

Why? Whether you want to assume a uniform dielectric beyond the cutoff has nothing to do with whether you're modelling a periodic system.

[jchodera]

Because the reaction field approximation is derived assuming that the entirety of the region within the cutoff is filled with explicit particles. Near the boundaries of a nonperiodic system, the reaction field form you use is no longer an accurate model of anything.

[jchodera]

I think we may have previously discussed leaving this in (during a discussion about GB with cutoffs) because it was sort of like a switching function, but now that we have a real switching function that is twice continuously differentiable, that need is no longer present.

[peastman]

The switching function is only applied to LJ, not to Coulomb. As for what's realistic, that depends on what you're trying to model. If you're doing explicit solvent without periodic boundary conditions, that won't be realistic at the edge of the system no matter what you do. But towards the interior (which is what you presumably care about in that case), RF will be more realistic.

If you're using implicit solvent, then things are completely different. We turn off RF in that case. We discussed that in lots of detail some months ago.

[jchodera]

The switching function is only applied to LJ, not to Coulomb. As for what's realistic, that depends on what you're trying to model. If you're doing explicit solvent without periodic boundary conditions, that won't be realistic at the edge of the system no matter what you do. But towards the interior (which is what you presumably care about in that case), RF will be more realistic.

Use of reaction field for nonperiodic cutoff simulations is highly nonstandard. I think it violates the idea that the default behavior should not be surprising to the user. I was just trying to simulate a little cluster of particles or molecules in the gas phase, and having RF on by default for CutoffNonPeriodic was very much surprising.

If you're using implicit solvent, then things are completely different. We turn off RF in that case. We discussed that in lots of detail some months ago.

This is also an argument for turning it off by default, since it could be surprising to the user that this is surprisingly that RF is turned off whenever a specific implicit solvent forces are added.

If there was a way to keep the RF functionality for CutoffNonPerdiodic by allow it to be switched off, and perhaps have it turned off by default, that would be fine! I don't think there's currently a way to do this via an on/off switch, but we can set the solvent dielectric to be 1, which is a bit nonintuitive.

An example use case that may apply to the forcefield parameterization work we've been discussing is the computation of potential energies of clusters of water molecules for comparison to QM data.

[peastman]

This is also an argument for turning it off by default, since it could be surprising to the user that this is surprisingly that RF is turned off whenever a specific implicit solvent forces are added.

But this is completely independent of periodic boundary conditions. GB can be used with or without periodic boundary conditions. RF can be used with or without periodic boundary conditions. Using RF and GB together is wrong, independent of whether or not you have periodic boundary conditions. Using the switching function would be equally wrong in that case. For explicit solvent, RF is a better approximation than the switching function, independent of whether or not you have periodic boundary conditions.

[jchodera]

Using the switching function would be equally wrong in that case. For explicit solvent, RF is a better approximation than the switching function, independent of whether or not you have periodic boundary conditions.

I think we're confusing the switching function (which is applied to Lennard-Jones interactions only) with various electrostatics models (e.g. RF and GB).

[peastman]

I think we're confusing the switching function (which is applied to Lennard-Jones interactions only) with various electrostatics models (e.g. RF and GB).

Isn't that exactly what you're proposing? That for CutoffNonPeriodic we should turn off RF and instead use the switching function for Coulomb? Or were you proposing just using a sharp cutoff with no sort of smoothing at all?

If that's not what you're proposing, then yes, I'm confused about what we're discussing.

[jchodera]

This is my fault for having started this thread while still in a confused state!

Now that I think I've got my head on right, here's the suggestion:

• We should have reaction-field electrostatics turned off for CutoffNonPeriodic by default, but we should still allow the user to use reaction field if desired. This would give the least surprising behavior, since (1) this is what most people expect when a simple cutoff is in use, and (2) would be consistent with the cutoff being disabled when GB forces are in use.
• We can allow the user to turn back on reaction-field for CutoffNonPeriodic if they wish.
• I think we should also add the capability to turn on a switching function (the same one used for Lennard-Jones interactions) for both CutoffNonPeriodic and CutoffPeriodic. This would allow the energy to be twice continuously differentiable at the cutoff. We may also allow this option for PME and Ewald, since there are some cases where the tiny truncation introduced by the sharp cutoff of the exp(-alpha*r)/r could still be noticeable.

Does that make more sense? Again, apologies about having things mixed up earlier!

[peastman]

I don't think we should default to using a sharp cutoff for Coulomb. That's just a really, really bad thing to do, and no one should ever do it. It produces all sorts of errors. The switching function is a little better, but still pretty bad: it tends to produce artifacts, such as artificially creating local minima in the potential.

If someone really really wants a sharp cutoff, they can do it by setting the reaction field dielectric to 1, or by using a CustomNonbondedForce. But we shouldn't make it easy for them, because they almost certainly will get a bad result.

[jchodera]

I agree that a sharp cutoff is bad, which is why a switching function is a good idea. I'd be happy with switching functions to always be on by default, actually, and the user switching them off if they think they can live with the consequences.

For a correctly-chosen cutoff, the artifacts introduced by switching should not introduce any big problems for the potential. If the switch is creating problematic artificial minima, then the cutoff is way too short.

I still think that a reaction field form for cutoff electrostatics is very, very weird. Have you seen anyone publish this method anywhere?

[peastman]

Have you seen anyone publish this method anywhere?

Of course. Reaction field is one of the standard methods used for this, and is generally preferred over all other cutoff methods when you aren't using PME. In contrast, switching functions are well known to produce serious artifacts. I can dig up some references for you if you really want. But let me just quote what the Gromacs manual says for the option coulombtype=switch:

"Analogous to Switch for vdwtype. Switching the Coulomb potential can lead to serious artifacts, advice: use Reaction-Field-zero instead."

[jchodera]

I sorry! I was specifically referring to the use of reaction field without periodic boundary conditions.

[peastman]

I sorry! I was specifically referring to the use of reaction field without periodic boundary conditions.

So was I. Or rather, I was nonspecifically referring to RF with or without periodic boundary conditions. :) The switching function is a bad choice either way.

I'm not clear exactly what problem you're trying to solve. Is it just that you want C2 continuity at the cutoff? If so, that's easy to do by setting the reaction field dielectric to infinity.

[jchodera]

So was I. Or rather, I was nonspecifically referring to RF with or without periodic boundary conditions. :) The switching function is a bad choice either way.

There are few points here:

• I do not believe reaction-field was derived to give anything sensible in the CutoffNonPeriodic case. I am requesting a copy of the original article to verify this.
• I don't believe anyone else has published a paper using reaction-field electrostatics for a CutoffNonPeriodic equivalent situation.
• Many other codes implement switching function options as standard way of dealing with achieving C2 continuity in the CutoffNonPeriodic case.
• We turn off the reaction-field electrostatics when certain Force objects are present in the system.

Given these four points, I am arguing that it violates good software principles that the default behavior for CutoffNonPeriodic is to use reaction-field electrostatics with a dielectric near 80.

I'm not clear exactly what problem you're trying to solve. Is it just that you want C2 continuity at the cutoff? If so, that's easy to do by setting the reaction field dielectric to infinity.

I think this is a good alternative to support in addition to switching, especially if you can actually set the dielectric to infinity. But I don't believe it should be the main or default choice, since there is very little literature support for the use of reaction field in a CutoffNonPeriodic case to achieve C2 continuity.

If we had time to assemble a paper characterizing the use of reaction-field in a CutoffNonPeriodic situation as a superior (or at least comparable) alternative to switching functions for achieving C2 continuity, that would be a totally different story. We could just cite that paper and be done with it. But I have many reservations about our current choice of default without this.

[peastman]

I really don't understand what you're getting at. The reaction field approximation is based on one simple assumption: that everything beyond the cutoff is filled with a uniform dielectric. Whether that's a good assumption depends on what physical system you're trying to model. But in many cases (including many non-periodic cases), it is.

In contrast, the switching function is not derived from any physical approximation, and is not a reasonable model of any physical system. Furthermore, every paper I've seen on the subject has concluded it's a terrible way to treat Coulomb interactions and no one should ever use it. Yes, some programs support it anyway, but that's mostly for historical reasons. They now advise people not to use it. Like the bit I quoted from the Gromacs manual. It recommends using reaction field instead, without reference to whether or not you have periodic boundary conditions.

We turn off the reaction-field electrostatics when certain Force objects are present in the system.

That's a completely different issue. We turn off reaction field when using GB because GB was not designed to work with reaction field. This is true independent of whether or not you have periodic boundary conditions. (Yes, GB is sometimes used with periodic boundary conditions: think of crystals, membranes, etc.) And as I said, we already turn it off in that case.

[jchodera]

I really don't understand what you're getting at. The reaction field approximation is based on one simple assumption: that everything beyond the cutoff is filled with a uniform dielectric. Whether that's a good assumption depends on what physical system you're trying to model. But in many cases (including many non-periodic cases), it is.

How is this a good assumption in non-periodic cases? Using the reaction-field approximation in this case will effectively assume there is a layer of empty vacuum the width of the cutoff surrounding your entire molecule, with uniform dielectric outside. What kind of physical situation does that correspond to?

In contrast, the switching function is not derived from any physical approximation, and is not a reasonable model of any physical system. Furthermore, every paper I've seen on the subject has concluded it's a terrible way to treat Coulomb interactions and no one should ever use it. Yes, some programs support it anyway, but that's mostly for historical reasons. They now advise people not to use it. Like the bit I quoted from the Gromacs manual. It recommends using reaction field instead, without reference to whether or not you have periodic boundary conditions.

Unless we can point to a paper showing the deficiencies in the switching function and superiority of the reaction-field for non-periodic cases, we really can't defensibly use this scheme! Also, there is plenty of literature showing why C2 continuity is important for Verlet-style integrators. Reaction field does not provide that, unless we support infinite dielectric, but then the deviation from the desired Coulomb potential grows rapidly with decreasing cutoff.

That's a completely different issue. We turn off reaction field when using GB because GB was not designed to work with reaction field. This is true independent of whether or not you have periodic boundary conditions. (Yes, GB is sometimes used with periodic boundary conditions: think of crystals, membranes, etc.) And as I said, we already turn it off in that case.

I'm just saying it is surprising behavior that adding some Force objects will affect the NonbondedForce objects, while other Force objects will not.

[jchodera]

Unless we can point to a paper showing the deficiencies in the switching function and superiority of the reaction-field for non-periodic cases, we really can't defensibly use this scheme!

I'll add that I would be very happy to work on such a paper together, if there was interest! We could come up with some tests to compare the various commonly-used switching and truncation schemes and compare with the RF scheme with infinite dielectric for a few systems in our testsystems module, for example, and compare to no-cutoff schemes. But I do think we need to be able to point to a piece of literature that shows the comparison.

[peastman]

See, for example,

http://www.unc.edu/~perera/PAPERS/JChemPhys_1995_102_450.pdf

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1301047/pdf/10969015.pdf

I'm sure I can dig up some others if you really want. The switching function is just a really bad way to treat Coulomb interactions, and no one should ever use it. We absolutely should not add it as an option, and we certainly shouldn't make it the default! If someone really really wants it, they can always implement it with a custom force. But I'm not going to make it easy for them, because like I said, no one should ever use it.

Using the reaction-field approximation in this case will effectively assume there is a layer of empty vacuum the width of the cutoff surrounding your entire molecule

That's not correct. When simulating explicit solvent without periodic boundary conditions, you typically embed the molecule of interest inside a droplet of water. You don't expect the behavior at the surface of the droplet to be accurate, but that's not the point of the simulation. The atoms in your molecule of interest do not "see" any sort of vacuum anywhere. They see explicit solvent out to the cutoff distance, and a constant dielectric beyond that. This is a good approximation.

Contrast that with the switching function, which assumes everything outside your explicit system is vacuum, and throws in some weird cutoff artifacts to boot.

By the way, here's a recent article that discusses several of the more modern methods for cutting off Coulomb interactions, including reaction field:

http://europepmc.org/articles/PMC3428531

[jchodera]

See, for example,
http://www.unc.edu/~perera/PAPERS/JChemPhys_1995_102_450.pdf
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1301047/pdf/10969015.pdf

Just so I know what to look for, are these supposed to be examples of non-periodic systems using reaction-field, or examples of deficiencies in switching functions?

(Just not sure what to be looking for specifically!)

[rmcgibbo]

I've been lurking on this thread, and learning a lot. (From that perspective, it's really awesome to have this conversation happen on github where I can lurk, as opposed to on email or something where it would be private).

With that said, the Berkowitz paper doesn't make any sense to me. They say:

Finally in the Ewald summation technique we take full advantage of the periodic boundary conditions by periodically replicating our simulation box. To achieve convergence of the Coulomb sum we build up the periodic system in roughly spherical layers. Outside the sphere we have to specify the nature of the dielectric surrounding the sphere, i.e., specify the dieletic constant of the continuum.... We run our Ewald simulation with these tin foil boundary conditions.

This makes absolutely no sense to me.

[swails]

@rmcgibbo, the statement you quoted is basically a paraphrasing of a description of the Ewald technique taken from either Allen and Tildesley or the Tuckerman book (both are great resources, IMO). I agree that in the given context they are rather opaque, so I would definitely suggest looking at the Ewald derivations in one of (or both of) those texts.

To achieve convergence of the Coulomb sum we build up the periodic system in roughly spherical layers.

As I understand it, the Ewald sum is only conditionally convergent. Since the net charge is formally 0, you have a sum that is more or less like 1 - 1/r + 1/r^2 - 1/r^3 + 1/r^4 --- an alternating harmonic series in the energy. We know that the harmonic series diverges (slowly, but it diverges), but the convergence criteria for alternating series is significantly less restrictive. Basically, your terms have to tend toward 0 in the infinite limit, which the alternating harmonic series does. One of the weird properties of an infinite, conditionally convergent sum like the one above is that its value actually depends on the order in which the summation is performed. Hence the term "roughly spherical layers" is defining the summation order by adding in particles in an order defined by their distance from the central particle (more or less -- closer periodic images are added before farther ones).

The boundary conditions (what is outside this infinitely periodic system that really has no boundary) adds a constant term to the Ewald sum based on the dielectric of this (non-existent) surrounding medium. If the dielectric is infinite (i.e., "tin-foil boundary"), no correction term need be added to the computed energy. If the dielectric is 1 (i.e., an infinite system embedded in a vacuum), there is some constant energy term you need to add to 'correct' for these boundary conditions. It does not affect forces at all, just the total energy (and not even the pressure tensor since it's not volume-dependent, to my understanding). It definitely helps to see the explanation in Allen and Tildesley.

[swails]

Hence the term "roughly spherical layers" is defining the summation order by adding in particles in an order defined by their distance from the central particle (more or less -- closer periodic images are added before farther ones).

As a clarification, this convention is used in the derivation of the Ewald summations. The Ewald summation itself is a set of two quickly converging series -- one in direct space using the minimum image convention with every charge damped by a neutralizing charge distribution and the other in reciprocal space computing the convolution of the cancelling charge distributions. This formulation is only possible, however, because an 'order' was defined for summing up the original conditionally convergent series. Basically, the Ewald sum only matches a specific choice for the order of summation of the direct particle-particle interactions over all periodic cells.

I can't do full justice to the Ewald sum in a Github post, though. Allen and Tildesley is to date the most helpful resource I've found (closely followed by Tuckerman). I have a somewhat more relaxed discussion of electrostatic summation methods (including Ewald, PME, cutoffs, switches, shifting, GRF, and IPS) in my dissertation ~p 55 [http://jswails.wikidot.com/about at the bottom]

[swails]

Back to the original discussion: I agree with both @peastman and @jchodera. Coulomb switching functions just plain suck. Between two oppositely charged particles I believe the switching region introduces a repulsive force, whereas two particles of the same charge introduces an attractive force where the switch is turned on. Yes the switch may be (twice!) differentiable, but this is just physically wrong.

The use cases for a RF approach in a non-periodic system seem far-fetched to me. There's an easy solution, though: use NoCutoff :). No switch, no RF. All of my GB or gas phase energy calculations employ an infinite cutoff so the RF is never toggled.

Unless somebody has a good reason for wanting to use a cutoff in a non-periodic system, they should use NoCutoff and this issue becomes moot. Besides, for the numbers I've seen NoCutoff is actually faster than CutoffNonPeriodic on the CUDA and OpenCL platforms, isn't it? And in systems so large that cutoffs become cheaper, wouldn't it be even cheaper to add some solvent molecules and switch to a O(N logN) method, anyway?

[jchodera]

Unless somebody has a good reason for wanting to use a cutoff in a non-periodic system, they should use NoCutoff and this issue becomes moot.

Well, I can confess to having played around with simulating a few copies of GroEL in a non-perodic system using implicit solvent models. NoCutoff really isn't an option there.

Besides, for the numbers I've seen NoCutoff is actually faster than CutoffNonPeriodic on the CUDA and OpenCL platforms, isn't it?

This is highly dependent on system size

[jchodera]

But I certainly see with your point, @swails. Doesn't this argue for simply disabling the CutoffNonPeriodic option?

[swails]

This is highly dependent on system size

Of course, I wasn't that clear in my statement, I don't think. At some point as the size of the system/solute increases, there's a "crossing point" where even omitting the solvent particles in a GB calculation isn't enough to make up for the O(N^2) scaling of the GB/gas-phase algorithm with no cutoff, and biting the bullet to add solvent molecules and run with PME where you can use a small cutoff with a O(N logN) method for long-range electrostatics is actually faster than the corresponding GB or gas-phase calculation.

If you use a cutoff with non-periodic systems this isn't true, I suppose (unless you're not using a neighbor list, but those can get huge for a very big cutoff).

[swails]

But I certainly see with your point, @swails. Doesn't this argue for simply disabling the CutoffNonPeriodic option?

To quote Doug Gwyn:

UNIX was not designed to stop its users from doing stupid things, as that would also stop them from doing clever things.

Substitute in OpenMM and I think that's a strong argument for keeping CutoffNonPeriodic around. Yes it can be abused (the same as someone may be tempted to run a minimal basis set HF calculation because MP2 with one of the Dunning sets is so much more expensive), but who knows how it can be used in the hands of somebody clever... You've used it by your own admission, and I don't think there's much risk of you abusing the method whose limitations you know well.

[peastman]

The use cases for a RF approach in a non-periodic system seem far-fetched to me. There's an easy solution, though: use NoCutoff :)

Perhaps so, but think of this from a Bayesian standpoint. :) The user is using a cutoff, they are not using periodic boundary conditions, and they are not using implicit solvent. How frequently that combination occurs is not the question. Consider all of those as givens, and our job is to figure out the best thing to do. There are two main possibilities:

1. This is a vacuum simulation. I hope they've specified a really large cutoff, because unscreened Coulomb interactions remain significant out to several nm. Perhaps they really shouldn't be using a cutoff, but maybe they're simulating a really huge system in vacuum where even a 5-10 nm cutoff is useful. Anyway, the best we can do is apply the cutoff in a way that has less tendency to produce artifacts (i.e. not a switching function).
2. This is an explicit solvent simulation. This does happen sometimes (I did some like this back in grad school), and the simulations are basically done in the way I described earlier. And RF is a pretty reasonable physical model in this situation.

[swails]

Perhaps so, but think of this from a Bayesian standpoint. :) The user is using a cutoff, they are not using periodic boundary conditions, and they are not using implicit solvent. How frequently that combination occurs is not the question. Consider all of those as givens, and our job is to figure out the best thing to do.

Indeed -- I was arguing for as much. I don't have any use-cases for cutoffs and RF in a non-periodic system, but we ought not to punish people that want to do something clever with this functionality just because somebody else might do something foolish :).

My summary from this thread:

Relating to CutoffNonPeriodic:

1. Using a hard cutoff is (almost) always foolish.
2. Using a switching function is almost always (maybe slightly less) foolish.
3. Using RF is usually (a little less) foolish.

The average user should avoid all 3, but the clever user should not be unnecessarily blocked.

Addressing @jchodera's concerns, it would probably be best if the defaults in the application layer steered people away from CutoffNonPeriodic, especially for users of the Amber, Gromacs, and Desmond (and hopefully soon CHARMM) wrappers.

[rmcgibbo]

@swails: Thanks for the info in Ewald!

[peastman]

I believe the conclusion of this discussion was not to change the behavior, so I'm closing it.

[jchodera]

I don't think this is quite settled. There have also been numerous reports that GBSA with a cutoff behaves very, very oddly, and part of that problem may be that CutoffNonPeriodic uses a reaction-field treatment of Coulomb electrostatics while the GBSAOBCForce is not adapted to use a reaction-field treatment for CutoffNonPeriodic.

[swails]

If we were picking the defaults today, I think that RF with CutoffNonPeriodic is marginally worse than no RF since I think it's more common to use GB in that scenario than no implicit solvent model at all.

The main reason IMO for keeping the current default is that it's been this way for a long time already. And of all nonbonded schemes to use, CutoffNonPeriodic is probably the stupidest choice for NonbondedForce. Combine that with the fact that all of the existing parsers already do the right thing, and I think this whole decision is rather unimportant. So I'm solidly +/-0 on changing this (I don't think it matters, and in that case inertia probably wins).

[jchodera]

Combine that with the fact that all of the existing parsers already do the right thing

Ah! I had forgotten about that line:
https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/internal/amber_file_parser.py#L1095

I'll keep looking for the origin of the issues with GB cutoffs then...

[peastman]

ForceField also disables RF if you add a GBSAOBCForce to the system:

https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/forcefield.py#L1267-L1272

[swails]

I think all System creators do. I checked Amber, Gromacs, and CHARMM (Desmond doesn't offer implicit solvent options at all), and as @peastman pointed out, so does ForceField.

[JoaoRodrigues]

Not to resurrect a dead issue, but I just wanted to add a potential use case where this (having the automatic RF correction) leads to very confusing results.

I wanted to calculate the pairwise energy of a particular pair of atoms and I found that setting force.setReactionFieldDielectric(1.0) still applies a correction to the energy. In CpuNonbondedForce.cpp, the total energy in case we specify CutoffNonPeriodic and a small cutoff, say 1 nm, is given by energy += (float) (chargeProd*(inverseR+krf*r2-crf));, where krf is indeed 0.0 but not crf. This is quite different from, for example, running the pairwise command on cpptraj, which gives the 'proper' energy without any correction.

As OpenMM is a library, it should not, in my opinion as a user, do this sort of decisions under the hood. I agree with @jchodera that there should be a way to absolutely disable the RF correction, even if it's the stupidest thing in the universe. I can use a CustomNonbondedForce to go around this, or just use a NoCutoff, but it means I'll suffer quite a bit with performance, either from the unoptimized kernels or because I'll have to go about calling setParticleParameters and updateParametersInContext for pair of atoms...

Sorry if this sounds a bit like a rant, but I spent two days figuring out why the NonbondedForce of a two atom system was not matching what I had on paper and what I wrote from a CustomNonbondedForce.

Edit. This also leads to this slightly mind-boggling effect: change the cutoff and see your energy change ...

If you want the code I was using, which is the simplest I could come up with, here it is.

import simtk.openmm as mm

system = mm.System()
system.addParticle(1.0)
system.addParticle(1.0)

force = mm.NonbondedForce()
force.addParticle(0.4478, 1, 0)
force.addParticle(-0.8188, 1, 0)
force.setCutoffDistance(2.9)  # anything lower should give 0 energy.
force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
force.setReactionFieldDielectric(1.0) 
system.addForce(force)

platform = mm.Platform.getPlatformByName('Reference')
context = mm.Context(system, mm.VerletIntegrator(0), platform)

positions = [(28.089, 41.809, -2.480), (30.005, 43.060, -4.204)]
context.setPositions(positions)

state = context.getState(getEnergy=True)
print(state.getPotentialEnergy())

[proteneer]

@JoaoRodrigues I'm literally debugging this right now and ran into the same behavior. I'm +1 on disabling the rf correction.

[proteneer]

What's more insane is that we do turn off the crf offset when computing the 1-4 exceptions

[peastman]

What's more insane is that we do turn off the crf offset when computing the 1-4 exceptions

It's not just the offset. Reaction field isn't applied to 1-4 interactions. Nor are cutoffs, nor periodic boundary conditions. That's because they're fundamentally bonded interactions, not nonbonded. Likewise, if you use PME, the exception is only applied to the two particular bonded atoms. The interactions of each one with all the infinite periodic images of the other are computed with the standard parameter values.