LAMMPS reports slightly different torsion energies

[mattwthompson]

          After #897 #914 I'm getting the following results:

              Bond       Angle     Torsion  Electrostatics        vdW  RBTorsion
OpenMM   20.836174  310.672706  896.221085    -3244.337363  14.406467        NaN
Amber    20.836320  310.672878  896.221168    -3244.227158  14.088783        NaN
GROMACS  20.836302  310.673157  896.221071    -3244.741211  14.430883        0.0
LAMMPS   20.836172  310.672698  896.334673    -3244.337702  13.744850        NaN

My two concerns are the slight difference in LAMMPS torsions and vdW energies (which should match Amber, since neither appear to have a SMIRNOFF-complatible switching function implemented, and so just have a hard cut-off in these code paths).

cc: @timbernat particularly the torsions - 0.113 / 896 kJ/mol isn't much, but I can't ignore it either

Originally posted by @mattwthompson in https://github.com/openforcefield/openff-interchange/issues/894#issuecomment-1976971873

[timbernat]

@mattwthompson Let me get back to you on this, I have ~1,300 polymer structures that I've successfully exported to both LAMMPS and OpenMM using Interchange 0.3.29 for a collaboration. Once I switch my workflow over to the latest NAGL and sort out some weirdness w/ OpenMM integrating out 1 step (KE for some reason isn't 0 for OpenMM), I can get an energy evaluation on the same set of structures to give a much more comprehensive benchmark of the scope of the issue (albeit without the GMX and Amber comparison)

[mattwthompson]

Wonderful, much appreciated

[timbernat]

A crucial detail I've run into while getting this ready; OpenMM's VerletIntegrator (and in general, any leapfrog algorithm integrator in OpenMM) will report non-zero kinetic energies when evaluating State, even if all velocities are zero. Taken straight from the horse's mouth: "Note that [KE] may be different from simply mv/2 summed over all particles. For example, a leapfrog integrator will store velocities offset by half a step, so they must be adjusted before computing the kinetic energy."

Some quick testing shows that for static structures exported from Interchange to OpenMM, VerletIntegrator and LangevinIntegrator both yield nonzero KEs, while LangevinMiddleIntegrator and BrownianIntegrator DO yield the expected 0 KE. I haven't yet checked the full list of supported Integrators exhaustively.

This becomes relevant to Interchange when noting that the drivers.openmm layer does use VerletIntegrator in its evaluation, and does not check 0 KE as far as I can tell. This has been introducing errors into my current data, and is likely doing the same for yours @mattwthompson

[mrshirts]

How important is it to match the kinetic energy? Usually when checking for things being equal, I just look at the potential energy.

[timbernat]

Well because that kinetic energy represents a piece of the potential energy accounted for in the wrong place; a small fraction of the potential is incorrectly reported as belonging to kinetic instead

[mrshirts]

a small fraction of the potential is incorrectly reported as belonging to kinetic instead

Ah, yea, that's a problem.

[mattwthompson]

Excellent find @timbernat! I'm splitting that direction of work off to #1096. It seems like it's worth changing but not not likely to be the reason torsions are coming out a little different