This is a computational portfolio example for a reader well acquainted with the theory of kinetic isotope effects (KIE; Refs. 1-3). He/she is well aware of the possibility to calculate these effects from first principles, and knows exactly why he/she needs that calculation. Below we provide a scheme and a case study to get started with this type of calculations. 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 is as following (the license for the software used differs depending on the application):
1. Molecular modeling package to draw structures and get rough 3D-placement of atoms.
2. Open Babel – a package to convert between plethora of molecule formats.
3. Packages to minimize the structures and compute vibration frequencies: GAMESS, Mopac2016, Gaussian.
4. Packages to visualize the calculated vibrations: Molekel (the interface to show vibrations can save the movements as avi or image sequence files); Molden (very convenient visualization of vibrations; displays Raman and IR spectra); ChemCraft.
5. ISOEFF98 (Ref.5) – a successor of BEBOVIB (Ref. 6) to calculate kinetic and equilibrium isotope effects. Available 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), we need 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 for 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 this 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 (back in 2008), which was a substantial gain.
So, the structures of GS and TS should look like this:
Now, if we want to save some more computing time, it is advisable to adjust the interatomic distances and valence angles according to the available structural data (if there is any). In this particular case Ref. 4 provides plethora of structural information both on GS and TS, which we can use. The step above is optional, but, in some instances, can be crucial for the success of conformation search.
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.
Preparation of the input files for GAMESS, Mopac, or Gaussian and running the calculation
The crucial point is to settle on which particular program you are going to use. We are limited to these three packages, because the current version of ISOEFF can read output only from these programs. The example inputs will be prepared for Gaussian 98.
First, we create draft input files using Babel package, 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
-----------------------------------------------------------------
A wrong choice of keywords will result in wasted computing time, abnormal termination of the calculation, and other misbehavior.
The structures above contain transition metals. Therefore, the only option is to employ the density functional theory (DFT) for calculation (and we are limited to choice between GAMESS and Gaussian in this case, because Mopac does not support DFT). Ref. 7 and 8 are good introductory texts for DFT theory. For usual organic molecules one could use more common ab initio or semi-empirical methods available in Mopac too.
For the input file for the GS below, the following keywords are required:
| %mem=100MW | Dynamic memory allocated for Gaussian is 100 megawords |
| %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. |
| opt(maxcycle=500) | Perform the structure optimization for maximum of 5000 cycles. |
| scf(maxcycle=2000) | Force the SCF procedure to converge at maximum of 2000 cycles. |
| 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. |
| 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 (see Ref.9 and Ref. 10 for more thorough DFT study of that reaction). |
The charge of molecule is obviously 0, and the multiplicity normally would be a singlet (1). But this system (see Ref. 4,9,10) in 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. It 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.
The transition state is trickier, and below are comments on special keywords 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; 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 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 possible isotopes using these two output log files.
Analysis of vibrational frequencies
The aim of the calculation is 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 KIE is experimentally measured (Ref. 4), so we can use it to verify the calculation scheme.
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.
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 2000 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.
Calculation of the KIE
This step is performed by ISOEFF98 program, for which a separate input file is needed together with both Gaussian output files in the same directory. For general instructions it’s better to consult with ISOEFF manual. Below the example of the input file is provided, for calculation of KIE using all possible vibration frequencies (but neglects trivial ones, namely, translations and rotations).
| $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 | A comment |
| 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. Carefully check your files for the atom numbering. When the transition state was built, the central nitrogens accidentally got the numbers 1 and 2 instead of 19 and 20 in GS (thanks that ISOEFF is 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 issue a command
>isoeff KIE.inp
and obtain the output file named KIE.inp.iso. The very last strings of this file have 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 anticipated, the computed KIE decreases with an increase in temperature, and the match with the experiment is very good:
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 assumption 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.