"""
Utilities for swapping fragments.
Copyright Schrodinger, LLC. All rights reserved.
"""
from schrodinger.application.matsci import reaction_workflow_utils as rxnwfu
from schrodinger.application.matsci import msutils
from schrodinger.infra import structure as infrastructure
from schrodinger.structutils import measure
from schrodinger.structutils import rmsd
KEEP_ATOM_PROP = rxnwfu.KEEP_ATOM_KEY
SUPERPOSITION_ATOM_PROP = rxnwfu.SUPERPOSITION_ATOM_KEY
DISTANCE_CELL_LEN = 0.75
[docs]def get_keep_idxs(st):
"""
Return a list of indices of keep atoms
in the given structure.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: list
:return: contains indices of keep atoms
"""
return rxnwfu.get_idxs_marked_atoms(st, KEEP_ATOM_PROP)
def _get_remaining_atom_idxs(struct, idxs):
"""
Return remaining atom indices in the given structure
that are not in the given indices.
:type struct: schrodinger.structure.Structure
:param struct: the structure
:type idxs: list
:param idxs: the indices
:rtype: list
:return: the remaining indices
"""
all_idxs = set(range(1, struct.atom_total + 1))
return sorted(all_idxs.difference(idxs))
[docs]def get_superposition_idxs(st):
"""
Return a list of indices of superposition atoms
in the given structure.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: list
:return: contains indices of superposition atoms
"""
idxs = rxnwfu.get_idxs_marked_atoms(st, SUPERPOSITION_ATOM_PROP)
getter = lambda x: st.atom[x].property[SUPERPOSITION_ATOM_PROP]
return sorted(idxs, key=getter)
[docs]def get_extracted_and_maps(st, idxs):
"""
Extract and return a structure from the given indices as
well as old-to-new and new-to-old atom index maps.
:type st: schrodinger.structure.Structure
:param st: the structure
:type idxs: list
:param idxs: the indices
:rtype: schrodinger.structure.Structure, dict, dict
:return: the extracted structure and old-to-new and
new-to-old index maps
"""
est = st.extract(idxs, copy_props=True)
old_to_new_map = dict(zip(idxs, range(1, len(idxs) + 1)))
new_to_old_map = dict(zip(range(1, len(idxs) + 1), idxs))
return est, old_to_new_map, new_to_old_map
[docs]def get_cut_bonds(st, keep_idxs):
"""
Return a list of (keep, replace, order) tuples, where
keep is an index that will be kept, replace is an index
that will be replaced, and order is the bond order, for
bonds involving atoms specified in the given keep indices.
:type st: schrodinger.structure.Structure
:param st: the structure
:type keep_idxs: list
:param keep_idxs: the keep indices
:rtype: list
:return: contains (keep, replace, order) tuples
"""
bonds = []
for bond in st.bond:
idx_1 = bond.atom1.index
idx_2 = bond.atom2.index
keep_1 = idx_1 in keep_idxs
keep_2 = idx_2 in keep_idxs
if keep_1 and not keep_2:
bonds.append((idx_1, idx_2, bond.order))
elif not keep_1 and keep_2:
bonds.append((idx_2, idx_1, bond.order))
return bonds
[docs]def get_closest_atom(nov_st, ref_st, nov_idx, ref_cell):
"""
Return the reference atom index closest to the given novel
atom index.
:type nov_st: schrodinger.structure.Structure
:param nov_st: first structure, called novel
:type ref_st: schrodinger.structure.Structure
:param ref_st: second structure, called reference
:type nov_idx: int
:param nov_idx: the novel index
:type ref_cell: infrastructure.DistanceCell
:param ref_cell: distance cell for the reference structure
:rtype: int
:return: the reference index
"""
nov_point = nov_st.atom[nov_idx].xyz
ref_idxs = [match.getIndex() for match in ref_cell.query_atoms(*nov_point)]
ref_idx_distance_pairs = []
for ref_idx in ref_idxs:
ref_point = ref_st.atom[ref_idx].xyz
pair = (ref_idx, measure.measure_distance(nov_point, ref_point))
ref_idx_distance_pairs.append(pair)
if ref_idx_distance_pairs:
return sorted(ref_idx_distance_pairs, key=lambda x: x[1])[0][0]
def _get_bonds_for_assembling(nov_st,
ref_st,
nov_keep_idxs,
ref_keep_idxs,
require_identical_bonds=True):
"""
Return a list of (novel atom index, reference atom index, order)
tuples of bonds that should be created upon assembly of a new
structure containing the specified keep indices from the specified
structures.
:type nov_st: schrodinger.structure.Structure
:param nov_st: first structure, called novel
:type ref_st: schrodinger.structure.Structure
:param ref_st: second structure, called reference
:type nov_keep_idxs: list
:param nov_keep_idxs: novel keep indices
:type ref_keep_idxs: list
:param ref_keep_idxs: reference keep indices
:type require_identical_bonds: bool
:param require_identical_bonds: whether to require that bonds to be
created must exist in both novel and reference structures and
be of the same bond order
:rtype: list
:return: contains (novel atom index, reference atom index, order)
tuples of bonds to be created
"""
nov_bonds = get_cut_bonds(nov_st, nov_keep_idxs)
ref_bonds = get_cut_bonds(ref_st, ref_keep_idxs)
if not nov_bonds or not ref_bonds:
return []
ref_cell = infrastructure.DistanceCell(ref_st, DISTANCE_CELL_LEN)
bonds_to_make = []
for nov_keep, nov_replace, order in nov_bonds:
ref_replace = get_closest_atom(nov_st, ref_st, nov_keep, ref_cell)
ref_keep = get_closest_atom(nov_st, ref_st, nov_replace, ref_cell)
if (ref_keep, ref_replace, order) in ref_bonds or \
(not require_identical_bonds and ref_keep in ref_keep_idxs):
bonds_to_make.append((nov_keep, ref_keep, order))
return bonds_to_make
def _set_property_on_assembled(st, ref_st, key):
"""
Set the given property on the assembled structure.
:type st: schrodinger.structure.Structure
:param st: the assembled structure
:type ref_st: schrodinger.structure.Structure
:param ref_st: the reference structure
:type key: str
:param key: the structure property key
"""
prop = st.property.get(key)
if prop is None:
prop = ref_st.property.get(key)
if prop is not None:
st.property[key] = prop
def _update_properties(st,
nov_st,
ref_st,
nov_old_to_new,
ref_old_to_new,
offset,
title=None):
"""
Update the properties of the given assembled structure.
:type st: schrodinger.structure.Structure
:param st: the assembled structure
:type nov_st: schrodinger.structure.Structure
:param nov_st: first structure, called novel
:type ref_st: schrodinger.structure.Structure
:param ref_st: second structure, called reference
:type nov_old_to_new: dict
:param nov_old_to_new: the old-to-new atom index map for the
extracted novel structure
:type ref_old_to_new: dict
:param ref_old_to_new: the old-to-new atom index map for the
extracted reference structure
:type offset: int
:param offset: an offset for reference indices given that the
assembled structure is a copy of part of the novel extended
by part of the reference
:type title: str
:param title: the title to be given to the assembled structure
"""
if not title:
title = '_'.join([nov_st.title, ref_st.title])
st.title = title
rxn_key = rxnwfu.REACTION_WF_STRUCTURE_KEY
if not nov_st.property.get(rxn_key) and not ref_st.property.get(rxn_key):
return
# the assembled structure has inherited structure properties
# from the novel and inherited atom properties from both the
# novel and reference
_set_property_on_assembled(st, ref_st, rxnwfu.CHARGE_KEY)
_set_property_on_assembled(st, ref_st, rxnwfu.MULTIPLICITY_KEY)
if ref_st.property.get(rxn_key):
conf_key = rxnwfu.CONFORMERS_GROUP_KEY
sibling_group_key = rxnwfu.SIBLING_GROUP_KEY
parent_sibling_groups_key = rxnwfu.PARENT_SIBLING_GROUPS_KEY
st.property[rxn_key] = True
st.property[conf_key] = ref_st.property[conf_key]
# allow missing sibling group key for backwards compatibility
st.property.pop(sibling_group_key, None)
sibling = ref_st.property.get(sibling_group_key)
if sibling:
st.property[sibling_group_key] = sibling
st.property.pop(parent_sibling_groups_key, None)
parents = ref_st.property.get(parent_sibling_groups_key)
if parents:
st.property[parent_sibling_groups_key] = parents
hierarchy = msutils.get_project_group_hierarchy(st=ref_st)
msutils.set_project_group_hierarchy(st, hierarchy)
ts_key = rxnwfu.TRANSITION_STATE_STRUCTURE_KEY
if nov_st.property.get(ts_key) or ref_st.property.get(ts_key):
st.property[ts_key] = True
# keep, superpose, and restrained atom indices are boolean atom
# properties and so they are automatically inherited
# handle restrain distances, angles, and dihedrals
keys = (rxnwfu.RESTRAINED_DISTANCES_KEY, rxnwfu.RESTRAINED_ANGLES_KEY,
rxnwfu.RESTRAINED_DIHEDRALS_KEY)
for key in keys:
st.property.pop(key, None)
idxs = rxnwfu._get_new_restrain_group_idxs(
nov_st, key, nov_old_to_new, offset=0)
idxs += rxnwfu._get_new_restrain_group_idxs(
ref_st, key, ref_old_to_new, offset=offset)
text = rxnwfu.get_idx_groups_str(idxs)
if text:
st.property[key] = text
[docs]class SwapFragmentsException(Exception):
pass
[docs]def get_assembled_structure(nov_st,
ref_st,
nov_superposition_idxs,
ref_superposition_idxs,
nov_keep_idxs,
ref_keep_idxs,
title=None,
require_identical_bonds=True):
"""
Return an assembled structure from the given structures
using superposition followed by extraction.
:type nov_st: schrodinger.structure.Structure
:param nov_st: first structure, called novel
:type ref_st: schrodinger.structure.Structure
:param ref_st: second structure, called reference
:type nov_superposition_idxs: list
:param nov_superposition_idxs: novel superposition indices
:type ref_superposition_idxs: list
:param ref_superposition_idxs: reference superposition indices
:type nov_keep_idxs: list
:param nov_keep_idxs: novel keep indices
:type ref_keep_idxs: list
:param ref_keep_idxs: reference keep indices
:type title: str
:param title: the title to be given to the assembled structure
:type require_identical_bonds: bool
:param require_identical_bonds: whether to require that bonds to be
created must exist in both novel and reference structures and
be of the same bond order
:raise SwapFragmentsException: if there is an issue
:rtype: schrodinger.structure.Structure
:return: the assembled structure
"""
if len(nov_superposition_idxs) != len(ref_superposition_idxs):
msg = ('The number of reference and novel superposition '
'atoms must be equivalent.')
raise SwapFragmentsException(msg)
rmsd.superimpose(
nov_st,
nov_superposition_idxs,
ref_st,
ref_superposition_idxs,
use_symmetry=False,
move_which=rmsd.CT)
nov_est, nov_old_to_new, nov_new_to_old = get_extracted_and_maps(
nov_st, nov_keep_idxs)
ref_est, ref_old_to_new, ref_new_to_old = get_extracted_and_maps(
ref_st, ref_keep_idxs)
old_idx_bonds = _get_bonds_for_assembling(
nov_st,
ref_st,
nov_keep_idxs,
ref_keep_idxs,
require_identical_bonds=require_identical_bonds)
new_idx_bonds = []
for old_nov_keep, old_ref_keep, order in old_idx_bonds:
new_nov_keep = nov_old_to_new[old_nov_keep]
new_ref_keep = ref_old_to_new[old_ref_keep] + nov_est.atom_total
new_idx_bonds.append((new_nov_keep, new_ref_keep, order))
st = nov_est.copy()
st.extend(ref_est)
offset = nov_est.atom_total
_update_properties(
st, nov_st, ref_st, nov_old_to_new, ref_old_to_new, offset, title=title)
for new_idx_bond in new_idx_bonds:
st.addBond(*new_idx_bond)
return st
[docs]def get_idxs_str(idxs, sort=True):
"""
Get a string representation of the given indices.
:type idxs: list
:param idxs: the idxs
:type sort: bool
:param sort: whether to sort
:rtype: str
:return: the string
"""
if sort:
jdxs = sorted(idxs)
else:
jdxs = list(idxs)
return rxnwfu.INDEX_SEPARATOR.join([str(j) for j in jdxs])
[docs]def get_idxs(le):
"""
Get indices from the given QLineEdit.
:type le: QtWidgets.QLineEdit
:param le: the line edit
:rtype: list
:return: the indices
"""
text = le.text().strip()
if not text:
return []
else:
return [int(s) for s in text.split(rxnwfu.INDEX_SEPARATOR) if s]