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, Topology
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 295.554000 295.516000
1 353.824000 353.775000
2 291.270000 291.208000
3 346.323000 346.268000
4 298.830000 298.783000
5 313.993000 313.938000
6 413.678000 413.611000
7 286.413000 286.387000
8 391.320000 391.257000
9 353.257000 353.254000
10 340.867000 340.783000
11 387.785000 387.695000
12 311.338000 311.306000
13 332.746000 332.684000
14 334.692000 334.597000
15 409.744000 409.675000
16 394.676000 394.620000
17 427.032000 426.973000
18 339.771000 339.658000
19 359.076000 358.989000