import schrodinger.utils.log as log
from schrodinger import structure
from schrodinger.infra import phase
logger = log.get_output_logger(name="fragment_snapcore")
[docs]def get_aligned_coordinates_from_fragments(ct1, ct2, match1, match2,
broken_bond_ct1, broken_bond_ct2):
"""
Send each connected fragment pairs for snapcore coordinates.
If the fragments have less than three atoms and contain no
dummy atoms, generate snapcore coordinates. If case of stereo failure
or phase exeception, return None for coordinate array.
ct1: Structure from first molecule
ct2: Structure from second molecule
matched1: list of matched atom indices in ct1
matched2: list of matched atom indices in ct2
broken_bond_ct1: list of bond tuples
broken_bond_ct2: list of bond tuples
return: tuple of snapped core coordinates for ct1 and ct2
if no failure, otherwise None for failed snapcore
"""
ORIG_ATOM_INDEX = "i_fep_original_atom_index"
atom_map = {a: b for (a, b) in zip(match1, match2)}
st1 = ct1.copy()
st2 = ct2.copy()
for ct in [st1, st2]:
for a in ct.atom:
a.property[ORIG_ATOM_INDEX] = a.index
for ct, bonds in zip((st1, st2), (broken_bond_ct1, broken_bond_ct2)):
for ai, aj in bonds:
ct.deleteBond(ai, aj)
frag1 = [mol.extractStructure(copy_props=True) for mol in st1.molecule]
all_wt_orig2frag = {}
for f_idx, frag in enumerate(frag1):
all_wt_orig2frag[f_idx] = {
a.property[ORIG_ATOM_INDEX]: a.index
for a in frag.atom
}
frag_used = set()
coord1 = st1.getXYZ()
coord2 = st2.getXYZ()
snap_success1 = False
snap_success2 = False
for fmut in [mol.extractStructure(copy_props=True) for mol in st2.molecule]:
mut_orig2frag = {
a.property[ORIG_ATOM_INDEX]: a.index
for a in fmut.atom
}
matched1 = []
matched2 = []
fwt = None
wt_orig2frag = None
for fwt_idx in set(range(len(frag1))) - frag_used:
wt_map = all_wt_orig2frag[fwt_idx]
for a in wt_map:
if a in atom_map and atom_map[a] in mut_orig2frag:
matched2.append(mut_orig2frag[atom_map[a]])
matched1.append(wt_map[a])
if len(matched1) > 0:
frag_used.add(fwt_idx)
fwt = frag1[fwt_idx]
wt_orig2frag = wt_map
break
if len(matched1) == 0:
continue
if len(matched1) > 2:
snap_success2 = get_aligned_fragment_coords(
fwt, fmut, matched1, matched2, mut_orig2frag, coord2)
snap_success1 = get_aligned_fragment_coords(
fmut, fwt, matched2, matched1, wt_orig2frag, coord1)
elif (len(matched1) == fwt.atom_total) and (
len(matched2) == fmut.atom_total):
# just copy all coordinates
for o_idx in wt_orig2frag:
coord1[o_idx - 1] = ct2.atom[atom_map[o_idx]].xyz
coord2[atom_map[o_idx] - 1] = ct1.atom[o_idx].xyz
else:
logger.info(
"skipping some fragments for snapcore because there are less than 3 atoms"
)
ret_val1 = None
ret_val2 = None
if snap_success1:
ret_val1 = coord1
if snap_success2:
ret_val2 = coord2
return (ret_val1, ret_val2)
[docs]def get_aligned_fragment_coords(fwt, fmut, matched1, matched2, map_orig2frag,
coord):
"""
fwt: fragment Structure from first molecule
fmut: fragment Structure from second molecule
matched1: list of matched atom indices in fwt
matched2: list of matched atom indices in fmut
map_orig2frag: dictionary mapping original atom indices to fragment indices
coord: original coordinates returned by mut_ct.getXYZ(), if snapcore is possible
aligned coordinates of fmut atoms are copied
return: False if alignment failed otherwise True
"""
try:
aligned_coord = get_aligned_coordinates(fwt, fmut, matched1, matched2)
except:
#In case of PhpException and RuntimeError for stereo failure
#don't update coordinates
return False
inv_dict = {v: k for (k, v) in map_orig2frag.items()}
for a in fmut.atom:
idx = inv_dict[a.index] - 1
for ii in range(3):
coord[idx][ii] = aligned_coord[(a.index - 1)][ii]
return True
[docs]def get_aligned_coordinates(ct1,
ct2,
matched_atom1,
matched_atom2,
broken_bond_ct2=[]): # noqa: M511
"""
wrapper for PhpAlignCore
ct1: first Structure
ct2: second Structure
matched_atom1: list of mateched atom in ct1
matched_atom2: list of mateched atom in ct2
broken_bond_ct2: list of broken bonds in ct2
return: numpy array of aligned coordinates for ct2 atoms
if successful, otherwise None
"""
align_options = phase.PhpAlignCoreOptions()
align_options.stereo_change_action = phase.PhpStereoChangeAction_SAVE_MAPPINGS
align_options.stereo_change = phase.PhpStereoChange_INVERSION_ONLY
aligned = phase.PhpAlignCore(ct1, ct2, matched_atom1, matched_atom2,
broken_bond_ct2, align_options)
aligned_ct_coord = None
if not aligned.stereoFailure():
aligned_ct = structure.Structure(aligned.getAlignedB()[0])
return aligned_ct.getXYZ()
else:
raise RuntimeError(
'WARNING: Call to phase_align_core failed due to stereoFailure: template ct: %s, other ct:%s'
% (ct1.title, ct2.title))