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