Switching function in LAMMPS export

[mattwthompson]

See #158 (comment)

[mrshirts]

If there are issues in the nonbonded matching, it's probably due to differences in the way they are handled. Is there a link to the input files/python and lammps files to look at?

[mrshirts]

OK, tracing back, I would agree that the switches are an issue.

Have you considered just using an abrupt cutoff to validate conversion? InterMol should have the example files that implement these: for example for GROMACS, https://github.com/shirtsgroup/InterMol/blob/master/intermol/tests/gromacs/grompp_vacuum.mdp, and for LAMMPS, https://github.com/shirtsgroup/InterMol/blob/master/intermol/tests/lammps/unit_tests/atom_style-full_vacuum/atom_style-full-data_vacuum.input

[mattwthompson]

Using a hard cutoff would be okay in validating that the per-particle parameters are properly set in each export, but it would not fully validate the exports from SMIRNOFF force fields since they specify that a switching function should be used.

Of course, there is confusion around what actually should be used, something that still needs input in order to be resolved in a clear way.

[mrshirts]

but it would not fully validate the exports from SMIRNOFF force fields since they specify that a switching function should be used.

A specific switching function or any switching function? Any switching function is ambiguous and a specific switching function may be unimplemented anywhere but one package (for any choice of switching function) - they generally differ a fair amount.

I assume this is for LJ? Ewald treatments can have larger differences (though I think mostly in the energy offset, not as much in the forces.

I see from the linked posts that this is definitely still being discussed. I would argue (and I can expand on it in openforcefield/openff-toolkit#896 as needed), that there are two separate components that are somewhat orthogonal - what the parameters are, and how the long range interactions are approximated. And trying to validate the combination simultaneously is almost impossible, given that there is no consensus on the right way to do the long range interaction approximated (though LJ or P3M for both LJ and Coulomb might be creeping towards being the "right" way).

[mattwthompson]

A specific switching function or any switching function? Any switching function is ambiguous and a specific switching function may be unimplemented anywhere but one package (for any choice of switching function) - they generally differ a fair amount.

This is the heart of the issue, before working out the implementation details of how each engine handles it. It's hard to get more information out of a force field than "apply a switching function starting at this distance," and your criticism of the ambiguity is spot-on. This is one of the many details of the SMIRNOFF spec that needs improvement. It's a part of the contents of a force field and therefore must be considered.

A good path forward from the current situation was suggested by @j-wags [1], and @SimonBoothroyd also suggested a solution in which the spec unambiguously specifies the form of the switching function that should be used. This would not solve the issue of different engines implementing different switching functions (no matter how different) but it would fix the current issue that OpenMM's implementation is silently and implicitly considered the default. I don't recall where this was posted, possibly just in meeting notes somewhere, but it should be considered as a spec change.

@jchodera in that thread also suggested that thread that the goal should just be making sure each export uses some sort of switching function, which would ensure continuity of the derivatives and should get good-but-imperfect agreement between engines' implementations. This bends the rules of what role a specification serves, but is an appealing option for implementation since it's something that can actually be done - compared to i.e. making sure 4 different engines use precisely the same polynomial.

I assume this is for LJ?

The context here is LJ/vdW only, yes, but there are also parallel discussions [2] around how electrostatic methods should be specified in general, including whether or not they should at all. Our mainline force fields use <Electrostatics ... switch_width="0.0 * angstrom" ...> which dodges this bullet for now.

The way the toolkit handles non-bonded details is generally OpenMM-centric [3] and I'm trying to avoid that wherever possible, since it leads to a whole host of ambiguities that are difficult to implement [4-6] - including this one.

1. Fix #882 openff-toolkit#896 (comment)
2. SMIRNOFF: Specifying "PME" as a non-bonded method is somewhat ambiguous standards#15
3. How should we expose keywords and implement logic to handle simulation periodicity and phase? openff-toolkit#219 (comment)
4. Cutoff Not Set for Topology Without Box Vectors openff-toolkit#735
5. ForceField with no Electrostatics tag can be used to make a system openff-toolkit#716
6. Parametrizing non-periodic topologies with mainline force fields changes non-bonded settings openff-toolkit#1084

[mrshirts]

@jchodera in that thread also suggested that thread that the goal should just be making sure each export uses some sort of switching function, which would ensure continuity of the derivatives and should get good-but-imperfect agreement between engines' implementations.

But would this include LJ-vdw, which is demonstrably more physical than an isotropic switching function (see, for example, the energy at the liquid/vacuum interface).

I would still argue for separating the specification: specification of parameters and specification of non-bonded treatment, such that the first can be validated without the second, which are essentially orthogonal. My playing around with InterMol seems to suggest that the 2nd part is much, much harder and may be impossible to do precisely. Whereas it's more straightforward to make sure that different engines are using compatible sigmas, epsilons, and partial charges (pick a single fixed cutoff, if energies matches, the parameters are consistent, and the total nonbonded energy WOULD be right IF the same method was used across implementations for engines).

I'll see if I can dive into the above issues later today, though if there is something that is approaching a decision point, that is where I'm happy to discuss things.

[mattwthompson]

This sounds like something that could be a proposal to update the SMIRNOFF spec (using OFF-EP 0), though I'm not quite sure how they'd go about splitting it up into two specs.

[mattwthompson]

It looks like this is implemented in a fairly esoteric pair potential ... it may be possible to force this through by setting the Gaussian widths to zero, but that probably produces more problems (computational cost, possibly numerical stability) than it solves (only the use of a - not necessarily the correct - switching function).

A band-aid which I plan to use for the next set of validation tests is to just turn this off as part of #626. This helps for our testing but doesn't really solve the problem of LAMMPS not clearly supporting a switching function that mainline OpenFF force fields use.

[mrshirts]

I think that we can only interchange what can be interchanged, and just report when it's not possible!

[mattwthompson]

Agreed. I came up with three paths forward here:

1. Implement a switching function in LAMMPS
2. Error out if asked to write LAMMPS files if the Interchange object is meant to be used with a switching function
3. Remove the use of a switching function when writing to LAMMPS, loudly warning the user

The only workable option here is 3. I also ran this by @j-wags who agreed.

[mrshirts]

One other possibility is switch to LJ/long styles (with loud warnings), as that is the way things are going (better supports surface tension, etc).

[mattwthompson]

That looks interesting

If flag_lj is set to long, no cutoff is used on the LJ 1/r^6 dispersion term. The long-range portion can be calculated by using the kspace_style ewald/disp or pppm/disp commands. The specified LJ cutoff then determines which portion of the LJ interactions are computed directly by the pair potential versus which part is computed in reciprocal space via the Kspace style. If flag_lj is set to cut, the LJ interactions are simply cutoff, as with pair_style lj/cut.

Not having read In‘t Veld's paper I'm not sure if this is more or less of a deviation than excluding the switching function altogether. I wonder if this is something the science team would consider for the mainline force fields, or if it's better to just make the jump to full LJPME.

[mrshirts]

pppm is essentially equal to LJ-PME - just a different mathematical formalism to get the infinite sum. I tested the rough numerical equivalency (i.e. the differences between pppm and different flavors of PME) with InterMol.

[mattwthompson]

Surely using LJPME (or a near equivalent) with parameters fit to cutoff + switching function produces more of a deviation than just not using the switching function with the same parameters?