9.8: Numerical QM Methods

Introduction

Quantum mechanical computations can be performed to estimate many properties of molecules: shape, charge distribution, dipole moments, spectroscopic properties, ionization energies, reactivity and much more. Realistic computations require significant computing power combined with sophisticated numerical methods that have been undergoing continuous development since the 1960s. The current state of the art ab initio (from first principles) methods come in two flavors Hartree-Fock self consistent field methods (HF-SCF) and Density functional theory methods (DFT). These methods are ab initio in the sense that they do not require any empirical (measured values) other than basic physical constants. This contrasts with semi-empirical methods, many of which are extensions of the Hückel method beyond \(\pi\) orbitals. In semi-empirical methods many of the necessary integrals are replaced with constants that have been adjusted to get results (e.g. molecular geometries) that match those of a particular set of molecules. Some of these semi-empirical methods (e.g. the extended Hückel) are used to make initial guesses at the electron distribution to begin an ab initio computation.

Both methods utilize the Born-Oppenheimer approximation to separate the nuclear and electronic parts of the calculation. They also depend on the variational principle, the provable assertion that the best guess at the ground state wavefunction is the lowest energy that can be achieved by adjusting parameters in a guessed wavefunction. Thus an optimized wavefunction provides an upper bound to the energy of the actual wavefunction.

Although it is now possible for anybody who can draw a Lewis structure to purchase or download software to do these types of computations, there are myriad pitfalls and ways to get results that are very far from what is observed experimentally. Care must be taken when doing these types of calculations and should generally be undertaken and interpreted with the assistance of an expert in quantum computations.

Hartree-Fock self consistent field methods (HF-SCF)¹

This method is based on the molecular orbital approximation, where each electron is assumed to be represented by a single electron wavefunction consisting of a linear combination of atomic orbitals (LCAO). So a single electron wavefunction is of the form:

\[\phi_j = c_{j1}\psi_1 + c_{j2}\psi_2 + c_{j3}\psi_3...\label{LCAO}\]

where \(\phi_j\) is the single electron wavefunction, the \(\psi_k\) are the atomic orbitals on all the atoms and the c_jk are the coefficients that scale the contribution of each of the atomic orbitals to the overall wavefunction. To make the computations numerically tractable, the \(\psi_k\) are built from sums of 3-D gaussians rather than using the functional forms of the s, p, d, etc. from the exact solution of the one electron atom. The overall wavefunction is then taken as a product of these wavefunctions in the same way as for a multi-electron atom.

\[\Psi_{molec} = \phi_1 \phi_2 \phi_3 \phi_4 ...\label{psimolec}\]

The Hartree-Fock method accounts for the interaction of the electrons in the different molecular orbitals, \(\phi_j\), by adding a potential term to the Hamiltonian that is the average electrostatic potential experienced by an electron in orbital j because of the probability distribution of the all the other electrons. This is then used to optimize the coefficients in equation \(\ref{LCAO}\). In order to start a calculation you must have an initial guess for the orbitals and their occupations. Most HF-SCF use an extended version of the Hückel method for this initial guess. Glossing over the details of setting up the complex secular equations for these large systems, the HF-SCF procedure can be outlined as follows:

1. Guess an initial geometry.
2. Choose a basis set (the \(\psi_q\) in equation \(\ref{LCAO}\)).
3. Use a method such as extended Hückel to generate an initial set of orbital energies and occupations.
4. Start the SCF process:
   1. Choose a wavefunction \(\phi_k\) to focus on.
   2. Calculate the average electrostatic potential experienced by an electron in orbital \(\phi_k\) using the current wavefunctions for all the other occupied orbitals and the nuclear positions.
   3. Adjust \(\phi_k\) to minimize its energy in the electrostatic field.
   4. Choose a new \(\phi_k\) to concentrate on and repeat (b) and (c) until the total energy has reached a minimum within the precision requested.
   5. This process may iterate through all the \(\phi\) multiple times.
5. Move the nuclei and restart SCF (step 4) until the total energy is minimized within the precision requested.

The SCF part of this process keeps adjusting the wavefunctions until they stop changing. This generally works, but there is no guarantee that the SCF process will converge to an optimal set of \(\phi_k\). In large systems and with large basis sets the SCF process can start oscillating. Significant effort has gone into developing techniques to improve convergence (a couple of common techniques are Damping and DIIS). Another issue with this method (and any simple MO calculation) is that it ignores electron correlation to the extent that the average electrostatic potential does not account for the fact that any chosen pair of electrons repels each other, so will have some correlation of their movement. The result is that HF-SCF energies are expected to be higher than the true energy, where electron correlation is taken into account. Extensions to the HF-SCF that include variable contributions from alternative (excited) electron configurations can correct for this but are computationally very expensive. There are many versions of these configuration interaction extensions (e.g. CI, MCSCF, CASSCF, CC). Some of these configuration interaction methods come in specialized flavors designed to calculate a specific molecular property such as ionization energies (IP-EOMCCSD), that are sensitive to electron correlations.

Density functional theory methods (DFT)²

The basic idea behind density functional theory methods is to calculate the electron density at a grid of points in the system rather than the quantum states (orbitals) of each of the electrons. Rather than having to account for 3 coordinates for each electron, you only have to account for the electron density at each point in 3-D space. It has been proven that describing the electron density is adequate to determine the overall system energy and in the exact form of DFT the variational principle applies (for a chosen basis set the minimum energy is the best estimate of the actual value). Because only the overall electron density is calculated rather than individual wavefunctions, this method is generally much less computationally intensive.

DFT methods were originally developed for solid state computations where the number of nuclei and electrons are very large. All practical implementations require an empirical function (the "functionals") that is guessed at to perform the calculation. It was not until the early 1990s that people began to guess at functionals able to yield reasonable agreement with properties of small molecules. Since then a large amount of work has been directed at improving these guessed functionals. In theory these functionals have to account for electron correlations and the non-classical behavior of electrons, but as of yet only guesses as to the form of these functionals are available.

The outline of the procedure for DFT calculations is similar to HF-SCF except that the SCF process optimizes the electron density rather than each orbital.

DFT is rapidly evolving and can be a good method for calculating molecular parameters. However, the quality of the results is highly dependent on the choice of functional (an early popular one was B3YLP, which is no longer recommended for most computations). Different functionals give better results for different kinds of computations. Another issue with DFT computations is that they are sensitive to the spacing chosen between points on the 3-D grid at which the electron densities are computed.

Examples of output from quantum computations

Modern quantum software produces lots of information about a molecular species. The most obvious are the shapes of the molecules, their orbitals and the orbital energies. Combined with modern graphics software this can help chemists visualize what molecules look like and how they might interact. The figure  below shows some images for methyl fluoride.

a)

b)

c)

d)

Figure \(\PageIndex{1}\): HF-SCF calculations using an ACCD basis set on fluoromethane: a) geometry; b) electrostatic potential, partial charges and vector representing the dipole moment; c) & d) the two degenerate HOMO orbitals.

We can compare the computed bond lengths in figure \(\PageIndex{1}\)a with experimental values listed in the NIST computational chemistry comparison database. The listed bond lengths are r_CH = 0.1087 nm, r_CF = 0.1383 nm and bond angles are a_HCH = 110.2˚ and a_HCF = 108.73˚. There are slight discrepencies from the experimental geometry, but they are very small, suggesting the ACCD basis set gives a reasonable approximation to the structure of the molecule. The dipole moment and electrostatic potential map (figure \(\PageIndex{1}\)b) are consistent with what we would expect from the much higher electronegativity of F compared to the other atoms in the molecule. In figure \(\PageIndex{1}\) c & d we can visually see that the two degenerate HOMO orbitals involve primarily p-type orbitals on the F atom, but on the methyl end the orbitals involve combinations of p-type orbitals on the C and s-type orbitals on the hydrogens. These fractional contributions to the MOs can be understood more quantitatively by examining the numerical coefficients multiplying the AO contributions to the MOs.

The figure below is a calculation on a fragment of a cold virus showing where it is in close contact (bullseyes) with the cell surface receptor molecule it uses to enter the cell.

Figure \(\PageIndex{2}\): The binding site on the rhinovirus capsid in contact with the receptor it uses to enter the cell.

MO calculations can also be used to generate estimates of UV-Vis, NMR, IR and Raman spectra as well. Unless great care is taken the generated spectra will have peaks significantly shifted from the experimental values. Even with these shifts they can often be used to assign the origin of peaks in the spectra to particular spectroscopic transitions. See the figure below for a comparison of an experimental and simulated IR spectrum for methylfluoride.

 

a)

b)

Figure \(\PageIndex{3}\): a) Experimental IR spectrum (from NIST Chemistry Webbook); b) Simulated spectrum generated from a GAMESS computation on methylfluoride using the ACCD basis set. Notice that the C-H stretching wavenumbers seen at about 3000 cm^-1 in the experimental spectrum are much too high in the simulated spectrum. This reflects the fact that the simulation assumes a harmonic vibrational potential, while the C-H bond stretch is extremely anharmonic leading to more closely spaced quantum states and a lower energy vibrational transition.

References

1. C. J. Cramer Essentials of Computational Chemistry, 2nd Ed. (John Wiley & Sons, Ltd, Hoboken, NJ, 2004) Ch. 4.

2. C. J. Cramer Essentials of Computational Chemistry, 2nd Ed. (John Wiley & Sons, Ltd, Hoboken, NJ, 2004) Ch. 8.

Contributors

• Jonathan Gutow (UW Oshkosh)