Download Notebook View in GitHub Open in Google Colab
Compute conformer energies for a small molecule with Interchange
This example generates conformers for a small molecule and then evaluates their energies in OpenMM and GROMACS.
import pandas as pd
from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.units import unit
from openff.interchange import Interchange
from openff.interchange.drivers import get_gromacs_energies, get_openmm_energies
Choose the molecule we want by specifying a SMILES:
SMILES = "c1n(CCO)c(C(F)(F)(F))cc1CNCCl"
Generate conformers with the OpenFF Toolkit. Note that these conformers are not energy-minimized with any force field, QM method, or other tool; their purpose is to generate diverse starting points for further study.
molecule = Molecule.from_smiles(SMILES)
(molecule.generate_conformers(n_conformers=20, rms_cutoff=0.1 * unit.angstrom),)
topology = molecule.to_topology()
Load the Sage 2.1.0 force field:
sage = ForceField("openff-2.0.0.offxml")
Apply the force field to the topology to generate an Interchange. This step might be slow as it needs to compute partial charges:
interchange = Interchange.from_smirnoff(force_field=sage, topology=topology)
Modern versions of GROMACS do not support open boundary conditions, so we’ll choose a huge box to mimic the gas phase for the single-point energy calculation:
interchange.box = unit.Quantity([4, 4, 4], unit.nanometer)
Now we go over each conformer, set the Interchange’s positions, perform the point energy calculation and record the results in a dataframe:
summary = pd.DataFrame()
kj_mol = unit.kilojoule / unit.mol
for idx, conformer in enumerate(molecule.conformers):
interchange.positions = conformer
openmm_energies = get_openmm_energies(interchange).total_energy.m_as(kj_mol)
gromacs_energies = get_gromacs_energies(interchange).total_energy.m_as(kj_mol)
summary = pd.concat(
[
summary,
pd.DataFrame.from_dict(
{
"Conformer No.": [idx],
"Interchange -> OpenMM (kJ/mol)": [round(openmm_energies, 3)],
"Interchange -> GROMACS": [round(gromacs_energies, 3)],
}
),
]
)
Finally, we can look at the results as stored in the dataframe:
summary.style.hide(axis="index")
Conformer No. | Interchange -> OpenMM (kJ/mol) | Interchange -> GROMACS |
---|---|---|
0 | 284.781000 | 284.767000 |
1 | 268.268000 | 268.243000 |
2 | 288.840000 | 288.838000 |
3 | 323.597000 | 323.590000 |
4 | 303.745000 | 303.736000 |
5 | 310.541000 | 310.521000 |
6 | 311.973000 | 311.943000 |
7 | 291.782000 | 291.761000 |
8 | 298.616000 | 298.588000 |
9 | 301.894000 | 301.869000 |
10 | 328.402000 | 328.387000 |
11 | 306.161000 | 306.146000 |
12 | 291.642000 | 291.600000 |
13 | 296.292000 | 296.258000 |
14 | 324.311000 | 324.291000 |
15 | 307.564000 | 307.563000 |
16 | 302.106000 | 302.095000 |
17 | 309.095000 | 309.088000 |
18 | 335.689000 | 335.666000 |