search

openforcefield / openff-interchange

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

Angle energy mismatch between GROMACS and Amber #249

Open mattwthompson opened 3 years ago

mattwthompson commented 3 years ago

I'm getting different angle energies for a simple single-molecule system.

Reproduction

Unzip these files: repro.zip

In: 

$ python repro.py
$ gmx grompp -f auto_generated.mdp -c internal.gro -p internal.top -o internal.tpr
$ gmx mdrun -s internal.tpr
$ echo 2 | gmx energy -f ener.edr
$ sander -i auto_generated.in -c parmed.inpcrd -p parmed.prmtop -O
$ cat mdinfo

Out (just the end):

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Angle                    7.04423e-05         --          0          0  (kJ/mol)

GROMACS reminds you: "And You Will Know That My Name is the Lord When I Lay My Vengeance Upon Thee." (Pulp Fiction)

$ cat mdinfo

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1       8.0523E-01     4.6024E+00     8.4576E+00     C2          3

 BOND    =        0.0991  ANGLE   =        0.2822  DIHED      =        0.0000
 VDWAALS =        0.0000  EEL     =       -0.0001  HBOND      =        0.0000
 1-4 VDW =       -0.0715  1-4 EEL =        0.4955  RESTRAINT  =        0.0000

GROMACS reports 7.04423e-05 kJ/mol for the angle component and Amber reports 0.2822 kcal/mol (I think). Other values (bond, 1-4 vdW, 1-4 Coulomb) match closel. The bond component matching closely indicates to me that the positions (and bond parameters) are being processed the same by each engine. In other cases, I've found that angle energies match but bond energies are wildly off, which strongly indicated that positions were not being written correctly.

Other troubleshooting I've done

mattwthompson commented 3 years ago

Okay, Amber defines harmonic bonds without the 1/2 in front; I was confused because (3) here does have the 1/2.

j-wags commented 3 years ago

Huh, I'm reproducing this using a code path that doesn't use Interchange at all (just openff-toolkit, openmm, and parmed)

from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField
import parmed
mol = Molecule.from_smiles('ClC#CCl')
mol.generate_unique_atom_names()
mol.generate_conformers()

ff = ForceField('openff_unconstrained-1.0.0.offxml')
sys = ff.create_openmm_system(mol.to_topology())
omm_top = mol.to_topology().to_openmm()
print(dir(omm_top))
print(dir([*omm_top.residues()][0]))
[*omm_top.residues()][0].name = 'ABC'
structure = parmed.openmm.load_topology(mol.to_topology().to_openmm(),
                                                   sys,
                                                   xyz=mol.conformers[0],
                                                   #condense_atom_types=False
                                                  )
structure.residues[0].name = 'XYZ'
structure.save('ligand.prmtop', overwrite=True)
structure.save('ligand.inpcrd', overwrite=True)

!echo "parm ligand.prmtop" > parmed.in
!echo "loadCoordinates ligand.inpcrd" >> parmed.in
!echo "energy" >> parmed.in
!echo "quit" >> parmed.in
!parmed -i parmed.in
Adding prmtop parmed.prmtop to parm list. parmed.prmtop is the active parm.
Adding coordinates to parmed.prmtop from parmed.inpcrd
Computing a single-point energy for parmed.prmtop
Bond     =            0.0991069     Angle    =            0.2821805
Dihedral =            0.0000000     1-4 vdW  =           -0.0715118
1-4 Elec =            0.4955206     vdWaals  =           -0.0015985
Elec.    =           -0.0001400
TOTAL    =            0.8035579
Done!

@swails, do you have any ideas here? The basic problem is that we're seeing a nonzero angle energy for a linear molecule. Here's the prmtop and inpcrd that come out of my code above. parmed.inpcrd.gz parmed.prmtop.gz

The molecule + geometry:

Screen Shot 2021-07-21 at 3 22 16 PM

dichloroethyne.sdf.gz

mattwthompson commented 3 years ago

We didn't get to the bottom of the original issue, but stripping out some complexity and dropping in a methane does not reproduce this issue.

In:

import numpy as np
import parmed as pmd
from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField

from openff.interchange.components.interchange import Interchange

molecule = Molecule.from_smiles("C")
molecule.generate_conformers(n_conformers=1)
molecule.name = "FOO"
topology = molecule.to_topology()
parsley = ForceField("openff-1.0.0.offxml")
out = Interchange.from_smirnoff(force_field=parsley, topology=topology)
out.box = [4, 4, 4] * np.eye(3)
out.positions = np.random.rand(15).reshape((5, 3))  # Repalce with valid data

out.to_gro("internal.gro")
out.to_top("internal.top")

gro_file = pmd.load_file("internal.gro")
top_file = pmd.load_file("internal.top")

gro_file.save("parmed.inpcrd")
top_file.save("parmed.prmtop")

"""
gmx grompp -f auto_generated.mdp -c internal.gro -p internal.top -o internal.tpr -maxwarn 1
gmx mdrun -s internal.tpr
echo 2 | gmx energy -f ener.edr
"""

"""
sander -i auto_generated.in -c parmed.inpcrd -p parmed.prmtop -O
"""

"""
cat mdinfo
"""

Out:

<snipping most of the output>

Back Off! I just backed up energy.xvg to ./#energy.xvg.7#
Last energy frame read 0 time    0.000

Statistics over 1 steps [ 0.0000 through 0.0000 ps ], 1 data sets
All statistics are over 1 points

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Angle                       1097.64         --          0          0  (kJ/mol)

GROMACS reminds you: "Wedged as we are between two eternities of idleness, there is no excuse for being idle now" (Anthony Burgess)

(openff-system) [repro] cat mdinfo                                                                                                                                        dc855e1  ✱

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1       1.2166E+05     6.7073E+03     1.8632E+04     C1          1

 BOND    =   121401.0796  ANGLE   =      262.3411  DIHED      =        0.0000
 VDWAALS =        0.0000  EEL     =       -0.0081  HBOND      =        0.0000
 1-4 VDW =        0.0000  1-4 EEL =        0.0000  RESTRAINT  =        0.0000

262.3411 kcal = 1097.64 kJ/mol (to reported precision)

mattwthompson commented 2 years ago

This is hard to reproduce, but work since this issue was opened (#247 being the main thread) has removed ParmEd as an intermediary for doing these comparisons. It's hard for me to go back and sift through all of this to diagnose the possibility that some part of OpenFF -> OpenMM -> ParmEd -> Engines did not work. Closing for now

mattwthompson commented 1 year ago

I can still reproduce this ... is there some gotcha with linear angles?

molecule = Molecule.from_smiles("ClC#CCl")
molecule.generate_conformers(n_conformers=1)
topology = molecule.to_topology()
topology.box_vectors = [4, 4, 4] * unit.nanometer
parsley = ForceField("openff-1.0.0.offxml")
out = Interchange.from_smirnoff(force_field=parsley, topology=topology)

get_summary_data(out)
             Bond     Angle  Torsion  Electrostatics       vdW  RBTorsion
OpenMM   0.416595  0.000066      0.0        2.073311 -0.304784        NaN
Amber    0.416726  1.180725      0.0        2.072754 -0.299156        NaN
GROMACS  0.416672  0.000070      0.0        2.067999 -0.299192        0.0
LAMMPS   0.416596  0.000066      0.0        2.056103 -0.308576        NaN

Here's the file:


%VERSION  VERSION_STAMP = V0001.000  DATE = 05/09/23  14:35:19
%FLAG TITLE
%FORMAT(20a4)

%FLAG POINTERS
%FORMAT(10I8)
       4       2       0       3       0       2       0       1       0       0
       7       1       3       2       1       2       1       1       1       0
       0       0       0       0       0       0       0       1       4       0
       0       0
%FLAG ATOM_NAME
%FORMAT(20a4)
  Cl   C   C  Cl
%FLAG CHARGE
%FORMAT(5E16.8)
  1.67280712E+00 -1.67298938E+00 -1.67298938E+00  1.67317165E+00
%FLAG ATOMIC_NUMBER
%FORMAT(10I8)
      17       6       6      17
%FLAG MASS
%FORMAT(5E16.8)
  3.54532000E+01  1.20107800E+01  1.20107800E+01  3.54532000E+01
%FLAG ATOM_TYPE_INDEX
%FORMAT(10I8)
       1       2       2       1
%FLAG NUMBER_EXCLUDED_ATOMS
%FORMAT(10I8)
       3       2       1       1
%FLAG NONBONDED_PARM_INDEX
%FORMAT(10I8)
       1       2       2       3
%FLAG RESIDUE_LABEL
%FORMAT(20a4)
RES
%FLAG RESIDUE_POINTER
%FORMAT(10I8)
       1
%FLAG BOND_FORCE_CONSTANT
%FORMAT(5E16.8)
  1.98830634E+02  3.64357562E+02
%FLAG BOND_EQUIL_VALUE
%FORMAT(5E16.8)
  1.75561946E+00  1.21628943E+00
%FLAG ANGLE_FORCE_CONSTANT
%FORMAT(5E16.8)
  7.05333893E+01
%FLAG ANGLE_EQUIL_VALUE
%FORMAT(5E16.8)
  3.14159265E+00
%FLAG DIHEDRAL_FORCE_CONSTANT
%FORMAT(5E16.8)
  0.00000000E+00
%FLAG DIHEDRAL_PERIODICITY
%FORMAT(5E16.8)
  1.00000000E+00
%FLAG DIHEDRAL_PHASE
%FORMAT(5E16.8)
  0.00000000E+00
%FLAG SCEE_SCALE_FACTOR
%FORMAT(5E16.8)
  1.20000000E+00
%FLAG SCNB_SCALE_FACTOR
%FORMAT(5E16.8)
  2.00000000E+00
%FLAG SOLTY
%FORMAT(5E16.8)
  0.00000000E+00
%FLAG LENNARD_JONES_ACOEF
%FORMAT(5E16.8)
  3.24095688E+06  2.54904139E+06  2.00225638E+06
%FLAG LENNARD_JONES_BCOEF
%FORMAT(5E16.8)
  1.85348706E+03  1.55090334E+03  1.29687909E+03
%FLAG BONDS_INC_HYDROGEN
%FORMAT(10I8)

%FLAG BONDS_WITHOUT_HYDROGEN
%FORMAT(10I8)
       0       3       1       3       6       2       6       9       1
%FLAG ANGLES_INC_HYDROGEN
%FORMAT(10I8)

%FLAG ANGLES_WITHOUT_HYDROGEN
%FORMAT(10I8)
       0       3       6       1       3       6       9       1
%FLAG DIHEDRALS_INC_HYDROGEN
%FORMAT(10I8)

%FLAG DIHEDRALS_WITHOUT_HYDROGEN
%FORMAT(10I8)
       0       3       6       9       1
%FLAG EXCLUDED_ATOMS_LIST
%FORMAT(10I8)
       2       3       4       3       4       4       0
%FLAG HBOND_ACOEF
%FORMAT(5E16.8)

%FLAG HBOND_BCOEF
%FORMAT(5E16.8)

%FLAG HBCUT
%FORMAT(5E16.8)

%FLAG AMBER_ATOM_TYPE
%FORMAT(20a4)
Cl1 C1  C2  Cl2
%FLAG TREE_CHAIN_CLASSIFICATION
%FORMAT(20a4)
BLA BLA BLA BLA
%FLAG JOIN_ARRAY
%FORMAT(10I8)
       0       0       0       0
%FLAG IROTAT
%FORMAT(10I8)
       0       0       0       0
%FLAG SOLVENT_POINTERS
%FORMAT(3I8)
       1       1       2
%FLAG ATOMS_PER_MOLECULE
%FORMAT(10I8)
       4
%FLAG BOX_DIMENSIONS
%FORMAT(5E16.8)
  9.00000000E+01  4.00000000E+01  4.00000000E+01  4.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
%FLAG SCREEN
%FORMAT(5E16.8)
  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
%FLAG IPOL
%FORMAT(1I8)
       0
mattwthompson commented 1 year ago

With some confidence, I can isolate the issue to the angles about adjacent linear bonds. After #735:

In [10]: molecule = Molecule.from_smiles("CCCCCC")
    ...: molecule.generate_conformers(n_conformers=1)
    ...: topology = molecule.to_topology()
    ...: topology.box_vectors = [4, 4, 4] * unit.nanometer
    ...: parsley = ForceField("openff-1.0.0.offxml")
    ...: out = Interchange.from_smirnoff(force_field=parsley, topology=topology)
    ...:
    ...: get_summary_data(out).std()["Angle"]
/Users/mattthompson/software/openff-interchange/openff/interchange/smirnoff/_valence.py:50: UserWarning: Automatically up-converting BondHandler from version 0.3 to 0.4. Consider manually upgrading this BondHandler (or <Bonds> section in an OFFXML file) to 0.4 or newer. For more details, see https://openforcefield.github.io/standards/standards/smirnoff/#bonds.
  warnings.warn(
Out[10]: 8.851925206475709e-05