Electronic Structure Theory
A. Electronic structure theory describes the motions of the electrons and produces energy surfaces
The shapes and geometries of molecules, their electronic, vibrational, and rotational energy levels and wavefunctions, as well as the interactions of these states with electromagnetic fields lie within the realm of structure theory.1. The Underlying Theoretical Basis
In the Born-Oppenheimer model of molecular structure, it is assumed that the electrons move so quickly that they can adjust their motions essentially instantaneously with respect to any movements of the heavier and slower moving atomic nuclei. This assumption motivates us to view the electrons moving in electronic wave functions (orbitals within the simplest and most commonly employed theories) that smoothly "ride" the molecule's atomic framework. These electronic functions are found by solving a Schrödinger equation whose Hamiltonian He contains the kinetic energy Te of the electrons, the Coulomb repulsions among all the molecule's electrons Vee, the Coulomb attractions Ven among the electrons and all of the molecule's nuclei treated with these nuclei held clamped, and the Coulomb repulsions Vnn among all of these nuclei. The electronic wave functions \(y_k\) and energies \(E_k\) that result from solving the electronic Schrödinger equation \(H_e \;y_k \;=\; E_k\; y_k\) thus depend on the locations \(Q_i\) at which the nuclei are sitting. That is, the \(E_k \)and \(y_k \) are parametric functions of the coordinates of the nuclei. These electronic energies' dependence on the positions of the atomic centers cause them to be referred to as electronic energy surfaces such as that depicted below for a diatomic molecule where the energy depends only on one interatomic distance R. For non-linear polyatomic molecules having N atoms, the energy surfaces depend on 3N-6 internal coordinates and thus can be very difficult to visualize. A slice through such a surface (i.e., a plot of the energy as a function of two out of 3N-6 coordinates) is shown below and various features of such a surface are detailed. The Born-Oppenheimer theory of molecular structure is soundly based in that it can be derived from a starting point consisting of a Schrödinger equation describing the kinetic energies of all electrons and of all N nuclei plus the Coulomb potential energies of interaction among all electrons and nuclei. By expanding the wavefunction Y that is an eigenfunction of this full Schrödinger equation in the complete set of functions \(y_k\) and then neglecting all terms that involve derivatives of any \(y_k\) with respect to the nuclear positions \(Q_i\), one can separate variables such that:
1. The electronic wavefunctions and energies must obey
\(H_e \;y_k \;=\; E_k\; y_k\)
2. The nuclear motion (i.e., vibration/rotation) wavefunctions must obey
\( (T_N + E_k ) c_{k,L} = E_{k,L} c_{k,L} \) ,
where TN is the kinetic energy operator for movement of all nuclei. That is, each and every electronic energy state, labeled k, has a set, labeled L, of vibration/rotation energy levels \(E_{k,L}\) and wavefunctions \(c_{k,L}\). Because the Born-Oppenheimer model is obtained from the full Schrödinger equation by making approximations (e.g., neglecting certain terms that are called non-adiabatic coupling terms), it is not exact. Thus, in certain circumstances it becomes necessary to correct the predictions of the Born-Oppenheimer theory (i.e., by including the effects of the neglected non-adiabatic coupling terms using perturbation theory).
For example, when developing a theoretical model to interpret the rate at which electrons are ejected from rotationally/vibrationally hot \(\text{NH}^-\) ions, we had to consider coupling between two states having the same total energy:
1. \(^2\)P NH\(^-\) in its v=1 vibrational level and in a high rotational level (e.g., J >30) prepared by laser excitation of vibrationally "cold" NH\(^-\) in v=0 having high J (due to natural Boltzmann populations); see the figure below, and
2. \(^3\)S NH neutral plus an ejected electron in which the NH is in its v=0 vibrational level (no higher level is energetically accessible) and in various rotational levels (labeled N).
Because NH has an electron affinity of 0.4 eV, the total energies of the above two states can be equal only if the kinetic energy KE carried away by the ejected electron obeys: \[\text{KE} = \text{E}_{vib/rot} (\text{NH}^- (v=1, J)) - \text{E}_{vib/rot} (\text{NH} (v=0, N)) - 0.4 eV\]
In the absence of any non-adiabatic coupling terms, these two isoenergetic states would not be coupled and no electron detachment would occur. It is only by the anion converting some of its vibration/rotation energy and angular momentum into electronic energy that the electron that occupies a bound N2p orbital in \(\text{NH}^-\) can gain enough energy to be ejected. Energies of \(\text{NH}^-\) and of NH pertinent to the autodetachment of v=1, J levels of \(\text{NH}^-\) formed by laser excitation of v=0 J'' \(\text{NH}^-\) .
2. What is Learned from an Electronic Structure Calculation? The knowledge gained via structure theory is great. The electronic energies Ek (Q) allow one to determine (see my book Energetic Principles of Chemical Reactions) the geometries and relative energies of various isomers that a molecule can assume by finding those geometries \(Q_i \) at which the energy surface \(E_k\) has minima \(\frac{\partial\text{E}_k}{\partial\text{Q}_i} = 0\), with all directions having positive curvature (this is monitored by considering the so-called Hessian matrix \(H_{i,j}=\frac{\partial^2\text{E}_k}{\partial_i Q_j}\) if none of its eigenvalues are negative, all directions have positive curvature). Such geometries describe stable isomers, and the energy at each such isomer geometry gives the relative energy of that isomer. Professor Berny Schlegel at Wayne State has been one of the leading figures whose efforts are devoted to using gradient and Hessian information to locate stable structures and transition states. Professor Peter Pulay has done as much as anyone to develop the theory that allows us to compute the gradients and Hessians appropriate to the most commonly used electronic structure methods. His group has also pioneered the development of so-called local correlation methods which focus on using localized orbitals to compute correlation energies in a manner that scales less severely with system size than when delocalized canonical molecular orbitals are employed. There may be other geometries on the \(E_k\) energy surface at which all "slopes" vanish \(\frac{\partial\text{E}_k}{\partial\text{Q}_i} = 0\), but at which not all directions possess positive curvature. If the Hessian matrix has only one negative eigenvalue there is only one direction leading downhill away from the point \(Q_i\) of zero force; all the remaining directions lead uphill from this point. Such a geometry describes that of a transition state, and its energy plays a central role in determining the rates of reactions which pass through this transition state. At any geometry \(Q_i\), the gradient or slope vector having components \(\frac{\text{E}_k}{\text{Q}_i}\) provides the forces (\(F_i = -\frac{\text{E}_k}{\text{Q}_i}\) ) along each of the coordinates \(Q_i \). These forces are used in molecular dynamics simulations (see the following section) which solve the Newton F = ma equations and in molecular mechanics studies which are aimed at locating those geometries where the F vector vanishes (i.e., the stable isomers and transition states discussed above). Also produced in electronic structure simulations are the electronic wavefunctions \(y_k\) and energies \(E_k\) of each of the electronic states. The separation in energies can be used to make predictions about the spectroscopy of the system. The wavefunctions can be used to evaluate properties of the system that depend on the spatial distribution of the electrons in the system. For example, the \(z^-\) component of the dipole moment of a molecule mz can be computed by integrating the probability density for finding an electron at position r multiplied by the z- coordinate of the electron and the electron's charge e: mz = Ú e \(y_k\)* \(y_k\) z dr . The average kinetic energy of an electron can also be computed by carrying out such an average-value integral: Ú \(y_k\)* (-\(\frac{h^2}{2m_e^2}\) ) \(y_k\) dr. The rules for computing the average value of any physical observable are developed and illustrated in popular undergraduate text books on physical chemistry (e.g., Atkins text or Berry, Rice, and Ross) and in graduate level texts (e.g., Levine, McQuarrie, Simons and Nichols).
Not only can electronic wavefunctions tell us about the average values of all physical properties for any particular state (i.e., \(y_k\) above), but they also allow us to tell how a specific "perturbation" (e.g., an electric field in the Stark effect, a magnetic field in the Zeeman effect, light's electromagnetic fields in spectroscopy) can alter the specific state of interest. For example, the perturbation arising from the electric field of a photon interacting with the electrons in a molecule is given within the so-called electric dipole approximation (see, for example, Simons and Nichols, Ch. 14) by: \[H_{pert} =S_j e^2 r_j\cdot \textbf{E} (t)\]
where E is the electric field vector of the light, which depends on time t in an oscillatory manner, and ri is the vector denoting the spatial coordinates of the ith electron. This perturbation, Hpert acting on an electronic state \(y_k\) can induce transitions to other states \(y^{\prime}_k\) with probabilities that are proportional to the square of the integral \[U^\prime y^{\prime}_k* H_{pert} y_k d\textbf{r} \] So, if this integral were to vanish, transitions between \(y_k\) and \(y^{\prime}_k\) would not occur, and would be referred to as "forbidden". Whether such integrals vanish or not often is determined by symmetry. For example, if \(y_k\) were of odd symmetry under a plane of symmetry sv of the molecule, while \(y^{\prime}_k\) were even under sv , then the integral would vanish unless one or more of the three cartesian components of the dot product \(r_j\cdot\)E were odd under sv The general idea is that for the integral to not vanish the direct product of the symmetries of \(y_k\) and of \(y^{\prime}_k\) must match the symmetry of at least one of the symmetry components present in Hpert . Professor Poul Jørgensen has been especially involved in developing such so-called response theories for perturbations that may be time dependent (e.g., as in the interaction of light's electromagnetic radiation).
3. Present Day Challenges in Structure Theory
a. Orbitals form the starting point; what are the orbitals?
The full N-electron Schrödinger equation governing the movement of the electrons in a molecule is \begin{align*} -\bigg[\frac{h^2}{2m_e}S_{i=1} \nabla^2_i - \frac{S_a S_i Z_a e^2}{r_{ia}} + \frac{\text{S}_{i,j}{e^2}}{r_{ij}} \bigg] y = E y \end{align*} In this equation, i and j label the electrons and a labels the nuclei. Even at the time this material was written, this equation had been solved only for the case N=1 (i.e., for H, He+ , Li2+ , Be3+ , etc.). What makes the problem difficult to solve for other cases is the fact that the Coulomb potential e2/rij acting between pairs of electrons depends upon the coordinates of the two electrons ri and rj in a way that does not allow the separations of variables to be used to decompose this single 3N dimensional second-order differential equation into N separate 3-dimensional equations. However, by approximating the full electron-electron Coulomb potential Si,j e2 /rij by a sum of terms, each depending on the coordinates of only one electron Si, V(ri ), one arrives at a Schrödinger equation [-h2 /2me Si=1 i2 - Sa Si Za e2 /ria + Si V(ri)] y = E y which is separable. That is, by assuming that y (r1 , r2 , ... rN ) = f1 (r1 ) f2 (r2) ... fN (rN), and inserting this ansatz into the approximate Schrödinger equation, one obtains N separate Schrödinger equations: \begin{equation*} \big[-\frac{\hbar^2}{2m_e} \nabla_i^2 -\frac{S_a Z_a e^2}{r_{ia}} + V\left(r_i\right)\big] f_i = E_i f_i \end{equation*} one for each of the N so-called orbitals fi whose energies Ei are called orbital energies. It turns out that much of the effort going on in the electronic structure area of theoretical chemistry has to do with how one can find the "best" effective potential V(r); that is, the V(r), which depends only on the coordinates r of one electron, that can best approximate the true pairwise additive Coulomb potential experienced by an electron due to the other electrons. The simplest and most commonly used approximation for V(r) is the so-called Hartree-Fock (HF) potential: V(r) fi (r) = Sj [ Ú|fj (r')|2 e2 /|r-r'| dr' fi (r) - Úfj*(r') fj (r') e2 /|r-r'| dr' fj (r) ]. This potential, when acting on the orbital fi , can be viewed as multiplying fi by a sum of potential energy terms (which is what makes it one-electron additive), each of which consists of two parts: a. An average Coulomb repulsion Ú |fj (r')|2 e2 /|r-r'| dr' between the electron in fi with another electron whose spatial charge distribution is given by the probability of finding this electron at location r' if it resides in orbital fj : |fj (r')|2 . b. A so-called exchange interaction between the electron in fi with the other electron that resides in fj.. The sum shown above runs over all of the orbitals that are occupied in the atom or molecule. For example, in a Carbon atom, the indices i and j run over the two 1s orbitals, the two 2s orbitals and the two 2p orbitals that have electrons in them (say 2px and 2py ) The potential felt by one of the 2s orbitals is obtained by setting fi = 2s, and summing j over j=1s, 1s, 2s, 2s, 2px , 2py . The term Ú |1s(r')|2 e2 /|r-r'| dr' 2s(r) gives the average Coulomb repulsion between an electron in the 2s orbital and one of the two 1s electrons; Ú |2px(r')|2 e2 /|r-r'| dr' 2s(r) gives the average repulsion between the electron in the 2px orbital and an electron in the 2s orbital; and Ú |2s(r')|2 e2 /|r-r'| dr' 2s(r) describes the Coulomb repulsion between one electron in the 2s orbital and the other electron in the 2s orbital. The exchange interactions, which arise because electrons are Fermion particles whose indistinguishability must be accounted for, have analogous interpretations. For example, Ú 1s*(r') 2s(r') e2 /|r-r'| dr' 1s(r) is the exchange interaction between an electron in the 1s orbital and the 2s electron; Ú 2px*(r') 2s(r') e2 /|r-r'| dr' 2px(r) is the exchange interaction between an electron in the 2px orbital and the 2s electron; and Ú 2s*(r') 2s(r') e2 /|r-r'| dr' 2s(r) is the exchange interaction between a 2s orbital and itself (note that this interaction exactly cancels the corresponding Coulomb repulsion Ú |2s(r')|2 e2 /|r-r'| dr' 2s(r), so one electron does not repel itself in the Hartree-Fock model). There are two primary deficiencies with the Hartree-Fock approximation: a. Even if the electrons were perfectly described by a wavefunction y (r1 , r2 , ... rN ) = f1 (r1 ) f2 (r2) ... fN (rN) in which each electron occupied an independent orbital and remained in that orbital for all time, the true interactions among the electrons Si,j e2 /rij are not perfectly represented by the sum of the average interactions. b. The electrons in a real atom or molecule do not exist in regions of space (this is what orbitals describe) for all time; there are times during which they must move away from the regions of space they occupy most of the time in order to avoid collisions with other electrons. For this reason, we say that the motions of the electrons are correlated (i.e., where one electron is at one instant of time depends on where the other electrons are at that same time). Let us consider the implications of each of these two deficiencies. b. The imperfections in the orbital-level picture are substantial To examine the difference between the true Coulomb repulsion between electrons and the Hartree-Fock potential between these same electrons, the figure shown below is useful. In this figure, which pertains to two 1s electrons in a Be atom, the nucleus is at the origin, and one of the electrons is placed at a distance from the nucleus equal to the maximum of the 1s orbital's radial probability density (near 0.13 Å). The radial coordinate of the second is plotted along the abscissa; this second electron is arbitrarily constrained to lie on the line connecting the nucleus and the first electron (along this direction, the inter-electronic interactions are largest). On the ordinate, there are two quantities plotted: (i) the Hartree-Fock (sometimes called the self-consistent field (SCF) potential) Ú |1s(r')|2 e2 /|r-r'| dr', and (ii) the so-called fluctuation potential (F), which is the true coulombic e2/|r-r'| interaction potential minus the SCF potential. As a function of the inter-electron distance, the fluctuation potential decays to zero more rapidly than does the SCF potential. However, the magnitude of F is quite large and remains so over an appreciable range of inter-electron distances. Hence, corrections to the HF-SCF picture are quite large when measured in kcal/mole. For example, the differences DE between the true (state-of-the-art quantum chemical calculation) energies of interaction among the four electrons in Be (a and b denote the spin states of the electrons) and the HF estimates of these interactions are given in the table shown below in eV (1 eV = 23.06 kcal/mole). Orb. Pair 1sa1sb 1sa2sa 1sa2sb 1sb2sa 1sb2sb 2sa2sb DE in eV 1.126 0.022 0.058 0.058 0.022 1.234 These errors inherent to the HF model must be compared to the total (kinetic plus potential) energies for the Be electrons. The average value of the kinetic energy plus the Coulomb attraction to the Be nucleus plus the HF interaction potential for one of the 2s orbitals of Be with the three remaining electrons -15.4 eV; the corresponding value for the 1s orbital is (negative and) of even larger magnitude. The HF average Coulomb interaction between the two 2s orbitals of 1s22s2 Be is 5.95 eV. This data clearly shows that corrections to the HF model represent significant fractions of the inter-electron interaction energies (e.g., 1.234 eV compared to 5.95- 1.234 = 4.72 eV for the two 2s electrons of Be), and that the inter-electron interaction energies, in turn, constitute significant fractions of the total energy of each orbital (e.g., 5.95 -1.234 eV = 4.72 eV out of -15.4 eV for a 2s orbital of Be). The task of describing the electronic states of atoms and molecules from first principles and in a chemically accurate manner (± 1 kcal/mole) is clearly quite formidable. The HF potential takes care of "most" of the interactions among the N electrons (which interact via long-range coulomb forces and whose dynamics requires the application of quantum physics and permutational symmetry). However, the residual fluctuation potential is large enough to cause significant corrections to the HF picture. This, in turn, necessitates the use of more sophisticated and computationally taxing techniques to reach the desired chemical accuracy. c. Going beyond the simplest orbital model is sometimes essential What about the second deficiency of the HF orbital-based model? Electrons in atoms and molecules undergo dynamical motions in which their coulomb repulsions cause them to "avoid" one another at every instant of time, not only in the average-repulsion manner that the mean-field models embody. The inclusion of instantaneous spatial correlations among electrons is necessary to achieve a more accurate description of atomic and molecular electronic structure. Some idea of how large the effects of electron correlation are and how difficult they are to treat using even the most up-to-date quantum chemistry computer codes was given above. Another way to see the problem is offered in the figure shown below. Here we have displayed on the ordinate, for Helium's 1S (1s2) state, the probability of finding an electron whose distance from the He nucleus is 0.13Å (the peak of the 1s orbital's density) and whose angular coordinate relative to that of the other electron is plotted on the absissa. The He nucleus is at the origin and the second electron also has a radial coordinate of 0.13 Å. As the relative angular coordinate varies away from 0 deg, the electrons move apart; near 0 deg, the electrons approach one another. Since both electrons have the same spin in this state, their mutual Coulomb repulsion alone acts to keep them apart. What this graph shows is that, for a highly accurate wavefunction (one constructed using so-called Hylleraas functions that depend explicitly on the coordinates of the two electrons as well as on their interparticle distance coordinate) that is not of the simple orbital product type, one finds a "cusp" in the probability density for finding one electron in the neighborhood of another electron with the same spin. The probability plot for the Hylleraas function is the lower dark line in the above figure. In contrast, this same probability density, when evaluated for an orbital-product wavefunction (e.g., for the Hartree-Fock function) has no such cusp because the probability density for finding one electron at r, q, f is independent of where the other electron is (due to the product nature of the wavefunction). The Hartree-Fock probability, which is not even displayed above, would thus, if plotted, be flat as a function of the angle shown above. Finally, the graph shown above that lies above the Hylleraas plot and that has no "sharp" cusp was extracted from a configuration interaction wavefunction for He obtained using a rather large correlation consistent polarized valence quadruple atomic basis set. Even for such a sophisticated wavefunction (of the type used in many state of the art ab initio calculations), the cusp in the relative probability distribution is clearly not well represented. d. For realistic accuracy, improvements to the orbital picture are required Although highly accurate methods do exist for handling the correlated motions of electrons (e.g., the Hylleraas method mentioned above), they have not proven to be sufficiently computationally practical to be of use on atoms and molecules containing more than a few electrons. Hence, it is common to find other methods employed in most chemical studies in which so-called correlated wavefunctions are used. By far, the most common and widely used class of such wavefunctions involve using linear combinations of orbital product functions (actually, one must use so-called antisymmetrized orbital products to properly account for the fact that Fermion wavefunctions such as those describing electrons are odd under permutations of the electrons' labels): Y = S J CJ| fJ1 fJ2 fJ3 ...fJ(N-!) fJN |, with the indices J1, J2, ..., JN labeling the spin-orbitals and the coefficients CJ telling how much of each particular orbital product to include. As an example, one could use Y = C1 |1sa 1sb | - C2 [|2pza2pzb | - |2pxa2pxb | -2pya2pyb |] as a wavefunction for the 1S state of He (the last three orbital products are combined to produce a state that is spherically symmetric and thus has L = 0 electronic angular momentum just as the |1sa1sb| state does). Using a little algebra, and employing the fact that the orbital products |f1 f2 | = (2)-1/2 [ f1 f2 - f2 f1 ] are really antisymmetric products, one can show that the above He wavefunction can be rewritten as follows: Y = C1/3 {|fz a f'z b| - |fx a f'xb| - |fya f'y b| }, where fz = 1s + (3C2/C1)1/2 2pz and f'z = 1s - (3C2/C1)1/2 2pz , with analogous definitions for fx , f'x , fy , and f'y . The physical interpretation of the three terms ({|fz a f'z b| , |fx a f'x b| , and |fy a f'y b| ) is that |fz a f'z b| describes a contribution to Y in which one electron of a spin resides in a region of space described by fz while the other electron of b spin is in a region of space described by f'z , and analogously for |fxa f'x b| and |fy a f'y b|. Such a wavefunction thus allows the two electrons to occupy different regions of space since each orbital f in a pair is different from its partner f'. The extent to which the orbital differ depends on the C2/C1 ratio which, in turn, is governed by how strong the mutual repulsions between the two electrons are. Such a pair of so-called polarized orbitals is shown in the figure below. e. Density Functional Theory (DFT) These approaches provide alternatives to the conventional tools of quantum chemistry. The CI, MCSCF, MPPT/MBPT, and CC methods move beyond the single-configuration picture by adding to the wave function more configurations whose amplitudes they each determine in their own way. This can lead to a very large number of CSFs in the correlated wave function, and, as a result, a need for extraordinary computer resources. The density functional approaches are different. Here one solves a set of orbital-level equations \[\bigg[ -\frac{h^2}{2m_e} - \frac{S_A Z_A e^2}{\big|r-R_A\big|} + \frac{\tilde{U^\prime}\left(r^\prime\right)e^2}{\big|r-r^\prime\big|} dr^\prime + U\left(r\right) \bigg] f_i = e_i f_i\] in which the orbitals {fi} 'feel' potentials due to the nuclear centers (having charges ZA), Coulombic interaction with the total electron density r(r'), and a so-called exchange-correlation potential denoted U(r'). The particular electronic state for which the calculation is being performed is specified by forming a corresponding density r(r'). Before going further in describing how DFT calculations are carried out, let us examine the origins underlying this theory. The so-called Hohenberg-Kohn theorem states that the ground-state electron density r(r) describing an N-electron system uniquely determines the potential V(r) in the Hamiltonian H = Sj {-h2/2mej2 + V(rj) + e2/2 Skj 1/rj,k }, and, because H determines the ground-state energy and wave function of the system, the ground-state density r(r) determines the ground-state properties of the system. The proof of this theorem proceeds as follows:
a. r(r) determines N because Ú r(r) d3r = N.
b. Assume that there are two distinct potentials (aside from an additive constant that simply shifts the zero of total energy) V(r) and V'(r) which, when used in H and H', respectively, to solve for a ground state produce E0, Y (r) and E0', Y'(r) that have the same one-electron density: Ú |Y|2 dr2 dr3 ... drN = r(r)= Ú |Y'|2 dr2 dr3 ... drN .
c. If we think of Y' as trial variational wave function for the Hamiltonian H, we know that \begin{equation*} E_0 \left
d. Similarly, taking Y as a trial function for the H' Hamiltonian, one finds that E0' < E0 + Ú r(r) [V'(r) - V(r)] d3r. e. Adding the equations in c and d gives E0 + E0' < E0 + E0', a clear contradiction. Hence, there cannot be two distinct potentials V and V' that give the same ground-state r(r). So, the ground-state density r(r) uniquely determines N and V, and thus H, and therefore Y and E0. Furthermore, because Y determines all properties of the ground state, then r(r), in principle, determines all such properties. This means that even the kinetic energy and the electron-electron interaction energy of the ground-state are determined by r(r). It is easy to see that Ú r(r) V(r) d3r = V[r] gives the average value of the electron-nuclear (plus any additional one-electron additive potential) interaction in terms of the ground-state density r(r), but how are the kinetic energy T[r] and the electron-electron interaction Vee[r] energy expressed in terms of r? The main difficulty with DFT is that the Hohenberg-Kohn theorem shows that the ground-state values of T, Vee , V, etc. are all unique functionals of the ground-state r (i.e., that they can, in principle, be determined once r is given), but it does not tell us what these functional relations are. To see how it might make sense that a property such as the kinetic energy, whose operator -h2 /2me 2 involves derivatives, can be related to the electron density, consider a simple system of N non-interacting electrons moving in a three-dimensional cubic "box" potential. The energy states of such electrons are known to be E = (h2/2meL2) (nx2 + ny2 +nz2 ), where L is the length of the box along the three axes, and nx , ny , and nz are the quantum numbers describing the state. We can view nx2 + ny2 +nz2 = R2 as defining the squared radius of a sphere in three dimensions, and we realize that the density of quantum states in this space is one state per unit volume in the nx , ny , nz space. Because nx , ny , and nz must be positive integers, the volume covering all states with energy less than or equal to a specified energy E = (h2/2meL2) R2 is 1/8 the volume of the sphere of radius R: F(E) = 1/8 (4p/3) R3 = (p/6) (8meL2E/h2)3/2 . Since there is one state per unit of such volume, F(E) is also the number of states with energy less than or equal to E, and is called the integrated density of states. The number of states g(E) dE with energy between E and E+dE, the density of states, is the derivative of F: g(E) = dF/dE = (p/4) (8meL2/h2)3/2 E1/2 . If we calculate the total energy for N electrons, with the states having energies up to the so-called Fermi energy (i.e., the energy of the highest occupied molecular orbital HOMO) doubly occupied, we obtain the ground-state energy: = (8p/5) (2me/h2)3/2 L3 EF5/2. The total number of electrons N can be expressed as N = 2 g(e)dE= (8p/3) (2me/h2)3/2 L3 EF3/2, which can be solved for EF in terms of N to then express E0 in terms of N instead of EF: E0 = (3h2/10me) (3/8p)2/3 L3 (N/L3)5/3 . This gives the total energy, which is also the kinetic energy in this case because the potential energy is zero within the "box", in terms of the electron density r (x,y,z) = (N/L3). It therefore may be plausible to express kinetic energies in terms of electron densities r(r), but it is by no means clear how to do so for "real" atoms and molecules with electron-nuclear and electron-electron interactions operative. In one of the earliest DFT models, the Thomas-Fermi theory, the kinetic energy of an atom or molecule is approximated using the above kind of treatment on a "local" level. That is, for each volume element in r space, one assumes the expression given above to be valid, and then one integrates over all r to compute the total kinetic energy: TTF[r] = Ú (3h2/10me) (3/8p)2/3 [r(r)]5/3 d3r = CF Ú [r(r)]5/3 d3r , where the last equality simply defines the CF constant (which is 2.8712 in atomic units). Ignoring the correlation and exchange contributions to the total energy, this T is combined with the electron-nuclear V and Coulombic electron-electron potential energies to give the Thomas-Fermi total energy: E0,TF [r] = CF Ú [r(r)]5/3 d3r + Ú V(r) r(r) d3r + e2/2 Ú r(r) r(r')/|r-r'| d3r d3r', This expression is an example of how E0 is given as a local density functional approximation (LDA). The term local means that the energy is given as a functional (i.e., a function of r) which depends only on r(r) at points in space but not on r(r) at more than one point in space. Unfortunately, the Thomas-Fermi energy functional does not produce results that are of sufficiently high accuracy to be of great use in chemistry. What is missing in this theory are a. the exchange energy and b. the correlation energy; moreover, the kinetic energy is treated only in the approximate manner described. In the book by Parr and Yang, it is shown how Dirac was able to address the exchange energy for the 'uniform electron gas' (N Coulomb interacting electrons moving in a uniform positive background charge whose magnitude balances the charge of the N electrons). If the exact expression for the exchange energy of the uniform electron gas is applied on a local level, one obtains the commonly used Dirac local density approximation to the exchange energy: Eex,Dirac[r] = - Cx Ú [r(r)]4/3 d3r, with Cx = (3/4) (3/p)1/3 = 0.7386 in atomic units. Adding this exchange energy to the Thomas-Fermi total energy E0,TF [r] gives the so-called Thomas-Fermi-Dirac (TFD) energy functional. Because electron densities vary rather strongly spatially near the nuclei, corrections to the above approximations to T[r] and Eex.Dirac are needed. One of the more commonly used so-called gradient-corrected approximations is that invented by Becke, and referred to as the Becke88 exchange functional: Eex(Becke88) = Eex,Dirac[r] -g Ú x2 r4/3 (1+6 g x si\(\text{NH}^-\)1(x))-1 dr, where x =r-4/3 |r|, and g is a parameter chosen so that the above exchange energy can best reproduce the known exchange energies of specific electronic states of the inert gas atoms (Becke finds g to equal 0.0042). A common gradient correction to the earlier T[r] is called the Weizsacker correction and is given by dTWeizsacker = (1/72)(h/me) Ú |r(r)|2/r(r) dr. Although the above discussion suggests how one might compute the ground-state energy once the ground-state density r(r) is given, one still needs to know how to obtain r. Kohn and Sham (KS) introduced a set of so-called KS orbitals obeying the following equation: {-1/22 + V(r) + e2/2 Ú r(r')/|r-r'| dr' + Uxc(r) }fj = ej fj , where the so-called exchange-correlation potential Uxc (r) = dExc[r]/dr(r) could be obtained by functional differentiation if the exchange-correlation energy functional Exc[r] were known. KS also showed that the KS orbitals {fj} could be used to compute the density r by simply adding up the orbital densities multiplied by orbital occupancies nj : r(r) = Sj nj |fj(r)|2. (here nj =0,1, or 2 is the occupation number of the orbital fj in the state being studied) and that the kinetic energy should be calculated as T = Sj nj