-
Notifications
You must be signed in to change notification settings - Fork 24
Description
Description
I was trying to calculate the potential energy of a ligand using openff-interchange and openmm. I expected the potential energy calculation to be consistent across different versions of openff-interchange when using the same ligand and force field. However, I observed a significant difference in the calculated potential energy between openff-interchange
version 0.3.18 and 0.4.2.
Reproduction
I am using the following test script, that is executed in two separate conda environments (see: #Software versions) with different versions of openff-interchange
.
import rdkit
from rdkit import Chem
import parmed
import openmm
import openff
from openff.interchange import Interchange
from openff.toolkit import ForceField, Molecule
from openff.units import unit
import os
print(f"rdkit version {rdkit.__version__}")
print(f"openmm version {openmm.__version__}")
print(f"parmed version {parmed.__version__}")
print(f"openff-interchange version {openff.interchange.__version__}")
print(f"openff-toolkit version {openff.toolkit.__version__}")
print(f"openff-units version {openff.units.__version__}")
conda_env_name = os.path.basename(os.environ["CONDA_PREFIX"])
print(f"Conda env name {conda_env_name}")
def load_ligand_from_rdkit_mol(
rdmol: rdkit.Chem.rdchem.Mol | None,
residue_name: str | None = None,
forcefield: str | None = "parsley",
) -> parmed.amber.AmberParm:
"""Load a ligand from an RDKit Mol object.
Args:
rdmol: An RDKit Mol object for the ligand
residue_name: The (optional) residue name to assign to the ligand.
Returns:
The loaded ligand
"""
mol = Molecule.from_rdkit(rdmol, allow_undefined_stereo=True)
if forcefield == "sage":
ff = ForceField("openff-2.2.0.offxml")
elif forcefield == "parsley":
ff = ForceField("openff-1.3.1.offxml")
else:
raise RuntimeError("Force-field is not available.")
interchange = Interchange.from_smirnoff(topology=[mol], force_field=ff, box=None)
system = interchange.to_openmm_system(add_constrained_forces=True)
structure = parmed.openmm.load_topology(
mol.to_topology().to_openmm(), system, interchange.positions.m_as(unit.angstrom)
)
ligand = parmed.amber.AmberParm.from_structure(structure)
for residue in ligand.residues:
residue.name = residue_name
ligand.parm_data["RESIDUE_LABEL"] = [residue_name]
return ligand
ligands_file = "./ligands.sdf"
supplier = Chem.SDMolSupplier(ligands_file, removeHs=False)
mol = supplier[0]
ligand = load_ligand_from_rdkit_mol(
rdmol=mol, residue_name="L1", forcefield="parsley"
)
create_system_args = {
"nonbondedMethod": openmm.app.NoCutoff,
"constraints": openmm.app.HBonds,
}
system = ligand.createSystem(**create_system_args)
with open(f"{conda_env_name}-system.xml", 'w') as f:
f.write(openmm.XmlSerializer.serialize(system))
platform = openmm.Platform.getPlatformByName('Reference')
integrator = openmm.VerletIntegrator(0.01)
context = openmm.openmm.Context(system, integrator, platform)
context.setPositions(ligand.positions)
print(f"Potential Energy {context.getState(getEnergy=True).getPotentialEnergy()}")
Given below is an example ligands.sdf
file
1H1Q
3D
Schrodinger Suite 2022-2.
45 48 0 0 1 0 999 V2000
7.8210 43.7020 49.3990 C 0 0 0 0 0 0
10.7390 46.1550 51.6470 C 0 0 0 0 0 0
5.1160 45.2740 51.6170 C 0 0 0 0 0 0
4.1070 45.7240 52.7030 C 0 0 0 0 0 0
2.7660 45.9860 51.9860 C 0 0 0 0 0 0
1.5620 46.0860 52.9390 C 0 0 0 0 0 0
1.4480 44.9690 53.9870 C 0 0 0 0 0 0
2.7870 44.7400 54.7250 C 0 0 0 0 0 0
4.0090 44.6090 53.7920 C 0 0 0 0 0 0
4.6670 40.9170 47.4030 C 0 0 0 0 0 0
4.0590 40.6480 48.6330 C 0 0 0 0 0 0
4.6450 41.0930 49.8150 C 0 0 0 0 0 0
5.8520 41.7890 49.8080 C 0 0 0 0 0 0
6.8700 44.4890 50.1810 N 0 0 0 0 0 0
7.2370 45.3850 51.0370 C 0 0 0 0 0 0
6.3200 46.0740 51.7360 O 0 0 0 0 0 0
8.6590 45.6250 51.1990 C 0 0 0 0 0 0
9.5490 46.4520 51.9630 N 0 0 0 0 0 0
10.7730 45.2130 50.7360 N 0 0 0 0 0 0
9.5340 44.9120 50.4780 C 0 0 0 0 0 0
9.0610 43.9200 49.5710 N 0 0 0 0 0 0
7.6000 42.7620 48.5320 N 0 0 0 0 0 0
6.4450 42.0610 48.5810 C 0 0 0 0 0 0
5.8550 41.6320 47.3960 C 0 0 0 0 0 0
4.6775 45.4225 50.6302 H 0 0 0 0 0 0
5.3628 44.2225 51.7641 H 0 0 0 0 0 0
4.4581 46.6481 53.1622 H 0 0 0 0 0 0
2.5854 45.1961 51.2570 H 0 0 0 0 0 0
2.8416 46.9029 51.4015 H 0 0 0 0 0 0
0.6428 46.1264 52.3546 H 0 0 0 0 0 0
1.5815 47.0521 53.4434 H 0 0 0 0 0 0
1.1399 44.0442 53.4993 H 0 0 0 0 0 0
0.6748 45.2279 54.7104 H 0 0 0 0 0 0
2.7092 43.8475 55.3459 H 0 0 0 0 0 0
2.9556 45.5551 55.4288 H 0 0 0 0 0 0
4.9208 44.6063 54.3892 H 0 0 0 0 0 0
3.9911 43.6323 53.3084 H 0 0 0 0 0 0
4.2108 40.5699 46.4877 H 0 0 0 0 0 0
3.1327 40.0939 48.6677 H 0 0 0 0 0 0
4.1634 40.9006 50.7623 H 0 0 0 0 0 0
6.3008 42.1012 50.7394 H 0 0 0 0 0 0
8.2898 42.5561 47.8235 H 0 0 0 0 0 0
6.3393 41.8654 46.4593 H 0 0 0 0 0 0
11.5200 44.7341 50.2536 H 0 0 0 0 0 0
11.5662 46.6643 52.1189 H 0 0 0 0 0 0
1 14 2 0 0 0
1 21 1 0 0 0
1 22 1 0 0 0
2 18 2 0 0 0
2 19 1 0 0 0
2 45 1 0 0 0
3 4 1 0 0 0
3 16 1 0 0 0
3 25 1 0 0 0
3 26 1 0 0 0
4 5 1 0 0 0
4 9 1 0 0 0
4 27 1 0 0 0
5 6 1 0 0 0
5 28 1 0 0 0
5 29 1 0 0 0
6 7 1 0 0 0
6 30 1 0 0 0
6 31 1 0 0 0
7 8 1 0 0 0
7 32 1 0 0 0
7 33 1 0 0 0
8 9 1 0 0 0
8 34 1 0 0 0
8 35 1 0 0 0
9 36 1 0 0 0
9 37 1 0 0 0
10 11 2 0 0 0
10 24 1 0 0 0
10 38 1 0 0 0
11 12 1 0 0 0
11 39 1 0 0 0
12 13 2 0 0 0
12 40 1 0 0 0
13 23 1 0 0 0
13 41 1 0 0 0
14 15 1 0 0 0
15 16 1 0 0 0
15 17 2 0 0 0
17 18 1 0 0 0
17 20 1 0 0 0
19 20 1 0 0 0
19 44 1 0 0 0
20 21 2 0 0 0
22 23 1 0 0 0
22 42 1 0 0 0
23 24 2 0 0 0
24 43 1 0 0 0
M END
> <s_m_entry_id>
9
> <s_m_entry_name>
1H1Q.1
> <s_m_Source_Path>
/Users/sheenam/Downloads
> <s_m_Source_File>
1H1Q.mol2
> <i_m_Source_File_Index>
1
$$$$
Output
Output when running the above script in new
conda environment.
/scratch/micromamba/envs/new/lib/python3.10/site-packages/smirnoff99frosst/smirnoff99frosst.py:11: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
from pkg_resources import resource_filename
rdkit version 2024.03.5
openmm version 8.2
parmed version 4.3.0
openff-interchange version 0.4.2
openff-toolkit version 0.16.9
openff-units version 0.2.1
Conda env name new
Potential Energy -385.61852136708023 kJ/mol
Output produced when running the above script in the old
conda environment.
/scratch/micromamba/envs/old/lib/python3.10/site-packages/smirnoff99frosst/smirnoff99frosst.py:11: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
from pkg_resources import resource_filename
rdkit version 2024.03.6
openmm version 8.2
parmed version 4.3.0
openff-interchange version 0.3.18
openff-toolkit version 0.14.3
openff-units version 0.2.1
Conda env name old
Potential Energy -390.3161788158834 kJ/mol
Software versions
Please find below two minimal env.yaml
I used to conda environments (via micromamba env create -f file.yaml
)
name: old
channels:
- conda-forge
- nvidia
dependencies:
- python=3.10
- ambertools=23.6
- pip
- numpy <2
- scipy
- pandas
- gputil
- jax ~=0.4
- jaxlib ~=0.4
- ambertools=23.6
- cudatoolkit=11.8.0
- openmm=8.2.0
- openff-toolkit=0.14.3
- openff-units=0.2.1
- rdkit
- parmed
name: new
channels:
- conda-forge
- nvidia
dependencies:
- python=3.10
- ambertools=23.6
- pip
- numpy <2
- scipy
- pandas
- gputil
- jax ~=0.4
- jaxlib ~=0.4
- ambertools=23.6
- cudatoolkit=11.8.0
- openmm=8.2.0
- openff-toolkit=0.16.9
- openff-units=0.2.1
- rdkit
- parmed