I posted about openforcefield(OpenFF) before. You know, old version of openff supports only OpenEyeTK but current version supports RDKit too.
It is worth to know that we can use openff with open source tool kit. I really appreciate developer’s work! It is great. Today I use the package and ipymol which can control pymol in ipython session. It means that you can control pymol from jupyter notebook! My example code is below.
import openforcefield as off from rdkit import Chem from rdkit.Chem import AllChem from simtk import openmm, unit from simtk.openmm import app from openforcefield.topology import Topology from openforcefield.topology import Molecule from openforcefield.typing.engines.smirnoff import ForceField from rdkit.Chem.Draw import IPythonConsole
Load SMIRNOF99FROSS forcefield.
ff = ForceField('test_forcefields/smirnoff99Frosst.offxml')
Then define get_energy function.
def get_energy(system, positions): integrator = openmm.VerletIntegrator(1.0 * unit.femtoseconds) context = openmm.Context(system, integrator) context.setPositions(positions) state = context.getState(getEnergy=True) energy = state.getPotentialEnergy().in_units_of(unit.kilocalories_per_mole) return energy
Make sample mol object with rdkit and convert openff mol object. It is easy, just call from_rdkit!.
rdmol = Chem.MolFromSmiles('c1c(c2sccc2)c(c3c[nH]cc3)oc1') rdmol = Chem.AddHs(rdmol) rdmol ofmol = Molecule.from_rdkit(rdmol)
Then generate conformer with openff function. And get topology.
ofmol.generate_conformers(n_conformers=10) topology = ofmol.to_topology()
Next, create openMM system with generated topology and get position of one conformer.
org_system = ff.create_openmm_system(topology) pos = ofmol.conformers[0]
At first, calculate energy with default conformation.
get_energy(org_system, pos) > Quantity(value=80.93719044789302, unit=kilocalorie/mole)
Next I tried to minimize energy with openMM method.
new_system = ff.create_openmm_system(topology) new_energy = get_energy(new_system, pos) from sys import stdout def minimizeOpenMM(Topology, System, Positions): integrator = openmm.LangevinIntegrator( 300.0 * unit.kelvin, 1.0 / unit.picosecond, 2.0 * unit.femtosecond) #.002 * unit.picoseconds) simulation = app.Simulation(Topology, System, integrator) simulation.context.setPositions(Positions) simulation.minimizeEnergy(tolerance=5.0E-9, maxIterations=2000) state = simulation.context.getState(getPositions=True, getEnergy=True) positions =state.getPositions(asNumpy=True) energy = state.getPotentialEnergy().in_units_of(unit.kilocalories_per_mole) simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True)) simulation.step(30000) print(energy) positions = positions / unit.angstroms coordlist = list() for atom_coords in positions: coordlist += [i for i in atom_coords] return coordlist, positions cl, pos=minimizeOpenMM(topology, org_system, pos)
Now I could get minimized position by calling minimizeOpenMM.
Next I updated atom position with optimized atom geometries.
from rdkit.Chem import rdGeometry from rdkit.Chem.rdchem import Conformer AllChem.EmbedMolecule(rdmol, useExpTorsionAnglePrefs = True , useBasicKnowledge = True) conf = rdmol.GetConformer() w=Chem.SDWriter("test1.sdf") w.write(rdmol) w.close() for i in range(rdmol.GetNumAtoms()): conf.SetAtomPosition(i, rdGeometry.Point3D(pos[i][0], pos[i][1],pos[i][2],)) w=Chem.SDWriter("test2.sdf") w.write(rdmol) w.close()
I made two SDF one is default conformer and another is minimized conformer.
OK, I use ipymol to communicate pymol and visualize in jupyter notebook.
from ipymol import viewer viewer.start() # this method launches pymol viewer.load("test1.sdf") viewer.load("test2.sdf") viewer.align("test1","test2")
from ipymol import viewer viewer.start() # this method launches pymol viewer.load("test1.sdf") viewer.load("test2.sdf") viewer.align("test1","test2") viewer.ray() viewer.display()

viewer.deleteAll() viewer.load("test2.sdf") viewer.ray() viewer.display()

viewer.deleteAll() viewer.load("test1.sdf") viewer.ray() viewer.display()

The result seems different to original article. https://www.biorxiv.org/content/early/2018/07/13/286542.full.pdf
Hmm why???
BTW, OpenFF is very attractive package for chemoinformatics I think.
My code can access below.
https://nbviewer.jupyter.org/github/iwatobipen/playground/blob/master/openff.ipynb