schrodinger.application.matsci.rxn_path module

Classes and functions for generating reaction paths.

Copyright Schrodinger, LLC. All rights reserved.

class schrodinger.application.matsci.rxn_path.ParserWrapper(scriptname, description)

Bases: 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.0, 0.25, 0.5, 0.75, 1.0]
FVAL_DICT = {'midway': 0.5, 'product': 1.0, 'product_like': 0.75, 'reactant': 0.0, 'reactant_like': 0.25}
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
__init__(scriptname, description)

Create a ParserWrapper object and process it.

Parameters:
  • scriptname (str) – name of this script
  • description (str) – description of this script
loadIt()

Load ParserWrapper with options.

parseArgs(args)

Parse the command line arguments.

Parameters:args (tuple) – command line arguments
class schrodinger.application.matsci.rxn_path.CheckInput

Bases: object

Check user input.

PAIRS = [('[]', ''), ('][', '-'), ('[', ''), (']', ''), (';', ''), (' ', '_')]
COMBIGLD_REPLACEMENTS = {' ': '_', ';': '', '[': '', '[]': '', ']': '', '][': '-'}
TITLEKEY = 's_m_title'
ENTRYNAMEKEY = 's_m_entry_name'
checkJobName(job_name, logger=None)

Check job_name option.

Parameters:
  • job_name (str) – name of job
  • logger (logging.getLogger) – output logger
Return type:

str

Returns:

job_name, name of job

checkInputFiles(inputfiles, logger)

Check input files.

Parameters:
  • inputfiles (list of str) – all provided input files
  • logger (logging.getLogger) – output logger
checkStructures(logger=None, *allstructures)

Check structures.

Parameters:
  • logger (logging.getLogger) – output logger
  • allstructures (tuple of schrodinger.structure.Structure) – all provided structures
Return type:

list of schrodinger.structure.Structure

Returns:

structures, updated list of structures

checkStructurePairs(reorder, logger=None, *allstructures)

Check reactant and product pairs of structures.

Parameters:
  • reorder (boolean) – normalize the atomic ordering
  • logger (logging.getLogger) – output logger
  • allstructures (tuple of schrodinger.structure.Structure) – all provided structure
Return type:

list of schrodinger.structure.Structure

Returns:

structures, updated list of structures

checkSample(sample, logger=None)

Check sample option.

Parameters:
  • sample (list of float) – sample points
  • logger (logging.getLogger) – output logger
Return type:

list of floats

Returns:

sample, sample points

checkPresumedTs(presumed_ts, logger=None)

Check the location of the presumed ts.

Parameters:
  • presumed_ts (str) – gives the location of the presumed ts
  • logger (logging.getLogger) – output logger
Return type:

str

Returns:

presumed_ts, the location of the presumed ts

checkInterpolation(interpolation, logger=None)

Check interpolation option.

Parameters:
  • interpolation (str) – coordinate system used for interpolating
  • logger (logging.getLogger) – output logger
Return type:

str

Returns:

interpolation, coordinate system used for interpolating

checkMixPrevious(mixprevious, logger=None)

Check mixprevious option.

Parameters:
  • mixprevious (float) – mixing weight of solution from previous path point
  • logger (logging.getLogger) – output logger
Return type:

float

Returns:

mixprevious, mixing weight of solution from previous path point

checkGuess(guess, logger=None)

Check guess option.

Parameters:
  • guess (str) – type of guess solution
  • logger (logging.getLogger) – output logger
Return type:

str

Returns:

guess, type of guess solution

checkConnectivity(connectivity, logger=None)

Check connectivity option.

Parameters:
  • connectivity (str) – type of connectivity
  • logger (logging.getLogger) – output logger
Return type:

str

Returns:

connectivity, type of connectivity

checkNoRxnComplex(norxncomplex, vdwscale, logger=None)

Check norxncomplex and vdwscale options.

Parameters:
  • norxncomplex (boolean) – disable preprocessing into a reaction complex
  • vdwscale (float) – scales the intermolecular distance
  • logger (logging.getLogger) – output logger
Return type:

boolean, float

Returns:

norxncomplex, vdwscale, disable preprocessing into a reaction complex and scales the intermolecular distance

checkReorder(reorder, logger=None)

Check reorder option.

Parameters:
  • reorder (boolean) – normalize the atomic ordering
  • logger (logging.getLogger) – output logger
Return type:

boolean

Returns:

reorder, normalize the atomic ordering

checkReverseInterpolation(reverse_interpolation, logger=None)

Check reverse interpolation option.

Parameters:
  • reverse_interpolation (boolean) – interpolate the reaction path in reverse
  • logger (logging.getLogger) – output logger
Return type:

boolean

Returns:

reverse_interpolation, interpolate the reaction path in reverse

checkWeights(bondweight, angleweight, dihedralweight, cartesianweight, penaltyweight, logger=None)

Check weights, i.e. bondweight, angleweight, dihedralweight, cartesianweight, and penaltyweight.

Parameters:
  • bondweight (float) – weight of the bond term
  • angleweight (float) – weight of the angle term
  • dihedralweight (float) – weight of the dihedral term
  • cartesianweight (float) – weight of the Cartesian term
  • penaltyweight (float) – weight of the bond penalty term
  • logger (logging.getLogger) – output logger

:rtype float, float, float, float, float :return: bondweight, angleweight, dihedralweight, cartesianweight,

penaltyweight, weights of the bond, angle, dihedral, Cartesian, and penalty terms
__init__

Initialize self. See help(type(self)) for accurate signature.

class schrodinger.application.matsci.rxn_path.Coord(indicies, names, value)

Bases: object

Manage the properties of an internal coordinate.

BOND = 'bond'
ANGLE = 'angle'
DIHEDRAL = 'dihedral'
__init__(indicies, names, value)

Create a Coord instance.

Parameters:
  • indiciees – atomic indicies
  • names (list of str) – atomic names
  • value (list of float) – value(s) of the internal coordinate in units of Angstrom or degree
class schrodinger.application.matsci.rxn_path.InternalCoords

Bases: object

Manage the internal coordinates of a structure.

ZMATNUMBER = 1
__init__()

Create an InternalCoords instance.

getZmatrix(astructure)

Get the Z-matrix for the structure.

Parameters:astructure (schrodinger.structure.Structure) – the structure
getDmatrix(astructure)

Get the distance matrix for the structure.

Parameters:astructure (schrodinger.structure.Structure) – the structure
printInternals(headermsg, maxindexwidth, logger)

Formatted print of header followed by the internal coordinates.

Parameters:
  • headermsg (str) – header
  • maxindexwidth (int) – number of characters in the largest atom index
  • logger (logging.getLogger) – output logger
schrodinger.application.matsci.rxn_path.max_pair_vdw_distance(astructure)

Find the largest atom-atom VDW distance in a structure.

Parameters:astructure (schrodinger.structure.Structure) – the structure
Return type:int, int, float
Returns: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.
schrodinger.application.matsci.rxn_path.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.

Parameters:
  • astructure (schrodinger.structure.Structure) – the structure containing the atom to which a hydrogen will be added
  • index (int) – the index of the atom to which to add the hydrogen
Return type:

int

Returns:

the index of the added temporary hydrogen

class schrodinger.application.matsci.rxn_path.ReactionCoords

Bases: object

Manage reaction coordinates.

REACTIONBONDTHRESH = 0.05
REACTIONANGLETHRESH = 1.0
REACTIONDIHEDRALTHRESH = 1.0
REVOLUTION = 360
HALFREVOLUTION = 180
REVOLUTIONTHRESH = 10
REACTANTPREFIX = 'pre-'
PRODUCTPREFIX = 'post-'
__init__()

Create a ReactionCoords instance.

getNormalOrdering(reactant, product, logger=None)

Attempt to normalize the atomic ordering between reactants and products.

Parameters:
Return type:

schrodinger.structure.Structure, schrodinger.structure.Structure

Returns:

newreactant, newproduct, if defined specifies the reordered reactant and product structures

makeRxnComplex(reactant, product, vdwscale, logger=None)

For certain bimolecular reactions preprocess reactants and products into reaction complexes.

Parameters:
Return type:

schrodinger.structure.Structure, schrodinger.structure.Structure

Returns:

newreactant, newproduct, if defined specifies the reactant and product structures in the created reaction complex

collectInternals(reactant, product, rinternals, pinternals, logger=None)

Find the redundant internal coordinates by merging the coordinates defined in the reactant and product.

Parameters:
getReactionInternals(rinternals, pinternals, reactioninternals)

Determine the reactive redundant internal coordinates.

Parameters:
runSuperposition(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.

Parameters:
Return type:

list of ints, float

Returns:

tosuperpose, armsd, atom indicies used to superpose and the final RMSD

getCartesianCoords(astructure)

Get the Cartesian coordinates of a structure as a 3N dimensional list of floats.

Parameters:astructure (schrodinger.structure.Structure) – structure
Return type:list of float
Returns:cartesians, the 3N Cartesian coordinates ordered like [x1, y1, z1, x2, …, zN]
prepare(reactant, product, interpolation, norxncomplex, vdwscale, samplepoints, logger=None)

Prepare reaction coordinates.

Parameters:
  • reactant (schrodinger.structure.Structure) – reactant
  • product (schrodinger.structure.Structure) – product
  • interpolation (str) – coordinate system used for interpolating
  • norxncomplex (boolean) – disable preprocessing into a reaction complex
  • vdwscale (float) – scales the intermolecular distance
  • samplepoints (list of float) – reaction path sample points
  • logger (logging.getLogger) – output logger
class schrodinger.application.matsci.rxn_path.Point(index, fval, name, astructure, internals, cartesians)

Bases: object

Collect properties of reaction path points.

__init__(index, fval, name, astructure, internals, cartesians)

Create a Point instance.

Parameters:
  • index (int) – path point index
  • fval (float) – path point value
  • name (str) – path point name
  • astructure (schrodinger.structure.Structure) – structure
  • internals (InternalCoords) – path point internal coordinates
  • cartesians (list of float) – the 3N Cartesian coordinates ordered like [x1, y1, z1, x2, …, zN]
class schrodinger.application.matsci.rxn_path.ReactionPath

Bases: object

Generate reaction path.

FVALINCREMENT = 0.001
NORMTHRESH = 1e-12
BONDPENALTYTHRESH = 1000000000.0
DIFFLOWTHRESH = 10.0
DIFFHIGHVAL = 1000000000.0
RXNINDEX = 'i_matsci_RXN_Index'
RXNCOORD = 'r_matsci_RXN_Coord'
REACTIVEATOM = 'b_matsci_Reactive_Atom'
__init__()

Create a ReactionPath instance.

getSamplePoints(sample, densearound, presumed_ts)

Determine the final set of sampling points.

Parameters:
  • sample (list of float) – points to be interpolated
  • densearound (bool) – include additional sampling points at specific regions
  • presumed_ts (str) – location of presumed ts
Return type:

list of floats

Returns:

samplepoints, the list of points to be sampled.

getReactiveAtoms(reactant, product)

Determine the reactive atoms, i.e. those which have changed Cartesian positions in the superposed reactant/product pair.

Parameters:
Return type:

list of int

Returns:

reactiveatoms, reactive atoms

interpolateReactionCoords(pointindex, fval, connectivity, rpoint, ppoint, reactiveatoms, logger=None)

Interpolate reaction coordinates between the reactant and product for this sample point.

Parameters:
  • pointindex (int) – sample point index
  • fval (float) – interpolated reaction path point
  • connectivity (str) – specifies the type of connectivity
  • rpoint (Point) – reactant information
  • ppoint (Point) – product information
  • reactiveatoms (list of ints) – reactive atoms
  • logger (logging.getLogger) – output logger
Return type:

Point

Returns:

ipoint, interpolated point information

getInitialGuess(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.

Parameters:
  • cartesians (list of floats) – all interpolated Cartesian coordinates
  • reactiveatoms (list of ints) – atom indicies of reactive atoms
Return type:

list of floats

Returns:

guessparams, interpolated Cartesian coordinates of the reactive atoms

doNonLinearFit(ipoint, guessparams, reactiveatoms, mixprevious, logger=None)

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.

Parameters:
  • ipoint (Point) – interpolated point information
  • guessparams (list of floats) – initial parameters, i.e. Cartesian coordinates of the reactive atoms
  • reactiveatoms (list of ints) – atomic indicies of reactive atoms
  • mixprevious (float) – 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.
  • logger (logging.getLogger) – output logger
Return type:

list

Returns:

optcartesians, non-linear-optimized Cartesian coordinates for this sample point.

getInterpolatedStructure(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.

Parameters:
  • ipoint (Point) – interpolated point information
  • optcartesians (list of floats) – optimized Cartesian coordinates for the reactive atoms
  • reactiveatoms (list of ints) – atom indicies of reactive atoms.
  • tosuperpose (list of ints) – contains the atom indicies of the atoms used in the superposition.
  • fval (float) – The interpolated reaction path point.
  • connectivity (str) – specifies the type of connectivity
  • presumed_ts (str) – specifies the location of the presumed ts
  • rpoint (Point) – reactant information
  • ppoint (Point) – product information
Return type:

list of floats, InternalCoords object

Returns:

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.

runIt(reactant, product, job_name='rxnpath', sample=[10.0], presumed_ts='midway', densearound=False, bondweight=1000.0, angleweight=1000.0, dihedralweight=1000.0, cartesianweight=1000.0, penaltyweight=1.0, interpolation='cartesian', mixprevious=0.5, guess='beforesuperposition', connectivity='ts', norxncomplex=False, vdwscale=1.0, reorder=False, reverse_interpolation=False, logger=None)

Function to orchestrate calculation of the reaction path.

Parameters:
  • reactant (schrodinger.structure.Structure) – reactant
  • product (schrodinger.structure.Structure) – product
  • job_name (str) – name of job
  • sample (list of floats) – contains either the list of sample points or the number of points to sample as a “decimal-less” float.
  • presumed_ts (str) – gives the location of the presumed ts.
  • densearound (bool) – Specifies if additional sampling points should be included in the interpolation.
  • bondweight (float) – Specifies the weight of the bonding term in the objective function which is minimized using non-linear least squares.
  • angleweight (float) – Specifies the weight of the angle term in the objective function which is minimized using non-linear least squares.
  • dihedralweight (float) – Specifies the weight of the dihedral term in the objective function which is minimized using non-linear least squares.
  • cartesianweight (float) – Specifies the weight of the Cartesian term in the objective function which is minimized using non-linear least squares.
  • penaltyweight (float) – Specifies the weight of the bond penalty term in the objective function which is minimized using non-linear least squares.
  • interpolation (str) – specifies the coordinate system in which the reaction path points are interpolated.
  • mixprevious (float) – 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.
  • guess (str) – specifies the type of solution guess generated from the optimized Cartesian coordinates of the previous reaction path point.
  • connectivity (str) – specifies the type of connectivity to use in defining the structure objects of points along the reaction path.
  • norxncomplex (boolean) – Disables the preprocessing of reactants and products into reaction complexes for certain bimolecular reactions.
  • vdwscale (float) – Scales the inter-molecular VDW distance used to separate reactant or product structures when forming reaction complexes for certain bimolecular reactions.
  • reorder (boolean) – Specifies to run the protocol to normalize the atom ordering in the reactants and products.
  • reverse_interpolation (boolean) – interpolate the reaction in reverse
  • logger (a logging.getLogger object) – The output logger for this script.