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 shutil
import numpy as np
from rdkit.Chem import rdMolAlign
try:
import openmm
from openmm import unit
except ImportError:
from simtk import openmm, unit
from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.utils import RDKitToolkitWrapper, get_data_file_path
Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing
# If using a OFF Toolkit version before 0.7.0, loading SDFs through RDKit and OpenEye may provide
# different behavior in some cases. So, here we force loading through RDKit to ensure the correct behavior
rdktkw = RDKitToolkitWrapper()
# Locate molecule in OpenFF Toolkit package data and copy to local directory
orig_path = get_data_file_path("molecules/ruxolitinib_conformers.sdf")
shutil.copy(orig_path, ".")
# Load in the molecule and its conformers.
# Note that all conformers of the same molecule are loaded as separate Molecule objects
loaded_molecules = Molecule.from_file(
"ruxolitinib_conformers.sdf",
toolkit_registry=rdktkw,
)
# The logic below only works for lists of molecules, so if a
# single molecule was loaded, cast it to list
if type(loaded_molecules) is not list:
loaded_molecules = [loaded_molecules]
# Collatate all conformers of the same molecule
# NOTE: This isn't necessary if you have already loaded or created multi-conformer molecules;
# it is just needed because our SDF reader does not automatically collapse conformers.
molecules = [loaded_molecules[0]]
for molecule in loaded_molecules[1:]:
if molecule == molecules[-1]:
for conformer in molecule.conformers:
molecules[-1].add_conformer(conformer)
else:
molecules.append(molecule)
n_molecules = len(molecules)
n_conformers = sum(mol.n_conformers for mol in molecules)
print(f"{n_molecules} unique molecule(s) loaded, with {n_conformers} total conformers")
1 unique molecule(s) loaded, with 10 total conformers
# Load the openff-2.0.0 force field appropriate for vacuum calculations (without constraints)
from openff.toolkit.typing.engines.smirnoff import ForceField
forcefield = ForceField("openff_unconstrained-2.0.0.offxml")
# Loop over molecules and minimize each conformer
for molecule in molecules:
# If the molecule doesn't have a name, set mol.name to be the hill formula
if molecule.name == "":
molecule.name = Topology._networkx_to_hill_formula(molecule.to_networkx())
print("%s : %d conformers" % (molecule.name, molecule.n_conformers))
# Make a temporary copy of the molecule that we can update for each minimization
mol_copy = Molecule(molecule)
# Make an OpenFF Topology so we can parameterize the system
off_top = molecule.to_topology()
print(f"Parametrizing {molecule.name} (may take a moment to calculate charges)")
system = forcefield.create_openmm_system(off_top)
# Use OpenMM to compute initial and minimized energy for all conformers
integrator = openmm.VerletIntegrator(1 * unit.femtoseconds)
platform = openmm.Platform.getPlatformByName("Reference")
omm_top = off_top.to_openmm()
simulation = openmm.app.Simulation(omm_top, system, integrator, platform)
# Print text header
print(
"Conformer Initial PE Minimized PE RMS between initial and minimized conformer"
)
output = [
[
"Conformer",
"Initial PE (kcal/mol)",
"Minimized PE (kcal/mol)",
"RMS between initial and minimized conformer (Angstrom)",
]
]
for conformer_index, conformer in enumerate(molecule.conformers):
simulation.context.setPositions(conformer)
orig_potential = simulation.context.getState(
getEnergy=True
).getPotentialEnergy()
simulation.minimizeEnergy()
min_state = simulation.context.getState(getEnergy=True, getPositions=True)
min_potential = min_state.getPotentialEnergy()
# Calculate the RMSD between the initial and minimized conformer
min_coords = min_state.getPositions()
min_coords = (
np.array([[atom.x, atom.y, atom.z] for atom in min_coords]) * unit.nanometer
)
mol_copy._conformers = None
mol_copy.add_conformer(conformer)
mol_copy.add_conformer(min_coords)
rdmol = mol_copy.to_rdkit()
rmslist = []
rdMolAlign.AlignMolConformers(rdmol, RMSlist=rmslist)
minimization_rms = rmslist[0]
# Save the minimized conformer to file
mol_copy._conformers = None
mol_copy.add_conformer(min_coords)
mol_copy.to_file(
f"{molecule.name}_conf{conformer_index+1}_minimized.sdf",
file_format="sdf",
)
print(
"%5d / %5d : %8.3f kcal/mol %8.3f kcal/mol %8.3f Angstroms"
% (
conformer_index + 1,
molecule.n_conformers,
orig_potential / unit.kilocalories_per_mole,
min_potential / unit.kilocalories_per_mole,
minimization_rms,
)
)
output.append(
[
str(conformer_index + 1),
f"{orig_potential/unit.kilocalories_per_mole:.3f}",
f"{min_potential/unit.kilocalories_per_mole:.3f}",
f"{minimization_rms:.3f}",
]
)
# Write the results out to CSV
with open(f"{molecule.name}.csv", "w") as of:
for line in output:
of.write(",".join(line) + "\n")
# Clean up OpenMM Simulation
del simulation, integrator
ruxolitinib : 10 conformers Parametrizing ruxolitinib (may take a moment to calculate charges) Conformer Initial PE Minimized PE RMS between initial and minimized conformer 1 / 10 : 17.890 kcal/mol -13.860 kcal/mol 0.606 Angstroms 2 / 10 : 24.274 kcal/mol -11.793 kcal/mol 0.547 Angstroms 3 / 10 : 59.493 kcal/mol -11.919 kcal/mol 0.963 Angstroms 4 / 10 : 25.481 kcal/mol -12.201 kcal/mol 0.733 Angstroms 5 / 10 : 21.963 kcal/mol -10.797 kcal/mol 0.919 Angstroms 6 / 10 : 21.315 kcal/mol -13.756 kcal/mol 0.588 Angstroms 7 / 10 : 24.900 kcal/mol -10.122 kcal/mol 0.693 Angstroms 8 / 10 : 14.714 kcal/mol -12.981 kcal/mol 0.715 Angstroms 9 / 10 : 18.114 kcal/mol -13.708 kcal/mol 0.925 Angstroms 10 / 10 : 21.270 kcal/mol -14.699 kcal/mol 0.776 Angstroms