Download Notebook View in GitHub Open in Google Colab

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.

from openff.interchange.drivers.openmm import get_openmm_energies
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.2.1 force field appropriate for vacuum calculations (without constraints)
forcefield = ForceField("openff_unconstrained-2.2.1.offxml")

We’ll now apply the Sage force field to the molecule to produce an Interchange, which stores the result of parametrizing with this force field

print(f"Parametrizing {molecule.name} (may take a moment to calculate charges)...")
interchange = forcefield.create_interchange(molecule.to_topology())
print("Done.")

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:
    # Use this conformer to update the positions of the Interchange object
    interchange.positions = conformer

    # Get the (total) initial energy from this conformer and store it
    initial_energies.append(get_openmm_energies(interchange).total_energy)

    # Minimize using Interchange.minimize, which wraps OpenMM
    interchange.minimize(engine="openmm")

    # Record minimized energy and positions
    minimized_energies.append(get_openmm_energies(interchange).total_energy)
    minimized_molecule.add_conformer(interchange.positions.to("angstrom"))

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.m_as("kilocalories_per_mole"),
            min_energy.m_as("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
/    10 :   -62.931 kcal/mol  -98.059 kcal/mol    0.494 Å
/    10 :   -57.099 kcal/mol  -94.665 kcal/mol    0.449 Å
/    10 :   -17.500 kcal/mol  -94.886 kcal/mol    0.929 Å
/    10 :   -55.874 kcal/mol  -96.397 kcal/mol    0.721 Å
/    10 :   -56.941 kcal/mol  -93.608 kcal/mol    0.667 Å
/    10 :   -59.210 kcal/mol  -96.819 kcal/mol    0.574 Å
/    10 :   -55.683 kcal/mol  -94.680 kcal/mol    1.027 Å
/    10 :   -64.933 kcal/mol  -96.754 kcal/mol    0.362 Å
/    10 :   -62.790 kcal/mol  -97.910 kcal/mol    0.749 Å
/    10 :   -60.172 kcal/mol  -97.586 kcal/mol    0.613 Å