search

openforcefield / openff-interchange

A project (and object) for storing, manipulating, and converting molecular mechanics data.
https://docs.openforcefield.org/projects/interchange
MIT License
71 stars 23 forks source link

Incorrect nonbonded energies from Amber if improper torsions are present #686

Open mattwthompson opened 1 year ago

mattwthompson commented 1 year ago

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
mattwthompson commented 1 year ago

Something's pretty significantly off here.

image
mattwthompson commented 1 year ago

688 is related, I think

dscerutti commented 1 year ago

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.

mattwthompson commented 1 year ago

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.