Scan dihedral angles using Open Babel / Pybel

scan-dihedrals.py

#!/usr/bin/env python

import sys
import os
import math

from openbabel import pybel

# read in a file and scan the torsions
# generate a bunch of XYZ files with varied phi / psi angles

if __name__ == "__main__":

    filename = sys.argv[1]
    extension = os.path.splitext(filename)[1]
    # read the molecule from the supplied file
    mol = next(pybel.readfile(extension[1:], filename))

    print("Molecule read. Number of atoms: ", len(mol.atoms))

    # example peptide bond (e.g., phi-psi)
    smarts = pybel.Smarts("O=CC([#1])[NH]C")
    matches = smarts.findall(mol)
    match = matches[0]
    if len(match) < 4:
        print("error, need at least 4 atoms for %d" % len(match))

    # get the atom indexes from the SMARTS match
    a = (match[0])
    b = (match[1])
    c = (match[2])
    d = (match[3])

    for theta in range(-180, 190, 20): # make sure we get +180
        mol.OBMol.SetTorsion(a, b, c, d, theta)

        # write the current coordinates to a file
        # (xyz here, but change as you want)
        outFileName = "{}_{}.xyz".format(os.path.splitext(filename)[0], theta)
        mol.write("xyz", outFileName, overwrite=True)