"""
AutoTS keywords input validation and specialized Exceptions
"""
# Contributors: Mark A. Watson
import collections
import schrodinger.application.jaguar.autots_input_constants as constants
import schrodinger.application.jaguar.utils as utils
import schrodinger.application.jaguar.workflow_validation as wv
from schrodinger.application.jaguar.autots_bonding import align_frozen_atoms
from schrodinger.application.jaguar.autots_bonding import coord_in_ring
from schrodinger.application.jaguar.autots_constants import AUTOTS_ATOM_CLASS
from schrodinger.application.jaguar.autots_exceptions import UnsupportedReaction
from schrodinger.application.jaguar.autots_rmsd import \
molecule_lists_are_conformers
from schrodinger.application.jaguar.exceptions import JaguarUserFacingException
from schrodinger.infra import mm
#------------------------------------------------------------------------------
[docs]class AutoTSAtomClassConflict(JaguarUserFacingException):
"""
Runtime Error indicating there is a problem with the atom
class definition
"""
[docs]class AutoTSInvalidConstraintError(JaguarUserFacingException):
"""
An invalid constraint
"""
[docs]def check_conflicts(kwd, all_keywords):
"""
Raise Exception if keyword value is inconsistent with the
other keywords. This is done in an adhoc case-by-case way.
:type kwd: AutoTSKeyword
:param kwd: reference AutoTS keyword
:type all_keywords: dict of AutoTSKeyword instances indexed by name
:param all_keywords: all the other AutoTS keywords set
:raise WorkflowKeywordConflictError if conflicting values found.
"""
# (1) skip_internal_mapping keyword conflict
if kwd.name == 'skip_internal_mapping' and kwd.value:
key = 'complex_formation'
val = constants.FORM_TYPE_SUPERIMPOSE
if all_keywords[key].value != val:
raise wv.WorkflowKeywordConflictError(kwd.name, key, val)
# (2) ts_acceptance_strictness of sloppy requires ts vetting of some kind
if kwd.name == 'ts_acceptance_strictness' and kwd.value == constants.STRICTNESS_SLOPPY:
key = 'ts_vetting_type'
val = constants.VET_NONE
if all_keywords[key].value == val:
raise wv.WorkflowKeywordConflictError(kwd.name, key, val)
# (3) jobtype JOB_TYPE_ENERGY_CORRECTION requires a full path file
if kwd.name == 'jobtype' and kwd.value == constants.JOB_TYPE_ENERGY_CORRECTION:
key = 'full_path_file'
val = '<name>'
if not all_keywords[key].value:
raise wv.WorkflowKeywordConflictError(kwd.name, key, val)
# (4) jobtype JOB_TYPE_CONF_SEARCH requires a full path file
if kwd.name == 'jobtype' and kwd.value == constants.JOB_TYPE_CONF_SEARCH:
key = 'full_path_file'
val = '<name>'
if not all_keywords[key].value:
raise wv.WorkflowKeywordConflictError(kwd.name, key, val)
# (5) add other known conflicts here...
if kwd.name == 'xxx' and kwd.value == 'yyy':
pass
return True
[docs]def reaction_is_supported(reactants, products):
"""
Check that the input reaction is supported, if not raise
an UnsupportedReaction.
"""
if not isinstance(reactants, list):
reactants = [reactants]
products = [products]
# must be a chemical change
if molecule_lists_are_conformers(reactants,
products,
same_molecularity=False):
msg = "This chemical reaction is not supported by AutoTS.\n"
msg += "The reactants and products are conformers.\n"
raise UnsupportedReaction(msg)
[docs]def matter_is_conserved(reactants, products):
"""
Check the input structures to ensure that matter is conserved.
If matter is not conserved a AutoTSConservationError is raised.
:type reactants: Structure object or iterable of Structure objects
:param reactants: the reactants or reactant complex
:type products: Structure object or iterable of Structure objects
:param products: the reactants or reactant complex
:rtype: bool
:return: whether or not matter is conserved
"""
rxn_stoich = ''
space = ' '
r_elmnts = []
if hasattr(reactants, "__iter__"):
for st in reactants:
elmnts = [at.element for at in st.atom if at.element]
r_elmnts += elmnts
rxn_stoich += utils.get_stoichiometry_string(elmnts) + space
else:
elmnts = [at.element for at in reactants.atom if at.element]
r_elmnts += elmnts
rxn_stoich += utils.get_stoichiometry_string(elmnts) + space
rxn_stoich += "-> "
p_elmnts = []
if hasattr(reactants, "__iter__"):
for st in products:
elmnts = [at.element for at in st.atom if at.element]
p_elmnts += elmnts
rxn_stoich += utils.get_stoichiometry_string(elmnts) + space
else:
elmnts = [at.element for at in products.atom if at.element]
p_elmnts += elmnts
rxn_stoich += utils.get_stoichiometry_string(elmnts) + space
isomers = True
r_elmnts.sort()
p_elmnts.sort()
if r_elmnts != p_elmnts:
msg = "Matter is not conserved:\n"
msg += "The number and/or type of atoms are not the same in the reactants and products.\n"
msg += "Please check the input structures and re-submit job.\n"
msg += "Input Reaction stoichiometry: %s\n" % rxn_stoich
raise wv.WorkflowConservationError(msg)
[docs]def validate_atom_classes(reactants, products):
"""
Check the input structures to ensure that if atom classes are defined
they are defined for an equal number of atoms in each class on each
side of the reaction. This is required to be able to compute an
atom map.
:type reactants: Structure object or iterable of Structure objects
:param reactants: the reactants or reactant complex
:type products: Structure object or iterable of Structure objects
:param products: the reactants or reactant complex
:rtype: bool
:return: whether or not atom classes are valid
"""
def _count_classes(st, counter):
for at in st.atom:
clss = at.property.get(AUTOTS_ATOM_CLASS, None)
if clss is not None:
counter[clss] += 1
class_counters = [collections.Counter(), collections.Counter()]
for counter, sts in zip(class_counters, [reactants, products]):
if hasattr(sts, "__iter__"):
for st in sts:
_count_classes(st, counter)
else:
_count_classes(sts, counter)
if class_counters[0] != class_counters[1]:
raise AutoTSAtomClassConflict(
"Atom classes are inconsistent between reactants and products")
[docs]def charge_is_conserved(reactants, products):
"""
Check the input structures to ensure that charge is conserved.
If charge is not conserved a AutoTSConservationError is raised.
:type reactants: Structure object or iterable of Structure objects
:param reactants: the reactants or reactant complex
:type products: Structure object or iterable of Structure objects
:param products: the reactants or reactant complex
:rtype: bool
:return: whether or not matter is conserved
"""
r_num_electrons = wv._calculate_num_electrons(reactants)
p_num_electrons = wv._calculate_num_electrons(products)
if r_num_electrons != p_num_electrons:
msg = "Charge is not conserved:\n"
msg += "The number of electrons is not the same in the reactants and products.\n\n"
msg += "Reactants: \n"
msg += get_stoich_charge(reactants)
msg += "Products: \n"
msg += get_stoich_charge(products)
msg += '\n'
msg += "Please check the input structures and re-submit job.\n"
raise wv.WorkflowConservationError(msg)
[docs]def get_stoich_charge(sts):
"""
Return the stoichiometry and charge of a list of structures or structure as a string
:type sts: list of Structures or Structure instance
:param sts: structures to print stoichiometry and charges of
:rtype: string
:return: the stoichiometry and charge of the input structure(s)
"""
output = ''
if hasattr(sts, "__iter__"):
for st in sts:
elmnts = [at.element for at in st.atom if at.element]
output += f"{utils.get_stoichiometry_string(elmnts)}, charge: {utils.get_total_charge(st)}\n"
else:
elmnts = [at.element for at in sts.atom if at.element]
output += f"{utils.get_stoichiometry_string(elmnts)}, charge: {utils.get_total_charge(sts)}\n"
return output
[docs]def validate_jag_keys(jaguar_keywords):
"""
raise an exception if tmpini, press, or pbf_after_pcm are keywords
:type jaguar_keywords: dict
:param jaguar_keywords: key: value dict of jaguar keywords
"""
lower_keys = [strng.casefold() for strng in jaguar_keywords]
if mm.MMJAG_RKEY_PRESS.casefold(
) in lower_keys or mm.MMJAG_RKEY_TMPINI.casefold() in lower_keys:
msg = "AutoTS requires that pressure and temperature are not explicitly set"
raise wv.JaguarKeywordConflict(msg)
if mm.MMJAG_IKEY_PBF_AFTER_PCM.casefold() in lower_keys:
msg = "pbf_after_pcm keyword is incompatible with AutoTS. Please remove."
raise wv.JaguarKeywordConflict(msg)
[docs]def validate_constraints(rinp):
"""
Validate any constraints that have been read.
This ensures that the constraints can be mapped
onto the input reactant/product molecules.
Specifically the properties that are required are:
1. relax_path = True cannot be set if we have constraints
2. The input reactants and products must have unique titles
3. Each constraint must map to some titled molecule or
COMPLEX_STRUCTURE_CONSTRAINT if both reactant and product complexes
are defined
4. Atom indexes of the constraint must be between 1 and the number of
atoms in that molecule.
5. constraints cannot be in rings
6. Atom indexes for complex constraints cannot span multiple molecules.
Unless it's an atomic constraint (frozen atom) in which case it's the wild west
"""
have_internal_constraints = rinp.have_constraints()
have_frozen_atoms = rinp.have_frozen_atoms()
if have_internal_constraints or have_frozen_atoms:
if have_internal_constraints and have_frozen_atoms:
msg = "Cannot have both internal constraints and frozen atoms"
raise AutoTSInvalidConstraintError(msg)
# cannot relax path with constraints
if rinp.getValue("relax_path"):
raise wv.WorkflowKeywordConflictError("constraints", "relax_path",
"False")
# cannot do conf search with constraints
if rinp.getValue("conf_search"):
raise wv.WorkflowKeywordConflictError("constraints", "conf_search",
"False")
# cannot do jobtype = conf_search with constraints
if rinp.getValue("jobtype") == constants.JOB_TYPE_CONF_SEARCH:
choices = list(constants.JOB_TYPE_CHOICES).remove(
constants.JOB_TYPE_CONF_SEARCH)
raise wv.WorkflowKeywordConflictError("constraints", "jobtype",
choices)
# all structures must have unique titles
all_strs = rinp.getReactants() + rinp.getProducts()
strs = {st.title: st for st in all_strs}
reactant_complex = rinp.getReactantComplex()
if len(strs) != len(all_strs) or None in strs:
msg = "All structures must have unique titles when using constraints"
raise JaguarUserFacingException(msg)
all_titles = list(rinp._constraints.keys()) + list(
rinp._frozen_atoms.keys())
for title in all_titles:
# each constraint should map to an input structure
# or to the constant COMPLEX_STRUCTURE_CONSTRAINT iff
# the reactant and product complexes are defined.
st = strs.get(title)
if st is None:
if reactant_complex is not None and \
title == constants.COMPLEX_STRUCTURE_CONSTRAINTS:
st = rinp.getReactantComplex()
else:
msg = "title for constraint (%s) does not correspond " \
"to an input reactant or product." % title
raise AutoTSInvalidConstraintError(msg)
for constraint in rinp._constraints.get(
title, []) + rinp._frozen_atoms.get(title, []):
# atom indexes of constraints should be in range
for indx in constraint.indices:
if indx < 1 or indx > st.atom_total:
msg = "Atom indexes for constraint are out of range:\n" \
" %s" % (title + str(constraint))
raise AutoTSInvalidConstraintError(msg)
if len(constraint.indices) > 1 and coord_in_ring(
st, *constraint.indices):
msg = "Atom indexes of constraint can not be in a ring:\n" \
" %s" % (title + str(constraint))
raise AutoTSInvalidConstraintError(msg)
# ensure constraints don't span several molecules if complexes are defined
complex_constraints = rinp._constraints.get(
constants.COMPLEX_STRUCTURE_CONSTRAINTS, [])
have_complexes = None not in [
rinp.getReactantComplex(),
rinp.getProductComplex()
]
if have_complexes and rinp.values.use_complex_inputs and have_internal_constraints and not complex_constraints:
raise AutoTSInvalidConstraintError(
"if using reactant and product complexes, complex constraints must be specified"
)
if complex_constraints:
r_mol_atoms = [
mol.getAtomIndices() for mol in rinp.getReactantComplex().molecule
]
p_mol_atoms = [
mol.getAtomIndices() for mol in rinp.getProductComplex().molecule
]
r = rinp.getReactantComplex()
p = rinp.getProductComplex()
for constraint in complex_constraints:
r_ats, p_ats = set(), set()
for i in constraint.indices:
r_ats.update(idx for idx, molats in enumerate(r_mol_atoms)
if i in molats)
p_ats.update(idx for idx, molats in enumerate(p_mol_atoms)
if i in molats)
if len(r_ats) > 1 or len(p_ats) > 1:
raise AutoTSInvalidConstraintError(
"Constraints cannot span multiple molecules")
complex_frozen_atoms = rinp.getFrozenAtoms().get(
constants.COMPLEX_STRUCTURE_CONSTRAINTS, [])
if have_complexes and rinp.values.use_complex_inputs and have_frozen_atoms and not complex_frozen_atoms:
raise AutoTSInvalidConstraintError(
"if using reactant and product complexes, complex constraints must be specified"
)
# complexes must be defined if using complex frozen atoms
if complex_frozen_atoms or complex_constraints:
if rinp.getReactantComplex() is None or rinp.getProductComplex(
) is None:
raise AutoTSInvalidConstraintError(
"Must supply a reactant and product complex if using complexes to define frozen atom constraints"
)
# must have single entries if not using complex frozen atoms
if have_frozen_atoms and not complex_frozen_atoms:
nr = len(rinp.getReactants())
np = len(rinp.getProducts())
if nr != 1 and np != 1:
raise AutoTSInvalidConstraintError(
"Frozen atom constraints require either complexes or single entry reactant and products"
)
# set reference values and check initial satisfaction for complex frozen atoms
# for non-complex frozen atoms this is done in renumber.py
if complex_frozen_atoms:
r = rinp.getReactantComplex()
p = rinp.getProductComplex()
for coord in complex_frozen_atoms:
coord.setValue(r)
rmsd_p = align_frozen_atoms(st, complex_frozen_atoms)
if rmsd_p > rinp.getValue("frozen_atoms_rmsd_thresh"):
raise JaguarUserFacingException(
"Frozen atom RMSD (%.4f) is not below a tolerance of %.4f Angstrom"
% (rmsd_p, rinp.getValue("frozen_atoms_rmsd_thresh")))
elif rmsd_p > constants.FROZEN_ATOM_WARN_THRESH:
print(
f"WARNING: Frozen atom RMS deviation ({rmsd_p:.3f}) is non-trivial"
)
[docs]def validate_structures(rinp):
"""
Perform a check to ensure that matter is conserved and that
charge/multiplicity are consistent with structures.
A WorkflowConservationError is raised if any test fails.
"""
charge = rinp.getValue("charge")
mult = rinp.getValue("multiplicity")
have_reaction = False
have_r_and_p = True
for r, p in [(rinp.getReactants(), rinp.getProducts()),
(rinp.getReactantComplex(), rinp.getProductComplex())]:
if not have_reaction and not (r and p):
have_r_and_p = False
if r and p:
_validate_reaction(r, p, charge, mult)
validate_atom_classes(r, p)
have_reaction = True
if not have_r_and_p:
if rinp.getReactants():
msg = 'Products have not been specified. ' \
'Please set "product" keyword in the AutoTS .in file'
elif rinp.getProducts():
msg = 'Reactants have not been specified. ' \
'Please set "reactant" keyword in the AutoTS .in file'
else:
msg = 'Neither reactants nor products have been specified. ' \
'Please set "reactant" and "product" keywords in the AutoTS .in file'
raise wv.WorkflowConservationError(msg)
# DSL: I don't believe it is possible to hit this because an exception will be
# raised above in any situation in which have_reaction is not True.
if not have_reaction:
msg = "The reaction must be defined by supplying reactant and product molecules"
raise wv.WorkflowConservationError(msg)
def _validate_reaction(r, p, charge, mult):
"""
Validate a reaction specification.
will raise a WorkflowConservationError if the specification is invalid.
"""
err_msg = ''
try:
matter_is_conserved(r, p)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + '\n'
try:
charge_is_conserved(r, p)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + '\n'
try:
reaction_is_supported(r, p)
except UnsupportedReaction as e:
err_msg += e.args[0] + '\n'
try:
wv.estate_is_physical(r, charge, mult)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + ' in the reactants\n'
try:
wv.charge_is_consistent(r, charge)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + ' in the reactants\n'
try:
wv.estate_is_physical(p, charge, mult)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + ' in the products\n'
try:
wv.charge_is_consistent(p, charge)
except wv.WorkflowConservationError as e:
err_msg += e.args[0] + ' in the products\n'
if err_msg:
raise wv.WorkflowConservationError(err_msg)