rdtorsion.py

#!/usr/bin/env python

# Thanks to Peter Schmidtke for the guts of this script
# https://pschmidtke.github.io/blog/rdkit/crystallography/small%20molecule%20xray/xray/database/2021/01/25/cod-and-torsion-angles.html

import pandas as pd
from rdkit.Chem import rdMolTransforms

import numpy as np

suppl = Chem.SDMolSupplier('out.sdf',removeHs=False)

torsions=pd.read_table("list_torsion_patterns.txt",header=None,usecols=[1])
patterns=torsions[1][:3]

for torsionSmarts in patterns:
    print(torsionSmarts)
    angles = np.zeros(360)
    total_matches = 0

    torsionQuery = Chem.MolFromSmarts(torsionSmarts)

    # these SMARTS have atom maps, so convert them
    # http://www.rdkit.org/docs/GettingStartedInPython.html#atom-map-indices-in-smarts
    index_map = {}
    for atom in torsionQuery.GetAtoms() :
        map_num = atom.GetAtomMapNum()
        if map_num:
            index_map[map_num-1] = atom.GetIdx()

    map_list = [index_map[x] for x in sorted(index_map)]

    # loop through the molecules
    for mol in suppl:
        
        if mol is None: continue

        # get the 3D geometry
        conf = mol.GetConformer(0)
        matches = mol.GetSubstructMatches(torsionQuery)
        for match in matches:
            total_matches += 1 # to normalize

            # get the atom maps from the SMARTS match
            mapped = [match[x] for x in map_list]
            angle = rdMolTransforms.GetDihedralDeg(conf, mapped[0],mapped[1],mapped[2],mappped[3])
            if (angle < 0.0):
                angle += 360.0

            angles[angle] += 1
    
    angles = angles / total_matches
    print(angles)