Compute conformer energies for a small molecule

This notebook illustrates reading conformers of a molecule from an SDF file and computation of vacuum conformer energies using a SMIRNOFF force field.

Note that absolute vacuum potential energies can be sensitive to small changes in partial charge, for example due to using OpenEye or AmberTools to generate AM1-BCC charges. However, in our experience, relative conformer energies are fairly consistent between AM1-BCC implementations.

Note also that the Open Force Field Toolkit produces deterministic charges that do not depend on the input conformation of parameterized molecules. See the FAQ for more information.

import openmm
from openff.interchange import Interchange
from openff.units.openmm import from_openmm
from rdkit.Chem import rdMolAlign

from openff.toolkit import ForceField, Molecule
from openff.toolkit.utils import get_data_file_path

First, load conformers from an SDF file. SDF files with multiple molecules are always interpreted as a list of molecules rather than as a single molecule with multiple conformers, so we also need to collapse them all into a single molecule.

loaded_molecules = Molecule.from_file(
    get_data_file_path("molecules/ruxolitinib_conformers.sdf"),
)

# Normalize to list
try:
    loaded_molecules = [*loaded_molecules]
except TypeError:
    loaded_molecules = [loaded_molecules]

# from_file loads each entry in the SDF into its own molecule,
# so collapse conformers into the same molecule
molecule = loaded_molecules.pop(0)
for next_molecule in loaded_molecules:
    if next_molecule == molecule:
        for conformer in next_molecule.conformers:
            molecule.add_conformer(conformer)
    else:
        # We're assuming the SDF just has multiple conformers of the
        # same molecule, so raise an error if that's not the case
        raise ValueError("Multiple chemical species loaded")

# Make sure the molecule has a name
if not molecule.name:
    molecule.name = molecule.to_hill_formula()

print(
    f"Loaded {molecule.n_conformers} conformers"
    + f" of {molecule.to_smiles(explicit_hydrogens=False)!r}"
    + f" ({molecule.name})"
)
molecule

Loaded 10 conformers of 'N#CC[C@H](C1CCCC1)n1cc(-c2ncnc3[nH]ccc23)cn1' (ruxolitinib)

Next, load the Sage force field. We’ll use the version without constraints, as it’s more appropriate for energy minimization - see the FAQ for more information about constraints.

# Load the openff-2.1.0 force field appropriate for vacuum calculations (without constraints)
forcefield = ForceField("openff_unconstrained-2.1.0.offxml")

We’ll now apply the Sage force field to the molecule to produce an Interchange, and then an OpenMM Simulation object. Interchange can produce inputs for a whole host of other MD engines though!

print(f"Parametrizing {molecule.name} (may take a moment to calculate charges)...")
interchange = Interchange.from_smirnoff(forcefield, [molecule])
print("Done.")

integrator = openmm.VerletIntegrator(1 * openmm.unit.femtoseconds)
simulation = interchange.to_openmm_simulation(integrator)

Parametrizing ruxolitinib (may take a moment to calculate charges)...

Done.

Now we’ll perform the actual minimizations, and store the outputs in a few lists and a second molecule:

# We'll store energies in two lists
initial_energies = []
minimized_energies = []

# And minimized conformers in a second molecule
minimized_molecule = Molecule(molecule)
minimized_molecule.conformers.clear()

for conformer in molecule.conformers:
    # Tell the OpenMM Simulation the positions of this conformer
    simulation.context.setPositions(conformer.to_openmm())

    # Keep a record of the initial energy
    initial_energies.append(
        simulation.context.getState(getEnergy=True).getPotentialEnergy()
    )

    # Perform the minimization
    simulation.minimizeEnergy()

    # Record minimized energy and positions
    min_state = simulation.context.getState(getEnergy=True, getPositions=True)

    minimized_energies.append(min_state.getPotentialEnergy())
    minimized_molecule.add_conformer(from_openmm(min_state.getPositions()))

We can visualize the minimised Molecule:

minimized_molecule

And we can write the results of the minimisation out to files as a permanent record:

# Get a shortcut to the number of conformers
n_confs = molecule.n_conformers
print(f"{molecule.name}: {n_confs} conformers")

# Create a copy of the molecule so we can work on it
working_mol = Molecule(molecule)

# Print text header
print("Conformer         Initial PE        Minimized PE        RMSD")
output = [
    [
        "Conformer",
        "Initial PE (kcal/mol)",
        "Minimized PE (kcal/mol)",
        "RMSD between initial and minimized conformer (Angstrom)",
    ]
]

for i, (init_energy, init_coords, min_energy, min_coords) in enumerate(
    zip(
        initial_energies,
        molecule.conformers,
        minimized_energies,
        minimized_molecule.conformers,
    )
):
    # Clear the conformers from the working molecule
    working_mol.conformers.clear()

    # Save the minimized conformer to file
    working_mol.add_conformer(min_coords)
    working_mol.to_file(
        f"{molecule.name}_conf{i+1}_minimized.sdf",
        file_format="sdf",
    )

    # Calculate the RMSD between the initial and minimized conformer
    working_mol.add_conformer(init_coords)
    rdmol = working_mol.to_rdkit()
    rmslist = []
    rdMolAlign.AlignMolConformers(rdmol, RMSlist=rmslist)
    minimization_rms = rmslist[0]

    # Record the results
    output.append(
        [
            i + 1,
            init_energy.value_in_unit(openmm.unit.kilocalories_per_mole),
            min_energy.value_in_unit(openmm.unit.kilocalories_per_mole),
            minimization_rms,
        ]
    )
    print(
        f"{{:5d}} / {n_confs:5d} :  {{:8.3f}} kcal/mol {{:8.3f}} kcal/mol {{:8.3f}} Å".format(
            *output[-1]
        )
    )

# Write the results out to CSV
with open(f"{molecule.name}.csv", "w") as of:
    of.write(", ".join(output.pop(0)) + "\n")
    for line in output:
        of.write("{}, {:.3f}, {:.3f}, {:.3f}".format(*line) + "\n")

ruxolitinib: 10 conformers
Conformer         Initial PE        Minimized PE        RMSD
    1 /    10 :   -36.905 kcal/mol  -69.365 kcal/mol    0.499 Å
    2 /    10 :   -29.678 kcal/mol  -66.365 kcal/mol    0.489 Å
    3 /    10 :     7.844 kcal/mol  -66.570 kcal/mol    0.907 Å
    4 /    10 :   -27.968 kcal/mol  -67.889 kcal/mol    0.732 Å
    5 /    10 :   -28.211 kcal/mol  -65.486 kcal/mol    0.874 Å
    6 /    10 :   -30.415 kcal/mol  -68.294 kcal/mol    0.606 Å
    7 /    10 :   -27.233 kcal/mol  -66.318 kcal/mol    1.027 Å
    8 /    10 :   -37.769 kcal/mol  -68.172 kcal/mol    0.370 Å
    9 /    10 :   -34.561 kcal/mol  -69.125 kcal/mol    0.768 Å
   10 /    10 :   -30.336 kcal/mol  -69.049 kcal/mol    0.628 Å