Conformer Minimization¶

This example will show how to optimize a conformer of paracetamol.

Load in a paracetamol molecule, generate a conformer for it, and perturb the conformer to ensure it needs minimization.

In [1]:

import openff.toolkit
import openff.units
import torch

molecule = openff.toolkit.Molecule.from_smiles("CC(=O)NC1=CC=C(C=C1)O")
molecule.generate_conformers(n_conformers=1)

conformer = torch.tensor(molecule.conformers[0].m_as(openff.units.unit.angstrom)) * 1.10
conformer.requires_grad = True

import openff.toolkit import openff.units import torch molecule = openff.toolkit.Molecule.from_smiles("CC(=O)NC1=CC=C(C=C1)O") molecule.generate_conformers(n_conformers=1) conformer = torch.tensor(molecule.conformers[0].m_as(openff.units.unit.angstrom)) * 1.10 conformer.requires_grad = True

We specify that the gradient of the conformer is required so that we can optimize it using PyTorch.

Parameterize the molecule using OpenFF Interchange and convert it into a PyTorch tensor representation.

In [2]:

import openff.interchange

interchange = openff.interchange.Interchange.from_smirnoff(
    openff.toolkit.ForceField("openff-2.1.0.offxml"),
    molecule.to_topology(),
)

import smee.converters

force_field, [topology] = smee.converters.convert_interchange(interchange)

import openff.interchange interchange = openff.interchange.Interchange.from_smirnoff( openff.toolkit.ForceField("openff-2.1.0.offxml"), molecule.to_topology(), ) import smee.converters force_field, [topology] = smee.converters.convert_interchange(interchange)

We can minimize the conformer using any of PyTorch's optimizers.

In [3]:

import smee

optimizer = torch.optim.Adam([conformer], lr=0.02)

for epoch in range(75):
    energy = smee.compute_energy(topology, force_field, conformer)
    energy.backward()

    optimizer.step()
    optimizer.zero_grad()

    if epoch % 5 == 0 or epoch == 74:
        print(f"Epoch {epoch}: E={energy.item()} kcal / mol")

import smee optimizer = torch.optim.Adam([conformer], lr=0.02) for epoch in range(75): energy = smee.compute_energy(topology, force_field, conformer) energy.backward() optimizer.step() optimizer.zero_grad() if epoch % 5 == 0 or epoch == 74: print(f"Epoch {epoch}: E={energy.item()} kcal / mol")

Epoch 0: E=102.10968017578125 kcal / mol
Epoch 5: E=7.088213920593262 kcal / mol
Epoch 10: E=-18.331130981445312 kcal / mol
Epoch 15: E=-22.182296752929688 kcal / mol
Epoch 20: E=-30.369152069091797 kcal / mol
Epoch 25: E=-36.81045150756836 kcal / mol
Epoch 30: E=-38.517852783203125 kcal / mol
Epoch 35: E=-40.50505828857422 kcal / mol
Epoch 40: E=-42.08476257324219 kcal / mol
Epoch 45: E=-42.19199752807617 kcal / mol
Epoch 50: E=-42.37827682495117 kcal / mol
Epoch 55: E=-42.6767692565918 kcal / mol
Epoch 60: E=-42.799903869628906 kcal / mol
Epoch 65: E=-42.94251251220703 kcal / mol
Epoch 70: E=-43.037200927734375 kcal / mol
Epoch 74: E=-43.084136962890625 kcal / mol

We can then re-store the optimized conformer back into the molecule. Here we add the conformer to the molecule's conformer list, but we could also replace the original conformer.

In [4]:

molecule.add_conformer(conformer.detach().numpy() * openff.units.unit.angstrom)
molecule.visualize(backend="nglview")

molecule.add_conformer(conformer.detach().numpy() * openff.units.unit.angstrom) molecule.visualize(backend="nglview")

NGLWidget(max_frame=1)