===============================================================

VRAI-selectivity Manual

Version 1.0

Copyright (c) 2020 Sanha Lee, Jonathan M. Goodman
University of Cambridge

This documents is adapted from 'VRAI-Selectivity: better 
calculation of selectivity than transition state theory' 
Supporting Information

===============================================================

CONTENTS
1) Requirements and Setup
2) Usage
3) Example Usage

===============================================================

1) REQUIREMENTS AND SETUP

The script is currently set up to run on Python version 2.7. The following Python modules must be installed before the script can be run:
- numpy
- argparse
- os
- math
- sys
- rdkit
- bisect
- pprint

The python script is designed to analyse Gaussian output files from Gaussian 16 only.

The script will only recognise the Gaussian output file if the frequency calculation is run independently. Do not run the calculation with 'opt' and 'freq'  keywords at the same time. The following frequency keyword must be included in the calculation 'freq=(hpmodes,intmodes)'.  

The python script requires the Gaussian output files to be in .mol file format.
It is possible to convert the frequency output files to .mol file using Gaussview via File -> Convert Files.


2) CORRECT USAGE

The .mol files and .out files must be inputed in the following order.

python VRAI-selectivity.py TS1molfile.mol INTmolfile.mol P1molfile.mol P2molfile.mol TS1freq.out TS2freq.out [OPTIONAL TS2Afreq.out] [OPTIONAL TS2Bfreq.out] [OPTIONAL -I]

TS1molfile.mol - geometry of the first transition state in .mol file format
INTmolfile.mol - geometry of the intermediate in .mol file format 
P1molfile.mol - geometry of the first product in .mol file format
P2molfile.mol - geometry of the second product in .mol file format
TS1freq.out - Gaussian16 frequency calculation output file for the first transition state
INTfreq.out - Gaussian16 frequency calculation output file for the intermediate
TS2Afreq.out - Gaussian16 frequency calculation output file for TS2A (optional)
TS2Bfreq.out - Gaussian16 frequency calculation output file for TS2B (optional)
-I - optional argument for intermediate activation
-S - optional argument for TST calculation

The analysis of reactions proceeding on the potential energy surface with an intermediate should use '-I' for intermediate activation.
- When -I option is used, the TS2A and TS2B Gaussian frequency output files must be provided (the algorithm will return an error message if these are not provided)
- When -I option is not used, do not provide TS2A or TS2B Gaussian frequency output files (the algorithm will return an error message if these are provided)
- The TS2A must correspond to the P1 forming TS, connecting INT with P1
- The TS2B must correspond to the P2 forming TS, connecting INT with P2
The algorithm will read in the TS2 frequency output files and extract the free energy of the molecule. The free energy difference between TS1 and TS2 will be used to decide whether the algorithm uses TST or dynamic analysis for the selectivity calculation.

When '-S' activation is used, the algorithm will simply output the product distribution calculated from TST. This use only requires the TS2A and TS2B Gaussian 16 output files as the input. The correct execution of the code on the command line would be:

python VRAI-selectivity.py TS2Afreq.out TS2Bfreq.out -S

For bifurcation analysis, the INTmolfile should be replaced by TS2molfile which is the geometry of the second transition state. Do not use intermediate activation.

The key output will be in the following format:
- files corresponding to the major and minor product are stated
- the two μ and λ values which decides the alignment type are given 
- the length of vector |g| and the angle φ (angle between the reaction vector and vector g_ is stated

The algorithm will create a log file and print out
- whether the rdkit object was successfully created
- vectors p1_, p2_ and the imaginary eigenvector
- φ: angle between the imaginary eigenvector, a_, and vector g_ (displacement of TS2 from TS1)
- dot products, angles and magnitude of the key vectors
- uncommon bonds the algorithm has used for the selectivity prediction
- returns a warning sign for large |g_| and when the predicted selectivity is low


3) EXAMPLE USAGE

The following files are provided for this section
- testP1freq.mol
- testP2freq.mol
- testTS1freq.mol
- testTS1freq.out
- testINTfreq.mol
- testINTfreq.out
- testTS2A.out
- testTS2B.out

Test the code on the example Gaussian output files given: 
python VRAI-selectivity_v2.py testTS1freq.mol testINTfreq.mol testP1freq.mol testP2freq.mol testTS1freq.out testINTfreq.out testTS2Afreq.out testTS2Bfreq.out -I

The Gaussian output files were generated by optimising the reaction published by Singleton et al in 2006 (Ussing, B. R.; Hang, C.; Singleton, D. A., J. Am. Chem. Soc. 2006, 128, 7594–7607). The correct execution should print out the following output.

**** Analysis Completed ****
Major product is testP2freq.mol
Minor product is testP1freq.mol

mu1_ = -0.0441868111245
mu2_ = 0.0484860277906
lambda1_ = 0.432580161859
lambda2_ = 0.438546804793
|g_| = 0.483733089154
phi = 6.62625149713

The algorith will now proceed to estimate the major and minor product ratios

Product Ratio Calculation Completed:
Major Product : Minor Product ratio
82.9 : 30.1

****************************

The algorithm should also create a log file called 'VRAI-selectivity_testTS1freq.log' with all the key information.