Applications of Density Functional Theory to Theoretical Organic Chemistry

An overview of applications of density functional theory (KS/DZVP-GGA, demon2k, or B3LYP/6-311G(d,p), Gaussian 09) to a wide range of problems in theoretical organic chemistry with examples is given, namely thermodynamic properties, geometries (bullvalene), charges (glycine cation), dipole moments, electrostatic potential (acetyl chloride and acetamide), spectroscopy (IR/Raman: acetaldehyde, UV/Vis: polyenes, NMR: thujone and EPR: Koelsch’s radical), gas phase acidities and pKa values (substituted benzoic acids), supramolecular chemistry (quinhydrone complex), and reaction pathways. numbers. Therefore it is advisable to work only with differences of DFT energies, e.g., in reaction pathways or in isodesmic transformations (in which the bonds formed have the same type as the bonds broken). If the focus is on enthalpies of formation, the use of semi empirical methods (such as MNDO, AM1, or PM3) might be more appropriate, the more so since they yield the standard enthalpies of formation directly. Thermodynamic properties such as enthalpies, entropies, Gibbs free energies or heat capacities are also available from Gaussian 09 or demon2k calculations, provided that vibrational frequencies have been calculated. For demon2k, there is a useful helper program called THERMO [7].


Introduction
With the improvement of computer technology and the development of efficient computational methods [1], it is nowadays possible to make detailed quantum-chemical computations even for medium-sized organic molecules (such as strychnine or cholesterol), at least in the gas phase, on personal computers. Basically every molecular property can be calculated: energies and thermodynamic data, geometries, charge distributions, vibrations, magnetic resonance parameters, and reaction pathways. The following text will focus on the application of routine computations based on density functional theory for the estimation of molecular properties. It is a kind of review but also a work-in-progress report. Most examples have been computed by the author, and details are available from our website, www.chemie. fu-berlin.de/dfttoc.
In the calculations presented here, the following software was used: demon2k [2] and Gaussian 09 [3] for quantum-chemical calculations, molden [4] for visualization and preparation of input files, an old version of Spartan (spartan5) [5] for the visualization of reaction pathways and babel [6] for some coordinate conversions.
An introduction to computational chemistry focusing on demon2k is given in Ref. [7], and amongst the relevant textbooks dealing with computational organic chemistry the following shall be cited [8][9][10][11].

Energies and thermodynamic properties
The result of an ab initio or DFT calculation is an energy value, usually given in atomic units (hartree, 1 hartree=2625.50 kJ/mol). This "total energy" refers to the formation of the atom or the molecule from atomic nuclei and electrons. It is straightforward to calculate the formation energy of a molecule from the constituting atoms in the gas phase. However, in physical chemistry standard enthalpies of formation (∆ f H o ) are used which refer to the elements in their standard states, e.g., carbon as graphite. Using a set of empirical parameters, an approximate interconversion is possible, at least for the formation energy in the gas phase. (The numerical difference between energy and enthalpy can usually be neglected, but it can also be corrected for, see below).
It turns out that the agreement between calculated (DFT) and experimental values for the enthalpies of formation is rather poor, with typical errors ranging from about 50 to 200 kJ/mol. It stands to reason that the calculated values are relatively small differences of large numbers. Therefore it is advisable to work only with differences of DFT energies, e.g., in reaction pathways or in isodesmic transformations (in which the bonds formed have the same type as the bonds broken). If the focus is on enthalpies of formation, the use of semi empirical methods (such as MNDO, AM1, or PM3) might be more appropriate, the more so since they yield the standard enthalpies of formation directly.
Thermodynamic properties such as enthalpies, entropies, Gibbs free energies or heat capacities are also available from Gaussian 09 or demon2k calculations, provided that vibrational frequencies have been calculated. For demon2k, there is a useful helper program called THERMO [7].

Geometries
The first step in a quantum-chemical calculation on organic molecules is usually geometry optimization. If successful, this will end up in the nearest local energy minimum of the multidimensional potential energy surface. In the case of flexible molecules, this is frequently not the global minimum. Techniques for conformational search usually make use of molecular mechanics or molecular dynamics and are off-topic here. Some molecules exhibit pronouncedly different conformations in the crystalline state, which is amenable to X-ray crystallography, and in the gas phase. A well-known example is biphenyl, which is planar in the solid state but has a torsional angle of about 35° between the phenyl rings in the gas phase. Some molecules even have fluctuating structure. The simplest example is the ammonia molecule with a double-minimum potential; the nitrogen atom oscillates through the plane of the three hydrogen atoms. The most famous example in organic chemistry is bullvalene (Figure 1), which oscillates between 10!/3=1.2 million equivalent structures with different connectivities.

Charges, dipole moments and electrostatic potentials
Since electron density is distributed over the whole molecule, there is no unique way to assign charges to individual atoms. Usually the procedure suggested by Mulliken is adopted. The calculation of dipole moments from the charge distribution is straightforward. Calculated dipole moments (KS/DZVP-GGA, demon2k) are, e.g., 4.44 debye for calicene and 0.96 debye for azulene, in good agreement with data found in the literature. Formal charges used in structural formulas often do not correspond to molecular reality. Thus, in ammonium groups, written with a formal positive charge on nitrogen, the positive partial charges actually reside on the hydrogen atoms, whereas the partial charge on nitrogen is even negative, see, e.g., the glycine cation ( Figure 2).
In the benzene diazonium ion, Ph-N + ≡N, the nitrogen atom with the formal positive charge carries only a partial charge of 0.16, whereas the maximum positive partial charge (0.29) is located at the neighboring carbon atom, and the remainder distributed among the hydrogen atoms.
In Figure 3, the electrostatic potentials of acetyl chloride and of acetamide are compared. The polarization of the carbonyl group in acetamide is more pronounced, in accordance with the picture given by mesomeric structures (i.e., the donation of an electron pair by the amino group).

Vibrational spectra (Infrared and Raman)
Vibrational spectra, i.e., normal modes of vibrations and their infrared (IR) or Raman intensities, can be calculated both with demon2k and with Gaussian 09. For acetaldehyde, the best agreement with experimental data [12] was obtained with the KS/DZVP-GGA method (demon2k). A part of the results is shown in Table 1. The most intense band (exp.: 1727 cm -1 ) is due to the carbonyl valence vibration.

Ionization energies
For ab initio Hartree-Fock calculations, Koopmans theorem is valid: The ionization energy is equal to the negative value of the orbital energy. However, this theorem is not valid for the Kohn-Sham orbitals used in DFT. Here, approximate ionization energies can be calculated as energy difference between the ground state molecule and its radical cation.

Electronic spectra (UV/Vis)
Only for the hydrogen atom or hydrogen-like ionized atoms with a single electron, the electronic excitation energy is equal to the energy difference between the respective orbitals (the HOMO-LUMO gap for the lowest-energy transition). In molecules, the electron-electron correlation cannot be neglected, and energy differences between electronic states have to be calculated. Especially in DFT calculations, the HOMO-LUMO gap is usually a poor approximation for the lowest-energy electronic transition for reasons related to the failure of Koopmans theorem mentioned above. An accepted method is the use of the respective time-dependent method (TD keyword in Gaussian 09).
As an example, results for polyenes (CH 3 -(CH=CH) n -CH 3 ) with n=2-7 conjugated double bonds are shown in Figure 4, experimental data have been taken from Ref. [13]. Here, experimental results lie in between calculations with the HF (Hartree-Fock) and the B3LYP functional, respectively; actually the better fit is obtained with the HF calculations in this case. Regarding a series of dyes, B3LYP gives usually better results than HF, but the agreement with experimental data is not always satisfactory.

NMR spectra
Methods for the calculation of EPR and NMR parameters have been collected in a recently published monograph [14]. Isotropic chemical    shielding values (σ) and shielding tensors can be obtained by means of Gaussian 09 (using the keyword NMR), either with the Hartree-Fock functional (HF) or the B3LYP hybrid method. For comparison with experiment, conversion to chemical shift values (δ in ppm) is necessary (δ=σ 0 -σ). Since tetramethylsilane (TMS) is the reference standard (δ=0) for the nuclei 1 H, 13 C, and 29 Si, a calculation for this molecule will provide the required reference values σ 0 . For the B3LYP/6-311G(d,p) method, σ 0 ( 13 C)=179.7024 ppm and σ 0 ( 1 H)=31.3919 ppm was obtained. Spin-spin couplings can also be calculated with Gaussian 09 (using the keyword NMR=SpinSpin), but the computational time is significantly increased. In this case, the B3LYP method should be used, because the HF method does not yield sensible results (at least for pyridine, e.g.,). Alternatively, chemical shifts can be estimated by means of the HOSE code approach, which uses base values and increments obtained from data bases of experimentally determined shifts. (NMR prediction is implemented, e.g., in ChemDraw Ultra or in MestreNova.) On the average, the error of both methods (prediction and DFT calculation) is similar. However, the quantum-chemical approach makes use of the three-dimensional molecular structure, whereas the prediction method does not take stereochemistry into account. Also, the latter method only works satisfactorily if a similar structure is found in the data base.
As an example, Figure 5 shows the calculated structure of the most stable conformer of α-thujone; note that this molecule has some conformational flexibility due to the rotation of the isopropyl group. (Thujone is found in a number of plants and is an ingredient of absinthe.) The calculated (and, in parentheses, experimental) 13 C and 1 H chemical shifts of this conformer are collected in Figure 6. For 13 C, the agreement is good to satisfactory, whereas the calculated 1 H-NMR shifts are systematically too small by about 0.7 ppm.

EPR/ENDOR spectra (free radicals)
Regarding free radicals, the spin distribution is of interest to organic chemists, i.e., the probability of finding the unpaired electron at a particular position (the spin population, imprecisely called spin density). The experimental methods of choice are electron paramagnetic resonance (abbreviated as EPR or ESR) [15] and ENDOR spectroscopy [16]. The experimental data of interest are the hyperfine coupling constants (HFC) which are related to the spin populations yet not directly proportional, but these details shall not be considered here [17]. A satisfactory calculation of HFC and hyperfine tensors can be achieved by Gaussian 09 using the methods UB3LYP/6-311G(d,p) (sufficient for proton HFC) or UB3LYP/EPR-III (preferable for 13 C and 14 N, though more time consuming).
As an example, Figure 7 shows the calculated structure of Koelsch's radical [18], and in Figure 8 (left), the calculated (UB3LYP/6-311G(d,p), Gaussian 09) and (in parentheses) the experimental proton HFC (in MHz) [19] are given. The right-hand side of Figure 8 shows the calculated spin densities at the carbon atoms. Experimental and calculated HFC agree well for the fluorenyl moieties, but the calculated HFC of the phenyl group are too large (in magnitude), hence the twist angle of the phenyl ring is probably larger than suggested by the calculation. The spin-density distribution can best be understood in terms of a substituted allyl radical, with the largest spin populations (0.371) at the allyl-1,3-positions and a negative spin population (-0.190) at the allyl-2 position.

Gas phase acidities/basicities
The calculation of gas phase acidities or basicities by DFT methods is straightforward, they are given by the energy differences (or rather   the differences in Gibbs free energies, ∆G o ) of the conjugated acids and bases (pK a =-lg K a =∆G o /(2.303 RT)). Gas phase acidities are conventionally reported as enthalpies, ∆H o . In Figure 9, calculated gas phase acidities (∆E, KS/DZVP-GGA, demon2k) are plotted against experimental data [20]. The correlation (y=1.184x -51, r 2 =0.984) is satisfactory, although it does not correspond exactly to theory.
However, the calculation of pK a values for aqueous solutions would have to take into account solvation energies, particularly the large contributions of hydration energies of ions (which are necessarily involved). Disregarding the influence of solvation for the moment, a reasonable correlation can be obtained nevertheless for a series of similar acids such as substituted benzoic acids, see Figure 10. In this case, a satisfactory correlation of pK a vs. ∆E (∆E in kJ/mol, y=0.0161x -19.4, r 2 =0.9694) is obtained, but the slope (0.0161) deviates from the theoretical value (0.17528) by an order of magnitude. Put differently, hydration effects mediate the strong variation of pK a values theoretically expected for the differences of calculated dissociation energies.
It stands to reason that the σ parameters of the Hammett equation could be estimated along the lines shown here (σ=pK 0 -pK with pK 0 =4.19 for unsubstituted benzoic acid).

Supramolecular chemistry and solvent effects
Since memory requirements and computational times for ab initio and DFT calculations increase strongly with molecular size (for ab initio calculations, computational time is proportional to the fourth power (!) of the number of basis functions), they are not well suited for the treatment of larger aggregates. As an example, the quinhydrone complex (p-benzoquinone+hydroquinone 1:1) shall be considered. The crystal structure is available from the Cambridge Structural Database (reference code QUIDON). In the crystal, the molecules are stacked alternately, but not like coins in a roll because of the repelling forces of the π orbitals. Despite claims in several textbooks, there are no hydrogen bonds between quinone and hydroquinone. The complex is a π complex or charge-transfer complex. In a DFT calculation, a dimeric structure taken from the X ray study was used as input. However, it turned out to be unstable towards the dissociation into the components. Solvation effects can be treated either as bulk effects (in Gaussian 09), or a few (!) solvent molecules can be included explicitly.

Reaction mechanisms (Reaction pathways)
A transition structure is a saddle point on the potential energy surface, i.e., the point of highest energy along the reaction pathway from reactant(s) to product(s). It is characterized by the fact that just one vibrational frequency in the normal mode analysis is imaginary, namely that of the vibration along the reaction coordinate. Quantumchemical calculations may find the transition structure (with the transition search option), provided that the starting structure has been chosen close enough. This type of calculation is particularly suitable for the investigation of pericyclic reactions. Reaction pathways can be followed approximately by keeping one geometrical constraint (atomic distance, bond angle, or dihedral angle) constant while optimizing all other internal coordinates.
Several examples with animated pictures can be found on our webserver (S N 2 reaction, Cope rearrangement, elimination of carbon monoxide).

Conclusion
DFT calculations have become a popular tool in theoretical organic chemistry because of their efficiency. Energies and thermodynamic properties can be calculated. The results for enthalpies of formation are frequently not satisfactory, therefore it is advisable to consider only energy or enthalpy differences of similar structures, e.g., in reaction pathways. All kinds of spectroscopic data can be calculated, at least for the gas phase (IR, UV/VIS, NMR, EPR). Regarding NMR, the quality of calculated chemical shifts is comparable to that obtained from chemical shift prediction (based on databases and increment systems), but they   are based on the full three-dimensional structure and not simply on connectivities. Regarding acid or base strengths, the calculation of gas phase acidities is straightforward, whereas pK a values in solution strongly depend on solvation effects which cannot be handled easily. Reaction mechanisms and reaction pathways can be studied. Transition structures are characterized by one imaginary frequency, but it is often difficult to spot them exactly in the calculations.