AutomaticReactionDiscovery

Using quantum chemical computation to find important reactions without requiring human intuition.

Description

This repository contains Python code for automatically discovering chemical reactions using a freezing string method with a subsequent exact transition state search and reaction path verification by intrinsic reaction coordinate calculation.

There are two parts required in order to run the whole program. First, initial geometries and transition state jobs have to be generated by executing python ard.py <infile> from the command line. This will generate a reactions folder containing input files for each job. Each transition state job can then be run by python -n <nproc> -m <mem> tssearch.py <infile>.

It is also possible to only run a freezing string method (without exact TS search and IRC calculation) by executing python sm.py <infile>.

Optional arguments for tssearch.py and sm.py are -n <nproc> and -m <mem>, where <mem> is specified in megabytes.

Several arguments can be specified in the input file in the format arg value. The possible arguments are:

• reac_smi - Valid SMILES string describing the reactant structure
• nbreak - Maximum number of bonds that may be broken
• nform - Maximum number of bonds that may be formed
• dH_cutoff - Heat of reaction cutoff (kcal/mol)
• forcefield - Force field for 3D geometry generation (mmff94 or uff)
• geometry - Reactant and product geometry (see main.py for details)
• nsteps - Number or gradient calculations per optimization step
• nnode - Number of nodes for calculation of interpolation distance
• lsf - Line search factor for Newton-Raphson optimization
• tol - Perpendicular gradient tolerance
• nLSTnodes - Number of high density LST nodes
• qprog - Program for quantum calculations (currently only gau)
• theory - Level of theory (e.g., m062x/cc-pvtz)
• name - Name of the log file for a TS search

Only reac_smi has to be specified for a full automatic discovery, all other arguments have default values or are set by running ard.py. Default values are:

• nbreak = 3
• nform = 3
• dH_cutoff = 20
• forcefield = mmff94
• nsteps = 4
• nnode = 15
• lsf = 0.7
• tol = 0.1
• nLSTnodes = 100
• qprog = gau
• theory = m062x/cc-pvtz
• name = 0000

The input file arguments can be specified in any order and comments can be added. An example is given in input.txt. If a freezing string method or single transition state search is run, geometry has to be specified.

Gaussian 03 or 09 has to be available at the command line and the appropriate environment variables have to be set.

Output

Running ard.py produces ARD.log, which is a log file containing all feasible reaction channels.

Each transition state search produces several files:

• output.####.log - Log file
• reac.####.out - Optimized reactant
• prod.####.out - Optimized product (intended or unintended)
• ts.####.out - Optimized transition state
• string.####.out - Nodes along FSM string
• irc.####.out - Points along IRC path
• bond_changes.out - Distance matrices during FSM steps
• fsmpath.####.png - Graph of FSM energies
• ircpath.####.png - Graph of IRC energies

There are also several Gaussian input, log, and checkpoint files that are written and rewritten over the course of the calculation.