"""
Classes and functions for generating reaction paths.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
import argparse
import logging
import math
import os
import sys
from collections import OrderedDict
from operator import itemgetter
from past.utils import old_div
import numpy
from scipy.optimize import leastsq
import schrodinger.application.jaguar.input as jaginp
import schrodinger.structure as structure
from schrodinger.application.matsci import parserutils
from schrodinger.application.matsci import rxn_channel
from schrodinger.infra import mm
from schrodinger.infra.mmcheck import MmException
from schrodinger.structutils import analyze
from schrodinger.structutils import build
from schrodinger.structutils import measure
from schrodinger.structutils import rmsd
from schrodinger.structutils import transform
from schrodinger.utils import fileutils
from .msutils import add_or_update_bond_type
_version = '$Revision 0.0 $'
[docs]class ParserWrapper(object):
"""
Manages the argparse module to parse user command line arguments.
"""
JOB_NAME = 'rxnpath'
REACTANT = 'reactant'
PRODUCT = 'product'
TS = 'ts'
REACTANT_LIKE = 'reactant_like'
MIDWAY = 'midway'
PRODUCT_LIKE = 'product_like'
PRESUMED_TS = MIDWAY
PRESUMED_TS_CHOICES = [REACTANT_LIKE, MIDWAY, PRODUCT_LIKE]
DENSEAROUND = False
SAMPLE_DEFAULT = [10.0]
SUPPORTEDINEXTS = ['.mae', '.mae.gz', '.maegz']
FVAL_KEYS = [REACTANT, REACTANT_LIKE, MIDWAY, PRODUCT_LIKE, \
PRODUCT]
FVAL_VALUES = [0.00, 0.25, 0.50, 0.75, 1.00]
FVAL_DICT = OrderedDict(list(zip(FVAL_KEYS, FVAL_VALUES)))
NUMDENSEPOINTS = 10
STEPDENSEPOINTS = 0.02
BONDWEIGHT = 1000.0
ANGLEWEIGHT = 1000.0
DIHEDRALWEIGHT = 1000.0
CARTESIANWEIGHT = 1000.0
PENALTYWEIGHT = 1.0
MIXPREVIOUS = 0.5
MIXPREVIOUSMIN = 0.0
MIXPREVIOUSMAX = 1.0
CARTESIAN = 'cartesian'
DISTANCE = 'distance'
INTERNAL = 'internal'
INTERPOLATIONCHOICES = [INTERNAL, DISTANCE, CARTESIAN]
BEFORESUPERPOSITION = 'beforesuperposition'
AFTERSUPERPOSITION = 'aftersuperposition'
GUESSCHOICES = [BEFORESUPERPOSITION, AFTERSUPERPOSITION]
CONNECTIVITYCHOICES = [REACTANT, TS, PRODUCT]
NORXNCOMPLEX = False
VDWSCALE = 1.0
REORDER = False
REVERSE_INTERPOLATION = False
[docs] def __init__(self, scriptname, description):
"""
Create a ParserWrapper object and process it.
:type scriptname: str
:param scriptname: name of this script
:type description: str
:param description: description of this script
"""
name = '$SCHRODINGER/run ' + scriptname
self.parserobj = parserutils.DriverParser(
prog=name,
description=description,
add_help=False,
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
[docs] def loadIt(self):
"""
Load ParserWrapper with options.
"""
job_name_help = """
Specify the job name of this calculation."""
self.parserobj.add_argument('-job_name',
type=str,
default=self.JOB_NAME,
help=job_name_help)
reorderhelp = """
Attempts to automatically determine the appropriate mapping of
reactant atoms to product atoms."""
self.parserobj.add_argument('-reorder',
action='store_true',
help=reorderhelp)
norxncomplexhelp = """
Disable the preprocessing of reactants and products into reaction
complexes for certain bimolecular reactions."""
self.parserobj.add_argument('-norxncomplex',
action='store_true',
help=norxncomplexhelp)
vdwscalehelp = """
Scales the inter-molecular VDW distance used to separate reactants
or products when forming reaction complexes. This option will not
used if the user chooses the option '-norxncomplex'."""
self.parserobj.add_argument('-vdwscale',
type=float,
default=self.VDWSCALE,
help=vdwscalehelp)
reverse_interpolation_help = """
Specify that the reaction path should be interpolated in reverse,
i.e. the reactant and product definitions do not change
however the interpolation is done from products to reactants rather
than from reactants to products."""
self.parserobj.add_argument('-reverse_interpolation',
action='store_true',
help=reverse_interpolation_help)
samplehelp = """
Specify the type of sampling used in generating the reaction
path. For uniform sampling between the reactant and product
pass N >= 1, where N is the number of points to sample. For
non-uniform sampling pass a whitespace delimited list of
numbers N, where %.1f (reactant) < N < %.1f (product), for
example %.1f might resemble the transition state."""
self.parserobj.add_argument(
'-sample',
nargs='+',
type=float,
default=self.SAMPLE_DEFAULT,
metavar='list of floats',
help=samplehelp %
(self.FVAL_DICT[self.REACTANT], self.FVAL_DICT[self.PRODUCT],
self.FVAL_DICT[self.MIDWAY]))
presumed_ts_help = """
Location of presumed transition state along the reaction path.
Choices are %s (%.2f), %s (%.2f), and %s (%.2f). Used only
in the '-densearound' and '-connectivity' options."""
fvals = []
for key, value in self.FVAL_DICT.items():
if key != self.REACTANT and key != self.PRODUCT:
fvals.extend([key, value])
self.parserobj.add_argument('-presumed_ts',
choices=self.PRESUMED_TS_CHOICES,
type=str,
default=self.PRESUMED_TS,
help=presumed_ts_help % tuple(fvals))
densearoundhelp = """
Increase the density of sampling points in the presumed transition
state region by adding %s sampling points, with a step size of %.2f,
centered on the transition state location. These points will override
any pre-existing sampling points in these regions, for example those
obtained from the '-sample' option."""
self.parserobj.add_argument('-densearound',
action='store_true',
help=densearoundhelp %
(self.NUMDENSEPOINTS, self.STEPDENSEPOINTS))
interpolationhelp = """
Specify the coordinate system in which to interpolate reaction
path points. Options include %s, %s, and %s. Option %s is the
fastest and works for certain reactions, usually simple reactions
involving a small number of reactive atoms. Option %s is more
computationally demanding but may have some advanages over %s,
for example it will not over-estimate bond lengths. Option %s
is also computationally demanding but works for the greatest
number of reactions, for example it can be used to avoid certain
types of atomic collisions along the reaction path."""
self.parserobj.add_argument(
'-interpolation',
choices=self.INTERPOLATIONCHOICES,
type=str,
default=self.CARTESIAN,
help=interpolationhelp %
(self.CARTESIAN, self.DISTANCE, self.INTERNAL, self.CARTESIAN,
self.DISTANCE, self.CARTESIAN, self.INTERNAL))
bondweighthelp = """
Specify the weight, w >= 0, of the bond terms used in fitting the
reaction path. Increase this value to increase the rigidity of
interpolated bonds or decrease it if one experiences difficulty
with convergence. Typically this value will not need to be
decreased."""
self.parserobj.add_argument('-bondweight',
type=float,
default=self.BONDWEIGHT,
help=bondweighthelp)
angleweighthelp = """
Specify the weight, w >= 0, of the angle terms used in fitting the
reaction path. Increase this value to increase the rigidity of
interpolated angles or decrease it if one experiences difficulty
with convergence. Typically this value will not need to be
decreased."""
self.parserobj.add_argument('-angleweight',
type=float,
default=self.ANGLEWEIGHT,
help=angleweighthelp)
dihedralweighthelp = """
Specify the weight, w >= 0, of the dihedral terms used in fitting the
reaction path. Increase this value to increase the rigidity of
interpolated dihedrals or decrease it if one experiences difficulty
with convergence. Typically this value will not need to be
decreased."""
self.parserobj.add_argument('-dihedralweight',
type=float,
default=self.DIHEDRALWEIGHT,
help=dihedralweighthelp)
cartesianweighthelp = """
Specify the weight, w >= 0, of the Cartesian terms used in fitting the
reaction path. Increase this value if one is confident that the
interpolated Cartesian coordinates well approximate the desired coordinates
of the reaction path point. Note that in such a situation the user might
prefer to use the option -interpolation cartesian. Decrease this value
if the interpolated Cartesian coordinates are a poor guess or if one
experiences difficulty with convergence. Increasing or decreasing this
value can typically remedy troublesome reaction paths more effectively
than by changing the other weighting terms."""
self.parserobj.add_argument('-cartesianweight',
type=float,
default=self.CARTESIANWEIGHT,
help=cartesianweighthelp)
penaltyweighthelp = """
Specify the weight, w >= 0, of the bond penalty terms used in fitting the
reaction path. Such weights prevent bond lengths from becoming too small.
Increase this value if the reaction path features atom-atom collisions or
decrease it if the user wants unconstrained bond lengths. Typically this
value will not need to be decreased."""
self.parserobj.add_argument('-penaltyweight',
type=float,
default=self.PENALTYWEIGHT,
help=penaltyweighthelp)
mixprevioushelp = """
Specify the extent to which the optimized Cartesian coordinates
from the previous reaction path point are mixed with the
interpolated Cartesian coordinates for the current reaction path
point. Typically, the interpolated Cartesian coordinates are used
as a guess solution to the optimizer. Mixing can help convergence
for certain types of reactions where these guesses can involve
overlapping atoms. Mixing additionally helps to create continuous
reaction paths. Accepted values are in [%s, %s]."""
self.parserobj.add_argument('-mixprevious',
type=float,
default=self.MIXPREVIOUS,
help=mixprevioushelp %
(self.MIXPREVIOUSMIN, self.MIXPREVIOUSMAX))
guesshelp = """
Specify the type of guess solution to use from the optimized
Cartesian coordinates of the previous reaction path point.
Note that if the '-mixprevious' option is zero then the
'-guess' option is irrelevant. Accepted values are %s and
%s. Option %s uses the optimized coordinates returned from
the optimizer while option %s will first superpose the
optimized structure on to the reactant."""
self.parserobj.add_argument(
'-guess',
choices=self.GUESSCHOICES,
type=str,
default=self.BEFORESUPERPOSITION,
help=guesshelp %
(self.BEFORESUPERPOSITION, self.AFTERSUPERPOSITION,
self.BEFORESUPERPOSITION, self.AFTERSUPERPOSITION))
connectivityhelp = """
Specify the type of connectivity, i.e. bond assignment, protocol
to use in defining the structures along the reaction path. Options
include %s, %s, or %s. Options %s and %s will use those bonds
defined in the reactant or product as the bonds in the structures
along the reaction path. Option %s specifies that the bonding of
the structures along the reaction path will change from that of
reactants, to something resembling both reactants and products
in the specified transition state region, and then finally to
products."""
self.parserobj.add_argument('-connectivity',
choices=self.CONNECTIVITYCHOICES,
type=str,
default=self.TS,
help=connectivityhelp %
(self.REACTANT, self.TS, self.PRODUCT,
self.REACTANT, self.PRODUCT, self.TS))
inputfileshelp = """
White space delimited list of input Maestro files where each
file may contain multiple reactant/product pairs of
structures ordered like (1) reactants for reaction 1, (2) products
for reaction 1, (3) reactants for reaction 2, (4) products for
reaction 2, etc."""
self.parserobj.add_argument('input_files',
nargs='+',
type=str,
help=inputfileshelp)
self.parserobj.add_argument('-verbose',
action='store_true',
help='Turn on verbose printing.')
self.parserobj.add_argument(
'-version',
'-v',
action='version',
default=argparse.SUPPRESS,
version=_version,
help="Show the script's version number and exit.")
driverhosthelp = """
Job Control option - if used with -HOST specifies the host on
which to run the script itself."""
self.parserobj.add_argument('-DRIVERHOST',
type=str,
default='localhost',
help=driverhosthelp)
hosthelp = """
Job Control option - specifies the host on which to run the
script or if used with -DRIVERHOST specifies the host to use
for any subjobs."""
self.parserobj.add_argument('-HOST',
type=str,
default='localhost',
help=hosthelp)
savehelp = """
Job Control option - copy the archived contents of the job
directory back to the submission directory after the job
finishes."""
self.parserobj.add_argument('-SAVE',
action='store_true',
default=False,
help=savehelp)
tmpdirhelp = """
Job Control option - specify the scratch directory for the
job."""
self.parserobj.add_argument('-TMPDIR',
type=str,
default='/scr',
help=tmpdirhelp)
localhelp = """
Job Control option - write temporary files in the submission
directory instead of the scratch directory."""
self.parserobj.add_argument('-LOCAL',
default=False,
action='store_true',
help=localhelp)
waithelp = """
Job Control option - wait for the job to finish before executing
another command."""
self.parserobj.add_argument('-WAIT',
default=False,
action='store_true',
help=waithelp)
debughelp = """
Job Control option - show the details of operation of the
top-level script."""
self.parserobj.add_argument('-DEBUG',
default=False,
action='store_true',
help=debughelp)
[docs] def parseArgs(self, args):
"""
Parse the command line arguments.
:type args: tuple
:param args: command line arguments
"""
# move input files, which will be the only positional arguments,
# to be first in the list of args
args = list(args)
copyofargs = list(args)
infilesfromargs = []
for arg in copyofargs:
if fileutils.splitext(arg)[1] in self.SUPPORTEDINEXTS:
args.remove(arg)
infilesfromargs.append(arg)
args = infilesfromargs + args
self.options = self.parserobj.parse_args(args)
[docs]class Coord(object):
"""
Manage the properties of an internal coordinate.
"""
BOND = 'bond'
ANGLE = 'angle'
DIHEDRAL = 'dihedral'
[docs] def __init__(self, indicies, names, value):
"""
Create a Coord instance.
:type indicies: list of int
:param indiciees: atomic indicies
:type names: list of str
:param names: atomic names
:type value: list of float
:param value: value(s) of the internal coordinate in units of
Angstrom or degree
"""
self.indicies = indicies
self.names = names
self.value = value
if len(self.indicies) == 2:
description = self.BOND
elif len(self.indicies) == 3:
description = self.ANGLE
elif len(self.indicies) == 4:
description = self.DIHEDRAL
self.description = description
[docs]class InternalCoords(object):
"""
Manage the internal coordinates of a structure.
"""
ZMATNUMBER = 1
[docs] def __init__(self):
"""
Create an InternalCoords instance.
"""
self.bonds = []
self.angles = []
self.dihedrals = []
[docs] def getZmatrix(self, astructure):
"""
Get the Z-matrix for the structure.
:raise ValueError: if there is a problem with the input
:type astructure: schrodinger.structure.Structure
:param astructure: the structure
"""
# from a Jaguar input object get a mmjag handle, convert it to
# the Z-matrix representation, and abstract all of the internal
# coordinate information for each atom
jaginpobj = jaginp.JaguarInput(structure=astructure)
mmjaghandle = jaginpobj.handle
try:
mm.mmjag_zmat_makeint(mmjaghandle, self.ZMATNUMBER - 1)
except MmException as err:
raise ValueError(
str(err) + ' This typically means that 3 or more'
' input atoms are in a line.')
numzmatatoms = mm.mmjag_zmat_atom_count(mmjaghandle,
self.ZMATNUMBER - 1)
for zmatatom in range(numzmatatoms):
atom1 = zmatatom + 1
(zmattype, atom2, bond, atom3, angle, atom4, dihedral) = \
mm.mmjag_zmat_atom_get(mmjaghandle, self.ZMATNUMBER - 1,
atom1)
# form arrays of indicies and names
bondindicies = [atom1, atom2]
angleindicies = [atom1, atom2, atom3]
dihedralindicies = [atom1, atom2, atom3, atom4]
atom1name = mm.mmjag_zmat_atom_get_name(mmjaghandle, atom1)
atom2name = mm.mmjag_zmat_atom_get_name(mmjaghandle, atom2)
atom3name = mm.mmjag_zmat_atom_get_name(mmjaghandle, atom3)
atom4name = mm.mmjag_zmat_atom_get_name(mmjaghandle, atom4)
bondnames = [atom1name, atom2name]
anglenames = [atom1name, atom2name, atom3name]
dihedralnames = [atom1name, atom2name, atom3name, atom4name]
# collect Coord instances for each type of internal coordinate
if atom2:
bondobj = Coord(bondindicies, bondnames, [bond])
self.bonds.append(bondobj)
if atom3:
angleobj = Coord(angleindicies, anglenames, [angle])
self.angles.append(angleobj)
if atom4:
dihedralobj = Coord(dihedralindicies, dihedralnames,
[dihedral])
self.dihedrals.append(dihedralobj)
[docs] def getDmatrix(self, astructure):
"""
Get the distance matrix for the structure.
:type astructure: schrodinger.structure.Structure
:param astructure: the structure
"""
for atom1 in astructure.atom:
for atom2 in astructure.atom:
if atom2.index > atom1.index:
bondindicies = [atom2.index, atom1.index]
# If somehow no Atom object 'name' properties are defined,
# either by not being defined in the input Maestro file
# for some reason, defined in the input Maestro file as
# "" for some atoms, or for whatever reason not defined upon
# structure object creation, then even though Maestro atom
# names are not simply the concatenation of the elemental
# symbol with the atom index we will use that nomenclature
# to reference them
if atom2.name and atom1.name:
bondnames = [atom2.name, atom1.name]
else:
atom2_name = atom2.element + str(atom2.index)
atom1_name = atom1.element + str(atom1.index)
bondnames = [atom2_name, atom1_name]
distance = measure.measure_distance(atom2, atom1)
bondobj = Coord(bondindicies, bondnames, [distance])
self.bonds.append(bondobj)
[docs] def printInternals(self, headermsg, maxindexwidth, logger):
"""
Formatted print of header followed by the internal coordinates.
:type headermsg: str
:param headermsg: header
:type maxindexwidth: int
:param maxindexwidth: number of characters in the largest atom index
:type logger: logging.getLogger
:param logger: output logger
"""
# find length of largest descriptor,
# find max length of entire index string,
# atom element names are no larger than two characters,
# find max length of entire name string,
# max length of formatted internal coord
maxdescwidth = len(Coord.DIHEDRAL)
dihedraldim = 4
maxindstrwidth = dihedraldim * (maxindexwidth + 1) - 1
maxelementwidth = 2
maxnamewidth = maxindexwidth + maxelementwidth
maxnamestrwidth = dihedraldim * (maxnamewidth + 1) - 1
maxvalue = len('-123.456')
# print header and number of internal coordinates or
# just return if there are no coordinates
if self.bonds or self.angles or self.dihedrals:
logger.info(headermsg)
logger.info('')
else:
return
if self.bonds:
logger.info('Number of bonds = %d', len(self.bonds))
if self.angles:
logger.info('Number of angles = %d', len(self.angles))
if self.dihedrals:
logger.info('Number of dihedrals = %d', len(self.dihedrals))
logger.info('')
# print formatted internal coordinates
allintcoords = [self.bonds, self.angles, self.dihedrals]
for intcoordtype in allintcoords:
for intcoord in intcoordtype:
description = intcoord.description.ljust(maxdescwidth)
# cat indicies
indexlabel = ''
for index in intcoord.indicies:
indexlabel += str(index).rjust(maxindexwidth) + '-'
indexlabel = indexlabel[:-1]
indexlabel = indexlabel.rjust(maxindstrwidth)
# cat names
namelabel = ''
for name in intcoord.names:
namelabel += name.rjust(maxnamewidth) + '-'
namelabel = namelabel[:-1]
namelabel = namelabel.rjust(maxnamestrwidth)
msg = '%s %s %s' % (description, indexlabel, namelabel)
# cat values
for value in intcoord.value:
valueform = '%.3f' % value
valueform = valueform.rjust(maxvalue)
msg += ' ' + valueform
logger.info(msg)
logger.info('')
[docs]def max_pair_vdw_distance(astructure):
"""
Find the largest atom-atom VDW distance in a structure.
:type astructure: schrodinger.structure.Structure
:param astructure: the structure
:rtype: int, int, float
:return: atom1, atom2, maxdistance, atom1 and atom2 are the first
and second atom indicies and maxdistance is the largest distance. If
input structure is a single atom then just return that atom index
twice followed by twice its VDW radius, i.e. the atomic diameter.
"""
if astructure.atom_total == 1:
atom = astructure.atom[1]
index = atom.index
diameter = 2 * atom.radius
return index, index, diameter
maxdistance = 0
for atoma in astructure.atom:
for atomb in astructure.atom:
if atomb.index > atoma.index:
distance = measure.measure_distance(atoma, atomb)
distance = distance + atoma.radius + atomb.radius
if distance > maxdistance:
atom1 = atoma.index
atom2 = atomb.index
maxdistance = distance
return atom1, atom2, maxdistance
[docs]def add_temp_hydrogen(astructure, index):
"""
To the given structure add a temporary hydrogen
to the atom with the given index. This function
is more robust than structutils.build.add_hydrogens.
:type astructure: schrodinger.structure.Structure
:param astructure: the structure containing the atom to
which a hydrogen will be added
:type index: int
:param index: the index of the atom to which to add the
hydrogen
:rtype: int
:return: the index of the added temporary hydrogen
"""
# the following class takes a schrodinger.structure._Molecule
# so just make one up and then hijack the class
molobj = rxn_channel.ReactantMolecule(astructure.molecule[1])
molobj.mol_struct = astructure
molobj.reactive_inds_new.append(index)
molobj.addTempHydrogens()
return molobj.tmp_hydrogen_index
[docs]class ReactionCoords(object):
"""
Manage reaction coordinates.
"""
REACTIONBONDTHRESH = 0.05
REACTIONANGLETHRESH = 1.00
REACTIONDIHEDRALTHRESH = 1.00
REVOLUTION = 360
HALFREVOLUTION = 180
REVOLUTIONTHRESH = 10
REACTANTPREFIX = 'pre-'
PRODUCTPREFIX = 'post-'
[docs] def __init__(self):
"""
Create a ReactionCoords instance.
"""
self.rcpypoint = None
self.rpoint = None
self.pcpypoint = None
self.ppoint = None
self.reactioninternals = None
self.tosuperpose = None
self.armsd = None
[docs] def getNormalOrdering(self, reactant, product, logger=None):
"""
Attempt to normalize the atomic ordering between reactants and
products.
:type reactant: schrodinger.structure.Structure
:param reactant: the reactant
:type product: schrodinger.structure.Structure
:param product: the product
:type logger: logging.getLogger
:param logger: output logger
:rtype: schrodinger.structure.Structure, schrodinger.structure.Structure
:return: newreactant, newproduct, if defined specifies the reordered
reactant and product structures
"""
origreactant = reactant.copy()
origproduct = product.copy()
newreactant = None
newproduct = None
rprops = origreactant.property
pprops = origproduct.property
rpobj = rmsd.ConformerRmsd(origreactant,
origproduct,
asl_expr='all',
in_place=True)
try:
armsd = rpobj.calculate()
newreactant = rpobj._working_ref_st.copy()
newreactant.property = rprops
newproduct = rpobj._working_test_st.copy()
newproduct.property = pprops
pordering = []
for sublist in rpobj.atom_index_map:
rorigindex, porigindex = sublist
for pnewatom in newproduct.atom:
poldindex = pnewatom.property[rpobj.orig_index_prop]
if poldindex == porigindex:
pordering.append(pnewatom.index)
break
newpordering = []
for pindex in range(1, len(pordering) + 1):
pmapindex = pordering.index(pindex)
newpordering.append(pmapindex + 1)
newproduct = build.reorder_atoms(newproduct, newpordering)
except RuntimeError:
pass
return newreactant, newproduct
[docs] def makeRxnComplex(self, reactant, product, vdwscale, logger=None):
"""
For certain bimolecular reactions preprocess reactants and products
into reaction complexes.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type vdwscale: float
:param vdwscale: scales the intermolecular distance
:type logger: logging.getLogger
:param logger: output logger
:rtype: schrodinger.structure.Structure, schrodinger.structure.Structure
:return: newreactant, newproduct, if defined specifies the reactant and
product structures in the created reaction complex
"""
def more_general_superposition(refst,
refatoms,
tarst,
taratoms,
logger=None):
"""
Superposes a list of atoms of a target structure on to a
list of atoms of a reference structure. It is more general
than the conventional superposition because it can handle the case
of superposing using two atoms subject to an ambiguous rotation
about an axis which turns out to be irrelevant for this application.
It can also handle the case of a single atom by a simple translation.
:type refst: schrodinger.structure.Structure
:param refst: reference
:type refatoms: list of ints
:param refatoms: reference atom indicies
:type tarst: schrodinger.structure.Structure
:param tarst: target
:type taratoms: list of ints
:param taratoms: target atom indicies
:type logger: logging.getLogger
:param logger: output logger
"""
if len(refatoms) == 1:
refatomobj = refst.atom[refatoms[0]]
taratomobj = tarst.atom[taratoms[0]]
refatomvec = numpy.array(
[refatomobj.x, refatomobj.y, refatomobj.z])
taratomvec = numpy.array(
[taratomobj.x, taratomobj.y, taratomobj.z])
atomatomvec = refatomvec - taratomvec
transform.translate_structure(tarst, atomatomvec[0],
atomatomvec[1], atomatomvec[2])
elif len(refatoms) == 2:
# first rotate the target bond vector so that it is parallel with
# the reference bond vector
refatom1obj = refst.atom[refatoms[0]]
refatom2obj = refst.atom[refatoms[1]]
taratom1obj = tarst.atom[taratoms[0]]
taratom2obj = tarst.atom[taratoms[1]]
refatom1vec = numpy.array(
[refatom1obj.x, refatom1obj.y, refatom1obj.z])
refatom2vec = numpy.array(
[refatom2obj.x, refatom2obj.y, refatom2obj.z])
taratom1vec = numpy.array(
[taratom1obj.x, taratom1obj.y, taratom1obj.z])
taratom2vec = numpy.array(
[taratom2obj.x, taratom2obj.y, taratom2obj.z])
refbond = refatom2vec - refatom1vec
tarbond = taratom2vec - taratom1vec
rotmatrix = transform.get_alignment_matrix(tarbond, refbond)
transform.transform_structure(tarst, rotmatrix)
# second translate the target so that the centroid of its bond vector
# coincides with that of the reference
refsuperposecentroid = transform.get_centroid(refst, refatoms)
tarsuperposecentroid = transform.get_centroid(tarst, taratoms)
intermolecularvec = refsuperposecentroid - tarsuperposecentroid
intermolecularvec = numpy.delete(intermolecularvec, 3)
intermolecularx = numpy.dot(intermolecularvec, transform.X_AXIS)
intermoleculary = numpy.dot(intermolecularvec, transform.Y_AXIS)
intermolecularz = numpy.dot(intermolecularvec, transform.Z_AXIS)
transform.translate_structure(tarst, intermolecularx,
intermoleculary, intermolecularz)
else:
# perform conventional superposition
armsd = rmsd.superimpose(refst,
refatoms,
tarst,
taratoms,
use_symmetry=False,
move_which=rmsd.CT)
def set_intermolecular_distance(refst, tarst, vdwscale):
"""
Set the intermolecular distance for the two input structures by translating
the target structure relative to the reference structure along the
centroid-to-centroid vector.
:type refst: schrodinger.structure.Structure
:param refst: reference
:type tarst: schrodinger.structure.Structure
:param tarst: target
:type vdwscale: float
:param vdwscale: scales the intermolecular distance
"""
# get vdw radius of each structure and sum them to get an intermolecular
# vdw distance, which when multiplied by vdwscale, will be the final
# intermolecular distance of the two structures
refatom1, refatom2, refmaxdist = max_pair_vdw_distance(refst)
refmaxrad = old_div(refmaxdist, 2)
taratom1, taratom2, tarmaxdist = max_pair_vdw_distance(tarst)
tarmaxrad = old_div(tarmaxdist, 2)
vdwdistance = refmaxrad + tarmaxrad
# get centroid-to-centroid vector along which we will translate the
# target structure
refcentroid = transform.get_centroid(refst)
tarcentroid = transform.get_centroid(tarst)
centroidvec = tarcentroid - refcentroid
centroiddistance = numpy.linalg.norm(centroidvec)
centroidunitvec = old_div(centroidvec, centroiddistance)
# determine the final length of the intermolecular vector, project it, and
# translate the target structure
intermolecularvec = centroidunitvec * (vdwdistance -
centroiddistance) * vdwscale
intermolecularvec = numpy.delete(intermolecularvec, 3)
intermolecularx = numpy.dot(intermolecularvec, transform.X_AXIS)
intermoleculary = numpy.dot(intermolecularvec, transform.Y_AXIS)
intermolecularz = numpy.dot(intermolecularvec, transform.Z_AXIS)
transform.translate_structure(tarst, intermolecularx,
intermoleculary, intermolecularz)
# given that there are many substructure manipulations here make sure
# that we do not change the original reactant and product passed in.
# not all reactions that reach this point can be handled and so initialize
# the reactant and product structures as None.
origreactant = reactant.copy()
origproduct = product.copy()
newreactant = None
newproduct = None
# catch any bimolecular reaction where one or more atoms, X, are transferred
# from a single reactant or product, A, to the other reactant or product, B.
# X can make up an entire molecule or be part of a molecule while A and B
# can have any number of atoms in addition to those from X (including the
# possibility of one of them having zero atoms), note that the total number
# of reactant and product atoms must still be greater than two, i.e. only
# polyatomic molecules are allowed. All conformations are free to
# change. Where N, N', etc. are the number of atoms):
#
# (A-X)(N+N') --> (A)(N) + (X)(N') (and reverse) (X must be at least a single
# atom and A must have at least two atoms so that we have a polyatomic
# molecule)
# (A-X)(N+N') + (B)(N'') --> (A)(N) + (B-X)(N'+N'') (and reverse) (X must be at
# least a single atoms and A and B can be one or more atoms)
#
# Figure out which atoms are transferred and for bimolecular reactions
# involving two reactants and two products obtain the mapping of
# molecule numbers.
numrmolecules = len(origreactant.molecule)
numpmolecules = len(origproduct.molecule)
if numrmolecules == 2 and numpmolecules <= 2 or \
numrmolecules <= 2 and numpmolecules == 2:
rreactmol, preactmol = 1, 1
if numrmolecules == 2 and numpmolecules == 2:
rmoloneatoms = origreactant.molecule[rreactmol].getAtomIndices()
pmoltwoatoms = origproduct.molecule[preactmol +
1].getAtomIndices()
if set(pmoltwoatoms).issubset(set(rmoloneatoms)) or \
set(rmoloneatoms).issubset(set(pmoltwoatoms)):
preactmol += 1
samereactant = origreactant.molecule[rreactmol].getAtomIndices()
sameproduct = origproduct.molecule[preactmol].getAtomIndices()
transferred = []
if set(sameproduct).issubset(set(samereactant)) or \
set(samereactant).issubset(set(sameproduct)):
transferred = set(samereactant).symmetric_difference(
set(sameproduct))
transferred = list(transferred)
# in the following note that there are three types of behaviors involving
# single atom entities:
# (1) the case where more than a single atom is transferred but where a
# reactant or product molecule may be an atom, for example SN2, i.e.
# Cl + H3C-F --> Cl-CH3 + F
# (2) the case where a single atoms is transferred but where neither reactant
# or product contains atomic species
# (3) the case where a single atoms is transferred and where the reactant
# and or product may contain an atomic species
# find out which reactant and product molecules are oxidized and reduced,
# i.e. those which have lost and gained the transferred set of atoms
rredmolnum = 0
predmolnum = 0
if transferred:
if logger:
logger.info(
'Building reaction complex for this bimolecular '
'reaction which features the transfer of atoms: %s' %
transferred)
logger.info('')
# handle the case of a single atom being transferred by getting
# a direction via a bond to a temporary hydrogen, the index of
# that hydrogen is the same for reactants and products
temp_index = None
if len(transferred) == 1:
temp_index = add_temp_hydrogen(origreactant, transferred[0])
temp_index = add_temp_hydrogen(origproduct, transferred[0])
transferred.append(temp_index)
for rmolecule in origreactant.molecule:
if set(rmolecule.getAtomIndices()).intersection(
set(transferred)):
roxmolnum = rmolecule.number
roxnumatoms = len(rmolecule.getAtomIndices())
roxmap = OrderedDict(
list(
zip(list(range(1, roxnumatoms + 1)),
rmolecule.getAtomIndices())))
else:
rredmolnum = rmolecule.number
rrednumatoms = len(rmolecule.getAtomIndices())
rredmap = OrderedDict(
list(
zip(list(range(1, rrednumatoms + 1)),
rmolecule.getAtomIndices())))
for pmolecule in origproduct.molecule:
if set(pmolecule.getAtomIndices()).intersection(
set(transferred)):
poxmolnum = pmolecule.number
poxnumatoms = len(pmolecule.getAtomIndices())
poxmap = OrderedDict(
list(
zip(list(range(1, poxnumatoms + 1)),
pmolecule.getAtomIndices())))
else:
predmolnum = pmolecule.number
prednumatoms = len(pmolecule.getAtomIndices())
predmap = OrderedDict(
list(
zip(list(range(1, prednumatoms + 1)),
pmolecule.getAtomIndices())))
# superpose reactant and product on the transferred atoms, handle the
# case of transferring one atom differently
temp_trans = list(transferred)
if temp_index:
temp_trans.reverse()
more_general_superposition(origreactant, transferred,
origproduct, temp_trans, logger)
# if the reactant has a molecule that is being reduced then superpose
# that reactant molecule alone on the product molecule that would be
# considered oxidized if the reaction were reversed. This way the two
# reactant molecules have a relative geometry with respect to each other
# that is favorable for reactivity, hence the calling of this arrangement
# of reactant molecules a reaction complex. Handle the case where the
# reduced reactant molecule contains one or two atoms.
roxmol = origreactant.molecule[roxmolnum].extractStructure(True)
rordering = list(roxmap.values())
if rredmolnum:
rredmol = origreactant.molecule[
rredmolnum].extractStructure(True)
more_general_superposition(origproduct,
list(rredmap.values()), rredmol,
list(rredmap), logger)
# set the intermolecular distance between the two reactants using
# VDW radii, then finally merge the structure objects and mappings
set_intermolecular_distance(roxmol, rredmol, vdwscale)
roxmol.extend(rredmol)
rordering += list(rredmap.values())
# if the product has a molecule that is being reduced then superpose
# that product molecule alone on the reactant molecule that would be
# considered oxidized if the reaction were reversed. This way the two
# product molecules have a relative geometry with respect to each other
# that is favorable for reactivity, hence the calling of this arrangement
# of product molecules a reaction complex. Handle the case where the
# reduced product molecule contains one or two atoms.
poxmol = origproduct.molecule[poxmolnum].extractStructure(True)
pordering = list(poxmap.values())
if predmolnum:
predmol = origproduct.molecule[predmolnum].extractStructure(
True)
more_general_superposition(origreactant,
list(predmap.values()), predmol,
list(predmap), logger)
# set the intermolecular distance between the two products using
# VDW radii, then finally merge the structure objects and mappings
set_intermolecular_distance(poxmol, predmol, vdwscale)
poxmol.extend(predmol)
pordering += list(predmap.values())
# restore original atom numbering to the new merged sets of reactant structures
# and product structures, which are now favorably oriented with respect to each
# other as they would be in a reaction complex.
newrordering = []
for rindex in range(1, len(rordering) + 1):
rmapindex = rordering.index(rindex)
newrordering.append(rmapindex + 1)
newreactant = build.reorder_atoms(roxmol, newrordering)
newpordering = []
for pindex in range(1, len(pordering) + 1):
pmapindex = pordering.index(pindex)
newpordering.append(pmapindex + 1)
newproduct = build.reorder_atoms(poxmol, newpordering)
# if a single atom was being transferred then delete the temporary
# hydrogens now
if temp_index:
newreactant.deleteAtoms([temp_index])
newproduct.deleteAtoms([temp_index])
return newreactant, newproduct
[docs] def collectInternals(self,
reactant,
product,
rinternals,
pinternals,
logger=None):
"""
Find the redundant internal coordinates by merging the
coordinates defined in the reactant and product.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type rinternals: InternalCoords
:param rinternals: reactant internal coordinates
:type pinternals: InternalCoords
:param pinternals: product internal coordinates
:type logger: logging.getLogger
:param logger: output logger
"""
def add_in_redundants(reactant, product, rcoords, pcoords, logger=None):
"""
Locate dissimilar internal coordinates between reactant
and product, define missing internal coordinates, and
update list of Coords.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type rcoords: list of Coords
:param rcoords: list for reactant
:type pcoords: list of Coords
:param pcoords: list for product
:type logger: logging.getLogger
:param logger: output logger
"""
# initialize list of Coords to be added, get insertion index and
# dissimilar reactant and product internals
coordstoadd = []
rpcoords = list(zip(rcoords, pcoords))
for index, rpcoord in enumerate(rpcoords):
rcoord, pcoord = rpcoord
if rcoord.indicies != pcoord.indicies:
# get reactant and product atoms for dissimilar internals and
# measure internals
ratoms = [
reactant.atom[atomindex]
for atomindex in pcoord.indicies
]
patoms = [
product.atom[atomindex] for atomindex in rcoord.indicies
]
if len(ratoms) == 2:
rvalue = measure.measure_distance(*ratoms)
pvalue = measure.measure_distance(*patoms)
if len(ratoms) == 3:
rvalue = measure.measure_bond_angle(*ratoms)
pvalue = measure.measure_bond_angle(*patoms)
if len(ratoms) == 4:
rvalue = measure.measure_dihedral_angle(*ratoms)
pvalue = measure.measure_dihedral_angle(*patoms)
# define the reactant Coord in the product and the product Coord
# in the reactant and append to list of Coords to be added
rnewcoord = Coord(pcoord.indicies, pcoord.names, [rvalue])
pnewcoord = Coord(rcoord.indicies, rcoord.names, [pvalue])
coordstoadd.append([index, rnewcoord, pnewcoord])
# add redundant Coords at the correct locations
for offset, sublist in enumerate(coordstoadd):
index, rcoord, pcoord = sublist
index += offset
rcoords.insert(index + 1, rcoord)
pcoords.insert(index, pcoord)
# do it for bonds, angles, and dihedrals
add_in_redundants(reactant, product, rinternals.bonds, pinternals.bonds)
add_in_redundants(reactant, product, rinternals.angles,
pinternals.angles)
add_in_redundants(reactant, product, rinternals.dihedrals,
pinternals.dihedrals)
[docs] def getReactionInternals(self, rinternals, pinternals, reactioninternals):
"""
Determine the reactive redundant internal coordinates.
:type rinternals: InternalCoords
:param rinternals: reactant internal coordinates
:type pinternals: InternalCoords
:param pinternals: product internal coordinates
:type reactioninternals: InternalCoords
:param reactioninternals: reaction internal coordinates
"""
def collect_reactive(rcoords, pcoords, reactioncoords, threshold):
"""
Collect the reactive coordinates of a given type.
:type rcoords: list of Coords
:param rcoords: Coords for reactant
:type pcoords: list of Coords
:param pcoords: Coords for product
:type reactioncoords: list of Coords
:param reactioncoords: Coords for reactions
:type threshold: float
:param threshold: cutoff for reactive internals
"""
# collect internals that change more than a threshold and store
# the initial and final values in a Coord object
for rcoord, pcoord in zip(rcoords, pcoords):
rvalue = rcoord.value[0]
pvalue = pcoord.value[0]
pphase = 0
if rvalue < 0 and pvalue > 0 or rvalue > 0 and pvalue < 0:
pphase = self.REVOLUTION
if rvalue < 0:
pphase = -1 * self.REVOLUTION
if abs(pvalue + pphase - rvalue) >= threshold:
value = [rvalue, pvalue]
coord = Coord(rcoord.indicies, rcoord.names, value)
reactioncoords.append(coord)
# collect reactive bonds, angles, and dihedrals. Ensure proper handling
# of cases where the dihedral angles differ in sign. Changes in dihedral
# angles are measured using the smaller of the clockwise and counter-clockwise
# rotations.
collect_reactive(rinternals.bonds, pinternals.bonds,
reactioninternals.bonds, self.REACTIONBONDTHRESH)
collect_reactive(rinternals.angles, pinternals.angles,
reactioninternals.angles, self.REACTIONANGLETHRESH)
collect_reactive(rinternals.dihedrals, pinternals.dihedrals,
reactioninternals.dihedrals,
self.REACTIONDIHEDRALTHRESH)
[docs] def runSuperposition(self,
reactant,
product,
reactioninternals,
logger=None):
"""
Superpose the product structure on to the reactant structure
using the non-reactive atoms, i.e. those that do not define
any reactive redundant internal coordinate.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type reactioninternals: InternalCoords
:param reactioninternals: reaction internal coordinates
:type logger: logging.getLogger
:param logger: output logger
:rtype: list of ints, float
:return: tosuperpose, armsd, atom indicies used to superpose and
the final RMSD
"""
# get set of atom indicies for atoms that belong to
# internal coordinates that define the reaction coordinate
reactionatoms = set()
allreactioninternals = [
reactioninternals.bonds, reactioninternals.angles,
reactioninternals.dihedrals
]
for coords in allreactioninternals:
for internal in coords:
for index in internal.indicies:
reactionatoms.add(index)
# get list of atoms to superpose
tosuperpose = set()
for atom in reactant.atom:
if atom.index not in reactionatoms:
tosuperpose.add(atom.index)
# if there are not enough atoms left to do the superposition
# then add some of the atoms that define the least reactive
# internal coordinates
rnatoms = reactant.atom_total
minnatoms = min(4, rnatoms)
if len(tosuperpose) < minnatoms:
allchanges = []
for coords in allreactioninternals:
for internal in coords:
rvalue = internal.value[0]
pvalue = internal.value[1]
pphase = 0
if rvalue < 0 and pvalue > 0 or rvalue > 0 and \
pvalue < 0:
pphase = self.REVOLUTION
if rvalue < 0:
pphase = -1 * self.REVOLUTION
change = abs(pvalue + pphase - rvalue)
allchanges.append([internal.indicies, change])
allchanges.sort(key=itemgetter(1))
for change in allchanges:
for index in change[0]:
if index not in tosuperpose:
tosuperpose.add(index)
if len(tosuperpose) >= minnatoms:
break
# superpose the product on to the reactant
tosuperpose = list(tosuperpose)
tosuperpose.sort()
armsd = rmsd.superimpose(reactant,
tosuperpose,
product,
tosuperpose,
use_symmetry=False,
move_which=rmsd.CT)
return tosuperpose, armsd
[docs] def getCartesianCoords(self, astructure):
"""
Get the Cartesian coordinates of a structure as a
3N dimensional list of floats.
:type astructure: schrodinger.structure.Structure
:param astructure: structure
:rtype: list of float
:return: cartesians, the 3N Cartesian coordinates
ordered like [x1, y1, z1, x2, ..., zN]
"""
cartesians = []
for atom in astructure.atom:
cartesians.append(atom.x)
cartesians.append(atom.y)
cartesians.append(atom.z)
return cartesians
[docs] def prepare(self,
reactant,
product,
interpolation,
norxncomplex,
vdwscale,
samplepoints,
logger=None):
"""
Prepare reaction coordinates.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type interpolation: str
:param interpolation: coordinate system used for interpolating
:type norxncomplex: boolean
:param norxncomplex: disable preprocessing into a reaction complex
:type vdwscale: float
:param vdwscale: scales the intermolecular distance
:type samplepoints: list of float
:param samplepoints: reaction path sample points
:type logger: logging.getLogger
:param logger: output logger
"""
# get the length of the largest atom index
maxindexwidth = len(str(reactant.atom_total))
# get level name of logger
if logger:
loglevel = logging.getLevelName(logger.getEffectiveLevel())
# for bimolecular reactions preprocess reactant and product into
# a reaction complex
reactantcpy = None
productcpy = None
if not norxncomplex:
newreactant, newproduct = \
self.makeRxnComplex(reactant, product, vdwscale, logger)
if newreactant:
reactantcpy = reactant.copy()
productcpy = product.copy()
reactant = newreactant.copy()
product = newproduct.copy()
# get Z-matrix coordinates for reactant and product
rinternals = InternalCoords()
pinternals = InternalCoords()
if interpolation == ParserWrapper.DISTANCE:
rinternals.getDmatrix(reactant)
pinternals.getDmatrix(product)
else:
rinternals.getZmatrix(reactant)
pinternals.getZmatrix(product)
if logger:
if loglevel == 'DEBUG':
rinternals.printInternals('These are the ' \
'Z-matrix coordinates for the reactant:',
maxindexwidth, logger)
pinternals.printInternals('These are the ' \
'Z-matrix coordinates for the product:',
maxindexwidth, logger)
# merge coordinates defined in Z-matricies
self.collectInternals(reactant, product, rinternals, pinternals, logger)
if logger:
if loglevel == 'DEBUG':
rinternals.printInternals('These are the ' \
'redundant internal coordinates for the reactant:',
maxindexwidth, logger)
pinternals.printInternals('These are the ' \
'redundant internal coordinates for the product:',
maxindexwidth, logger)
# determine reaction coordinates
reactioninternals = InternalCoords()
self.getReactionInternals(rinternals, pinternals, reactioninternals)
if logger:
reactioninternals.printInternals('These are the reaction ' \
'coordinates:', maxindexwidth, logger)
self.reactioninternals = reactioninternals
# superpose the product on to the reactant
self.tosuperpose, self.armsd = self.runSuperposition(
reactant, product, reactioninternals, logger)
if logger:
logger.debug('The following atoms of the reactant/product pair have ' \
'been used to perform this superposition: %s' % self.tosuperpose)
logger.debug('')
# get the Cartesian coordinates of the reactant and product
rcartesians = self.getCartesianCoords(reactant)
pcartesians = self.getCartesianCoords(product)
# collect Point objects for reactant
rfval = ParserWrapper.FVAL_DICT[ParserWrapper.REACTANT]
rindex = 1
rtitle = reactant.property[CheckInput.TITLEKEY]
rentry = reactant.property[CheckInput.ENTRYNAMEKEY]
if reactantcpy:
prertitle = self.REACTANTPREFIX + rtitle
prerentry = self.REACTANTPREFIX + rentry
reactantcpy.property[CheckInput.TITLEKEY] = prertitle
reactantcpy.property[CheckInput.ENTRYNAMEKEY] = prerentry
reactantcpy.property[ReactionPath.RXNINDEX] = rindex
reactantcpy.property[ReactionPath.RXNCOORD] = rfval
self.rcpypoint = Point(rindex, rfval, prertitle, reactantcpy, None,
None)
rindex += 1
reactant.property[ReactionPath.RXNINDEX] = rindex
reactant.property[ReactionPath.RXNCOORD] = rfval
self.rpoint = Point(rindex, rfval, rtitle, reactant, rinternals,
rcartesians)
# collect Point objects for product
pfval = ParserWrapper.FVAL_DICT[ParserWrapper.PRODUCT]
pindex = len(samplepoints) + rindex + 1
ptitle = product.property[CheckInput.TITLEKEY]
pentry = product.property[CheckInput.ENTRYNAMEKEY]
product.property[ReactionPath.RXNINDEX] = pindex
product.property[ReactionPath.RXNCOORD] = pfval
self.ppoint = Point(pindex, pfval, ptitle, product, pinternals,
pcartesians)
if reactantcpy:
pindex += 1
preptitle = self.PRODUCTPREFIX + ptitle
prepentry = self.PRODUCTPREFIX + pentry
productcpy.property[CheckInput.TITLEKEY] = preptitle
productcpy.property[CheckInput.ENTRYNAMEKEY] = prepentry
productcpy.property[ReactionPath.RXNINDEX] = pindex
productcpy.property[ReactionPath.RXNCOORD] = pfval
self.pcpypoint = Point(pindex, pfval, preptitle, productcpy, None,
None)
[docs]class Point(object):
"""
Collect properties of reaction path points.
"""
[docs] def __init__(self, index, fval, name, astructure, internals, cartesians):
"""
Create a Point instance.
:type index: int
:param index: path point index
:type fval: float
:param fval: path point value
:type name: str
:param name: path point name
:type astructure: schrodinger.structure.Structure
:param astructure: structure
:type internals: InternalCoords
:param internals: path point internal coordinates
:type cartesians: list of float
:param cartesians: the 3N Cartesian coordinates
ordered like [x1, y1, z1, x2, ..., zN]
"""
self.index = index
self.fval = fval
self.name = name
self.astructure = astructure
self.internals = internals
self.cartesians = cartesians
[docs]class ReactionPath(object):
"""
Generate reaction path.
"""
FVALINCREMENT = 0.001
NORMTHRESH = 0.000000000001
BONDPENALTYTHRESH = 1000000000.0
DIFFLOWTHRESH = 10.0
DIFFHIGHVAL = 1000000000.0
RXNINDEX = 'i_matsci_RXN_Index'
RXNCOORD = 'r_matsci_RXN_Coord'
REACTIVEATOM = 'b_matsci_Reactive_Atom'
[docs] def __init__(self):
"""
Create a ReactionPath instance.
"""
self.samplepoints = []
self.reactioninternals = None
self.tosuperpose = None
self.armsd = None
self.reactiveatoms = []
self.points = []
self.bondweight = None
self.angleweight = None
self.dihedralweight = None
self.cartesianweight = None
self.penaltyweight = None
[docs] def getSamplePoints(self, sample, densearound, presumed_ts):
"""
Determine the final set of sampling points.
:type sample: list of float
:param sample: points to be interpolated
:type densearound: bool
:param densearound: include additional sampling points at
specific regions
:type presumed_ts: str
:param presumed_ts: location of presumed ts
:rtype: list of floats
:return: samplepoints, the list of points to be sampled.
"""
# if the user supplied more than a single value then just use
# those values provided.
if len(sample) > 1:
self.samplepoints.extend(sample)
# get reactant and product fvals
rfval = ParserWrapper.FVAL_DICT[ParserWrapper.REACTANT]
pfval = ParserWrapper.FVAL_DICT[ParserWrapper.PRODUCT]
# if the user supplied a single value less than 1 then just use
# that value otherwise take that single greater-than-1 value
# and find that number of uniformly distributed sampling points
if len(sample) == 1:
if sample[0] < 1:
fval = sample[0]
self.samplepoints.append(fval)
else:
npoints = sample[0]
numerator = pfval - rfval
stepsize = old_div(numerator, (npoints + 1))
for point in range(int(npoints)):
fval = stepsize * (point + 1)
self.samplepoints.append(fval)
densepoints = []
if densearound:
center = ParserWrapper.FVAL_DICT[presumed_ts]
densepoints.append(center)
# half of the dense set of points go to the left of
# the center point and the other half to the right
numdensepoints = ParserWrapper.NUMDENSEPOINTS
stepsize = ParserWrapper.STEPDENSEPOINTS
for point in range(old_div(numdensepoints, 2)):
fval = stepsize * (point + 1)
densepoints.append(center - fval)
densepoints.append(center + fval)
densepoints.sort()
# override with the dense set any preexisting sample
# points in the range that is covered by the dense set
mindensefval = min(densepoints)
maxdensefval = max(densepoints)
self.samplepoints = [fval for fval in self.samplepoints \
if fval < mindensefval or fval > maxdensefval]
self.samplepoints.extend(densepoints)
self.samplepoints.sort()
return self.samplepoints
[docs] def getReactiveAtoms(self, reactant, product):
"""
Determine the reactive atoms, i.e. those which have changed
Cartesian positions in the superposed reactant/product pair.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:rtype: list of int
:return: reactiveatoms, reactive atoms
"""
for ratom, patom in zip(reactant.atom, product.atom):
distance = measure.measure_distance(ratom, patom)
if distance >= ReactionCoords.REACTIONBONDTHRESH:
self.reactiveatoms.append(ratom.index)
return self.reactiveatoms
[docs] def interpolateReactionCoords(self,
pointindex,
fval,
connectivity,
rpoint,
ppoint,
reactiveatoms,
logger=None):
"""
Interpolate reaction coordinates between the reactant and
product for this sample point.
:type pointindex: int
:param pointindex: sample point index
:type fval: float
:param fval: interpolated reaction path point
:type connectivity: str
:param connectivity: specifies the type of connectivity
:type rpoint: Point
:param rpoint: reactant information
:type ppoint: Point
:param ppoint: product information
:type reactiveatoms: list of ints
:param reactiveatoms: reactive atoms
:type logger: logging.getLogger
:param logger: output logger
:rtype: Point
:return: ipoint, interpolated point information
"""
def interpolate_internals(rcoords,
pcoords,
icoords,
reactiveatoms,
fval,
logger=None):
"""
Collect interpolated internals.
:type rcoords: list of Coords
:param rcoords: Coords for reactant
:type pcoords: list of Coords
:param pcoords: Coords for product
:type icoords: list of Coords
:param icoords: Coords for interpolated point
:type reactiveatoms: list of ints
:param reactiveatoms: indicies of reactive atoms
:type fval: float
:param fval: interpolated reaction path point
:type logger: logging.getLogger
:param logger: output logger
"""
for rcoord, pcoord in zip(rcoords, pcoords):
if set(rcoord.indicies).intersection(set(reactiveatoms)):
rvalue = rcoord.value[0]
pvalue = pcoord.value[0]
pphase = 0
if rvalue < 0 and pvalue > 0 or rvalue > 0 and pvalue < 0:
pphase = ReactionCoords.REVOLUTION
if rvalue < 0:
pphase = -1 * ReactionCoords.REVOLUTION
coord = (1 - fval) * rvalue + fval * (pvalue + pphase)
if pphase:
if logger:
logger.debug('Treating the phase of dihedral %s' %
rcoord.indicies)
logger.debug('reactant dihedral = %.3f' %
rvalue)
logger.debug('product dihedral w/o phase = %.3f' %
pvalue)
logger.debug('product dihedral w phase = %.3f' %
(pvalue + pphase))
logger.debug('interpolated dihedral = %.3f' %
coord)
if coord < -1 * ReactionCoords.HALFREVOLUTION:
coord += ReactionCoords.REVOLUTION
if coord > ReactionCoords.HALFREVOLUTION:
coord += -1 * ReactionCoords.REVOLUTION
if logger:
logger.debug('final interpolated dihedral = %.3f' %
coord)
logger.debug('')
icoord = Coord(rcoord.indicies, rcoord.names, [coord])
icoords.append(icoord)
# collect interpolated bonds, angles, and dihedrals. Ensure proper handling
# of dihedral angles that differ in sign. Changes in dihedral angles are measured
# using the smaller of the clockwise and counter-clockwise rotations.
internals = InternalCoords()
interpolate_internals(rpoint.internals.bonds, ppoint.internals.bonds,
internals.bonds, reactiveatoms, fval, logger)
interpolate_internals(rpoint.internals.angles, ppoint.internals.angles,
internals.angles, reactiveatoms, fval, logger)
interpolate_internals(rpoint.internals.dihedrals,
ppoint.internals.dihedrals, internals.dihedrals,
reactiveatoms, fval, logger)
# interpolate all of the Cartesian coordinates. The
# fixed Cartesians will keep their fixed values and
# will be used along with the optimizable Cartesians
# of the reactive atoms to compute the internal
# coordinates which are to reproduce in a non-linear
# least squares sense the interpolated internals.
# The increment here prevents atoms from having the
# same interpolated Cartesian coordinates, for example
# as I saw with HCN -> CNH. This is a numerical
# precision issue, meaning that it is only coordinates
# which are identical according to the machine
# are a problem. Coordinates that are just close
# to one another are alright. The shift applied
# is fairly arbitrary because the interpolated
# Cartesians are simply an initial guess
# to the optimizer.
tmpcartesians = []
rpcartesians = list(zip(rpoint.cartesians, ppoint.cartesians))
for index in range(0, len(rpcartesians), 3):
rx, px = rpcartesians[index:index + 3][0]
ry, py = rpcartesians[index:index + 3][1]
rz, pz = rpcartesians[index:index + 3][2]
ix = (1 - fval) * rx + fval * px
iy = (1 - fval) * ry + fval * py
iz = (1 - fval) * rz + fval * pz
ixyz = [ix, iy, iz]
if ixyz not in tmpcartesians:
tmpcartesians.append(ixyz)
else:
shiftedfval = fval + self.FVALINCREMENT
ix = (1 - shiftedfval) * rx + shiftedfval * px
iy = (1 - shiftedfval) * ry + shiftedfval * py
iz = (1 - shiftedfval) * rz + shiftedfval * pz
ixyz = [ix, iy, iz]
tmpcartesians.append(ixyz)
cartesians = []
for xyz in tmpcartesians:
cartesians.extend(xyz)
# for convenience create a Point object for this interpolated
# point now even though the structure, internals, and
# cartesians will need to be updated later
if connectivity != ParserWrapper.PRODUCT:
refstructure = rpoint.astructure
else:
refstructure = ppoint.astructure
astructure = refstructure.copy()
pointname = str(pointindex)
ipoint = Point(pointindex, fval, pointname, astructure, internals,
cartesians)
return ipoint
[docs] def getInitialGuess(self, cartesians, reactiveatoms):
"""
Obtain the initial guess Cartesians for the non-linear least
squares solver. The guess is the interpolated Cartesian
coordinates for the reactive atoms.
:type cartesians: list of floats
:param cartesians: all interpolated Cartesian coordinates
:type reactiveatoms: list of ints
:param reactiveatoms: atom indicies of reactive atoms
:rtype: list of floats
:return: guessparams, interpolated Cartesian coordinates
of the reactive atoms
"""
guessparams = []
for atomindex in reactiveatoms:
xyzstart = 3 * (atomindex - 1)
guessparams.extend(cartesians[xyzstart:xyzstart + 3])
return guessparams
[docs] def doNonLinearFit(self,
ipoint,
guessparams,
reactiveatoms,
mixprevious,
logger=None):
r"""
Using the interpolated redundant internal coordinates and
interpolated Cartesian coordinates obtain the final set of
Cartesian coordinates for this reaction path point by
minimizing a sum-of-squares error function using non-linear
least sqaures, i.e.::
error = \sum_{bonds} bondweight*(r(a,b) - r^{i}(a,b))**2
+ \sum_{angles} angleweight*(theta(a,b,c) - theta^{i}(a,b,c))**2
+ \sum_{dihedrals} dihedralweight*(tau(a,b,c,d) - tau^{i}(a,b,c,d))**2
+ \sum_{atoms} cartweight*[(x(a) - x^{i}(a))**2 + (y(a) - y^{i}(a))**2
+ (z(a) - z^{i}(a))**2]
where those variables marked with "^{i}" are the interpolated
quantities and where::
r(a,b) = r(x(a), y(a), z(a), x(b), y(b), z(b)) = norm(vec(a,b))
theta(a,b,c) = arccos[(vec(a,b) \dot vec(c,b))/(norm(vec(a,b))*norm(vec(c,b)))]
tau(a,b,c,d) = arccos[((vec(c,b) \cross vec(a,b)) \dot (vec(d,c) \cross vec(b,c)))
/ (norm((vec(c,b) \cross vec(a,b)))*norm((vec(d,c) \cross vec(b,c))))]
The 3N Cartesian coordinates, x(a), y(a), z(a), x(b), ..., z(N), are
choosen so as to minimize the error.
:type ipoint: Point
:param ipoint: interpolated point information
:type guessparams: list of floats
:param guessparams: initial parameters, i.e. Cartesian coordinates of the
reactive atoms
:type reactiveatoms: list of ints
:param reactiveatoms: atomic indicies of reactive atoms
:type mixprevious: float
:param mixprevious: specifies to what extent the optimized
Cartesian coordinates from the previous reaction path point
are mixed with the coordinates determined by interpolation
at the current point.
:type logger: logging.getLogger
:param logger: output logger
:rtype: list
:return: optcartesians, non-linear-optimized Cartesian
coordinates for this sample point.
"""
def get_vectors(params,
indicies,
cartesians,
reactiveatoms,
logger=None):
"""
Return a list of vectors pointing to each of the atoms in
indicies.
:type params: numpy.array
:param params: optimizable Cartesians
:type indicies: list of ints
:param indicies: indicies defining the internal coordinate
:type cartesians: list
:param cartesians: Cartesians
:type reactiveatoms: list of ints
:param reactiveatoms: reactive atom indicies
:type logger: logging.getLogger
:param logger: output logger
:rtype: list of numpy.array
:return: vecs, vectors pointing to atoms in indicies
"""
# for non-reactive and reactive atoms use the original and
# optimized Cartesians, respectively
vecs = []
for atomindex in indicies:
if atomindex in reactiveatoms:
tmpindex = reactiveatoms.index(atomindex)
tmpindex = 3 * tmpindex
values = list(params)
else:
tmpindex = 3 * (atomindex - 1)
values = list(cartesians)
vecs.append(numpy.array(values[tmpindex:tmpindex + 3]))
return vecs
def get_bond_distance(params,
indicies,
cartesians,
reactiveatoms,
logger=None):
"""
Return the bond distance.
:type params: numpy.array
:param params: parameters being optimized
:type indicies: list of ints
:param indicies: bond indicies
:type cartesians: list
:param cartesians: all Cartesian coordinates
:type reactiveatoms: list of ints
:param reactiveatoms: atom indicies of reactive atoms.
:type logger: logging.getLogger
:param logger: output logger
:rtype: float
:return: distance, bond distance in Ang.
"""
# form vectors and find distance
vecs = get_vectors(params, indicies, cartesians, reactiveatoms,
logger)
veca = vecs[0]
vecb = vecs[1]
bondvec = vecb - veca
distance = numpy.linalg.norm(bondvec)
return distance
def get_bond_weight(interpolatedval):
"""
Return the weight of the bonding term.
:type interpolatedval: float
:param interpolatedval: interpolated bond distance
:rtype: float
:return: weight, the weight of this bond term.
"""
return self.bondweight
def get_bond_angle(params,
indicies,
cartesians,
reactiveatoms,
logger=None):
"""
Return the bond angle.
:type params: numpy.array
:param params: parameters being optimized
:type indicies: list of ints
:param indicies: angle indicies
:type cartesians: list
:param cartesians: all Cartesian coordinates
:type reactiveatoms: list of ints
:param reactiveatoms: atom indicies of reactive atoms
:type logger: logging.getLogger
:param logger: output logger
:rtype: float
:return: angle, bond angle in degree
"""
# form vectors and find angle, note that math.acos
# is multivalued, the domain is [-1.0, 1.0] while
# by convention the range is [0.0, 180.0] (degrees).
# Here the argument to math.acos may suffer from
# numerical precision errors and so must be bounded
# in [-1.0, 1.0] in order to avoid a domain error.
# Bond angles are always defined in [0.0, 180.0] so
# we do not need to check the sign of the angle
# returned.
vecs = get_vectors(params, indicies, cartesians, reactiveatoms,
logger)
veca = vecs[0]
vecb = vecs[1]
vecc = vecs[2]
bondvecab = veca - vecb
bondveccb = vecc - vecb
normab = numpy.linalg.norm(bondvecab)
normcb = numpy.linalg.norm(bondveccb)
if normab < self.NORMTHRESH or normcb < self.NORMTHRESH \
and logger:
logger.warning('Vectors defining internal coordinate %s ' \
'are becoming vanishingly small.' % indicies)
acosarg = old_div(numpy.dot(bondvecab, bondveccb),
(normab * normcb))
acosarg = max(min(acosarg, 1), -1)
angle = math.acos(acosarg)
angle = math.degrees(angle)
return angle
def get_angle_weight(interpolatedval):
"""
Return the weight of the angle term.
:type interpolatedval: float
:param interpolatedval: interpolated bond angle
:rtype: float
:return: weight, the weight of this angle term.
"""
return self.angleweight
def get_dihedral_angle(params,
indicies,
cartesians,
reactiveatoms,
logger=None):
"""
Return the dihedral angle.
:type params: numpy.array
:param params: parameters being optimized
:type indicies: list of ints
:param indicies: dihedral indicies
:type cartesians: list
:param cartesians: all Cartesian coordinates
:type reactiveatoms: list of ints
:param reactiveatoms: atom indicies of reactive atoms.
:type logger: logging.getLogger
:param logger: output logger
:rtype: float
:return: dihedral, the dihedral angle in units of degree.
"""
# form vectors and find dihedral, note that math.acos
# is multivalued, the domain is [-1.0, 1.0] while
# by convention the range is [0.0, 180.0] (degrees).
# Here the argument to math.acos may suffer from
# numerical precision errors and so must be bounded
# in [-1.0, 1.0] in order to avoid a domain error.
# Dihedral angles are defined in [-180.0, 180.0] and
# so we need to determine the proper sign of the
# dihedral returned.
vecs = get_vectors(params, indicies, cartesians, reactiveatoms,
logger)
veca = vecs[0]
vecb = vecs[1]
vecc = vecs[2]
vecd = vecs[3]
bondvecab = veca - vecb
bondveccb = vecc - vecb
bondvecbc = vecb - vecc
bondvecdc = vecd - vecc
planevecabc = numpy.cross(bondveccb, bondvecab)
planevecbcd = numpy.cross(bondvecdc, bondvecbc)
normabc = numpy.linalg.norm(planevecabc)
normbcd = numpy.linalg.norm(planevecbcd)
if normabc < self.NORMTHRESH or normbcd < self.NORMTHRESH \
and logger:
logger.warning('Vectors defining internal coordinate %s ' \
'are becoming vanishingly small.' % indicies)
acosarg = old_div(numpy.dot(planevecabc, planevecbcd),
(normabc * normbcd))
acosarg = max(min(acosarg, 1), -1)
dihedral = math.acos(acosarg)
dihedral = math.degrees(dihedral)
sign = numpy.dot(bondvecab, planevecbcd)
if sign > 0:
dihedral = -1 * dihedral
return dihedral
def get_dihedral_weight(interpolatedval):
"""
Return the weight of the dihedral term.
:type interpolatedval: float
:param interpolatedval: interpolated dihedral angle
:rtype: float
:return: weight, the weight of this dihedral term.
"""
return self.dihedralweight
def get_cartesian_weight():
"""
Return the weight of the Cartesian term.
:rtype: float
:return: weight, the weight of this Cartesian term.
"""
return self.cartesianweight
def residuals(params,
bondinds,
bondvals,
angleinds,
anglevals,
dihedralinds,
dihedralvals,
cartesians,
reactiveatoms,
initialparams,
logger=None):
r"""
The main argument to scipy.optimize.leastsq. Determine
the residuals which are to be minimized by leastsq where
the residuals are defined as the argument to the square,
i.e.
error = \sum_{residuals} (residual)**2.
:type params: numpy.array
:param params: parameters being optimized
:type bondinds: list of lists
:param bondinds: all bond indicies
:type bondvals: list
:param bondvals: interpolated bond distances
:type angleinds: list of lists
:param angleinds: all angle indicies
:type anglevals: list
:param anglevals: interpolated bond angles
:type dihedralinds: list of lists
:param dihedralinds: all dihedral indicies
:type dihedralvals: list
:param dihedralvals: interpolated dihedral angles
:type cartesians: list
:param cartesians: all Cartesian coordinates
:type reactiveatoms: list
:param reactiveatoms: atom indicies for reactive atoms.
:type initialparams: list of floats
:param initialparams: list of interpolated Cartesian coordinates
for the reactive atoms.
:type logger: logging.getLogger
:param logger: output logger
:rtype: numpy.array
:return: residual, residuals ordered as bonds, angles, dihedrals,
and cartesians.
"""
residual = []
# for each bonding reaction coordinate form a weighted
# residual with the difference between the interpolated
# distance and the distance calculated from the Cartesian
# coordinates of the two atoms. At least one of the atoms
# will be a reactive atom. For bonds include the bond penalty
# terms which should guide the optimizer away from bonds
# of length zero.
for indicies, interpolatedval in zip(bondinds, bondvals):
distance = get_bond_distance(params, indicies, cartesians,
reactiveatoms, logger)
weight = get_bond_weight(interpolatedval)
diff = interpolatedval - distance
bondresidual = weight * diff / interpolatedval
if distance <= self.BONDPENALTYTHRESH:
bondpenalty = old_div(self.penaltyweight, distance)
bondresidual += bondpenalty
residual.append(bondresidual)
# for each angle reaction coordinate form a weighted
# residual with the difference between the interpolated
# angle and the angle calculated from the Cartesian
# coordinates of the three atoms. At least one of the atoms
# will be a reactive atom.
for indicies, interpolatedval in zip(angleinds, anglevals):
angle = get_bond_angle(params, indicies, cartesians,
reactiveatoms, logger)
weight = get_angle_weight(interpolatedval)
diff = interpolatedval - angle
try:
angleresidual = weight * diff / interpolatedval
except ZeroDivisionError:
if logger:
logger.debug('The interpolated angle %s ' \
'is zero.' % indicies)
logger.debug('')
if abs(diff) <= self.DIFFLOWTHRESH:
angleresidual = 0.0
else:
angleresidual = self.DIFFHIGHVAL
residual.append(angleresidual)
# for each dihedral reaction coordinate form a weighted
# residual with the difference between the interpolated
# dihedral and the dihedral calculated from the Cartesian
# coordinates of the four atoms. At least one of the atoms
# will be a reactive atom.
for indicies, interpolatedval in zip(dihedralinds, dihedralvals):
dihedral = get_dihedral_angle(params, indicies, cartesians,
reactiveatoms, logger)
weight = get_dihedral_weight(interpolatedval)
diff = interpolatedval - dihedral
if (ReactionCoords.REVOLUTION - abs(diff)) < \
ReactionCoords.REVOLUTIONTHRESH:
if logger:
logger.debug('Change in dihedral %s is close to a full '
'revolution: %s %s.' %
(indicies, interpolatedval, dihedral))
logger.debug('')
diff = abs(interpolatedval) - abs(dihedral)
try:
dihedralresidual = weight * diff / interpolatedval
except ZeroDivisionError:
if logger:
logger.debug('The interpolated dihedral angle %s ' \
'is zero.' % indicies)
logger.debug('')
if abs(diff) <= self.DIFFLOWTHRESH:
dihedralresidual = 0.0
else:
dihedralresidual = self.DIFFHIGHVAL
residual.append(dihedralresidual)
# for each Cartesian coordinate of each reactive atom
# add a term to the residual that is the weighted difference
# between the coordinate being optimized and its original
# guess value. This is supposed to ensure that the
# optimizer does not resort to a solution that is that
# different from the guess solution.
weight = get_cartesian_weight()
for param, initialparam in zip(params, initialparams):
cartresidual = weight * (initialparam - param)
residual.append(cartresidual)
residual = numpy.array(residual)
return residual
# isolate out indicies and values of internal coordinates so
# that we are not passing a large object in the args to
# leastsq
bonds = ipoint.internals.bonds
angles = ipoint.internals.angles
dihedrals = ipoint.internals.dihedrals
bondinds, bondvals = [], []
for bond in bonds:
bondinds.append(bond.indicies)
bondvals.append(bond.value[0])
angleinds, anglevals = [], []
for angle in angles:
angleinds.append(angle.indicies)
anglevals.append(angle.value[0])
dihedralinds, dihedralvals = [], []
for dihedral in dihedrals:
dihedralinds.append(dihedral.indicies)
dihedralvals.append(dihedral.value[0])
# the initial guess for the Cartesian coordinates of
# the reactive atoms will be a mixture of the coordinates
# interpolated from the reactant and product for this sample
# point and the optimized coordinates from the last reaction
# path point. Taking a mixture seems to help avoid discontinuous
# reaction paths. We will need two copies (1) initialparams
# which will not be updated and (2) guessparams which will
# be updated.
interpolatedguess = self.getInitialGuess(ipoint.cartesians,
reactiveatoms)
previousguess = list(guessparams)
avgguess = []
for iguess, pguess in zip(interpolatedguess, previousguess):
guess = (1 - mixprevious) * iguess + mixprevious * pguess
avgguess.append(guess)
initialparams = list(avgguess)
guessparams = numpy.array(avgguess)
# minimize the residuals vector by varying the Cartesian
# coordinates of the reactive atoms. It is necessary that
# residuals, guessparams, and optcartesians are numpy.array
# objects.
optcartesians, cov, infodict, mesg, ier = leastsq(
residuals,
guessparams,
args=(bondinds, bondvals, angleinds, anglevals, dihedralinds,
dihedralvals, ipoint.cartesians, reactiveatoms, initialparams,
logger),
full_output=True)
optcartesians = optcartesians.tolist()
# print some leastsq exit info
if logger:
logger.debug('Optimizer: number of function calls = %s' %
infodict['nfev'])
logger.debug('Optimizer: exit status = %s' % ier)
if ier not in [1, 2, 3, 4]:
logger.critical('Optimizer: failure message is %s' % mesg)
logger.debug('')
return optcartesians
[docs] def getInterpolatedStructure(self, ipoint, optcartesians, reactiveatoms,
tosuperpose, fval, connectivity, presumed_ts,
rpoint, ppoint):
"""
Build the schrodinger.structure.Structure object from
the optimized Cartesian coordinates for the interpolated
structure at this sample point and update the internal
and Cartesian coordinates in the Point object.
:type ipoint: Point
:param ipoint: interpolated point information
:type optcartesians: list of floats
:param optcartesians: optimized Cartesian coordinates for
the reactive atoms
:type reactiveatoms: list of ints
:param reactiveatoms: atom indicies of reactive atoms.
:type tosuperpose: list of ints
:param tosuperpose: contains the atom indicies of
the atoms used in the superposition.
:type fval: float
:param fval: The interpolated reaction path point.
:type connectivity: str
:param connectivity: specifies the type of connectivity
:type presumed_ts: str
:param presumed_ts: specifies the location of the presumed ts
:type rpoint: Point
:param rpoint: reactant information
:type ppoint: Point
:param ppoint: product information
:rtype: list of floats, InternalCoords object
:return: optcartesianssuperposed, reactiveinternals, the optimized Cartesian
coordinates for the reactive atoms after superposition on to the reactant
structure and an object containing the interpolated and calculated reactive
internal coordinates.
"""
# update the structure object, i.e. move reactive atoms,
# define connectivity, superpose it on to the reactant,
# define some properties, etc.
if connectivity == ParserWrapper.REACTANT or \
connectivity == ParserWrapper.PRODUCT:
for tmpindex, atomindex in enumerate(reactiveatoms):
iatom = ipoint.astructure.atom[atomindex]
index = 3 * (tmpindex)
iatom.x = optcartesians[index]
iatom.y = optcartesians[index + 1]
iatom.z = optcartesians[index + 2]
else:
center = ParserWrapper.FVAL_DICT[presumed_ts]
numdensepoints = ParserWrapper.NUMDENSEPOINTS
stepsize = ParserWrapper.STEPDENSEPOINTS
transitionstart = center - stepsize * numdensepoints / 2
transitionend = center + stepsize * numdensepoints / 2
if fval > transitionend:
ipoint.astructure = ppoint.astructure.copy()
for tmpindex, atomindex in enumerate(reactiveatoms):
iatom = ipoint.astructure.atom[atomindex]
index = 3 * (tmpindex)
iatom.x = optcartesians[index]
iatom.y = optcartesians[index + 1]
iatom.z = optcartesians[index + 2]
if fval >= transitionstart and fval <= transitionend:
for rbond in rpoint.astructure.bond:
if not ppoint.astructure.areBound(rbond.atom1, rbond.atom2):
add_or_update_bond_type(ipoint.astructure, rbond.atom1,
rbond.atom2,
structure.BondType.Zero)
for pbond in ppoint.astructure.bond:
if not rpoint.astructure.areBound(pbond.atom1, pbond.atom2):
add_or_update_bond_type(ipoint.astructure, pbond.atom1,
pbond.atom2,
structure.BondType.Zero)
ipoint.astructure.retype()
if not set(tosuperpose).issubset(set(reactiveatoms)):
armsd = rmsd.superimpose(rpoint.astructure,
tosuperpose,
ipoint.astructure,
tosuperpose,
use_symmetry=False,
move_which=rmsd.CT)
ipoint.astructure.property[CheckInput.TITLEKEY] = ipoint.name
ipoint.astructure.property[CheckInput.ENTRYNAMEKEY] = ipoint.name
ipoint.astructure.property[self.RXNINDEX] = ipoint.index
ipoint.astructure.property[self.RXNCOORD] = ipoint.fval
def get_final_internals(iinternals, istructure, refinternals):
"""
For every internal coordinate in the reference set define one
for the interpolated structure. While doing so grab, for the
reactive internals, both the initial interpolated internal as well
as the final fitted internal so that we can print them.
:type iinternals: list of Coord
:param iinternals: Coords for interpolated point
:type istructure: schrodinger.structure.Structure
:param istructure: structure for interpolated point
:type refinternals: list of Coord
:param refinternals: Coords for reference point
:rtype: list of Coord, list of Coord
:return: reactiveinternals, allinternals, contains reactive and all Coords
"""
# get dictionary mapping internal indicies to values for fitted internals
fittedinternals = {}
for iinternal in iinternals:
fittedinternals[tuple(iinternal.indicies)] = iinternal.value[0]
# initialize two lists, one for printing information about the fitted
# internals and one for updating the interpolated Point with the final
# set of internals, define internal coordinates for the interpolated
# point using the set of reference interal coordinates
reactiveinternals = []
allinternals = []
for refinternal in refinternals:
iatoms = [
istructure.atom[index] for index in refinternal.indicies
]
if len(iatoms) == 2:
value = measure.measure_distance(*iatoms)
if len(iatoms) == 3:
value = measure.measure_bond_angle(*iatoms)
if len(iatoms) == 4:
value = measure.measure_dihedral_angle(*iatoms)
icoord = Coord(refinternal.indicies, refinternal.names, [value])
allinternals.append(icoord)
# if this internal was fitted then collect both the interpolated and
# fitted values for printing, be sure to include cases where the
# interpolated coordinate is zero
interpolval = fittedinternals.get(tuple(refinternal.indicies))
if interpolval is not None:
icoord = Coord(refinternal.indicies, refinternal.names,
[interpolval, value])
reactiveinternals.append(icoord)
return reactiveinternals, allinternals
# update the internal coordinates in the Point object of the interpolated
# point, also for those reactive internal coordinates, i.e. the ones that
# are changing along the reaction path, collect the original interpolated
# values as well as the final values predicted by the non-linear fit this
# way we can do a nice print out of what has happened
reactiveibonds, allibonds = get_final_internals(ipoint.internals.bonds,
ipoint.astructure,
rpoint.internals.bonds)
reactiveiangles, alliangles = get_final_internals(
ipoint.internals.angles, ipoint.astructure, rpoint.internals.angles)
reactiveidihedrals, allidihedrals = get_final_internals(
ipoint.internals.dihedrals, ipoint.astructure,
rpoint.internals.dihedrals)
reactiveinternals = InternalCoords()
reactiveinternals.bonds = list(reactiveibonds)
reactiveinternals.angles = list(reactiveiangles)
reactiveinternals.dihedrals = list(reactiveidihedrals)
ipoint.internals.bonds = list(allibonds)
ipoint.internals.angles = list(alliangles)
ipoint.internals.dihedrals = list(allidihedrals)
# update the Cartesian coordinates in the Point object of the interpolated
# point
ipoint.cartesians = ReactionCoords().getCartesianCoords(
ipoint.astructure)
# return optcartesians of the interpolated structure after the
# superposition to help form the solution guess for the next
# reaction path point and while we are at it mark the reactive atoms
# in the interpolated path point structure object
optcartesianssuperposed = []
for atomindex in reactiveatoms:
iatom = ipoint.astructure.atom[atomindex]
optcartesianssuperposed.append(iatom.x)
optcartesianssuperposed.append(iatom.y)
optcartesianssuperposed.append(iatom.z)
ipoint.astructure.atom[atomindex].property[self.REACTIVEATOM] = True
return optcartesianssuperposed, reactiveinternals
[docs] def runIt(self,
reactant,
product,
job_name=ParserWrapper.JOB_NAME,
sample=ParserWrapper.SAMPLE_DEFAULT,
presumed_ts=ParserWrapper.PRESUMED_TS,
densearound=ParserWrapper.DENSEAROUND,
bondweight=ParserWrapper.BONDWEIGHT,
angleweight=ParserWrapper.ANGLEWEIGHT,
dihedralweight=ParserWrapper.DIHEDRALWEIGHT,
cartesianweight=ParserWrapper.CARTESIANWEIGHT,
penaltyweight=ParserWrapper.PENALTYWEIGHT,
interpolation=ParserWrapper.CARTESIAN,
mixprevious=ParserWrapper.MIXPREVIOUS,
guess=ParserWrapper.BEFORESUPERPOSITION,
connectivity=ParserWrapper.TS,
norxncomplex=ParserWrapper.NORXNCOMPLEX,
vdwscale=ParserWrapper.VDWSCALE,
reorder=ParserWrapper.REORDER,
reverse_interpolation=ParserWrapper.REVERSE_INTERPOLATION,
logger=None):
"""
Function to orchestrate calculation of the reaction path.
:type reactant: schrodinger.structure.Structure
:param reactant: reactant
:type product: schrodinger.structure.Structure
:param product: product
:type job_name: str
:param job_name: name of job
:type sample: list of floats
:param sample: contains either the list of sample points or
the number of points to sample as a "decimal-less" float.
:type presumed_ts: str
:param presumed_ts: gives the location of the presumed ts.
:type densearound: bool
:param densearound: Specifies if additional sampling points
should be included in the interpolation.
:type bondweight: float
:param bondweight: Specifies the weight of the bonding term
in the objective function which is minimized using non-linear
least squares.
:type angleweight: float
:param angleweight: Specifies the weight of the angle term
in the objective function which is minimized using non-linear
least squares.
:type dihedralweight: float
:param dihedralweight: Specifies the weight of the dihedral term
in the objective function which is minimized using non-linear
least squares.
:type cartesianweight: float
:param cartesianweight: Specifies the weight of the Cartesian term
in the objective function which is minimized using non-linear
least squares.
:type penaltyweight: float
:param penaltyweight: Specifies the weight of the bond penalty term
in the objective function which is minimized using non-linear
least squares.
:type interpolation: str
:param interpolation: specifies the coordinate system in which
the reaction path points are interpolated.
:type mixprevious: float
:param mixprevious: specifies to what extent the optimized
Cartesian coordinates from the previous reaction path point
are mixed with the coordinates determined by interpolation
at the current point.
:type guess: str
:param guess: specifies the type of solution guess generated
from the optimized Cartesian coordinates of the previous
reaction path point.
:type connectivity: str
:param connectivity: specifies the type of connectivity
to use in defining the structure objects of points along
the reaction path.
:type norxncomplex: boolean
:param norxncomplex: Disables the preprocessing of reactants and
products into reaction complexes for certain bimolecular reactions.
:type vdwscale: float
:param vdwscale: Scales the inter-molecular VDW distance used to
separate reactant or product structures when forming reaction
complexes for certain bimolecular reactions.
:type reorder: boolean
:param reorder: Specifies to run the protocol to normalize the atom
ordering in the reactants and products.
:type reverse_interpolation: boolean
:param reverse_interpolation: interpolate the reaction in reverse
:type logger: a logging.getLogger object
:param logger: The output logger for this script.
"""
if logger:
logger.info('Working on the %s --> %s reaction:' %
(reactant.property[CheckInput.TITLEKEY],
product.property[CheckInput.TITLEKEY]))
logger.info('')
if reverse_interpolation:
interpolation_msg = (
'The reaction path for the above reaction '
'will be determined by interpolating in reverse, i.e. from '
'the products to the reactants rather than from the reactants '
'to the products.')
logger.info(interpolation_msg)
logger.info('')
# check job name
job_name = CheckInput().checkJobName(job_name, logger)
# check reorder
reorder = CheckInput().checkReorder(reorder, logger)
# check structures and pairs of structures
reactant, product = CheckInput().checkStructures(
logger, reactant, product)
reactant, product = CheckInput().checkStructurePairs(
reorder, logger, reactant, product)
# check sample, no check needed for densearound
sample = CheckInput().checkSample(sample, logger)
# check presumed_ts
presumed_ts = CheckInput().checkPresumedTs(presumed_ts, logger)
# check interpolation
interpolation = CheckInput().checkInterpolation(interpolation, logger)
# check mixprevious
mixprevious = CheckInput().checkMixPrevious(mixprevious, logger)
# check guess
guess = CheckInput().checkGuess(guess, logger)
# check connectivity
connectivity = CheckInput().checkConnectivity(connectivity, logger)
# check norxncomplex and vdwscale
norxncomplex, vdwscale = CheckInput().checkNoRxnComplex(
norxncomplex, vdwscale, logger)
# check all weights, i.e. bondweight, angleweight, dihedralweight,
# cartesianweight, penaltyweight
bondweight, angleweight, dihedralweight, cartesianweight, penaltyweight = \
CheckInput().checkWeights(bondweight, angleweight, dihedralweight,
cartesianweight, penaltyweight, logger)
# check reverse
reverse_interpolation = \
CheckInput().checkReverseInterpolation(reverse_interpolation, logger)
# get sample points
samplepoints = self.getSamplePoints(sample, densearound, presumed_ts)
samplepointsmsg = """
This calculation will sample the following %s points: %s"""
samplepointsrounded = []
for spoint in samplepoints:
samplepointsrounded.append(round(spoint, 2))
if logger:
logger.info(samplepointsmsg %
(len(samplepointsrounded), samplepointsrounded))
logger.info('')
# get reaction coordinates
rxncoordsobj = ReactionCoords()
rxncoordsobj.prepare(reactant, product, interpolation, norxncomplex, \
vdwscale, samplepoints, logger)
rcpypoint = rxncoordsobj.rcpypoint
rpoint = rxncoordsobj.rpoint
pcpypoint = rxncoordsobj.pcpypoint
ppoint = rxncoordsobj.ppoint
reactioninternals = rxncoordsobj.reactioninternals
tosuperpose = rxncoordsobj.tosuperpose
armsd = rxncoordsobj.armsd
self.reactioninternals = reactioninternals
self.tosuperpose = tosuperpose
self.armsd = armsd
# determine the reactive atoms
reactiveatoms = self.getReactiveAtoms(reactant, product)
if logger:
logger.info('These are the reactive atoms: %s' % reactiveatoms)
logger.info('')
# mark reactive atoms in reactant and product structures
for atomindex in reactiveatoms:
rpoint.astructure.atom[atomindex].property[self.REACTIVEATOM] = True
ppoint.astructure.atom[atomindex].property[self.REACTIVEATOM] = True
if rcpypoint:
rcpypoint.astructure.atom[atomindex].property[
self.REACTIVEATOM] = True
pcpypoint.astructure.atom[atomindex].property[
self.REACTIVEATOM] = True
# initialize points list with reactant and product point objects
rindex = rpoint.index
pindex = ppoint.index
if rcpypoint:
pindex = pcpypoint.index
self.points.extend([''] * pindex)
self.points[rindex - 1] = rpoint
self.points[pindex - 1] = ppoint
if rcpypoint:
self.points[rindex - 2] = rcpypoint
self.points[rindex - 1] = rpoint
self.points[pindex - 2] = ppoint
self.points[pindex - 1] = pcpypoint
# return if there is no reaction but first update the points array
rxnyesno = reactioninternals.bonds or reactioninternals.angles or \
reactioninternals.dihedrals
if not rxnyesno:
if logger:
rtitle = rpoint.name
ptitle = ppoint.name
logger.warning(
'The reactants and products you have provided '
'in %s and %s are not for a reactive system and so no '
'reaction path will be computed.' % (rtitle, ptitle))
while '' in self.points:
self.points.remove('')
pindex = rindex + 1
self.points[pindex - 1].astructure.property[self.RXNINDEX] = pindex
if rcpypoint:
self.points[pindex].astructure.property[
self.RXNINDEX] = pindex + 1
return
# set weights
self.bondweight = bondweight
self.angleweight = angleweight
self.dihedralweight = dihedralweight
self.cartesianweight = cartesianweight
self.penaltyweight = penaltyweight
# get the length of the largest atom index
maxindexwidth = len(str(reactant.atom_total))
# get level name of logger
if logger:
loglevel = logging.getLevelName(logger.getEffectiveLevel())
# get parameters for forward or reverse interpolation
fval_iter = list(enumerate(samplepoints))
first_interp_point = 0
if reverse_interpolation:
fval_iter = fval_iter[::-1]
first_interp_point = len(fval_iter) - 1
# interpolate reaction path points
for pointindex, fval in fval_iter:
# get Point object for the interpolated point
ipoint = self.interpolateReactionCoords(pointindex + rindex + 1,
fval, connectivity, rpoint,
ppoint, reactiveatoms,
logger)
self.points[pointindex + rindex] = ipoint
if logger:
msg = 'Optimizing coordinates for reaction path point %s:' % ipoint.name
msglen = len(msg)
logger.info(msg)
banner = '=' * msglen
logger.info(banner)
logger.info('')
# handle initial guess for Cartesian coordinates of
# reactive atoms and run the optimizer
if interpolation == ParserWrapper.CARTESIAN:
optcartesians = self.getInitialGuess(ipoint.cartesians,
reactiveatoms)
else:
if pointindex == first_interp_point:
guessparams = self.getInitialGuess(ipoint.cartesians,
reactiveatoms)
optcartesians = self.doNonLinearFit(ipoint, guessparams,
reactiveatoms, mixprevious,
logger)
if guess == ParserWrapper.BEFORESUPERPOSITION:
guessparams = list(optcartesians)
# update Point object for the interpolated point
optcartesianssuperposed, comprxninternals = \
self.getInterpolatedStructure(ipoint, optcartesians, reactiveatoms,
tosuperpose, fval, connectivity, presumed_ts, rpoint, ppoint)
if guess == ParserWrapper.AFTERSUPERPOSITION:
guessparams = list(optcartesianssuperposed)
if logger:
comprxninternals.printInternals('Interpolated and optimized internal ' \
'coordinates for sample point %s:' % ipoint.name, maxindexwidth,
logger)