In [1]:
# import the necessary package elements
import numpy as np
from kie import KIE_Calculation
In [2]:
claisen_calculation = KIE_Calculation("../test/claisen_demo.config", "../test/claisen_gs.out", "../test/claisen_ts.out", style="g09")
To print the KIEs:
In [3]:
print "The claisen_calculation object:"
print claisen_calculation
We can also access the KIEs directly via the KIE dictionary belonging to our KIE_Calculation. This can be useful for automating KIE analyses over a large number of files. This functionality is further developed in the autoquiver routine (see below).
In [4]:
print "Iterating through the KIEs dictionary:"
for name,kie_object in claisen_calculation.KIES.iteritems():
# prettily print the underlying KIE object:
print kie_object
# or pull the name and value directly:
print type(kie_object)
print kie_object.name
print kie_object.value
break
System objectsWe didn't have to specify file paths as the targets of our KIE_Calculation constructor. Instead, we can work directly with System objects, which contain the geometry and Hessian as fields. Below, we load the Claisen ground state and transition state and print the position of the first atom in the ground state and print the first row of the transition state Hessian.
In [5]:
from quiver import System
gs = System("../test/claisen_gs.out")
ts = System("../test/claisen_ts.out")
print "Position of atom 0 in the ground state:"
print gs.positions[0]
print "First row of the Carteisan transition state Hessian:"
print ts.hessian[0]
In [6]:
claisen_calculation2 = KIE_Calculation("../test/claisen_demo.config", gs, ts)
print claisen_calculation2
IsotopologuesOnce we have access to the underlying System objects, it is easy to make substituted Isotopologues and perform frequency calculations ourselves. To make an Isotopologue we need to provide a name, a corresponding System, and a list of masses - one for each atom in the System.
Let's build the default light ground state Isotopologue:
In [7]:
from quiver import Isotopologue
from constants import DEFAULT_MASSES
# we build the default masses by using the DEFAULT_MASSES dictionary which maps from atomic numbers to masses
# default masses in DEFAULT_MASSES are the average atomic weight of each element
gs_masses = [DEFAULT_MASSES[z] for z in gs.atomic_numbers]
gs_light = Isotopologue("the light ground state", gs, gs_masses)
print gs_light
In [8]:
# make a copy of gs_masses
sub_masses = list(gs_masses)
# index 3 corresponds to atom number 4
sub_masses[3] = 5000.0
gs_heavy = Isotopologue("the super heavy ground state", gs, sub_masses)
KIEs are essentially ratios of reduced isotopic partition functions (RPFRs). To calculate these, we use the function calculate_rpfr from the kie module. This function takes a tuple of the form (light, heavy), a frequency threshold, a scaling factor, and a temperature. (All of these are discussed in detail in the README.)
calculate_rpfr returns four values in a tuple: the first value is the RPFR, the second value is a ratio of large imaginary frequencies (if present), the third value is the frequencies in the light isotopomer, and the fourth value is the frequencies in the heavy isotopomer.
To print the individual contributions to the RPFR, set quiver.DEBUG = True in quiver.py and restart the kernel.
In [9]:
from kie import calculate_rpfr
gs_rpfr, gs_imag_ratio, gs_light_freqs, gs_heavy_freqs = calculate_rpfr((gs_light, gs_heavy), 50.0, 1.0, 273)
print gs_rpfr
print gs_light_freqs
It's nice to be able to calculate RPFRs by hand, but there is also an object available to calculate KIEs automatically: the KIE class.
We need to create a KIE object by passing it a pair of ground states (light and heavy) and a pair of transition states (light and heavy) as well as some temperature information.
We will make the desired substitution at Carbon 4 in the transition state and use a KIE object to calculate the KIE.
In [12]:
ts_masses = [DEFAULT_MASSES[z] for z in gs.atomic_numbers]
ts_light = Isotopologue("the light transition state", ts, ts_masses)
print ts_light
print
ts_sub_masses = list(ts_masses)
ts_sub_masses[3] = 5000.0
ts_heavy = Isotopologue("the heavy transition state", ts, ts_sub_masses)
print ts_heavy
gs_tuple = (gs_light, gs_heavy)
ts_tuple = (ts_light, ts_heavy)
from kie import KIE
# we make a KIE object using the gs and ts tuples from above for a calculation at 273 degrees K, with a scaling factor of 1.0, and a frequency threshold of 50 cm^-1
carbon5000_kie = KIE("Carbon 5000 at C4 KIE", gs_tuple, ts_tuple, 273, 1.0, 50.0)
Now, just like before, we have access to a KIE object (earlier we pulled them from the KIES dictionary). These can be printed prettily or you can take the value directly. As expected, we have large normal isotope effects because we used a carbon-5000 nucleus! (The three numbers correspond to the uncorrected, Wigner, and Bell KIEs.)
In [14]:
print carbon5000_kie
print carbon5000_kie.value
autoquiver.pyA common problem when computing KIEs is the need to screen a large number of levels of theory with numerous ground and transition state files. So far, we have only seen how to make a KIE_Calculation object that corresponds to a single configuration, ground state, and transition state.
The autoquiver module implements functionality to perform screening over a large number of files. Its command line use is described in depth in the README. Here we will see how to leverage Python to have more sophisticated uses for autoquiver. Suppose we have the following files that we want to use as inputs to PyQuiver:
config.config
ground_state1.out ts-1.out
gs1rotomer.out tsrotomer1.out
GROUNDSTATE2.out transition_state2.out
gs3.out ts3.outWe want to run PyQuiver on pair of files above in the same row using the configuration file config.config. The command line functionality of autoquiver is great when the file names have a consistent pattern, but here it's not obvious how to generate the proper pairings or detect which files are ground or transition states.
The function autoquiver (which is the only object available in the autoquiver module) requires the following arguments:
filepath: this is the target directory.config_path: the location of the configuration file to use for all PyQuiver calculationsinput_extension: the file extension for ground and transition state files (default=.out)style: the style of the ground and transition state files (default=g09)In the README, we discussed the command:
python src/autoquiver.py -e .output auto/ auto/substitutions.config gs ts -The strings "gs" and "ts" were used to find ground and transition files. A ground state file needed to have "gs" as a substring and likewise for "gs" and transition state files. The character "-" was used as a field delimiter and "gs" and "ts" matched if all fields after the first were identical.
That is some sensible default behaviour, but for more complex cases, we can provide the following functions:
gs_p: a function that takes a filename and returns a boolean value. True if the filename corresponds to a ground state and False if it does not.ts_p: a function that takes a filename and returns a boolean value. True if the filename corresponds to a transition state and False if it does not.gs_ts_match_p: a function that takes a ground state filename and a transition state filename and returns a boolean value. True if the ground state and transition state match and Fase if they do not.gs_p and ts_pFor the previous example, we need to detect that all of the following are ground state files:
ground_state1.out
gs1_rotomer.out
GROUNDSTATE2.out
gs3.outWe can accomplish this with the following Python function:
In [16]:
def gs_p(filename):
# convert to lowercase
filename = filename.lower()
# remove the _ character
filename = filename.replace("_", "")
# replace groundstate with gs
filename = filename.replace("groundstate", "gs")
# check if the modified filename begins with the string "gs"
if filename[:2] == "gs":
return True
else:
return False
Let's test this on all the filenames in the example:
In [17]:
gs_filenames = ["ground_state1.out", "gs1_rotomer.out", "GROUNDSTATE2.out", "gs3.out"]
ts_filenames = ["ts-1.out", "tsrotomer1.out", "transition_state2.out", "ts3.out"]
print "Checking ground state matches:"
for n in gs_filenames:
print gs_p(n)
print "Checking transition state matches:"
for n in ts_filenames:
print gs_p(n)
Lo and behold! The function successfully separated the ground state and transition state files. We can similarly implement ts_p:
In [18]:
def ts_p(filename):
# convert to lowercase
filename = filename.lower()
# remove the _ character
filename = filename.replace("_", "")
# replace groundstate with gs
filename = filename.replace("transitionstate", "ts")
# check if the modified filename begins with the string "gs"
if filename[:2] == "ts":
return True
else:
return False
gs_ts_match_pNow all we have to do is write a function gs_ts_match_p that can take two filenames (one ground state and one transition state) and detect if they are matches. There are numerous ways to accomplish this. We will use a basic regular expression to highlight a common method for writing these functions.
In [19]:
gs_filenames = ["ground_state1.out", "gs1_rotomer.out", "GROUNDSTATE2.out", "gs3.out"]
ts_filenames = ["ts-1.out", "tsrotomer1.out", "transition_state2.out", "ts3.out"]
import re
def gs_ts_p(gs_name, ts_name):
# a regular expression that finds and extracts the first integer substring in a filename
gs_match = re.search(".*([0-9]+).*", gs_name)
gs_number = gs_match.group(1)
ts_match = re.search(".*([0-9]+).*", ts_name)
ts_number = ts_match.group(1)
if ts_number == gs_number:
return True
else:
return False
import itertools
for gs_n, ts_n in itertools.product(gs_filenames, ts_filenames):
if gs_ts_p(gs_n, ts_n):
print "Match: {0} {1}".format(gs_n, ts_n)
In [ ]:
# missing example of how to actually use autoquiver