Dehydration of substituted alcohols produces a mixture of isomeric alkenes. The product ratios in such elimination reaction are determined largely by the thermodynamic stability (heat of formation) of the carbocation intermediates and the final alkenes. The kinetic barriers (transition states) play a secondary role in determining the observed product ratios. The following tutorial illustrates how quantum chemical methods can be used to predict or rationalize outcomes of chemical reactions when the product ratios are determined by their thermodynamic stability.
Draw a reaction mechanism that would rationalize the formation of 1-methylcyclohexene, 3-methylcyclohexene, and methylenecyclohexane when 2-methylcyclohexanol is dehydrated in the presence of a strong acid.
Notice that the reactant can exist either as cis or trans 2-methylcyclohexanol. Furthermore, each of these two stereoisomers can exist in multiple conformation. For example, cis-2-methylcyclohexanol can exist in a form where the methyl group is axial and hydroxyl group is equatorial, or in a conformation where the methyl group is equatorial and the hydroxyl group is axial. Furthermore, the rotation of the hydroxyl hydrogen can generate different conformations. It is thus interesting to know what is the population of different conformers and if different conformers show different reactivities. For simplicity, consider two structures that differ in equatorial/axial placement of the two substituents. Prepare two structures of the cis-2-methylcyclohexanol using a molecular editor, such as MOLDEN. Save each input geometry in the XYZ format as a text file. (PROVIDE MOL2 FILES FOR TUTORIAL) Edit the molecular geometry such that each line contains the atom name, the nuclear charge (either as integer or floating point number) and the XYZ coordinates as three floating point numbers.
Next, the geometries of the two conformers must be optimized. Reasonable geometries of stable molecules can be obtained at the HF/6-31+G(d,p) level. Even though this approach suffers both from lack of description of electron correlation and limited basis set, the errors from these two sources cancel to a some degree at this level of theory.
Add appropriate keywords to perform HF/6-31+G(d,p) geometry optimization.
The input file for one of the conformers should look like this:
$contrl scftyp=rhf runtyp=optimize coord=unique nprint=-5 $end $system mwords=64 timlim=600 $end $basis gbasis=N31 ngauss=6 ndfunc=1 npfunc=1 diffsp=1 $end $statpt method=qa hess=calc hssend=.t. $end $data cis-2-methylcyclohexanol: axial Me; equatorial OH C1 C 6 -1.898291 0.276754 0.174432 H 1 -2.956336 0.394090 -0.129446 C 6 -1.004520 1.100869 -0.740572 H 1 -1.837561 0.667448 1.210409 C 6 0.484558 0.913653 -0.437693 H 1 -1.199372 0.824161 -1.796339 H 1 -1.270257 2.172824 -0.655905 C 6 0.813410 -0.593458 -0.527990 H 1 1.069068 1.452782 -1.222637 C 6 0.861269 1.494997 0.915978 C 6 -0.040886 -1.372410 0.485727 H 1 0.586330 -0.953511 -1.563125 O 8 2.201346 -0.750943 -0.303561 C 6 -1.513650 -1.194910 0.156713 H 1 0.227965 -2.445671 0.491971 H 1 0.170271 -0.985402 1.508367 H 1 -1.735349 -1.633550 -0.836865 H 1 -2.134064 -1.757705 0.880235 H 1 2.372905 -0.406381 0.562542 H 1 1.949140 1.597116 1.017876 H 1 0.417879 2.487376 1.066253 H 1 0.517283 0.840998 1.735584 $end
Submit your calculation. Under Linux, specify the path to pcgamess libraries with the -ex flag
pcgamess -i CycHex_axMe_eqOH.inp -o CycHex_axMe_eqOH.out -ex /usr/local/pcgamess &
PC GAMESS and US GAMESS offer several methods for geometry optimization. The default method, quadratic approximation (method=qa), starts with a guess hessian and updates the hessian during the optimization. This approach tends to converge slowly for when larger molecules are optimized in Cartesian coordinates, partially because the initial guess hessian (1/3 of the unit matrix) is too different from the true hessian. The evaluation of the accurate hessian can be requested with hess=calc in case of a new jobs or hess=read when the hessian was obtained as a part of a previous calculation. Even though the full evaluation of second derivative matrix can be expensive, the decrease in optimization steps may well pay off at the end. For example, default minimization of cis-2-methylcyclohexanol would lead to abortion of the calculation because the default number of steps (nstep=20) is insufficient to achieve the convergence. In this case, the initial calculation of hessian allows successful convergence in less than 15 steps.
The XYZ input used here is often not the best choice for geometry optimization. Optimization in correctly chosen internal coordinates (i.e. bond lengths, bond angles, and dihedrals) is typically much faster but the choice of correct internal coordinates is not trivial. The document http://phoenix.liu.edu/~nmatsuna/gamess/refs/howto.geom.html further describes how to best optimize molecular structures with GAMESS.
Quantum chemical calculations of molecules are necessarily approximate and thus inaccurate. In general, more accurate methods are more time-consuming. This makes highly accurate methods ill-suited for geometry optimizations. Furthermore, efficient geometry optimization requires analytical energy gradients, which may not be possible with more accurate methods. For example, PC GAMESS only supports analytical energy gradients at the Hartree-Fock level. Thus, geometry optimizations are often performed with modest basis sets ignoring dynamic electron correlation.
It is common to obtain a better estimate for the energy of the optimized structure by performing a single point calculation with a more accurate method using the structure optimized at the lower level. The single point calculation typically will include some treatment of electron correlation and may employ a larger basis. We will next perform MP3/6-31+G(d,p) calculation on the HF/6-31+G(d,p) optimized geometry to get a more accurate enery value for this structure.
PC GAMESS provides some facilities to speed up single-point energy calculations calculations after geometry optimization. Each PC GAMESS calculation will create a file named PUNCH which will include information about the final optimized geometry, molecular orbitals, and if the hession calculation was requested, the hessian and vibrational modes. A single point calculation can be started using the final geometry and molecular orbital coefficient from this file. To reuse the molecular orbital coefficients, the basis set of the optimization job and single point energy must be the same.
Similarly, a new higher-level optimization can be started from the final geometry and saved hessian matrix.
A new input file can be generated by copyig the converged geometry, molecular orbitals, and hessian from the existing PUNCH file. The process can be automated, for example a small Python script parse_punch.py will create relevant parts of the new input file based on the existing PUNCH file. The beginning of a input file that requests single point MP3/6-31+G(d,p) calculation based on stored MO's is shown below. Notice how the guess group keywords now specify that existing 222 molecular orbitals must be used to form the initial guess.
$contrl scftyp=rhf runtyp=energy mplevl=3 coord=unique nprint=-5 $end $system mwords=128 timlim=600 $end $basis gbasis=N31 ngauss=6 ndfunc=1 npfunc=1 diffsp=1 $end $guess guess=moread norb=222 $end $mp3 stpair=.t. $end $data cis-2-methylcyclohexanol: eq Me; ax OH: HF/6-31+G(d,p) minimum C1 0 C 6.0 1.4251583629 1.3393202875 0.0754414697 H 1.0 1.9809199407 1.9999049676 0.7360591014 C 6.0 -0.0711553425 1.4104605329 0.3947179793 H 1.0 1.5912694996 1.7018933046 -0.9383277786 C 6.0 -0.8929432610 0.4448332521 -0.4722016423 H 1.0 -0.2320285909 1.1868235236 1.4501351932 H 1.0 -0.4342271790 2.4246633509 0.2485467849 C 6.0 -0.3637471018 -0.9946740474 -0.3593930999 C 6.0 -2.3889442442 0.5341008250 -0.1709573484 H 1.0 -0.7433180984 0.7379242246 -1.5124187707 C 6.0 1.1341103924 -1.0605531021 -0.6684538429 O 8.0 -0.6586548737 -1.5786486999 0.8918324051 H 1.0 -0.8965210547 -1.6148631193 -1.0715983621 C 6.0 1.9557386823 -0.0921122570 0.1893658366 H 1.0 1.4777908454 -2.0814466141 -0.5302858492 H 1.0 1.2792849805 -0.8134003706 -1.7188263426 H 1.0 1.9350780463 -0.4056773949 1.2333300587 H 1.0 2.9994764596 -0.1326502060 -0.1099163488 H 1.0 -2.9573630890 -0.1329660813 -0.8130764416 H 1.0 -2.6011241635 0.2598398106 0.8563711637 H 1.0 -2.7533264086 1.5445317324 -0.3352890388 H 1.0 -0.1795548024 -1.1582519190 1.5858668733 $end $VEC 1 1-1.15049003E-05-6.75879423E-05-7.50616237E-06-1.12241772E-05 4.65277592E-06 1 2 6.08486041E-04-9.68350486E-05-1.75422051E-04 2.85571695E-05 5.65572418E-03
Prepare structures of the two carbocation intermediates using the Z-Matrix Editor in the program MOLDEN. The structure of the secondary carbocation can be obtained from the reactant by removing the axial hydroxyl group. The structure of the tertiary carbocation can be obtained from the secondary carbocation by deleting a hydrogen atom at the tertiary carbon and adding a hydrogen to the carbon that was the cationic site in the secondary carbocation. Your optimization will finish faster if you alter the bond distances and angles in the Z-matrix such that they are close to the actual structure of the secondary and tertiary carbocation (e.g. planar trigonal geometry at the cationic center). Save the structures as XYZ files naming them CycHex_sec_cation.inp and CycHex_ter_cation.inp, for example.
Prepare structures of the three alkenes that form in this reaction ...
(MP2/cc-pVTZ takes two hours)
This could look really cool with new Mathematica 6
PC GAMESS reports energy values in the Hartree units (1 Hartree = 627.51 kcal/mol). You can use the Unix command grep to display relevant lines: grep 'ENERGY= ' filename. Alternatively, you can use MOLDEN to visualize the structure and the progress of the optimization; holding your mouse over the point will display the Hartree-Fock SCF energy.
TO BE CONTINUED ...