"""
Answers to the question of 'are these two structures the same?'
using various definitions of 'the same'
"""
from collections import Counter
from schrodinger.comparison.atom_mapper import ConnectivityAtomMapper
from schrodinger.comparison.atom_mapper import LewisAtomMapper
from schrodinger.comparison.chirality import get_chiral_centers
from schrodinger.structutils import rmsd
[docs]def are_isomers(st1, st2):
"""
Determines if the two structures are isomers
(same number of each type of element)
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:return bool:
"""
isomers = st1.atom_total == st2.atom_total
if isomers:
isomers = st1.property.get('i_m_Molecular_charge',
st1.formal_charge) == st2.property.get(
'i_m_Molecular_charge',
st2.formal_charge)
if isomers:
cntr1 = Counter(at.atomic_number for at in st1.atom)
cntr2 = Counter(at.atomic_number for at in st2.atom)
isomers = cntr1 == cntr2
return isomers
[docs]def are_tautomers(st1, st2, use_stereo=False):
"""
Determines if the two structures are tautomers by inspecting if the
heavy atom connectivity is the same. Here, the definition of the same
uses connectivity only, it is not assumed bond orders are the same.
This does not make any assumptions about atom ordering being initially
correct and, if that is your case this may not be the appropriate tool.
We remove hydrogens that are connected to atoms that are not chiral centers.
This is to avoid hydrogen atoms, which are ignored in this procedure, from
causing symmetry-breaking at remote heavy atom sites. We leave H's directly
attached to chiral centers effectively as dummy atoms to ensure that chirality
of the remaining three substituents is computed.
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:return: bool
"""
tautomers = are_isomers(st1, st2)
if tautomers:
tautomers = are_heavy_atom_conformers(st1, st2, use_stereo=use_stereo)
return tautomers
[docs]def are_stereoisomers(st1, st2, use_lewis_structure=True):
"""
Determines if the two structures are stereoisomers (atom order independent).
A molecule is not a stereoisomer of itself, i.e. are_stereoisomers(st, st)==False
This does not make any assumptions about atom ordering being initially
correct and, if that is your case this may not be the appropriate tool.
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:type use_lewis_structure: bool
:param use_lewis_structure: if True, use atomic symbol, formal charge
and bond orders to define equivalence of atoms, else use only atomic symbol.
The former definition is more traditional whereas the latter
only requires consistent connectivity
:return: bool
"""
stereoisomers = are_isomers(st1, st2)
if stereoisomers:
stereoisomers = are_conformers(st1,
st2,
use_lewis_structure=use_lewis_structure,
use_stereo=False)
if stereoisomers:
stereoisomers = not are_conformers(
st1, st2, use_lewis_structure=use_lewis_structure, use_stereo=True)
return stereoisomers
[docs]def are_enantiomers(st1, st2, use_lewis_structure=True):
"""
Determines if the two structures are enantiomers (atom order independent).
This does not make any assumptions about atom ordering being initially
correct and, if that is your case this may not be the appropriate tool.
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:type use_lewis_structure: bool
:param use_lewis_structure: if True, use atomic symbol, formal charge
and bond orders to define equivalence of atoms, else use only atomic symbol.
The former definition is more traditional whereas the latter
only requires consistent connectivity
:return: bool
"""
enantiomers = are_isomers(st1, st2)
if enantiomers:
enantiomers = not are_conformers(
st1, st2, use_lewis_structure=use_lewis_structure, use_stereo=True)
if enantiomers:
if use_lewis_structure:
mapper = LewisAtomMapper(use_chirality=True)
else:
mapper = ConnectivityAtomMapper(use_chirality=True)
atlist = list(range(1, st1.atom_total + 1))
enantiomers = mapper.are_enantiomers(st1, atlist, st2, atlist)
return enantiomers
[docs]def are_diastereomers(st1, st2, use_lewis_structure=True):
"""
Determines if the two structures are diastereomers (atom order independent).
This does not make any assumptions about atom ordering being initially
correct and, if that is your case this may not be the appropriate tool.
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:type use_lewis_structure: bool
:param use_lewis_structure: if True, use atomic symbol, formal charge
and bond orders to define equivalence of atoms, else use only atomic symbol.
The former definition is more traditional whereas the latter
only requires consistent connectivity
:return: bool
"""
diastereomers = are_stereoisomers(st1,
st2,
use_lewis_structure=use_lewis_structure)
if diastereomers:
if use_lewis_structure:
mapper = LewisAtomMapper(use_chirality=True)
else:
mapper = ConnectivityAtomMapper(use_chirality=True)
atlist = list(range(1, st1.atom_total + 1))
diastereomers = not mapper.are_enantiomers(st1, atlist, st2, atlist)
return diastereomers
[docs]def are_same_geometry(st1, st2, rms_thresh=0.25):
"""
Determines if two structures are conformers with the same geometry.
Assumes that structures are already consistently numbered.
:type st1: Structure
:param st1: the first structure
:type st2: Structure
:param st2: the second structure
:type use_lewis_structure: bool
:param use_lewis_structure: if True, use atomic symbol, formal charge
and bond orders to define equivalence of atoms, else use only atomic symbol.
The former definition is more traditional whereas the latter
only requires consistent connectivity
:type rms_thesh: float
:param rms_thresh: Threshold for RMSD to be considered the same geometry
:return: bool
"""
are_same = are_isomers(st1, st2)
if are_same:
# Opt for permissiveness, rmsd will handle stereochemical and Lewis questions
mapper = ConnectivityAtomMapper(use_chirality=False)
atlist = list(range(1, st1.atom_total + 1))
are_same = mapper.are_consistently_ordered_conformers(st1, st2, atlist)
if are_same:
at_list = list(range(1, st1.atom_total + 1))
rms = rmsd.superimpose(st1.copy(), at_list, st2.copy(), at_list)
are_same = rms < rms_thresh
return are_same