Source code for schrodinger.analysis.substructure
"""
Functionality for generating connected molecule substructures.
Copyright Schrodinger, LLC. All rights reserved.
"""
[docs]class ConnectedMoleculeChain:
"""
Superpose st2 and st1 by moving each molecule independently.
This is accomplished by fixing the first molecule of st1 (with atoms in the atom
list), and aligning all molecules in st2 that contain atoms which map to the atoms
in this molecule of st1. These atoms are now fixed and all molecules in st1 which
contain atoms which map to the now-fixed atoms of st2 are aligned to the existing
fixed st1/st2 molecules. This fixing and aligning is continued until all molecules
in the ensemble connected to the first molecule in st1 is completed. If this did not
include all molecules in the st1/st2 system, the next, non-fixed st1 molecule is
used to start a new fixing/aligning chain. If there are still un-fixed molecules
when all the st1 molecules are fixed, then these atoms only appear in st2 and so
are fixed at the input st2 locations.
:param st1: first structure
:type st1: Structure
:param atlist1: list of atom indexes for structure 1
:type atlist1: list of ints
:param st2: second structure
:type st2: Structure
:param atlist2: list of atom indexes for structure 2
:type atlist2: list of ints
:return: rmsd, st2 and st1 rotated and translated to align
"""
[docs] def __init__(self, st1, atlist1, st2, atlist2):
self._st1 = st1
self._st2 = st2
self._atlist1 = atlist1
self._atlist2 = atlist2
self._atmap12 = dict(zip(atlist1, atlist2))
self._atmap21 = dict(zip(atlist2, atlist1))
self._n_mols_st1 = len(st1.molecule)
self._n_mols = self._n_mols_st1 + len(st2.molecule)
def _initializeChain(self):
"""
Begin the molecule connectivity chain.
Arbitrarily, we consider the molecule which contains the first
atom in the atom list in the first structure to be the first
fixed molecule and start the chain here.
"""
start_mol = self._st1.atom[self._atlist1[0]].molecule_number
self._completed_mols = [start_mol]
self._mol_offset = self._n_mols_st1
self._atmap_move_to_fixed = self._atmap12
self._st_fixed = self._st2
self._st_move = self._st1
self._st1_located_atoms = set(
at for at in self._st1.molecule[start_mol].getAtomIndices()
if at in self._atlist1)
self._connected_mols = self._getNewConnections(start_mol)
self._moving_st = 1
def _setMovingParts(self):
"""
Set the relevant parameters to propogate the molecule connectivity
chain.
This involves determining whether the next molecule in _connected_mols
is in Structure1 or Structure2 and setting the corresponding atom map
directions accordingly.
"""
self._st_fixed = self._st2
self._st_move = self._st1
mol_idx = self._connected_mols.pop()
self._completed_mols.append(mol_idx)
self._atlist = self._atlist1
self._atmap_move_to_fixed = self._atmap12
self._mol_offset = self._n_mols_st1
self._located_atoms = self._st1_located_atoms
self._moving_st = 0
if mol_idx > self._n_mols_st1:
self._st_fixed = self._st1
self._st_move = self._st2
mol_idx -= self._n_mols_st1
self._atlist = self._atlist2
self._atmap_move_to_fixed = self._atmap21
self._mol_offset = 0
self._located_atoms = set(
self._atmap12[idx] for idx in self._st1_located_atoms)
self._moving_st = 1
return mol_idx
def _getNewConnections(self, mol_idx):
"""
Get the molecules that have not already been moved in the fixed
structure that are connected to the moved molecule
A note about the variables used below:
st_fixed: fixed structure that was not just moved
st_move: structure that was just moved
atmap_move_to_fixed: map from the moved structure to the fixed structure
mol_offset: offset either 0 or # of molecules in structure 1, determined
by which molecule is being moved
mol_idx: index of just moved molecule (in the moved structure)
This function returns a set of unfixed molecules in the chain
connected to the just moved molecule.
"""
at_in_fixed_st = [
self._atmap_move_to_fixed[idx]
for idx in self._st_move.molecule[mol_idx].getAtomIndices()
if idx in self._atmap_move_to_fixed.keys()
]
new_cnxns = set(self._st_fixed.atom[idx].molecule_number +
self._mol_offset for idx in at_in_fixed_st)
new_cnxns = set(new_mols for new_mols in new_cnxns
if new_mols not in self._completed_mols)
return new_cnxns
[docs] def getConnectivityChain(self):
"""
Build a chain of which structure moves, which atoms are fixed, and the
corresponding moving atoms.
:param return: The chain has the format (st_moving, fixed_list, moving_list) where
st_moving is 0 for the first structure and 1 for the second, and the
lists are lists of ints.
:param type: 3-tuple of int, list of ints, list of ints
"""
chain = []
self._initializeChain()
while len(self._completed_mols) != self._n_mols:
while len(self._connected_mols) > 0:
mol_idx = self._setMovingParts()
move_list = [
idx
for idx in self._st_move.molecule[mol_idx].getAtomIndices()
if idx in self._atlist and idx in self._located_atoms
]
if move_list:
fixed_list = [
self._atmap_move_to_fixed[i] for i in move_list
]
chain.append((self._moving_st, fixed_list, move_list))
at_to_st1 = dict(zip(self._atlist, self._atlist1))
self._st1_located_atoms |= set(
at_to_st1[idx]
for idx in
self._st_move.molecule[mol_idx].getAtomIndices()
if idx in self._atlist)
self._connected_mols |= self._getNewConnections(mol_idx)
missing_ats = set(self._atlist1) - self._st1_located_atoms
if len(missing_ats) > 0:
new_start_mol = self._st1.atom[
missing_ats.pop()].molecule_number
self._connected_mols = set([new_start_mol])
self._st1_located_atoms |= set(
at for at in
self._st1.molecule[new_start_mol].getAtomIndices()
if at in self._atlist1)
else:
missing_mols = set(range(1, self._n_mols + 1)) - set(
self._completed_mols)
if len(missing_mols) > 0:
self._connected_mols = set([missing_mols.pop()])
return chain