Force field energies and gradients using Open Babel

ff-gradients.py

#!/usr/bin/env python

from __future__ import print_function

import sys
import os

import pybel


if __name__ == "__main__":
    # Syntax  script.py  filename.mol2
    # read molecule
    if (sys.argv) != 2:
        exit

    filename = sys.argv[1]
    extension = os.path.splitext(filename)[1]
    # read a single molecule from the supplied file
    # turn this into a while loop if you want all molecules in a file
    mol = next(pybel.readfile(extension[1:], filename))
    # global variable for the MMFF94 force field
    ff = pybel._forcefields["mmff94"]

    if ff.Setup(mol.OBMol) is False:
        print("Cannot set up force field, exiting")
        exit()

    print("Current energy {}".format(ff.Energy()) )

    # see about gradients
    for atom in mol.atoms:
        # vector objects have to be de-referenced individually (sigh)
        coords = atom.vector
        grad = ff.GetGradient(atom.OBAtom)
        print(grad.GetX(), grad.GetY(), grad.GetZ(), sep=',')

This needs openbabel/openbabel#1833 first.