openff-interchange
openff-interchange copied to clipboard
openforcefield
Incorrect nonbonded energies from Amber if improper torsions are present
Something is wrong with the interplay between setting the third atom negative (its 1-4 interaction was already accounted for) and the fourth atom negative (this torsion is an improper). It makes no sense that they'd ever both be negative but I also can't reason through what it should be.
Based on b29ed479acfbfc0df7cd733304b3807f67ae0998 in #685:
In [1]: import numpy
...: import parmed
...: from openff.toolkit import ForceField, Molecule
...: from openff.units import unit
...:
...: from openff.interchange.drivers.all import get_summary_data
...:
...: molecule = Molecule.from_iupac("imidazole")
...: molecule.generate_conformers(n_conformers=1)
...:
...: force_field = ForceField("openff_unconstrained-2.0.0.offxml")
...:
...: interchange = force_field.create_interchange(molecule.to_topology())
...:
...: interchange.box = [5, 5, 5]
...:
...: print(get_summary_data(interchange))
...:
...: interchange.collections.pop("ImproperTorsions")
...:
...: print(get_summary_data(interchange))
Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.
/Users/mattthompson/mambaforge/envs/openff-interchange-env/lib/python3.10/site-packages/jupyter_client/connect.py:20: DeprecationWarning: Jupyter is migrating its paths to use standard platformdirs
given by the platformdirs library. To remove this warning and
see the appropriate new directories, set the environment variable
`JUPYTER_PLATFORM_DIRS=1` and then run `jupyter --paths`.
The use of platformdirs will be the default in `jupyter_core` v6
from jupyter_core.paths import jupyter_data_dir, jupyter_runtime_dir, secure_write
Bond Angle Torsion Nonbonded vdW Electrostatics
OpenMM 18.619643 118.574111 0.000081 -13.240277 NaN NaN
Amber 18.619637 118.574142 0.000000 NaN 117810.943599 -50.193356
GROMACS 18.620449 118.574234 0.000079 NaN -0.594816 -12.712922
LAMMPS 18.619652 118.574104 0.000081 NaN -0.598532 -12.741790
Bond Angle Torsion Nonbonded vdW Electrostatics
OpenMM 18.619643 118.574111 0.000069 -13.240277 NaN NaN
Amber 18.619637 118.574142 0.000000 NaN -0.594546 -12.638190
GROMACS 18.620449 118.574234 0.000067 NaN -0.594816 -12.712922
LAMMPS 18.619652 118.574104 0.000069 NaN -0.598532 -12.741790
The file (prior to popping the impropers):
%VERSION VERSION_STAMP = V0001.000 DATE = 04/27/23 14:52:31
%FLAG TITLE
%FORMAT(20a4)
%FLAG POINTERS
%FORMAT(10I8)
9 5 4 5 8 5 23 5 0 0
35 1 5 5 5 6 5 8 1 0
0 0 0 0 0 0 0 1 9 0
0 0
%FLAG ATOM_NAME
%FORMAT(20a4)
N1 C1 N2 C2 C3 H H H H
%FLAG CHARGE
%FORMAT(5E16.8)
-5.88270499E+00 6.96219421E+00 -1.21333182E+01 5.29740484E+00 -4.75602047E+00
5.45921867E+00 1.08295124E+00 7.67705477E-01 3.20256923E+00
%FLAG ATOMIC_NUMBER
%FORMAT(10I8)
7 6 7 6 6 1 1 1 1
%FLAG MASS
%FORMAT(5E16.8)
1.40067200E+01 1.20107800E+01 1.40067200E+01 1.20107800E+01 1.20107800E+01
1.00794700E+00 1.00794700E+00 1.00794700E+00 1.00794700E+00
%FLAG ATOM_TYPE_INDEX
%FORMAT(10I8)
1 2 1 2 2 3 4 5 5
%FLAG NUMBER_EXCLUDED_ATOMS
%FORMAT(10I8)
8 7 6 5 4 2 1 1 1
%FLAG NONBONDED_PARM_INDEX
%FORMAT(10I8)
1 2 4 7 11 2 3 5 8 12
4 5 6 9 13 7 8 9 10 14
11 12 13 14 15
%FLAG RESIDUE_LABEL
%FORMAT(20a4)
RES
%FLAG RESIDUE_POINTER
%FORMAT(10I8)
1
%FLAG BOND_FORCE_CONSTANT
%FORMAT(5E16.8)
3.90035547E+02 5.05144496E+02 5.05481966E+02 3.97254579E+02 3.82545017E+02
3.99159295E+02
%FLAG BOND_EQUIL_VALUE
%FORMAT(5E16.8)
1.39016055E+00 1.01948187E+00 1.30482618E+00 1.08535850E+00 1.37924439E+00
1.37168856E+00
%FLAG ANGLE_FORCE_CONSTANT
%FORMAT(5E16.8)
5.78982894E+01 3.96281169E+01 5.54218215E+01 6.34442260E+01 5.50788255E+01
%FLAG ANGLE_EQUIL_VALUE
%FORMAT(5E16.8)
2.21318644E+00 2.26548502E+00 2.04808379E+00 2.01597978E+00 1.96997758E+00
%FLAG DIHEDRAL_FORCE_CONSTANT
%FORMAT(5E16.8)
6.73676248E+00 5.29484300E+00 1.57666787E+00 1.55051542E+00 1.76281822E+00
3.50000000E+00 3.66666667E-01 3.50000000E+00
%FLAG DIHEDRAL_PERIODICITY
%FORMAT(5E16.8)
2.00000000E+00 2.00000000E+00 2.00000000E+00 2.00000000E+00 2.00000000E+00
2.00000000E+00 2.00000000E+00 2.00000000E+00
%FLAG DIHEDRAL_PHASE
%FORMAT(5E16.8)
3.14159265E+00 3.14159265E+00 3.14159265E+00 3.14159265E+00 3.14159265E+00
3.14159265E+00 3.14159265E+00 3.14159265E+00
%FLAG SCEE_SCALE_FACTOR
%FORMAT(5E16.8)
1.20000000E+00 1.20000000E+00 1.20000000E+00 1.20000000E+00 1.20000000E+00
1.20000000E+00 1.20000000E+00 1.20000000E+00
%FLAG SCNB_SCALE_FACTOR
%FORMAT(5E16.8)
2.00000000E+00 2.00000000E+00 2.00000000E+00 2.00000000E+00 2.00000000E+00
2.00000000E+00 2.00000000E+00 2.00000000E+00
%FLAG SOLTY
%FORMAT(5E16.8)
0.00000000E+00
%FLAG LENNARD_JONES_ACOEF
%FORMAT(5E16.8)
7.93518630E+05 9.43185634E+05 1.09874263E+06 1.95217207E+03 2.94202507E+03
1.83621096E-01 5.13478436E+04 6.52189376E+04 5.56123804E+01 2.66090756E+03
4.95412655E+04 6.28522004E+04 5.44509573E+01 2.57728328E+03 2.49611586E+03
%FLAG LENNARD_JONES_BCOEF
%FORMAT(5E16.8)
7.29565189E+02 6.74816129E+02 6.17925587E+02 1.94827762E+01 2.02916006E+01
1.01732411E-01 1.01149923E+02 9.67147330E+01 1.79224117E+00 1.25498420E+01
9.63867036E+01 9.21075325E+01 1.72045213E+00 1.19821174E+01 1.14396830E+01
%FLAG BONDS_INC_HYDROGEN
%FORMAT(10I8)
0 15 2 3 18 4 9 21 4 12
24 4
%FLAG BONDS_WITHOUT_HYDROGEN
%FORMAT(10I8)
0 3 1 0 12 1 3 6 3 6
9 5 9 12 6
%FLAG ANGLES_INC_HYDROGEN
%FORMAT(10I8)
0 3 18 2 0 12 24 2 3 0
15 4 6 3 18 2 6 9 21 2
9 12 24 2 12 0 15 4 12 9
21 2
%FLAG ANGLES_WITHOUT_HYDROGEN
%FORMAT(10I8)
0 3 6 1 0 12 9 1 3 0
12 3 3 6 9 5 6 9 12 1
%FLAG DIHEDRALS_INC_HYDROGEN
%FORMAT(10I8)
0 12 9 21 2 3 0 12 24 3
3 6 9 21 5 15 0 3 6 3
6 9 12 24 2 9 6 3 18 1
15 0 12 9 3 12 0 3 18 3
15 0 3 18 3 15 0 12 24 3
21 9 12 24 2 3 0 6 -18 6
3 6 18 0 6 3 18 0 -6 6
12 0 9 -24 7 12 9 24 0 7
12 24 0 -9 7 0 3 12 -15 8
0 12 15 -3 8 0 15 3 -12 8
9 6 12 -21 7 9 12 21 -6 7
9 21 6 -12 7
%FLAG DIHEDRALS_WITHOUT_HYDROGEN
%FORMAT(10I8)
0 3 -6 9 1 0 12 -9 6 2
3 0 -12 9 3 3 6 -9 12 4
12 0 -3 6 3
%FLAG EXCLUDED_ATOMS_LIST
%FORMAT(10I8)
2 3 4 5 6 7 8 9 3 4
5 6 7 8 9 4 5 6 7 8
9 5 6 7 8 9 6 7 8 9
7 9 0 9 0
%FLAG HBOND_ACOEF
%FORMAT(5E16.8)
%FLAG HBOND_BCOEF
%FORMAT(5E16.8)
%FLAG HBCUT
%FORMAT(5E16.8)
%FLAG AMBER_ATOM_TYPE
%FORMAT(20a4)
N1 C1 N2 C2 C3 H1 H2 H3 H4
%FLAG TREE_CHAIN_CLASSIFICATION
%FORMAT(20a4)
BLA BLA BLA BLA BLA BLA BLA BLA BLA
%FLAG JOIN_ARRAY
%FORMAT(10I8)
0 0 0 0 0 0 0 0 0
%FLAG IROTAT
%FORMAT(10I8)
0 0 0 0 0 0 0 0 0
%FLAG SOLVENT_POINTERS
%FORMAT(3I8)
1 1 2
%FLAG ATOMS_PER_MOLECULE
%FORMAT(10I8)
9
%FLAG BOX_DIMENSIONS
%FORMAT(5E16.8)
9.00000000E+01 5.00000000E+01 5.00000000E+01 5.00000000E+01
%FLAG RADIUS_SET
%FORMAT(1a80)
0
%FLAG RADII
%FORMAT(5E16.8)
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
%FLAG SCREEN
%FORMAT(5E16.8)
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
%FLAG IPOL
%FORMAT(1I8)
0
Can you reprint what the Amber mdout provides? Amber differentiates between 1:4 vdW / electrostatics and "the rest" of nonbonded vdW / electrostatics. From the energies you reported, there is probably an interaction that should be 1:2 excluded masquerading as a 1:4 scaled non-bonded pair.
I can't reproduce this anymore, but I'm pretty sure I've seen similar behavior in other systems. I'll try to remember to repost them here.