""" Methods to handle structure reading and bonding in AutoTS. """
# Contributors: Leif D. Jacobson
import collections
import numpy as np
from schrodinger.application.jaguar import autots_constants as constants
from schrodinger.application.jaguar import utils
from schrodinger.application.matsci.buildcomplex import fix_metal_bond_orders
from schrodinger.infra import mm
from schrodinger.infra.mmerr import ErrorHandler
from schrodinger.structure import AtomsInRingError
from schrodinger.structure import StructureReader
from schrodinger.structure import create_new_structure
from schrodinger.structutils import analyze
from schrodinger.structutils import build
from schrodinger.structutils.rmsd import CT
from schrodinger.structutils.rmsd import superimpose
[docs]def clean_st(st, reset_bonding=True):
"""
Clean up a st via redefining bonding.
We also delete formal charges because they get in the way of the SMARTS
pattern based matching used in many places.
:type st: schrodinger.structure.Structure instance
:param st: structure to clean
:type reset_bonding: boolean
:param reset_bonding: recompute bonding with mmjag
:rtype st: schrodinger.structure.Structure instance
:return: the cleaned up structure
"""
# no need for bonding if there is only one atom
if st.atom_total == 1:
return st.copy()
# reset the connectivity and Lewis structure using mmjag
st_out = st.copy()
if reset_bonding:
utils.mmjag_reset_connectivity(st_out)
# remove ghost atoms and bond metals to non-metals connected by ghosts
st_out = delete_ghosts(st_out)
get_mmlewis_bonding(st_out)
zero_order_metal_bonds(st_out)
remove_agostic_bonds(st_out)
canonicalize_atom_names(st_out)
return st_out
[docs]def canonicalize_atom_names(st):
"""
Canonicalize the atom names
:type st: Structure
:param st: Structure containing atoms to name
"""
cnt = collections.Counter()
for at in st.atom:
cnt[at.element] += 1
at.name = at.element + str(cnt[at.element])
[docs]class AutoTSStructureReader():
"""
Local version of StructureReader which redefines bonding
when reading
"""
[docs] def __init__(self, *args, reset_bonding=True, **kwargs):
self.reset_bonding = reset_bonding
self.reader = StructureReader(*args, **kwargs)
def __iter__(self):
for x in self.reader:
yield clean_st(x, reset_bonding=self.reset_bonding)
def __next__(self):
return clean_st(next(self.reader), reset_bonding=self.reset_bonding)
[docs]def copy_autots_atom_properties(st1, st2):
"""
copy all known atom-level AutoTS-specific properties
from st1 to st2. Atoms in st1 and st2 must be in
the same order.
"""
for at in st1.atom:
for key in constants.AUTOTS_ATOM_PROPERTIES:
if key in list(at.property):
st2.atom[at.index].property[key] = at.property[key]
[docs]def copy_autots_st_properties(st1, st2):
"""
copy all known structure-level AutoTS-specific properties
from st1 to st2.
"""
for key in constants.AUTOTS_ST_PROPERTIES:
if key in list(st1.property):
st2.property[key] = st1.property[key]
# free energy is special as it has a
# temperature which is explicit in the property key
gibbs_data = utils.parse_gibbs_energies(st1, inf_sep=False)
for g in gibbs_data.values():
st2.property[g.property_key] = g.gibbs
inf_sep_gibbs_data = utils.parse_gibbs_energies(st1, inf_sep=True)
for g in inf_sep_gibbs_data.values():
st2.property[g.property_key] = g.gibbs
[docs]def copy_autots_bond_properties(st1, st2):
"""
copy all known bond-level AutoTS-specific properties.
Atoms in st1 and st2 must be in the same order but the
bonds do not. However, a bond with the same atom indexes
and bond order must exist in order to copy the properties.
"""
for bond in st1.bond:
if st2.areBound(bond.atom1, bond.atom2):
bond2 = st2.getBond(bond.atom1, bond.atom2)
for key in constants.AUTOTS_BOND_PROPERTIES:
if key in list(bond.property):
bond2.property[key] = bond.property[key]
[docs]def copy_autots_properties(st1, st2):
"""
copy all known atom-, bond- and st-level AutoTS-specific properties from
st1 to st2. Atoms in st1 and st2 must be in the same order.
"""
copy_autots_atom_properties(st1, st2)
copy_autots_bond_properties(st1, st2)
copy_autots_st_properties(st1, st2)
[docs]def clear_autots_atom_properties(st):
"""
clear all known atom-level AutoTS-specific properties.
"""
for at in st.atom:
for key in constants.AUTOTS_ATOM_PROPERTIES:
try:
at.property.pop(key)
except (ValueError, KeyError):
pass
[docs]def clear_autots_st_properties(st,
exceptions=(constants.PROPERTY_KEY_CHARGE,
constants.PROPERTY_KEY_MULT)):
"""
Copy all known structure-level AutoTS-specific properties.
Any properties in exceptions are not cleared (charge and mult by default)
"""
for key in constants.AUTOTS_ST_PROPERTIES:
if key not in exceptions:
try:
st.property.pop(key)
except (ValueError, KeyError):
pass
[docs]def clear_autots_bond_properties(st):
"""
clear all known bond-level AutoTS-specific properties.
"""
for bond in st.bond:
for key in constants.AUTOTS_BOND_PROPERTIES:
try:
bond.property.pop(key)
except (ValueError, KeyError):
pass
[docs]def clear_autots_properties(st):
"""
Remove all known atom-, bond- and st-level AutoTS specific properties.
"""
clear_autots_atom_properties(st)
clear_autots_bond_properties(st)
clear_autots_st_properties(st)
[docs]def delete_ghosts(st):
"""
Delete ghost/dummy atoms from the structure and
return the modified structure.
Ghosts are identified by the attribute _StructureAtom.atomic_number = 0.
The original structure is unmodified.
:type st: Structure instance
:param st: the structure with the ghost atoms
:rtype: Structure instance
:return: the structure without the ghost atoms
"""
st_out = st.copy()
nats = st.atom_total
reorder = False
key = 'i_j_delete_ghosts_index'
for at in st_out.atom:
at.property[key] = at.index
ghosts = []
for at in st_out.atom:
if at.atomic_number == 0:
ghosts.append(at.index)
if ghosts:
_replace_ghost_bonding(st_out, ghosts)
st_out.deleteAtoms(ghosts)
# we ensure the atoms are in the same order as they were input,
# just with the ghosts removed
pairs = [[at.property.get(key), at.index] for at in st_out.atom]
pairs.sort()
atlist = [j for i, j in pairs]
st_out = build.reorder_atoms(st_out, atlist)
# remove temporary property
for at in st_out.atom:
at.property.pop(key)
return st_out
def _replace_ghost_bonding(st, ghosts):
r"""
Replace the bonding to ghost atoms for any bond to transition metals (TM).
The idea is to replace::
C
/|
/ |
M-Gh |
\ |
\|
C
with::
C
/|
/ |
M |
\ |
\|
C
We only add bonds between the metal and the other atoms the ghost is bonded to.
We'll also add bonds if the Ghost is only bonded to two atoms, i.e.::
A-Gh-A -> A-A
All added bonds are zero order.
:type st: Structure instance
:param st: structure to modify bonding for
:type ghosts: list of ints
:param ghosts: list of indexes for ghost atoms
"""
for indx in ghosts:
gh = st.atom[indx]
bonded_atoms = [at for at in gh.bonded_atoms]
nats = len(bonded_atoms)
if nats == 2:
if not st.areBound(*bonded_atoms):
st.addBond(bonded_atoms[0], bonded_atoms[1], 0)
elif nats > 2:
metal_indices = analyze.evaluate_asl(st, 'metals')
metal_atoms = []
for at in bonded_atoms:
if at.index in metal_indices:
metal_atoms.append(at)
if len(metal_atoms) == 1:
metal_atom = metal_atoms[0]
remaining = list(set(bonded_atoms).difference(set(metal_atoms)))
# bond metal to all remaining
for at in remaining:
if not st.areBound(at, metal_atom):
st.addBond(at, metal_atom, 0)
[docs]def get_mmlewis_bonding(st, require_charge_conservation=True, debug=False):
"""
Get bonding from mmlewis. Do it molecule by molecule.
:type st: Structure instance
:param st: the structure to get bond orders for (must be connected)
:type require_charge_conservation: boolean
:param require_charge_conservation: if True we require the sum of formal
charges after running through mmlewis to equal the total charge.
(defined either by PROPERTY_KEY_CHARGE or by the sum of formal charges).
"""
st_copy = st.copy()
q = utils.get_total_charge(st)
for mol in st.molecule:
# extract structure and map atoms to original structure
mol_st = mol.extractStructure()
mol_map = mol.getAtomIndices()
# get bond orders
_get_mmlewis_bonding(mol_st, debug)
# copy over the bond orders
for bond in mol_st.bond:
i = mol_map[bond.atom1.index - 1]
j = mol_map[bond.atom2.index - 1]
st_bond = st.getBond(i, j)
st_bond.order = bond.order
# copy over formal charges
for at in mol_st.atom:
i = mol_map[at.index - 1]
st.atom[i].formal_charge = at.formal_charge
qnew = sum([at.formal_charge for at in st.atom])
if require_charge_conservation and q != qnew:
# revert to initial bond orders and formal charges
for bond in st_copy.bond:
st_bond = st.getBond(bond.atom1.index, bond.atom2.index)
st_bond.order = bond.order
for at in st_copy.atom:
st.atom[at.index].formal_charge = at.formal_charge
def _get_mmlewis_bonding(st, debug=False):
"""
get bonding from mmlewis
:type st: Structure instance
:param st: the structure to get bond orders for (must be connected)
"""
# hides output from mmlewis
mm.mmlewis_initialize(mm.MMERR_DEFAULT_HANDLER)
old_handler = mm.mmlewis_get_error_handler()
new_handler = ErrorHandler(queued=True, silent=True)
new_handler.push_level(6)
mm.mmlewis_set_error_handler(new_handler.handle)
# see function mmlewis_from_ct_Z in mmlewis.cxx for a description of IFO values
# 0 means define bond orders and charges regardless of validity of input.
# 5 means define bond orders and charges only if structure is invalid.
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_IFO, 5, 0.0, '')
# not verbose
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_VERBOSE, 0, 0.0, '')
# according to MATSCI-1838 this is better for nano-tubes
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_BOND_PLACEMENT, \
mm.MMLEWIS_BOND_PLACEMENT_BFS, 0.0, '')
# don't want dummys
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_MMLEWIS_APPLY_DELETE_DUMMIES, 1,
0.0, '')
try:
mm.mmlewis_apply(st)
except Exception as e:
if debug:
msg = "Warning: mmlewis issued error message: %s\n" % e
print(msg)
mm.mmlewis_set_error_handler(old_handler)
mm.mmlewis_terminate()
[docs]def simplify_structure(st):
"""
Make all bonds single bonds and remove all charges.
:type st: schrodinger.structure.Structure
:param st: a structure
"""
for bond in st.bond:
bond.order = 1
remove_formal_charges(st)
st.property.pop(constants.PROPERTY_KEY_CHARGE, None)
st.retype()
[docs]def active_reactant_atom_pairs(reactant, product):
"""
return active atom pairs in reactant structure
as a list of pairs of integers
"""
r_active_pairs = []
for bond in reactant.bond:
if not product.areBound(product.atom[bond.atom1.index],
product.atom[bond.atom2.index]):
r_active_pairs.append((bond.atom1.index, bond.atom2.index))
return r_active_pairs
[docs]def active_atom_pairs(reactant, product):
"""
Determine active bonds and return them as lists of pairs of atoms
"""
active_pairs = active_reactant_atom_pairs(reactant, product) \
+ active_reactant_atom_pairs(product, reactant)
return active_pairs
[docs]def copy_bonding(st1, st2):
"""
Impose the bonding and formal charges of st1 onto st2
The two structures must have the same
number of atoms or a ValueError is raised.
:type st1: Structure
:type st2: Structure
"""
build.delete_bonds(st2.bond)
bond_list = [
(bond.atom1.index, bond.atom2.index, bond.order) for bond in st1.bond
]
st2.addBonds(bond_list)
for iat, jat in zip(st1.atom, st2.atom):
jat.formal_charge = iat.formal_charge
[docs]class Coordinate(object):
"""
An internal coordinate.
The value and indexes are stored as data "value" and "indices".
"""
[docs] def __init__(self, value, *args):
"""
Create an internal coordinate with a value
:type value: float or np.ndarray
:param value: value of coordinate, if one index (atom) then
this must be an ndarray with leading dimension 3
:type args: tuple
:param args: atom indexes defining constraint
Example usage
torsion = Coordinate(91.2, 4, 5, 8, 12)
bond = Coordinate(1.2, 12, 14)
atom = Coordinate(st.atom[1].xyz, 1)
"""
self.value = value
self.indices = tuple(args)
if len(self.indices) == 1:
assert isinstance(self.value, np.ndarray)
assert self.value.shape[0] == 3
def __str__(self):
indices_str = [str(i) for i in self.indices]
strng = {
1: "Atomic position: ",
2: "Distance: ",
3: "Angle: ",
4: "Torsion: "
}[len(self.indices)]
if len(self.indices) == 1:
return strng + " %d" % self.indices[0]
else:
return strng + " ".join(indices_str) + " %.4f" % self.value
def __eq__(self, other):
same = False
if sorted(self.indices) == sorted(other.indices):
if len(self.indices) == 1:
same = all(self.value[i] == other.value[i] for i in range(3))
else:
same = self.value == other.value
return same
def __hash__(self):
return hash(self.indices) + hash(self.value)
[docs] def similar_coordinate(self, other):
"""
A similar coordinate is one that describes the
same degree of freedom. For example, the
torsion 1 2 3 4 and torsion 1 2 3 6 describe
rotation about the same bond.
:type other: Coordinate
:param other: the other coordinate for comparison
:return: True if the other coordinate is similar, else False
"""
similar = len(self.indices) == len(other.indices)
if similar:
# if torsion we only compare the indices of the rotatable bond
is_torsion = len(self.indices) == 4
my_index = self.indices[1:3] if is_torsion else self.indices
other_index = other.indices[1:3] if is_torsion else other.indices
similar = sorted(my_index) == sorted(other_index)
return similar
[docs] def adjust(self, st):
"""
adjust this coordinate for a structure
If coordinate is in a ring no adjustment is made
"""
if len(self.indices) == 1:
assert self.value.shape[0] == 3
st.atom[self.indices[0]].x = self.value[0]
st.atom[self.indices[0]].y = self.value[1]
st.atom[self.indices[0]].z = self.value[2]
else:
try:
st.adjust(self.value, *self.indices)
except AtomsInRingError:
pass
[docs] def setValue(self, st):
"""
Set value of coordinate using structure
"""
self.value = self.getValue(st)
[docs] def getValue(self, st):
"""
return current value of coordinate
"""
if len(self.indices) == 1:
value = np.array(st.atom[self.indices[0]].xyz)
else:
value = st.measure(*self.indices)
return value
[docs] def getDifference(self, st):
"""
return difference between value of coordinate and Structure value
"""
value = self.getValue(st)
if len(self.indices) == 1:
dv = value - self.value
err = np.sqrt(np.dot(dv, dv))
else:
err = abs(self.value - value)
return err
[docs]def examine_constraints(constraints, frozen_atoms, st, tol=0.01, enforce=True):
"""
Examine the satisfaction of constraints.
Returns a list of Coordinate instances representing
the error in any constraints which are not satisfied to
a specified tolerance.
:type constraints: list of Coordinate
:param constraints: the constraints
:type frozen_atoms: list of Coordinate
:param frozen_atoms: list of frozen atoms
:type st: Structure
:param st: structure which should satisfy constraints
:type tol: float
:param tol: tolerance
:type enforce: boolean
:param enforce: if True attempt to enforce constraints with Structure.adjust
"""
unsatisfied_constraints = []
# enforce all constraints
if enforce:
for constraint in constraints:
assert len(constraint.indices) > 1
try:
st.adjust(constraint.value, *constraint.indices)
except AtomsInRingError:
pass
if frozen_atoms:
_ = align_frozen_atoms(st, frozen_atoms)
for constraint in constraints:
err = constraint.getDifference(st)
if err > tol:
unsatisfied_constraints.append(Coordinate(err, *constraint.indices))
for constraint in frozen_atoms:
err = constraint.getDifference(st)
if err > tol * 0.01:
dx = constraint.getValue(st) - constraint.value
unsatisfied_constraints.append(Coordinate(dx, *constraint.indices))
return unsatisfied_constraints
[docs]def align_frozen_atoms(st, frozen_atoms):
"""
Rotate structure to align frozen atoms and return the RMSD
to the constraint values. After this analysis the atoms are moved to be
in exact agreement with the reference positions. It is assumed that
constraints of xyz positions of input are copied to the list of constraints.
This is done in renumber.map_constraints_to_complexes
:type st: Structure
:param st: Structure to enforce constraints on
:type frozen_atoms: list of Coordinate
:param frozen_atoms: the constraints
:return: float giving RMSD from constraints
"""
rmsd = 0.0
if frozen_atoms:
mock_atlist = []
atlist = []
mock_st = create_new_structure()
for icoord, coord in enumerate(frozen_atoms, start=1):
mock_st.addAtom(st.atom[coord.indices[0]].element, *coord.value)
atlist.append(coord.indices[0])
mock_atlist.append(icoord)
rmsd = superimpose(mock_st, mock_atlist, st, atlist, move_which=CT)
# explicitly enforce
for coord in frozen_atoms:
st.atom[coord.indices[0]].x = coord.value[0]
st.atom[coord.indices[0]].y = coord.value[1]
st.atom[coord.indices[0]].z = coord.value[2]
return rmsd
[docs]def get_static_constraints(constraints):
"""
Convert a list of constraints into static constraints.
This just means returning a list of internal coordinate constraints
with the value set to None
:type constraints: list
:param constraints: list of Coordinate instances
:return: a list of Coordinates
"""
constraints = [
Coordinate(None, *c.indices) for c in constraints if len(c.indices) > 1
]
return constraints
[docs]def atoms_in_ring(st, i, j):
"""
Determine if two atoms are in the same ring
:type st: Structure
:param st: structure
:type i: int
:param i: atom index
:type j: int
:param j: atom index
"""
return any(i in ring and j in ring for ring in st.find_rings())
[docs]def coord_in_ring(st, *args):
"""
Determine whether or not an internal coordinate is in a ring.
For a bond: returns True if both atoms are in the same ring.
For an angle: returns True if all three atoms are in the same ring.
For a torsion: returns True if the two central atoms are in the same ring.
:type st: Structure
:param st: Structure used to test if indexes are in the same ring
:type args: tuple/list
:param args: indexes of the coordinate.
"""
# for a bond both atoms are in the same structure
if len(args) == 2:
in_ring = atoms_in_ring(st, *args)
# for angles all three atoms must be in the same ring.
elif len(args) == 3:
in_ring = atoms_in_ring(st, *args[0:2]) and atoms_in_ring(
st, *args[1:3]) and atoms_in_ring(st, args[0], args[2])
# for torsions the central (rotatable) bond must be in the same ring
elif len(args) == 4:
in_ring = atoms_in_ring(st, *args[1:3])
else:
in_ring = False
return in_ring
[docs]def remove_agostic_bonds(st):
"""
Remove zero-order bonds from metals to valence saturated
C and H atom.
These agostic interactions are very weak and easily broken
and their presence or absence can mess up codes that use
connectivity.
:param st: Structure to remove agostic bonds
:type st: Structure instance
"""
metal_indices = set(analyze.evaluate_asl(st, 'metals'))
for bond in list(st.bond):
if bond.order > 0:
continue
bond_ats = {bond.atom1.index, bond.atom2.index}
if len(bond_ats.intersection(metal_indices)) == 1:
other_at = st.atom[(bond_ats - metal_indices).pop()]
if other_at.atomic_number in {1, 6}:
if bound_hydrogen(bond, metal_indices):
continue
n_neighbors = sum(1 for at in st.atom[other_at].bonded_atoms
if at.index not in metal_indices)
if other_at.atomic_number==1 and n_neighbors==1 or \
other_at.atomic_number==6 and n_neighbors==4:
st.deleteBond(bond.atom1.index, bond.atom2.index)
[docs]def bound_hydrogen(bond, metal_indices):
"""
Returns whether a hydrogen atom bound to a metal is a bound H2 molecule
:param bond: the H-metal bond of interest
:type bond: StructureBond
:param metal_indices: set of atom indexes which are metals
:type metal_indices: set of ints
:return param: True if H-metal bond is part of bound H2 molecule
:return type: boolean
"""
bound_H2 = False
if bond.atom1.atomic_number == 1:
the_h = bond.atom1
else:
the_h = bond.atom2
neighbors = [
atom for atom in the_h.bonded_atoms if atom.index not in metal_indices
]
if len(neighbors) == 1 and \
neighbors[0].atomic_number == 1 and \
len([atom for atom in neighbors[0].bonded_atoms if atom.index not in metal_indices]) == 1:
bound_H2 = True
return bound_H2