You can access this at http://bessie.che.uc.edu/tlb/teach
Density functional theory has gained a great deal of popularity as a tool for quantum chemistry in the last ten years. A clear indication of this was last year's Nobel Prize in Chemistry to John Pople and Walter Kohn. Kohn developed DFT from its very foundations, while Pople made several advances in quantum chemistry and initiated the writing of the Gaussian code which is the most widely used program for quantum chemistry. Solid state physicists have used DFT methods to calculate the band structures of solids ever since Kohn and Sham put DFT into a useful form in 1965. It has taken chemists a rather long time to latch onto this method, primarily because in its initial form, one could not get highly accurate binding energies for molecules. Now in the last ten years, people have come up with better ways to include the effects of electron correlation in DFT in a practical sense, and one can obtain highly accurate binding energies for molecules. This effort has been led by Perdew and Becke. Computationally the reason for the growth of DFT over Hartree-Fock and higher scaling correlated methods is that the method generally scales as N3 (N is the number of atoms or electrons) at worst, and includes correlation to varying degrees of accuracy, which Hartree-Fock doesn't do. The perturbation theory methods (MP2 etc.) scale with higher powers than N4 (the theoretical limit for Hartree-Fock), so DFT really gains something in relation to scaling. The scaling behavior determines how large a system you can study. Note the prefactor to the scaling relation is important too but it is relatively similar for HF and DFT (for the basis set implementations). People now routinely apply DFT methods to calculations on molecules with tens of atoms, and often treat all the electrons explicitly.
Some references that might be helpful are [1], [2], and [3]. Other prominent names of people doing large scale applications are: Siegbahn, Ziegler, Salahub, Andzelm, Baerends, Trickey, Scuseria, Frisch, White, Challacombe, Noodleman, and Case. We also have a manual for the Gaussian code in which you can read briefly about DFT methods. The standard reference on Hartree-Fock based approaches is [4].
So what's the big difference between HF-based and DFT approaches to quantum chemistry? In their practical implementation, they actually look fairly similar. However, they come from very different starting points, and I'd like to first discuss that. Let's start with the HF method. First, we assume the Born-Oppenheimer approximation. Then, we take the full Hamiltonian of a molecule which includes all electron-electron and electron-nucleus interactions (the nucleus-nucleus term can be added easily at the end). Next, we postulate that the wavefunction for all the electrons is a Slater determinant. This is the key assumption of the HF method. Really, the many electron wavefunction is:
| (1) |
where the coordinates for each electron include both the spatial and spin components. You could just add another coordinate like a and b for spin-up and spin-down if you wanted. The Slater determinant introduces from the start the notion of `orbital', which is purely a mathematical operation (although a good one). The Slater determinant is:
| (2) |
Here the subscript on the orbital labels the orbital number, and that on the coordinate labels the particle number. So we have N orbitals and N electrons. This would correspond to the unrestricted HF theory: each electron can fill a different spatial orbital. At this point, I'd like to point out that orbitals and basis functions aren't the same thing! This is a common misconception. The orbitals are the `physical' one electron functions that we postulated from the very beginning, and basis functions are just the mathematical entities (Gaussians, Slaters, etc.) that we use to `build up' the actual orbitals. Let's say the basis functions are a set of functions like zk. Then:
| (3) |
where the coefficient matrix of the expansion cik can be varied to get the best possible answer, meaning lowest possible total energy. So if we used a really good basis (that is lots of functions to build up orbitals), then the result we would get for the orbitals would be independent of the particular set of functions used. The orbitals are the actual physically occupied states under the assumption we can even represent the many electron wavefunction with an approximation so severe as that above.
Now we assume all the orbitals are normalized, and we construct the expectation value of the Hamiltonian H:
| (4) |
The next step is formally to minimize this total energy with respect to variations in one of the orbitals. You can imagine this as varying the values of the f functions (orbitals) in space and only keeping those variations which lower the total energy. That is, you perform a variational calculation. When you represent the orbitals with a linear combination of basis functions, like above, then you get the canonical HF coupled equations[4]. But note you don't have to introduce the linear combination of basis functions to get out the HF equations, which are expressed entirely in terms of the orbitals. The basis functions are simply a tool to solve the HF equations, which are independent of basis. So here are the equations (given for a diatomic molecule, that is two nuclei A and B:
| (5) |
[^F] is the called the Fock operator and the ei's are the energy levels for the orbitals. This is an integro-differential equation for the f orbitals. The first term on the rhs is the one electron part which has the kinetic energy and the interactions with the nuclei. The second term gives the potential (acting on the wavefunction fi) due to the other electrons. This is just the classical electrostatic potential from the electron distribution. The third term is the exchange operator, which has no classical counterpart. It is directly due to the Pauli exclusion principle. We can write a Poisson type equation which when solved will give the potentials for the last two terms which act on the orbitals:
| (6) |
Then the action of the Fock operator on the orbitals can be rewritten:
| (7) |
Notice that in the second term on the rhs, the wavefunction can come out front because its index is i. At any rate, the last term due to exchange is some kind of nonclassical potential resulting from the Pauli principle. To solve for the classical electrostatic potential (second term), all you have to do is to compute the electron density:
| (8) |
and solve the standard Poisson equation once using that density. To get the action of the exchange operator, you have to solve N Poisson-type equations. That shows why that term is trickier. In a standard basis set calculation, say Gaussians, then analytical forms for the integrals may be available and you can assemble the operators that way, rather than numerically. But there are lots of such integrals, and it takes a lot of time and storage to get them. The N4 scaling in a basis set HF calculation comes because you need four index integrals (to assemble the Hamiltonian) for the exchange term. A more numerical approach can reduce that overhead a lot.
The HF method is something called a mean field theory. This means physically the electrons in the HF calculation move in the average potential (with both classical electrostatic and quantum exchange contributions) of the other electrons.
The DFT approach starts from a very different direction. Namely, the most important quantity in the theory is the electron density r(r), and orbitals are only introduced later as a tool to get accurate kinetic energies. In one of the simplest and most beautiful theorems in physics, Hohenberg and Kohn (HK) proved that the electron density determines the external potential from the nuclei. The electron density is measureable by Xray scattering, so this is very important. Other implications of this theorem are that you can obtain the ground state wavefunction from this density, and in addition the expectation value of any operator, including the energy. Thus you can get the full Hamiltonian for the system, in principle. However, of course we don't have the recipe for getting out this energy, we only know in principle it can be done. If you do figure it out, book your flights to Stockholm! They also proved that any density which is not the ground state density gives an energy greater than the ground state energy, so they have a variational principle. These are very deep and powerful theorems, but essentially useless as they stand.
Kohn and Sham (KS) developed a computational approach in 1965 which made the theory useful. The HK theorem says we can calculate the total energy of a system as follows:
| (9) |
where T is the total kinetic energy. The electron-electron term includes the classical electrostatic interaction:
| (10) |
KS cleverly recast the problem by saying:
| (11) |
where Ts is the kinetic energy of a noninteracting set of electrons moving in an external potential.
We can rewrite this as:
| (12) |
The first term is the interaction of the electrons with the nuclei. Now make:
| (13) |
Exc includes both kinetic and potential terms. Then:
| (14) |
The derivative (functional derivative really) of Exc[r] with respect to r gives the important exchange correlation part of an effective potential:
| (15) |
One can write down a one electron Schrödinger type equation in the overall effective potential:
| (16) |
These are the KS equations for the electron density:
| (17) |
Again, the KS orbitals were introduced only as a tool to get out a more accurate kinetic energy. If we could somehow calculate the exact vxc, we would have solved everything, but of course that is not possible. We know there is one, we just don't know what it is. However, we can estimate it quite well by now.
How do we solve these equations? In principle it is quite simple (and a little harder in practise). You simply make an initial guess for the orbtals, construct the effective potential for that set of orbitals, and then solve the eigenvalue equation to get a new set of orbitals, which should be `better' in some sense. This in turn generates a new effective potential and you go round and round until the problem converges. You can look at convergence of the total energy, the individual eigenvalues, or the electron density, to beneath some tolerance.
If we take a look at the HF and KS equations, they are actually very similar. The one particle and classical electrostatic terms are exactly the same. They only differ in the last term. Quite a while ago, Dirac and then Slater came up with a simple expression for vxc which is a local function of the density and is proportional to r1/3. This approach is called the Xa method and is remarkably accurate for such a simple expression. Slater used this for many years in the study of atoms, molecules, and solids. The next level of refinement is the Local Density Approximation (LDA). Here Exc is calculated as:
| (18) |
where exc is the exact value of the energy for a uniform electron gas (gotten from Monte Carlo simulation). This expression includes both exchange and correlation, but it assumes uniformity, which is certainly not true in atoms or molecules. Since the representation is local, it is very easy to evaluate, and adds little to the cost of the algorithm. Typically, this method makes slight errors in bond lengths but overbinds the bonds by about 10-20%.
The great advances in quantum chemistry applications have happened in the last ten years due to improvements in vxc. This has occurred due to more accurate estimates of the exchange term and through the use of gradient expansions of the energy in terms of derivatives of the electron density. I won't go into details of this except to say that Perdew and Becke are the big players in this and through a lot of hard work they have been able to generate vxc's which can yield chemical accuracies. The modern version of vxc which are used are termed GGA, B-LYP, B3-LYP, etc. Plain old LDA is called S-VWN. (All these acronyms are after people's names.) I won't go into a critical analysis of all these potentials, but each one has its pros and cons. The B3-LYP seems to be the method of choice these days. This field is still under development. The main computational difference between HF and DFT is that there are no analytic integrals for the vxc terms, so the resulting integrals must be evaluated numerically on grids in space. This requires some extra overhead, but it is not horribly costly. The scaling improvement over HF appears to balance that extra cost. For equal cost or less, DFT far outperforms HF. The negative side is that it is not systematically refinable as are the perturbation expansions or the coupled cluster method. The Gaussian code includes all the modern versions of vxc.
In my lecture I'll discuss applications of modern DFT calculations to small molecules (mainly first row elements) for tests of accuracy. Then I will discuss calculations on transition metal complexes, where now DFT appears to be the method of choice (and necessity?). In the last few years, people have begun to apply DFT methods to systems of interest in biolgy, for example DNA and interactions of molecules with DNA. Also, studies have appeared which examine spin densities in bio-inorganic complexes. These are very challenging calculations involving up to hundreds of electrons. In about 1985, Car and Parrinello introduced a new method whereby one can solve for the electron density for a configuration of nuclei, then move the nuclei based on the resulting forces, re-solve the electronic structure problem, and so on. This means one can do real time simulations without any `made up' forcefields. This technique has been applied to many problems in chemistry and materials science. Examples are water and ions in water, the proton in water, silicon surfaces, chemical reactions, etc. In the last few years, a lot of work has gone into developing methods which scale linearly with system size. My group has been working on such a method which uses something called a multigrid technique. If time permits I will discuss this method, and my biased opinions as to where these calculations will head in the future.