In this morning I found great news on my tiwtter TL that openforcefield supports RDKit! Older version of openforcefiled is required openeyeTK which is commercial license for industry. I had interested in the package but could not install because I’m not academia and our company does not OETK license now.
I can’t wait to use it so I installed the package and check it immediately.
Following code is ran on google colab. If would like to run on local env, you can skip several steps.
# install conda !wget https://repo.continuum.io/miniconda/Miniconda3-4.4.10-Linux-x86_64.sh !chmod ./Miniconda3-4.4.10-Linux-x86_64.sh !time bash ./Miniconda3-4.4.10-Linux-x86_64.sh -b -f -p /usr/local # install rdkit via conda command !time conda install -c conda-forge rdkit -y
import sys import os sys.path.append('/usr/local/lib/python3.6/site-packages/') from rdkit import Chem from rdkit import rdBase print(rdBase.rdkitVersion)
Then install openforcefield. Openforcefield is required openmm. Cuda version of colab was 10.0 so I did not need change of openmm install version. So, just type install ‘conda install openforcefield’. ;)
! conda config --add channels omnia --add channels conda-forge ! conda install -y openforcefield
Now ready, let’s call openforcefield.
from openforcefield.topology import Molecule from openforcefield.topology import Topology from openforcefield.utils import RDKitToolkitWrapper from openforcefield.typing.engines.smirnoff import ForceField from simtk import openmm from simtk import unit
Following code was borrowed from original repo.
def get_energy(system, positions): """ Return the potential energy. Parameters ---------- system : simtk.openmm.System The system to check positions : simtk.unit.Quantity of dimension (natoms,3) with units of length The positions to use Returns --------- energy """ 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
Load molecule with openFF method and calculate energy.
mol = Molecule.from_smiles('CCO') mol.generate_conformers() positions = mol.conformers #load Force Field ff = ForceField('smirnoff99Frosst.offxml') topology = mol.to_topology() orginal_system = ff.create_openmm_system(topology) orig_energy = get_energy(orginal_system, positions) print(orig_energy) > 2.3399880921635527 kcal/mol
Next calculate energy from rdkit mol object.
mol = Chem.MolFromSmiles('c1ccccc1') rdw = RDKitToolkitWrapper() ofmol = rdw.from_rdkit(mol) ofmol.generate_conformers() positions = ofmol.conformers topology = ofmol.to_topology() orginal_system = ff.create_openmm_system(topology) orig_energy = get_energy(orginal_system, positions) print(orig_energy) >7.987402593200566 kcal/mol
This is very simple example of OpenFF usage with RDKit. The package seems powerful for CADD. I would like to learn the package more and molecular dynamics.
Whole code of the post can check from my repo.
And I really thank to developer of RDKit, Openforcefield and many open source chemoinformatics tools.