Calculation of Heavy Isotope Effects:
Case Study and Tutorial

Suppose the reader is an experimentalist (or a strayed theoretician) well acquainted with the theory of kinetic isotope effects (KIE; Refs. 1-3). He is well aware of the possibility to calculate these effects from first principles, and knows exactly why he needs that calculation. It will just take him to look into one comprehensive example to get started with his own calculations. The goal of this tutorial is to present such an example for calculation of heavy isotope effects. We are going to reproduce the observed 14N/15N KIE for the reaction of nitrogen cleavage by molybdenum complexes (Ref. 4):


The KIE of a reaction is determined by the vibration frequencies of the reagent state and transition state (highlighted red above). The general workflow for calculation will be as following:

The toolbox

The software toolbox to get through all those steps will be as following (if not stated otherwise, all software is free for academy):

1. Molecular modeling package to draw structures and get rough 3D-placement of atoms. Absolutely the best modeling tool, capable of building metal complexes, as well as small organic molecules, proteins, and nucleic acids is DS Visualizer from Accelrys. Here I assume there are no maniacs left willing to manually enter the atom coordinates in the input file.

2. Open Babel - a package to convert between plethora of molecule formats.

3. Packages to minimize the structures and compute vibration frequencies: GAMESS, Mopac2012 (free clone of Mopac; I personally never had a chance to play around its capabilities. It should be useful for those who never goes beyond standard elements of organic chemistry), Gaussian (not free).

4. Packages to visualize the calculated vibrations: Molekel (the interface to show vibrations is not very convenient, but it is able to save the movements as avi or image sequence files); Molden (very convenient visualization of vibrations; displays Raman and IR spectra, but has the ugliest UNIX-style interface and needs X-server from Cygwin to run); ChemCraft (the interface is good... when it works. Looks like it dislikes some output files (e.g. with imaginary frequencies), and it's not free).

5. ISOEFF98 (Ref.5) - a successor of BEBOVIB (Ref. 6) to calculate kinetic and equilibrium isotope effects. Available free from Prof. Piotr Paneth, Technical University of Lodz.

First steps

Once we know how the ground (GS) and transition (TS) state looks like (highlighted pink in the picture above), it's time to build draft structures of these molecules. Catch number 1: if we take a look at the structures above, there are a lot of atoms, which will have no significant effect on the KIE. These are aryl rings and isobutyl moieties. They must be present in the real molybdenum complex to provide it with the necessary stability and reactivity. In a model structure approximated by a certain force field, their presence is not necessary, so the truncation of these atoms will save us a whole lot of computer time. The extent of truncation is to be chosen mostly in accordance with available computational power. In my case, the substitution of N(C4H9)Ar with NH2 in the Mo complexes brought the computer time from 2-3 days down to 7 hours, which is a substantial gain.

So, let's build the structures of GS and TS in DS Visualizer. We should get something like this:



Now, we want to save some more computing time, so we adjust the interatomic distances and valence angles according to the available structural data (if there is any). Luckily, in this particular case Ref. 4 provides plethora of structural information both on GS and TS, which we can use. Anyways, this step is optional, but, in some instances, can be crucial for the success of conformation search.

Now it's time to save each molecule in separate file. You can use any molecular format (mol, xyz, pdb...) provided it can be read by Babel. By this time you should have two files containing the structures above, say, GS.mol and TS.mol.

 

Step 2: Preparation of the input files for GAMESS, Mopac, or Gaussian and running the calculation

Of course, the crucial point is to orient yourself which particular program you are going to use and to read the user manual, of course. We are limited to these three packages, because current version of ISOEFF can read output only from these programs (but this is ample for most purposes). Historically it happened so that I had to use Gaussian 98, therefore, the example inputs will be prepared for this program.

The most tedious part in the preparation of the input files is entering the coordinates. Luckily, we have such a beautiful program as Babel, which can take almost any molecular file and convert it into Gaussian input file:

>obabel.exe -imol GS.mol -ogau GS.inp 

>obabel.exe -imol TS.mol -ogau TS.inp 

What we are getting in the .inp files are actually stubs, which must be completed with necessary Gaussian keywords, correct molecular charge and multiplicity:

-------------------GS.inp: raw output from Babel----------------------------------
#Put Keywords Here, check Charge and Multiplicity.

YUPGUL10

0 1
N -1.78240 4.77160 -8.67180 
N -4.60070 4.62300 -6.76640 
N -1.93030 2.63840 -6.06130 
Mo -0.28440 6.96110 -3.79490 
N -0.92680 5.70500 -4.84550 
N -2.00950 5.49720 -5.92650 
Mo -2.65180 4.24120 -6.97710 
N -1.13010 6.41470 -2.09290 
N -0.99960 8.56890 -4.72690 
N 1.65520 6.57450 -3.99120 
H -0.71940 4.74480 -8.55160 
H -2.08630 5.76350 -8.93390 
H -4.80660 5.62080 -7.09340 
H -4.87530 4.52160 -5.73720 
H -2.30260 2.60200 -5.05880 
H -0.86160 2.68830 -6.04290 
H -0.81450 5.42430 -1.83910 
H -2.19450 6.43270 -2.20100 
H -2.06830 8.52200 -4.75090 
H -0.62180 8.59930 -5.72760 
H 1.93660 6.67220 -5.01900 
H 1.85590 5.57700 -3.66020 
------------------------------------------------------------------------------------

Now it's time to think and think well. Because wrong choice of keywords will result in wasted computing time, abnormal termination of the calculation, and abundance of other unpleasant things. However, no matter how carefully you think, several (tens) first runs will be wasted anyways - this is a law of Nature.

The structures above contain transition metals. Therefore, we have no other options, but to use density functional theory (DFT) for calculation (and we are limited to choice between GAMESS and Gaussian, because Mopac does not have DFT methods). It would be useful to become introduced to the theory. Ref. 7 and 8 are very good introductory texts for both theoreticians and laymen (like me). For "plain" organic molecules it could be better to use more common ab initio or semi-empirical methods available in Mopac.

For the input file for the GS below, I'm not going to reason how I got to that or another keyword, but I will just explain what each keyword means.

%mem=100MW 

Dynamic memory allocated for Gaussian is 100 megawords (I didn't really care, but this was enough)

%chk=gs.chk

Scratch file name and location.

b3lyp/gen

General Becke's 3 parameter functional - most commonly used method in DFT calculations.

pseudo(read)

Read pseudopotentials for the core electrons from the bottom section (needed if general basis set is declared with /gen keyword)

6d

Use 6d functions

10f

Use 10f functions. Not sure exactly what are these two, but they made the calculations to run OK.

opt(maxcycle=500)

Perform the structure optimization (and do that to the bitter end).

scf(maxcycle=2000)

Force the SCF procedure to converge at any cost (optional).

freq(raman)

Calculate the force constants and vibrational frequencies (make sure that Raman frequencies will be calculated as well). This is the very crucial keyword for further KIE calculation. The GS of our system includes very important N-N stretching Raman frequency, that's why we do care about raman keyword.

Bottom section (needed if general basis set is declared with /gen keyword)

N H 0
6-31g** 
****
Use 6-31g** basis for N and H atoms.
Mo 0
LANL2dz
f 1
0.4 1.0
****

Mo 0
LANL2dz
Use LANL2dz basis for Mo atoms and add diffuse function on Mo atoms (no sane chemist would ever want to know what is that, but for those concerned, please see Ref.9 and Ref. 10 for more thorough DFT study of the very same reaction).

The charge of molecule is obviously 0, and the multiplicity normally would be a singlet (1). But what we are doing is not normal human activity from the very beginning, and so is the system (see Ref. 4,9,10 for proof). The ground state is a paramagnetic colored complex in a triplet state (3).

----------------GS.inp: the completed Gaussian input for the ground state-------
%mem=100MW
%chk=gs.chk
# b3lyp/gen pseudo(read) 6d 10f opt(maxcycle=500) scf(maxcycle=2000) freq(raman)

Ground state optimization and frequencies

0 3
N -0.97550 1.68970 2.78920 
N -0.97550 -1.68970 2.78920 
N 0.97550 1.68970 -2.78920 
N -1.95110 0.00000 -2.78920 
N 0.97550 -1.68970 -2.78920 
N 1.95110 0.00000 2.78920 
H -1.59970 1.76700 4.06560 
H -0.73040 -2.26890 4.06560 
H 0.73040 2.26890 -4.06560 
H -2.33010 -0.50190 -4.06560 
H 1.59970 -1.76700 -4.06560 
H 2.33010 0.50190 4.06560 
H -1.17190 2.83240 1.83350 
H -1.86700 -2.43110 1.83350 
H 1.86700 2.43110 -1.83350 
H -3.03890 0.40130 -1.83350 
H 1.17190 -2.83240 -1.83350 
H 3.03890 -0.40130 1.83350 
N 0.00000 0.00000 -0.59630 
N 0.00000 0.00000 0.59630 
Mo 0.00000 0.00000 -2.47230 
Mo 0.00000 0.00000 2.47230 

N H 0
6-31g** 
****
Mo 0
LANL2dz
f 1
0.4 1.0
****

Mo 0
LANL2dz
------------------------------------------------------------------------------------

Now feed the input file above to Gaussian 98 and get the result containing optimized geometry and frequencies in GS.log file. What is important is to check that the located geometry is a true GS minimum (i.e. there are no imaginary frequencies). Look for "NImag=" string in the log file and make sure that NImag=0. Take a deep breath - we made it, to the ground state.

Now, the transition state is trickier, and below are comments on special keyword we use to calculate TS.

rb3lyp/gen

Restricted (closed shell) general Becke's 3 parameter functional. Because TS is a singlet, we can use restricted wavefunctions to speed up the calculation.

opt(ts,noeigentest,calcfc, maxcycle=500)

Perform the structure optimization for the transition state (ts); suppress the test for the curvature (not sure if it's really needed); calculate the force constants for the initial geometry.

symm(loose)

Makes it easier for the program to determine the symmetry of the molecule.

---------------------TS.inp: completed Gaussian input for the transition state ---------- %mem=20MW
%chk=TS.chk
# rb3lyp/gen pseudo(read) 6d 10f opt(ts,noeigentest,calcfc, maxcycle=500)
scf(maxcycle=2000) freq(raman) symm(loose)

Transition Structure Search and Freq Calc

0 1
N -0.161665 0.722336 0.000000
N 0.161665 -0.722336 0.000000
Mo -1.667031 1.662125 0.000000
Mo 1.667031 -1.662125 0.000000
N -1.667031 2.661908 1.701945
H -0.995335 2.524620 2.454380
H -2.401412 3.320757 1.965154
N -2.696191 0.013323 0.000000
H -3.707584 -0.109304 0.000000
H -2.177581 -0.877325 0.000000
N -1.667031 2.661908 -1.701945
H -0.995335 2.524620 -2.454380
H -2.401412 3.320757 -1.965154
N 2.696191 -0.013323 0.000000
H 3.707584 0.109304 0.000000
H 2.177581 0.877325 0.000000
N 1.667031 -2.661908 1.701945
H 0.995335 -2.524620 2.454380
H 2.401412 -3.320757 1.965154
N 1.667031 -2.661908 -1.701945
H 0.995335 -2.524620 -2.454380
H 2.401412 -3.320757 -1.965154

N H 0
6-31g** 
****
Mo 0
LANL2dz
f 1
0.4 1.0
****

Mo 0
LANL2dz
------------------------------------------------------------------------------------

After successfully getting the output log file from Gaussian, we need to check again that the TS structure was optimized correctly. Now, TS must have exactly one imaginary frequency ("NImag=1"). If this is so, we are all set to perform KIE calculations with any imaginable isotopes using these two output log files. Save them in several different places - they are precious.

Step 3. Analysis of vibrational frequencies

We are going to find out what will be the KIE in our reaction if one isotopomer has two 14N isotopes in the middle, and another - two 15N isotopes. This is important, because the molecule is a promising scaffold for designing a low molecular nitrogenase mimic, and even more important because we have an experimental measure of this KIE (Ref. 4), so we'll be able to compare our speculations with the actuality.

As long as we want the KIE on central nitrogen atoms, we should presume that the vibrations responsible for this effect will involve these atoms. There are two complementary ways to find these vibrations: the fun one and the boring one.

The fun way involves animation of the frequencies contained in the Gaussian output files. Below the two examples are shown (produced with Molekel program). You should start the vibration animation in the Molekel window (with the Gaussian file opened) and then just visualize each frequency in the Animation/Per Molecule Setings... dialog. Below two animations are shown: Raman N-N stretch at 1837 cm-1 for the ground state, and combination of 880 and 909 cm-1 modes for the transition state.


1837cm-1 Raman stretch of the GS

 


880 and 909 cm-1 modes for the TS

 

The boring way is to analyze Raman and IR spectra for both GS and TS looking for promising vibration modes. Below these spectra are depicted (created with Molden program from the same Gaussian outputs). Qualitatively, the most changed area of the spectrum lies between 900 and 200 cm-1 and we are going to check how much does it contribute to the KIE.

 

Raman (blue) and IR (red) spectra obtained from GS.log (upper panel) and TS.log (lower panel) files using Molden program. The inset on the upper panel is the experimental Raman stretching vibration of the N-N bond for actual molybdenum purple complex (with all aryl and isobutyl substituents in place) taken in toluene (from Ref. 4). a) 14N-14N compound; b) 15N-15N compound. Area shaded pink - range of molecular vibrations contributing the most to 14N/15N KIE.

Final step: Calculation of the KIE

Let's now prepare the input file for ISOEFF98. Additionally you will need both Gaussian output files in the same directory. For general instructions it's better to consult with ISOEF manual. Below I'm giving the example of the input file, which calculates KIE using all possible vibration frequencies (but neglects trivial ones - translation and rotation of the molecules as a whole).

$IODEF SSINP = GS.log TSINP = TS.log $END

Global definition of the data files for the ground and transition states (file names of Gaussian output)

30N/28N for Mo complexes

Just a comment for humans.

LMASS(19)=14,14 HMASS(19)=15,15 

Arrays of atomic masses for light (LMASS) and heavy (HMASS) isotopomers of the GS. In this examples atoms # 19 and 20 have mass of 14 (14N) in the light form and mass of 15 (15N) in the heavy form.
LMATS(1)=14,14 HMATS(1)=15,15
The same for the TS. Note the catch: carefully check your files for the atom numbering. When the transition state was build, the central nitrogens got the numbers 1 and 2 instead of 19 and 20 in GS (my fault, but ISOEFF is really flexible in this regard). 
RMTRIV = .t. 
Remove trivial frequencies (translations and rotations of a whole molecule) from the calculations (.t. stands for TRUE). 
KIE=.t. EIE=.f.
Calculate KIE, not the equilibrium isotope effect
TKELV(1)=300,320,340
Calculate KIE at 300, 320, and 340K

-------------- KIE.inp: input file for ISOEFF98 ------------------------------------

$IODEF SSINP = GS.log TSINP = TS.log $END

30N/28N for Mo complexes
$ISOEFF LMASS(19)=14,14 HMASS(19)=15,15 
LMATS(1)=14,14 HMATS(1)=15,15 
RMTRIV = .t. 
KIE=.t. EIE=.f.
TKELV(1)=300,320,340 $END
------------------------------------------------------------------------------

Now we issue a command >isoeff KIE.inp and in a couple of seconds obtain the output file named KIE.inp.iso. Peep in the very last strings of this file and find the following records:

ISOTOPE EFFECT(300.00) = 1.107572
ISOTOPE EFFECT(320.00) = 1.102131
ISOTOPE EFFECT(340.00) = 1.097315

This is our final result, which we can compare with the experimental data of Ref. 4 (below). As expected, the computed KIE goes to unity with an increase in temperature, but the match with the experiment is very good. 

Overlay of computed KIEs (red dots) at three different temperatures with experimental data from Ref. 4

 

The last stroke: let's now check how much does the highlighted spectrum range contributes to the KIE. The modified file below does exactly that: it neglects the frequencies below 900 and above 2000 cm-1.

-------------- KIE.inp: input file for ISOEFF98 with vibration cutoffs------------------
$IODEF SSINP = GS.log TSINP = TS.log $END

30N/28N for Mo complexes
$ISOEFF LMASS(19)=14,14 HMASS(19)=15,15 
LMATS(1)=14,14 HMATS(1)=15,15 
RMTRIV = .f. cutlss=900.0 cutlts=900.0 cutuss=2000.0 cututs=2000.0
KIE=.t. EIE=.f.
TKELV(1)=300,320,340 $END
-------------------------------------------------------------------------------

And the result

ISOTOPE EFFECT(300.00) = 1.119611
ISOTOPE EFFECT(320.00) = 1.113151
ISOTOPE EFFECT(340.00) = 1.107423

almost perfectly reproduces the KIE obtained using the whole range of vibrations, confirming our speculations that the direct vibrations of central nitrogens contribute the most to the KIE.

References

1. Melander, L. C. S.; Saunders, W. H., Reaction rates of isotopic molecules. Wiley: New York, 1980; xiv, 331 pp.
2. Enzyme Mechanism from Isotope Effects. Cook, P. F., Ed. CRC Press: 1991.
3. Isotope Effects in Chemistry and Biology. Kohen, A.; Limbach, H.-H., Eds. Taylor & Francis: 2006.

4. Laplaza, C.; Johnson, M.; Peters, J.; Odom, A.; Kim, E.; Cummins, C. C.; George, G.; Pickering, I., Dinitrogen Cleavage by Three-Coordinate Molybdenum(III) Complexes: Mechanistic and Structural Data. JACS 1996, 118, 8623-8638.

5. Anisimov, V.; Paneth, P., ISOEFF98. A program for studies of isotope effects using Hessian modifications. Journal of Mathematical Chemistry 1999, 26, 75-86.
6. L.B. Sims, G. Burton and D.E. Lewis, Q.C.P.E. 337 (1985).

7.Cramer, C. J., Essentials of computational chemistry : theories and models. 2nd ed.; Wiley: Chichester, West Sussex, England ; Hoboken, NJ, 2004; 596 pp.
8. Koch, W.; Holthausen, M. C., A chemist's guide to density functional theory. 2nd ed.; Wiley-VCH: Weinheim ; New York, 2001; xiii, 300 pp.

9. Christian, G. J.; Stranger, R.; Yates, B. F., Optimizing small molecule activation and cleavage in three-coordinate M[N(R)Ar]3 complexes. Inorg Chem 2006, 45, (17), 6851-9.
10. Stranger, R.; Yates, B. F., Mixing of electronic states in molybdenum complexes involved in nitrogen activation. Chemical Physics 2006, 324, (1), 202-209.

Contacts

 

Please send your questions to: