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

Return type

float, float, float, float, float

Returns

bondweight, angleweight, dihedralweight, cartesianweight, penaltyweight, weights of the bond, angle, dihedral, Cartesian, and penalty terms

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.

Raises

ValueError – if there is a problem with the input

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

• reactant (schrodinger.structure.Structure) – the reactant

• product (schrodinger.structure.Structure) – the product

• logger (logging.getLogger) – output logger

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

• reactant (schrodinger.structure.Structure) – reactant

• product (schrodinger.structure.Structure) – product

• vdwscale (float) – scales the intermolecular distance

• logger (logging.getLogger) – output logger

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

• reactant (schrodinger.structure.Structure) – reactant

• product (schrodinger.structure.Structure) – product

• rinternals (InternalCoords) – reactant internal coordinates

• pinternals (InternalCoords) – product internal coordinates

• logger (logging.getLogger) – output logger

getReactionInternals(rinternals, pinternals, reactioninternals)

Determine the reactive redundant internal coordinates.

Parameters

• rinternals (InternalCoords) – reactant internal coordinates

• pinternals (InternalCoords) – product internal coordinates

• reactioninternals (InternalCoords) – reaction internal coordinates

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

• reactant (schrodinger.structure.Structure) – reactant

• product (schrodinger.structure.Structure) – product

• reactioninternals (InternalCoords) – reaction internal coordinates

• logger (logging.getLogger) – output logger

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

• reactant (schrodinger.structure.Structure) – reactant

• product (schrodinger.structure.Structure) – product

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.