Test intramolecular interactions in polycyclic systems with Amber

[mattwthompson]

Description

#688

Checklist

• Add tests
• Lint
• Update docstrings

[codecov]

Codecov Report

Merging #900 (de61ec3) into main (6b9f8d2) will decrease coverage by 0.43%.
Report is 2 commits behind head on main.
The diff coverage is n/a.

Additional details and impacted files

[mattwthompson]

Working through a funky ring-like molecule, trying to keep hydrogens off to make the files almost human-legible:

In:

ForceField("openff-2.1.0.offxml").create_interchange(
    Molecule.from_smiles(
        "C1#CC#CC#C1",
    ).to_topology()
).to_prmtop("tmp")

Out:
https://gist.github.com/mattwthompson/fcba880ea1dab43bc7d721f8f3a29eb2

Just the relevant bit:

%FLAG NUMBER_EXCLUDED_ATOMS
%FORMAT(10I8)
       5       4       3       2       1       1
%FLAG DIHEDRALS_WITHOUT_HYDROGEN
%FORMAT(10I8)
       0       3       6       9       1       0      15     -12       9       2
       3       0      15      12       1       3       6      -9      12       2
      15       0       3       6       2       6       9     -12      15       1
%FLAG EXCLUDED_ATOMS_LIST
%FORMAT(10I8)
       2       3       4       5       6       3       4       5       6       4
       5       6       5       6       6       0

In the dihedrals section, there are a total of 6 dihedrals, covering the forward direction 0-1-2-3, 1-2-3-4, 2-3-4-5 and the backwards direction 0-5-4-3, 1-0-5-4, 5-0-1-2. The 1-4 pairs are 0-3, 1-4, 2-5; the second time a dihedral is listed, the third atom has a negative sign on it. (In this case no third atoms are index 0, conveniently or by rearrangement, I'm not sure). This all seems accurate with respect to what's documented. (To get from the internal to Amber index in this section, just multiply by three, these are indices into an array of dihedral parameters.)

Working through the exclusions, now: all intramolecular pairs should have their interactions zeroed out (1-2, 1-3) or scaled (1-4). Nothing in this bond graph is farther than three nodes away when considering minimum distances. If NUMBER_EXCLUDED_ATOMS is meant to be trusted, atom 0 has 5 exclusions not previously listed, atom 1 has 4, etc. until atoms 4 and 5 have 1 each that weren't already listed. Two confusing caveats here are that an atom whose exclusions are sufficiently specified earlier in the list is given an "number" of exclusions of 1 and only atom 0 (which does not exist) is listed. Anyway, indexing into EXCLUDED_ATOMS_LIST (to go from internal to Amber indices here, just add 1, since Amber starts counting atoms at 1 for some reason) this appears to give the following exclusions:

Atom	Number of excluded atoms	Excluded atoms (not appearing earlier)
0	5	1 2 3 4 5
1	4	2 3 4 5
2	3	3 4 5
3	2	4 5
4	1	5
5	"1"	N/A (see above)

Note that Amber only lists i-j exclusions once, i.e. j will be present as an exclusion for i but i need not be listed as an exclusion for j.

That all seems fine to me. For good measure, I ran it through OpenMM to ParmEd

In [1]: from openff.toolkit import *

In [2]: import parmed

In [3]: topology = Molecule.from_smiles(        "C1#CC#CC#C1",).to_topology()

In [4]: !rm parmed.prmtop

In [5]: parmed.openmm.load_topology(
   ...:     topology=topology.to_openmm(),
   ...:     system=ForceField("openff-2.1.0.offxml").create_openmm_system(topology),
   ...: ).save("parmed.prmtop")

and got a file that treats these sections identically: https://gist.github.com/mattwthompson/4cc00f9257a53a261322fcff416102f8

%FLAG NUMBER_EXCLUDED_ATOMS
%FORMAT(10I8)
       5       4       3       2       1       1
%FLAG DIHEDRALS_WITHOUT_HYDROGEN
%FORMAT(10I8)
       0       3       6       9       2       0      15     -12       9       1
       3       0      15      12       3       3       6      -9      12       1
      15       0       3       6       4       6       9     -12      15       1
%FLAG EXCLUDED_ATOMS_LIST
%FORMAT(10I8)
       2       3       4       5       6       3       4       5       6       4
       5       6       5       6       6       0

Something's quite wrong, still. Running a larger molecule through the machinery shows that vdW energies differ with larger polycyclic rings, something like 0.2 kJ/mol/ring. This would explain why @Yoshanuikabundi is seeing horrible results with NADP

In [11]: smiles = "C1=CC=C2C(=C1)C=CC3=C2C=CC4=CC=CC=C43"

In [12]:
    ...: molecule = MoleculeWithConformer.from_smiles(smiles)
    ...: topology = molecule.to_topology()
    ...: topology.box_vectors = Quantity([10, 10, 10], "nanometer")

In [13]: get_summary_data(
    ...:     sage_unconstrained.create_interchange(topology),
    ...:     _engines=["OpenMM", "Amber", "GROMACS"],
    ...: )
['\n', '\n', '   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER\n', '      1       3.1456E+01     5.0509E+00     2.2557E+01        C       17\n', '\n', ' BOND    =        3.5399  ANGLE   =        3.8041  DIHED      =        0.0001\n', ' VDWAALS =       -0.6553  EEL     =       -6.5118  HBOND      =        0.0000\n', ' 1-4 VDW =       20.6784  1-4 EEL =       10.6007  RESTRAINT  =        0.0000\n']
Out[13]:
              Bond      Angle   Torsion  Electrostatics        vdW  RBTorsion
OpenMM   14.810762  15.916498  0.000281       17.114335  83.155060        NaN
Amber    14.810942  15.916354  0.000418       17.107958  83.776650        NaN
GROMACS  14.811943  15.916394  0.000277       17.092522  83.155503        0.0

(The extra output in the middle is a raw print of the mdinfo file, just to verify that I'm not accidentally missing or double-counting the 1-4 contributions; the difference between engines is much less than what's reported as the 1-4 contribution).

As another check I thought it's be useful to look at a long alkane, maybe isolating if this is a problem with rings or not. It appears to just be a problem everywhere:

In [6]:
   ...:     molecule = MoleculeWithConformer.from_smiles(10 * "CC")
   ...:     topology = molecule.to_topology()
   ...:     topology.box_vectors = Quantity([10, 10, 10], "nanometer")

In [7]: get_summary_data(sage_unconstrained.create_interchange(topology))
['\n', '\n', '   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER\n', '      1       3.9987E+01     5.9529E+00     2.1936E+01        C       19\n', '\n', ' BOND    =        0.0530  ANGLE   =        5.4211  DIHED      =       29.9423\n', ' VDWAALS =       -2.5087  EEL     =        4.1413  HBOND      =        0.0000\n', ' 1-4 VDW =        5.7720  1-4 EEL =       -2.8337  RESTRAINT  =        0.0000\n']
Out[7]:
             Bond      Angle     Torsion  Electrostatics        vdW  RBTorsion
OpenMM   0.221711  22.681757  125.278745        5.470604  12.886736        NaN
Amber    0.221752  22.681882  125.278583        5.470998  13.653647        NaN
GROMACS  0.221702  22.683037  125.278786        5.458763  12.887473        0.0
LAMMPS   2.800538  22.681743  125.278742        5.556665  12.810040        NaN

[mattwthompson]

Here's a representative input file (not sure the jargon, it's what I pass through to sander -i):

single-point energy
&cntrl
imin=1,
maxcyc=0,
ntb=1,
fswitch=8.0,
ntc=1,
ntf=2,
cut=9.0,
/
&ewald
order=4
skinnb=1.0
//

[mattwthompson]

A helpful commenter points out that Amber uses the same exclusion policy for electrostatics and vdW interactions; so, given that the electrostatics energies match decently well, that's probably not the culprit.

On a whim I decided to see what happens when I turn off the switching function and, well, that seems to fix it. So #470 lives, and the switching function needs to be looked into more.