Source code for schrodinger.rdkit.alignment

from math import ceil

import numpy
from rdkit import Chem
from rdkit.Chem import Conformer
from rdkit.Chem import rdFMCS
from rdkit.Chem import rdMolAlign
from rdkit.Chem import rdMolTransforms
from rdkit.Chem.rdchem import Mol
from rdkit.Chem.rdCoordGen import AddCoords
from scipy.spatial.transform import Rotation


[docs]def transform_conformer(conformer: Conformer, matrix: numpy.ndarray):
    """
    Apply a 3x3 transformational matrix to the coordinates of the conformer.

    This is a convenience function that fits the 3x3 matrix into a 4x4 matrix to
    include the extra translational dimension required by `TransformConformer`.

    :param conformer: a conformer with nonzero coordinates
    :param matrix: a 3x3 matrix to transform the conformer coordinates
    """
    mat_with_translation = numpy.identity(4)
    for idx, col in enumerate(matrix):
        mat_with_translation[idx][:3] = col
    rdMolTransforms.TransformConformer(conformer, mat_with_translation)


[docs]def get_height(conformer: Conformer) -> float:
    """
    :param conformer: a conformer with nonzero coordinates
    :return: the maximum y-axis distance between positions of the conformer
    """
    pos = conformer.GetPositions()
    return pos[:, 1].max() - pos[:, 1].min()


[docs]def rotate_2D(conformer: Conformer, rot_angle: float):
    """
    Rotate a conformer clockwise around the z axis.

    :param conformer: a conformer with nonzero coordinates
    :param rot_angle: the angle at which to rotate the conformer (clockwise in
            around the z axis, in degrees)
    """
    rot = Rotation.from_euler(seq='z', angles=rot_angle, degrees=True)
    transform_conformer(conformer, rot.as_matrix())


[docs]def align_horizontally(conformer: Conformer):
    """
    Rotate the specified conformer to minimize its height.

    :param conformer: a conformer with nonzero coordinates
    :note: this optimization is performed over a limited number of possible
            angles in order to preserve aesthetically-pleasing bond angles
            generated for this conformer by `Chem.rdCoordGen`
    """
    degrees_per_rotation = 30
    max_rotation_count = ceil(90 / degrees_per_rotation) + 1
    conf_copy = Conformer(conformer)

    min_height = get_height(conf_copy)
    min_height_angle = 0
    for rotation_count in range(1, max_rotation_count):
        rotate_2D(conf_copy, rot_angle=degrees_per_rotation)
        height = get_height(conf_copy)
        if height < min_height:
            min_height = height
            min_height_angle = rotation_count * degrees_per_rotation

    rotate_2D(conformer, rot_angle=min_height_angle)


[docs]def generate_min_height_coords(mol: Mol):
    """
    Clear the existing conformers of the input `Mol` object and generate one
    such that its height (greatest y-axis distance between two points) is
    minimized to some extent.

    :param mol: input RDKit molecule object
    """

    AddCoords(mol)
    align_horizontally(mol.GetConformer())


[docs]def align_to_template(mol, template_mol):
    """
    Align a 2D `Mol` object against a similar template.

    :param mol: the `Mol` to be aligned
    :type mol: rdchem.Mol
    :param template_mol: the template
    :type template_mol: rdchem.Mol
    :return: the RMSD of the aligned structures
    :rtype: float
    """
    # Use the same method to find MCS as 2D Viewer. Strip explicit hydrogens
    # for a very significant speed up (fixes Maestro hangs):
    try:
        mols = [Chem.RemoveHs(template_mol), Chem.RemoveHs(mol)]
    except Chem.rdchem.AtomValenceException:
        # Handle https://github.com/rdkit/rdkit/issues/3400
        mols = [template_mol, mol]

    mcs = rdFMCS.FindMCS(mols,
                         ringMatchesRingOnly=True,
                         completeRingsOnly=False,
                         timeout=1)

    # Workaround for SHARED-8169: make sure rings are initialized
    Chem.FastFindRings(template_mol)
    Chem.FastFindRings(mol)

    # Now find the MCS query molecule in each input structure, and get
    # the original atom numbers for the matches:
    ref_match = template_mol.GetSubstructMatch(mcs.queryMol)
    match = mol.GetSubstructMatch(mcs.queryMol)
    atom_map = list(zip(match, ref_match))
    return rdMolAlign.AlignMol(mol, template_mol, atomMap=atom_map)