Introduction

What do we care about?

Things chemistry/materials-related:

What are the properties of atoms?
What molecules do they make? What other substances do they make?
What are the shapes of those molecules? Structures of those solids? Properties of them?
How do those substances react with each other?
What are the energies of those reactions?
What are the rates of those reactions?
What is the strongest substance?
How do we make a substance to do....?
add your own questions....

Things that relate to the chemical properties of substances.

How are we going to figure these out? With only a computer?

1926: Erwin Schrödinger equation: $\hat{H}Ψ=EΨ$
1929: Paul Dirac, British physicist

The fundamental laws necessary for the mathematical treatment of a large part of physics and the whole of chemistry are thus completely known, and the difficulty lies only in the fact that application of these laws leads to equations that are too complex to be solved.

It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation.

1930’s-1950’s: Elaboration, analytical applications
1950’s: Computers start to appear for technical applications
1960’s: Density functional theory emerges.
1960’s-1970’s: Numerical solutions of Schrödinger equation for atoms/molecules—expert users
1980’s: “Supercomputer” era—routine use of computational chemistry software becomes possible

1990’s: “Chemical accuracy” era—very precise solutions routinely available, for a cost! See Schneider, JPC *1994*.
1990’s: Density functional theory (DFT) allows applications to solids/surfaces/liquids to become common. See Hass, Science, *1998*
1990’s: Visualization moves to the desktop
2000’s: Computational “screening,” materials discovery (Gurkan, J. Phys. Chem. Lett., *2010*), materials genome (https://materialsproject.org/).
Today: Computational chemistry widely integrated into all aspects of chemical, materials, biological research

Computational chemistry is now so vast it is impossible to cover everything completely. We limit ourselves to quantum-mechanics-based calculations.

Our goals

Understand when it is appropriate to use quantum-mechanics methods.
Be able to state the basic theoretical, mathematical, and numerical concepts behind quantum mechanical “wavefunction theory” (WFT) and “density functional theory,” (DFT) calculations.
Understand the terminology and practical issues associated with doing quantum chemical simulations.
Get hands-on experience with these concepts using popular computational tools of today, including GAMESS for molecular systems and Vasp for condensed phase systems.
Learn how to set up, execute, and analyze results in a modern, Python notebook environment.
Learn how to apply the results of quantum chemical simulations to calculate things you care about.
Demonstrate an ability to formulate a problem and apply QM methods to it.
Develop the skills to understand a literature paper in the area.

Reading resources

These notes
Chris Cramer, Essentials of Computational Chemistry, Wiley, 2004
Martin, Electronic Structure, Cambridge, 2004
Sholl and Steckel, Density Functional Theory: A Practical Introduction, Wiley, 2009
Kitchin book, http://kitchingroup.cheme.cmu.edu/dft-book/

Software tools

Notebooks

jupyter/ipython
emacs/org-mode

\newpage

Refresher on Quantum Mechanics

Why quantum mechanics?

Want to describe “mechanics” (equations of motion) of atomic-scale things, like electrons in atoms and molecules

Why? These ultimately determine the energy, the shape, and all the properties of matter.

de Broglie wavelength (1924) \begin{equation} λ = h/p = h/mv \end{equation} \begin{equation} h = \SI{6.626e-34}{J.s} \text{(Planck’s constant)} \end{equation}

	Car	Electron
mass $m$	\SI{1000}{kg}	\SI{9.1e-31}{kg}
velocity $v$	\SI{100}{km/hr}	\SI{0.01}{c}
	typical value on the highway	typical value in an atom
momentum $p$	\SI{2.8e-4}{kg.m/s}	\SI{2.7e-24}{kg.m/s}
wavelength λ	\SI{2.4e-38}{m}	\SI{2.4e-10}{m}
	too small to detect. Classical!	Comparable to size of an atom.
		Must treat with QM!

How to describe wave properties of an electron? Schrödinger equation (1926)

Kinetic energy + Potential energy = Total Energy

Expressed as differential equation (Single particle, non-relativistic): \begin{equation} -\frac{\hbar^2}{2m}∇^2 Ψ(\mathbf{r},t) + V(\mathbf{r},t) Ψ(\mathbf{r},t) = -i \hbar \frac{∂}{∂ t} Ψ(\mathbf{r},t) \end{equation}

If the potential $V$ is time-invariant, can use separation of variables to show that the steady-state, time-independent solutions are characterized by an energy $E$ and described by: \begin{eqnarray} -\frac{\hbar^2}{2m}∇^2 ψ(\mathbf{r}) + V(\mathbf{r}) ψ(\mathbf{r}) = E ψ(\mathbf{r})
Ψ(\mathbf{r},t) = ψ(\mathbf{r})e^-iEt/\hbar \end{eqnarray}

Postulates of non-relativistic quantum mechanics

Notes on constants and units

Resource on physical constants: http://physics.nist.gov/cuu/Constants/ Resource for unit conversions: http://www.digitaldutch.com/unitconverter/

Unit converter available in Calc mode in Gnu emacs highly recommended

	Atomic unit	SI unit	Common unit
Charge	$e = 1$	\SI{1.6021e-19}{C}
Length	$a_0 = 1$ (bohr)	\SI{5.29177e-11}{m}	\SI{0.529177}{Å}
Mass	$m_e = 1$	\SI{9.10938e-31}{kg}
Angular momentum	$\hbar = 1$	\SI{1.054572e-34}{J.s}
Energy	$E_h = 1$ (hartree)	\SI{4.359744e-18}{J}	\SI{27.2114}{eV}
Electrostatic force	$1/(4π\epsilon_0) = 1$	\SI{8.987552e-9}{N.m^2/C^2}
Boltzmann constant		\SI{1.38065e-23}{J\per K}	\SI{8.61733e-5}{eV/K}

Energy units 1 eV = \SI{1.60218e-19}{J} = \SI{96.485}{kJ/mol} = \SI{8065.5}{cm^-1} = \SI{11064}{K.k_B}

Example: Energy states of a particle in a box

System defined by potential experienced by particle:

$V(\mathbf{r}) = 0,\qquad 0 < x,y,z < L$

$V(\mathbf{r}) = ∞,\qquad x,y,z \leq 0,\ x,y,z \geq L$

3D box → 3 degrees of freedom/coordinates

Schrödinger equation \begin{equation} -\frac{\hbar^2}{2m_e} \left ( \frac{∂^2 }{∂ x^2} + \frac{∂^2 }{∂ y^2} + \frac{∂^2 }{∂ z^2} \right ) ψ(x,y,z) = E ψ(x,y,z) \end{equation} \begin{equation} ψ(x,y,z) = 0, \quad x,y,z \leq 0,\ x,y,z \geq L \end{equation}

A second-order, linear, partial differential equation. Boundary value problem. Solve by separation of variables. Postulate $ψ(x,y,z) = X(x)Y(y)Z(z)$. Substituting and rearrange to get

\begin{equation} -\frac{\hbar^2}{2m_e} \left (\frac{1}{X(x)}\frac{∂^2 X(x)}{∂ x^2} + \frac{1}{Y(y)}\frac{∂^2 Y(y)}{∂ y^2} + \frac{1}{Z(z)}\frac{∂^2 Z(z)}{∂ z^2} \right ) = E \qquad 0 < x,y,z <L \end{equation}

ftn x + ftn y + ftn z = constant → each term must be constant.

Equation for each dimension \begin{equation} -\frac{\hbar^2}{2m_e}\frac{∂^2 X(x)}{∂ x^2} = E_x X(x), \qquad X(0)=X(L) = 0 \end{equation}

Seek function that twice differentiated returns itself and satisfies boundary conditions. \begin{equation} X(x) = sin\frac{n_xπ x}{L},\qquad n_x = 1,2,3,\ldots \end{equation}

\begin{equation} E_{n_x} = \frac{n_x^2π^2\hbar^2}{2 m_e L^2} \end{equation}

Solutions called eigenfunctions (or wavefunctions) and eigenvalues. Characterized by quantum numbers, one for each degree of freedom. These (and all QM) solutions have certain special properties, including that they orthonormal and form a complete set.

Normalization

Seek a constant such that the inner eigenfunction product is unity. \begin{eqnarray} C^2 ∫_0^L sin^2 \frac{n_xπ x}{L} dx = C^2 L/2 = 1 → C=±\sqrt{\frac{2}{L}}
X(x) = ±\sqrt{\frac{2}{L}}sin\frac{n_xπ x}{L},\qquad n_x = 1,2,3,\ldots \end{eqnarray}

Orthonormal \begin{equation} \langle X_{n_x} | X_{n^′_x} \rangle = δ_{n_{x},n_x^′}\qquad \text{Dirac notation} \end{equation}

Energy increase with number of nodes.
Is this real? See Ho, J. Phys. Chem. B 2005, 109, 20657. Where is the electron?

Complete set Any function on the same domain that satisfies the same boundary conditions can be represented as a linear combination of these solutions: \begin{eqnarray} f(x) &=& ∑_i \left \{∫_0^L X_i(x) f(x) dx \right \} X_i(x) = ∑_i C_i X_i(x)

f\rangle &= &∑_i

X_i\rangle\langle X_i

f\rangle\quad\text{Dirac notation}

\end{eqnarray} Illustrates idea of a basis set. These functions are the basis in “plane wave” supercell methods.

Three-dimensional solutions \begin{equation} ψ(x,y,z) = X(x)Y(y)Z(z) = \left ( \frac{2}{L} \right )^3/2 sin\frac{n_xπ x}{L}sin\frac{n_yπ y}{L}sin\frac{n_zπ z}{L},\qquad n_x,n_y,n_z=1,2,3,\ldots \end{equation} \begin{equation} \label{eq:2} E = E_x+E_y+E_z=\frac{(n_x²+n_y²+n_z²) π²\hbar²}{2 m L²} \end{equation}

Properties of solutions:

Symmetry of system introduces degeneracy in solutions
Energy depends on volume → pressure!

\newpage

Hydrogen atom: simplest chemical “thing”

Schrödinger equation

Place massive nucleus at origin and describe position of electron in spherical coordinates

\begin{eqnarray} \left\{-\frac{\hbar^2}{2m_e}∇^2 + V(r)\right\} ψ(\bm{r}) &=& E ψ(\bm{r})
V(r) &=& -\frac{e^2}{4π\epsilon_0}\frac{1}{|\bm{r}|} \end{eqnarray} Coulomb potential—our nemesis! Decays slowly with distance. Boundary conditions?

Analytical solutions

Separate: $ψ(r,θ,φ)=R(r)Θ(θ,φ)$
Angular equation $\hat{L}^2Θ(θ,φ) = E_LΘ(θ,φ)$
1. $Θ=Y_{lm_l}(θ,φ)$ are “spherical harmonics”, describe angular motion
2. Azimuthal quantum number $l=0,1,…,n-1$, correspond to $s$, $p$, $d$, \ldots orbital sub-shells; angular “shape,” number of angular nodes, angular momentum of electron
3. Magnetic quantum number $m_l=-l,-l+1,…,l$, \ldots orientation of orbital
Radial equation \[\left \{ -\frac{\hbar^2}{2m_e} \frac{d}{dr^2} + \frac{\hbar^2 l(l+1)}{2 m_e r^2} -\frac{e^2}{4π\epsilon_0}\frac{1}{r}\right \} r R(r) = E r R(r)\]

Solutions are a polynomial * exponential decay. Exponential part called a Slater function. Larger the exponent, faster the decay. Degree of polynomial determined by principle quantum number $n=1,2,\ldots$.

Energy expression, corresponds to our conventional H atom spectrum \[E_n = -\frac{1}{n^2}\left(\frac{e^2}{2 a_0}\right) = -\SI{13.6}{eV}⋅\frac{1}{n^2},\qquad n = 1,2,\ldots \]

Questions: Ionization energy of an H atom? $1s→ 2s$ energy? Thermal populations?

Integrate out angular components to get radial probability function $P_nl(r)=r^2 R_nl^2(r)$ \[\langle r\rangle = ∫ r P_nl(r) dr = \left(\frac{3}{2}n^2-l(l+1)\right)a_0\]

Note darn electron doesn’t want to stay in the Coulomb well! Wavefunction extends beyond the classical region defined by $E_n = V(r_\text{classical})$. This phenomenon is called tunneling, is a purely quantum mechanical effect, is pervasive in chemistry, leading for instance to chemical bonding.

Variational principle

What if we don’t know where to look to find the $R(r)$? Or an analytical solution doesn’t exist? Solve numerically.

$l=0$ case, in atomic units:

\[\left\{-\frac{1}{2}\frac{d^2}{dr^2} -\frac{1}{r}\right\} rR(r) = ErR(r),\quad 0<r∞ \]

Guess something. Must obey appropriate boundary conditions and be a well-behaved function. For example, a Gaussian: \[g_ξ(r) = e^{-ξ r^2}\] Let’s normalize: \[N = \left\{ ∫_0^∞ g_ξ^2(r) r^2 dr \right\}^-1/2 = 2 \left(\frac{8ξ^3}{π}\right)^1/4 \]

\[˜{g}_ξ(r) = N g_ξ(r)\] Now evaluate energy, say for $ξ=1$: \[\langle E \rangle = \langle ˜{g}_1|\hat{H}|˜{g}_1\rangle = \SI{-0.096}{Hartree} \] Hmmm, not very good, much higher in energy than true answer of \SI{-0.5}{Hartree}.

Let’s try adding two Gaussians together, with equal weight: \[b(r) = N_s \left(˜{g}_1(r) + ˜{g}_0.5(r) \right )\] Normalize: \begin{eqnarray} \langle b(r)|b(r)\rangle & = & N_s^2 \left(\langle\tilde{g}_1 | ˜{g}_1\rangle + \langle\tilde{g}_0.5 | ˜{g}_0.5\rangle + 2\langle\tilde{g}_1 | ˜{g}_0.5\rangle\right)
&=& N_s(1+1 + 2 S) 1 \\ N_s & &\frac{1}{\sqrt{2(1+S)}} \end{eqnarray}

Note appearance of “overlap integral” $S=\langle ˜{g}_1|˜{g}_0.5\rangle$, shows how similar or different $g_i$ are.

Re-evaluate energy \[\langle b(r) |\hat{H} |b(r) \rangle = \SI{-0.306}{Hartree} \] Much closer to the truth!

Could even weight the two Gaussians differently: \[c(r) = N_s^′ \left(˜{g}_1(r) + 1.5 ˜{g}_0.5(r) \right )\] \[\langle c(r) |\hat{H} |c(r) \rangle = \SI{-0.333}{Hartree} \] Better yet!

Could continue to add Gaussians of various exponents, and could vary weights, or could even add in any other functions that we want that are “well-behaved.” Would find that no matter what we do, the “model” energy would be greater than the “true” value. Basis of the variational principle:

For any system described by the Hamiltonian $\hat{H}$ and any satisfactory trial wavefunction $Ψ$, \[ E = \frac{\langle Ψ | \hat{H} | Ψ \rangle}{\langle Ψ | \rangle } ≥ E_0\] where $E_0$ is the true ground state energy of the system.

Consequence of the completeness of the solutions of the Schrödinger equation. Extremely important to us, because we can use the calculus of variations to seek energy-optimal $Ψ$.

Basis functions

Recognize that we are approximating

\[R₁₀(r) ≈ f(r) = ∑_i c_i φ_i(r)\]

$φ_i$ are basis functions and $c_i$ are variational parameters, or coefficients.

If we find $c_i$ that minimize energy, then we have an optimal approximation to $E_0$ within our basis, and we are sure that $E_0$ is an upper bound on the truth. Adding more basis functions /must/lower energy.

Common trade-offs:

Want to choose $φ_i$ that are good approximation to the “truth”
Want to choose $φ_i$ that are mathematically convenient

Slaters	Gaussians	Plane waves	Mixed basis
Accurate	Moderate accuracy	Poor accuracy
Expensive	modest cost	cheap!
`ADF`	`Gaussian`, `GAMESS`,	`Vasp`, `CPMD`,	`FLAPW`, `CP2k`
	`NWChem`, `Qchem`	`QuantumEspresso`

Virtually all quantum codes work on this principle, and the main differences are in details of implementation and ancillary functionality provided.

There are exceptions. GPAW for instance solves the QM equations by finite difference expressions on a numerical grid. Lends itself to parallelization and may be the future\ldots remains to be seen!

Secular equations

Apply variational principle to two basis functions for the H atom:

\[f(r) =c_1 φ_1(r) + c_2φ_2(r) \]

\[\langle E \rangle = \frac{\langle f(r)|\hat{H}|f(r)\rangle}{\langle f(r)|f(r)\rangle} \]

Substitute and solve $∂ \langle E \rangle/∂ c_1 = ∂ \langle E \rangle/∂ c_2 = 0$. Each gives a linear secular equation in $c_1$ and $c_2$:

\[ \left(\begin{array}{cc} H₁₁ - E & H₁₂-S₁₂E
H₁₂-S₁₂E & H₂₂ \end{array}\right)\left( \begin{array}{c} c_1 \c_2\end{array}\right) =0 \] where $H_ij = \langle\phi_i|\hat{H}|φ_j\rangle$ is a matrix element and $S_ij=\langle φ_i | \rangle φ_j \rangle$ is an overlap. If $S_ij = 0$, basis is orthogonal, problem simplifies. If $S_ij≈ 1$, basis is redundant, not efficient!

Evaluate secular determinant \[ \left| \begin{array}{cc} H₁₁ - E & H₁₂-S₁₂E \ H₁₂-S₁₂E & H₂₂ \end{array}\right| =0 \] Gives a quadratic in E. Yields two solutions, which would be approximations to the 1s and 2s orbital energies of hydrogen. Back substitution to get coeffients:

1s: $\langle E_1s\rangle > E_{1s,\mathrm{true}}$	$c^1s_1$ and $c^1s_2$
2s: $\langle E_2s\rangle > E_{2s,\mathrm{true}}$	$c^2s_1$ and $c^2s_2$

Note we always get one solution for each basis function. Secular matrix grows as the square of the number of basis functions, gets expensive to find roots.

If basis is not orthogonal, common to orthogonalize. Find linear transformation that makes $\langle S_ij = δ_ij$. Evaluate \[ \bm{c}^′ = \bm{S}^1/2 \bm{c} → \hat{\bm{H}}^′ = \bm{S}^-1/2\bm{H}\bm{S}^1/2\]

\[ \left(\begin{array}{cc} H^′₁₁ - E & H^′₁₂
H^′₁₂ & H^′₂₂ \end{array}\right)\left( \begin{array}{c} c^′_1 \c^′_2\end{array}\right) =0 \]

\[\bm{H}^′\bm{c}^′ = E \bm{c}^′ \]

Secular equations reduce to standard linear eigenvalue problem. All the tricks of linear algebra can be applied to find the eigenvalues (orbital energies) and eigenvectors (wavefunctions). Called diagonalizing the matrix. Tricks can be used to find the lowest energy roots only.

Same basic idea is used in virtually all calculations on atoms, molecules, \ldots. Basis of semi-empirical and first principles methods.

Spin

Can’t leave the H atom without mentioning electron spin. Non-relativistic QM gives us three quantum numbers. Relativity teaches that space and time are equivalent. Relativistic H atom solutions introduce a 4th degree of freedom that, under many circumstances, decouples from other three. Call it the electron spin, because it behaves like the electron has an intrinsic quantum angular momentum with magnitude $s = 1/2$.

\[m_s = +1/2, \quad\text{“spin up”},\quad\alpha\] \[m_s = -1/2, \quad\text{“spin down”},\quad\beta\]

$m_s$ specifies $z$ component of angular momentum, $s_z = m_s\hbar$.

To fully specify state of H atom, must specify all four quantum numbers. \newpage

(Two is too) many electrons

Helium: next (after hydrogen) simplest atom

In a sense, we “know” the answer\ldots $1s^2$. But is this same $1s$ as H? No! Different nuclear charge, interactions between the two electrons. This is an approximation and a very convenient shorthand!

Schrödinger equation for He

Wavefunction $Ψ(\mathbf{r}_1,\mathbf{r}_2)$, atom energy $E$.

Define 1-electron operator for each electron, in atomic units. Include kinetic energy of electron and its attraction to nucleus of charge $Z=2$: \[\hat{h}_i = -\frac{1}{2}∇^2_i -\frac{Z}{|\mathbf{r}_i|}\] Looks similar to hydrogen atom.

BUT, electrons also repel. Total Schrödinger equation for He: \[\left\{\hat{h}_1 + \hat{h}_2 + \frac{1}{|\mathbf{r}_2 - \mathbf{r}_1|\right\}Ψ(\mathbf{r}_1,\mathbf{r}_2) = EΨ(\mathbf{r}_1,\mathbf{r}_2) }\] Last term accounts for electron-electron electrostatic repulsion. Makes problem non-separable and really hard to solve. (How many solutions are there?)

Generalize to $n$-electron atom, in atomic units: \begin{eqnarray} \hat{H} & = & ∑_i \hat{h}_i + ∑_j>i+1\frac{1}{|\mathbf{r}_j - \mathbf{r}_i|}
\hat{H} Ψ(\mathbf{r}_1,\ldots,\mathbf{r}_n) & = & EΨ(\mathbf{r}_1,\ldots,\mathbf{r}_n) \end{eqnarray}

First summation over all electrons, second gets all electron pairs.

Solutions are many-dimensional functions of the coordinates of all the electrons. Cannot solve this analytically, although approaches exist (eg quantum Monte Carlo) that can in principle get very close. Thankfully, though, we can make approximations that work out really well. We’ll look at three historically important ones.

The Hartree atom

Simplest approach is to approximate $Ψ$. Douglas Hartree (1897-1958) writes: \[ Ψ(\mathbf{r}_1,\mathbf{r}_2) ≈ ψ_1(\mathbf{r}_1)⋅ ψ_2(\mathbf{r}_2)\] So-called Hartree product. Can’t be right. It gives the probability of two electrons being in the same place as some number $> 0$! Neglects electron correlation. How to apply?

Apply variational principle: What’s the best possible set of $ψ_i$? We’ll say best are the set that give the lowest expectation value of energy. \begin{eqnarray} \langle E \rangle & = & \langle Ψ | \hat{H} | Ψ \rangle /\langle Ψ | Ψ\rangle
\frac{ δ \langle E \rangle}{δ ψ_i} &= &0, ∀ i \end{eqnarray}
Lagrange multipliers to impose orthonormality constraint on $ψ_i$: \begin{eqnarray} \langle ψ_i | ψ_j \rangle &=& δ_ij
L & = & \langle E \rangle - ∑_i,jε_ij\left ( \langle ψ_i | ψ_j \rangle - δ_ij \right ) \ δ L &=& 0 \end{eqnarray}
Coupled, one-electron Hartree eigenvalue equations for energy-optimal $ψ_i$: \begin{eqnarray} \left \{ \hat{h}_i + \hat{v}^{\text{Hartree}}_i \right \} ψ_i(\mathbf{r}_1) &=& ε_i ψ_i(\mathbf{r}_1)
\hat{v}^{\text{Hartree}}_i(\mathbf{r}_i) & =& ∑_{j≠ i} ∫ \left | ψ_j(\mathbf{r}_2) \right |^2\frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \end{eqnarray}

Have to solve this for all n electrons of an atom/molecule. “Hartree potential” represents Coulomb repulsion between electron i and all other electrons, averaged over position of those electrons. Always positive. This is a mean field approximation. Note appearance of “one electron” energies, $ε_i$, kinetic energy plus repulsion of electron with all others. Total energy is sum of these $ε_i$ corrected to avoid overcounting repulsions: \[ \langle E \rangle = ∑_i ε_i - \frac{1}{2}∑_i \langle ψ_i| \hat{v}^{\text{Hartree}}_i|ψ_i\rangle \]

Presents an obvious difficulty. If we don’t know $ψ_j$ ahead of time, how can we even construct Hartree equations, let alone solve them? Hartree offered a numerical solution, in the 1930’s, called the /self-consistent field/ (SCF) approach:

Guess an initial set of $ψ_i$, one for each electron (he did this on a grid, and jumping ahead a bit, allowed each $ψ_i$ to represent two electrons)
Construct Hartree potential for each $ψ_i$
Solve the $n$ differential equations for $n$ new $ψ_i$
Compare new to old $ψ_i$
If the same within a desired tolerance, you are done!
If not, return to step 2, using new $ψ_i$, and repeat.

Hartree’s father did this by hand for all the atoms of the periodic table, tabulating wavefunctions and energies for all the electrons in each. See Hartree, Douglas R. The Calculation of Atomic Structures (1957). For instance, for He, he’d solve one equation, self-consistently, to get one $ψ_1$, and then combine to get $Ψ(1,2) = ψ_1(1)α(1)ψ_1(2)β_2$. Tedious! Qualitatively great, quantitatively not so hot. Mean-field approximation just not so hot.

Nonetheless, basic idea of representing many-body wavefunction in terms of “orbitals,” of setting up orbital equations, and solving using a self-consistent procedure, remain today at the heart of virtually all electronic structure calculations. Hurrah Hartree!

Note: It would be very cool to write a simple Python code to illustrate the SCF procedure for two electrons in an atom. Could be done on a grid or in a basis. See eg http://www.users.csbsju.edu/~frioux/scf/scf-he.pdf.

The Pauli principle

One big conceptual short-coming of the Hartree model is that it treats the electrons as if they were distinguishable. QM says electrons are indistinguishable. Furthermore, they have a quantized angular momentum, called a spin, that is either up or down, making them fermions.

Pauli principle: The wavefunction of a multi-particle fermion system must be anti-symmetric to coordinate exchange.

\[Ψ(\bm{x}_1, \bm{x}_2) = -Ψ(\bm{x}_2, \bm{x}_1) \]

Here the coordinate x includes both the position and the spin (up or down, α or β) of the electron.

Sorry Hartree. Can fix for He by writing \[Ψ(\bm{x}_1,\bm{x}_2) = ψ_1(\bm{r}_1)ψ_1(\bm{r}_2)\left(α(1)β(2) - β(1)α(2)\right) \] Hey, gentle reader, check, does this work? Yes! Exchanging the coordinates changes the sign but keeps everything else the same. Normalizing is easy if we take $ψ_1$ to be normalized and recall that spin functions are orthogonal: \[Ψ(\bm{x}_1,\bm{x}_2) = \frac{1}{\sqrt{2}}ψ_1(\bm{r}_1)ψ_1(\bm{r}_2)\left(α(1)β(2) - β(1)α(2)\right) \]

Note it is impossible to construct an antisymmetric wavefunction in which both electrons have the same spatial function and the same spin. Two electrons cannot have the same space and spin variables.

Slater determinants and Hartree-Fock

Slater determinant a general way to assure that a wavefunction satisfies Pauli principle: \[ Ψ = \frac{1}{\sqrt{n!}}\left | \begin{array}{cccc} ψ_1(1) & ψ_2(1) & \cdots & ψ_n(1)
ψ_1(2) & ψ_2(2) & \cdots & ψ_n(2) \ \vdots & \vdots & \ddots & \vdots \ ψ_1(n) & ψ_2(n) & \cdots & ψ_n(n) \end{array} \right | = |ψ_1ψ_2\cdots\psi_n\rangle \] Swapping rows swaps coordinates and, by rules of determinants, changes sign.

Let’s compare. Two spin-paired electrons in two different orbitals: \[ \left | \begin{array}{cc} ψ_1(1)α(1) & ψ_2(1) β(1)
ψ_1(2)α(2) & ψ_2(2) β(2)\end{array} \right | = ψ_1(1)ψ_2(2)α(1)β(2)-ψ_2(1)ψ_1(2)β(1)α(2) \] Antisymmetric? What happens when the two electrons have the same spatial coordinate?

Two spin-aligned electrons in two different orbitals: \[ \left | \begin{array}{cc} ψ_1(1)α(1) & ψ_2(1) α(1)
ψ_1(2)α(2) & ψ_2(2) α(2) \end{array} \right | = \left (ψ_1(1)ψ_2(2)-ψ_2(1)ψ_1(2)\right ) α(1)α(2) \] What happens now?

Exchange guarantees that two electrons of same spin cannot be in the same place! No such guarantee for electrons of opposite spin.

Hartree-Fock equation

Same song and dance:

Apply variational principle to Slater determinant. The best $ψ_i$ are those that minimize the expectation value of the energy.
Use method of Lagrange multipliers to keep $ψ_i$ orthogonal.

For simplicity, restrict ourself to cases with an even number of electrons N, all spin-paired, so two electrons in every orbital. Let index j run over all occupied orbitals. Arrive at restricted Hartree-Fock equations: \begin{eqnarray} \left \{ \hat{h}_i + \hat{v}^{\text{Hartree}}_i \right + \hat{v}^{\text{exchange}}\} ψ_i(\mathbf{r}_1) &=& ε_i ψ_i(\mathbf{r}_1)
\hat{v}^{\text{Hartree}}_i(\mathbf{r}_1) & =& 2 ∑_{j≠ i} ∫ \left | ψ_j(\mathbf{r}_2) \right |^2\frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \ \hat{v}^{\text{exchange}}_i(\mathbf{r}_1) ψ_i(\mathbf{r}_1) & = &-ψ_j(\mathbf{r}_1) ∑_{j≠ i} ∫ ψ_j (\mathbf{r}_2) ⋅ ψ_i (\mathbf{r}_2) \frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \end{eqnarray}

Yikes! Slater determinant wavefunction results in appearance of the “exchange” operator, which turns a $ψ_i$ into a $ψ_j$. Exchange operator is not a simple multiplication. Must be solved self-consistently, and is much harder to do than the simple Hartree expression.

Slight simplification possible, noting that $i=j$ terms cancel out and slightly redefining operators: \begin{eqnarray} \left \{ \hat{h}_i + \hat{v}^{\text{Hartree}}_i \right + \hat{v}^{\text{exchange}}\} ψ_i(\mathbf{r}_1) &=& ε_i ψ_i(\mathbf{r}_1)
\hat{v}^{\text{Hartree}}(\mathbf{r}_1) & =& 2 ∑_j ∫ \left | ψ_j(\mathbf{r}_2) \right |^2\frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \ \hat{v}^{\text{exchange}}_i(\mathbf{r}_1) ψ_i(\mathbf{r}_1) & = &-ψ_j(\mathbf{r}_1) ∑_j ∫ ψ_j (\mathbf{r}_2) ⋅ ψ_i (\mathbf{r}_2) \frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \end{eqnarray}

Now “Hartree potential” is the same for all orbitals/electrons. We can define the “charge density” to be \[ ρ(\mathbf{r}) = 2 ∑_j |ψ_j(\mathbf{r})^2|\] (units of charge/unit volume, multiply by e to get a charge). The Hartree potential can be written \[ \hat{v}^{\text{Hartree}}(\mathbf{r}_1) = ∫\frac{ρ(\mathbf{r}_2)}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \] This is the Coulomb repulsion of an electron with all electrons, including itself! Called Poisson equation, well known in classical physics. Because it involves a Coulomb repulsion, will see either $\hat{v}^{\text{Hartree}}$ or $\hat{v}^{\text{Coulomb}}$. I’ll often write $\hat{v}^{\text{Coulomb}}[ρ(\mathbf{r})]$, to emphasize that the Coulomb potential is a functional of the charge density.

(Just to make sure we are following units around, the Coulomb potential has units of energy/charge, eg in SI it would be J/C and would have $e/4π\epsilon_0$ in front.)

The “exchange potential” cancels out the “self-interaction” of an electron with itself, and insures that two electrons of the same spin cannot be in the same place, ie, the wavefunction vanishes whenever the spatial coordinates of two electrons are the same. It cannot be written simply in terms of the charge density.

Basis of wavefunction theory (WFT)

Hartree-Fock model is much better than Hartree alone, widely implemented in codes. Not particularly good by today’s standards. However, it is systematic and rigorous, by requiring exact adherence to the Pauli principle, and it can be systematically improved. It is the foundational basis of all wavefunction theory (WFT) models, all of which are characterized by exactly treating exchange. The only approach some people call ab initio.

Hartree-Fock-Slater

In 1951 John Slater introduced an approximation to the Hartree-Fock model that turned out to anticipate a whole new approach to solving the electronic structure problem, called density functional theory.

Rewrite exchange part as (and shorten “exchange” to “x”): \begin{eqnarray} \hat{v}^\text{x}_i(\mathbf{r}_1) ψ_i(\mathbf{r}_1) & = &-ψ_j(\mathbf{r}_1) ∑_j ∫ ψ_j (\mathbf{r}_2) ⋅ ψ_i (\mathbf{r}_2) \frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2
& = & -ψ_j(\mathbf{r}_1) ∑_j ∫ ψ_j (\mathbf{r}_2) ⋅ ψ_i (\mathbf{r}_2) \frac{1}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 ⋅ \frac{ψ_i(\mathbf{r_1})ψ_i(\mathbf{r_1})}{ψ_i(\mathbf{r_1})ψ_i(\mathbf{r_1})} \ & = & -\left [ ∫ \frac{ρ_i^\text{x}(\mathbf{r}_1;\mathbf{r}_2)}{|\mathbf{r}_2-\mathbf{r}_1|}d\mathbf{r}_2 \right ] ψ_i(\mathbf{r}_1) \ ρ_i^\text{x}(\mathbf{r}_1;\mathbf{r}_2)& =& ∑_j \frac{ψ_i(\mathbf{r_1})ψ_j(\mathbf{r_2})ψ_j(\mathbf{r_1})ψ_i(\mathbf{r_2})}{ψ_i(\mathbf{r_1})ψ_i(\mathbf{r_2})} \end{eqnarray}

This looks like the Coulomb expression, but the density thing is different for each orbital i. The “exchange density” does have units of charge density enters in the same way, but with minus sign, to the electron density. Suggests that exchange can be thought of as an electron “hole” around an electron. This exchange density has some special properties:

Every electron at any position $\mathbf{r}_1$ has an exchange hole around it equal to one electron of the same spin as itself: \[ ∫ ρ_i^\text{x}(\mathbf{r}_1;\mathbf{r}_2) d\mathbf{r}_2 = 1 \]
The exchange hole exactly cancels out all electrons of the same spin at the electron location. Ie, it exactly fixes self-interaction: \[ρ_i^\text{x}(\mathbf{r}_1;\mathbf{r}_1) = ∑_j |ψ_j(\mathbf{r})^2|\]

Thus, the Coulomb repulsion felt by an electron is diminished by an exchange hole that follows the electron around, exactly canceling out the charge at its current location. It’s not necessarily spherical and is not the same for all orbitals, but the fact that it has these general properties gives hope that it can be approximated somehow.

Hey, I have an idea! (Actually, Slater had an idea.) What if we had a homogeneous (density the same everywhere) gas of electrons, like electrons of a given density $ρ$ in an infinite box? By symmetry the exchange hole would be spherical, and if it must integrate to 1, then it must have a radius (factor of 2 comes from fact we are only including electrons of the same spin): \[R_\text{hole} = \left [\frac{4π\rho/2}{3} \right ]^1/3 \] The potential felt by an electron due to this spherical hole is \[\hat{v}^\ext{x} = -\frac{1}{2}∫_\text{sphere}\frac{ρ}{r}dr = - \left [\frac{9π\rho}{4}\right ]^1/3 \]

Now, let’s assume that an electron in a real system experiences an exchange hole potential at any point exactly like that of a homogeneous electron gas of the same density at that point. This is the basis of the Hartree-Fock-Slater model: \[\hat{v}^\ext{x,HFS}(\mathbf{r}_1) =-\frac{3}{2}\left[\frac{3ρ(\mathbf{r}_1)}{π} \right] -Cρ(\mathbf{r}_1)^1/3 \] Some ambiguity as to the right value of the constant C, so sometimes just taken as a parameter.

Can now write the Hartree-Fock-Slater equation: \[ \left\{ \hat{h} + \hat{v}^\text{Coulomb}[ρ] + \hat{v}^\text{x,HFS}[ρ]\right \} ψ_i(\mathbf{r}) = ε_iψ_i(\mathbf{r}) \] This is much simpler to solve than Hartree-Fock equation, because the left hand side is the same for all electrons given a total density $ρ(\mathbf{r})$. Still must be solved iteratively, using the self-consistent field.

Notes

Exchange potential scales with the total number of electrons around: more electrons (like near a nucleus) means a more compact, denser exchange hole, more electron exchange “screening,” a greater decrease in potential energy. Further from nucleus, more diffuse exchange hole, less screening.
Screening is not exact, though; does not exactly cancel self-interaction. Clearest in one-electron case: Coulomb and exchange potentials should then exactly cancel, which they evidently do not! HFS energy of an H atom is not exactly −0.5 au!
From a computational point of view, the exchange potential goes from being the hard thing to evaluate to being the easy thing. The Coulomb potential takes more effort to evaluate, and tricks are often implemented to simplify that, like fitting the density to an auxiliary basis set. On the other hand, the 1/3 power makes the exchange potential non-analytic, and solution of the HFS equation (and all DFT methods) involves some form of numerical quadrature.
How does the HFS model do? Pretty darn well, in particular for calculating the structures of things, and it works nicely for things like metals. Not a bad deal! Way to go, Slater!
Another aside: back in the day, the numerical implementations of Xα were very crude and sometimes gave unreasonable results (like linear water!). Slater still sold it very hard, which did not enamor him or all of DFT to the chemical community, although the physics community was far more accepting. For many years DFT was unaccepted by chemists, until solid numerical implementations in codes like Gaussian brought it to the mainstream.

Basis of density functional theory (DFT)

Slater’s arguments are not rigorous. However, as we will see later, they can be made rigorous. HFS is the very simplest example of a density functional theory model, because it is a model built entirely on charge density. Such approach is justifiable.

Implementations

`GAMESS`

Hartree-Fock method always paired with basis set methods and implemented in the codes available at http://webmo.net. Example GAMESS input for Hartree-Fock Ar:

 $CONTRL SCFTYP=RHF RUNTYP=ENERGY ISPHER=1
       ICHARG=0 MULT=1 COORD=CART $END
 $BASIS GBASIS=CCT $END
 $DATA
Ar
OH 1

Ar 18 0.00000000 0.00000000 0.00000000
 $END

And for Hartree-Fock-Slater Ar:

 $CONTRL SCFTYP= RHF RUNTYP=ENERGY DFTTYP=Slater ISPHER=1
       ICHARG=0 MULT=1 COORD=CART $END
 $BASIS GBASIS=CCT $END
 $DATA
Ar
OH 1

Ar 18 0.00000000 0.00000000 0.00000000
 $END

`FDA`

Hartree-Fock-Slater is intrinsically numerical. Historically interesting is the Herman-Skillman code, that solves the problem numerically on a grid. Available to us as the fda code, see https://www.chemsoft.ch/qc/fda.htm and ./Resources/00READ.ME.

Ar input:

300 0.0001 30.0
50 0.00001 0.10 0.50  0.682 0.0042
18.0 5
1 0 1.0 1.0
2 0 1.0 1.0
2 1 3.0 3.0
3 0 1.0 1.0
3 1 3.0 3.0

Output at ./Resources/Ar.out.

Performance

One metric is the ability to predict ionization energies.

Koopman’s theorem: The negative of the energy of an occupied orbital ($-ε_i)$) approximates the energy to extract an electron from that orbital, ie to ionize the system. The energy of a virtual orbital approximates the energy to add an additional electron to a system, i.e. the electron affinity. Assumes no relaxation of orbitals.

Correlation

If solved to reasonable precision, both the Hartree-Fock and Hartree-Fock-Slater models work pretty darn well for things like the shapes of molecules, structures of solids, charge distributions, vibrational frequencies, .... Don’t work so well for computing things that involve making and breaking bonds, like a reaction energy or an activation energy.

Why? All the models discussed here neglect electron correlation, the fact that the potential felt by an electron is a function of the instantaneous positions of all the other electrons. The contribution of correlation to absolute energies is not big by proportion, but it is very imporant to energy differences. Any “orbital” model cannot capture correlation. It can be introduced systematically and exactly into H-F models (at great computational expense) and systematically and approximately into DFT models (at much more modest expense). Hence the popularity of DFT!

\newpage

Practical electronic structure

Born-Oppenheimer approximation

In principle all nuclei and electrons should be described quantum mechanically. For \ce{H2}, for instance, true wavefunction would be a function of the positions of nuclei and electrons, $Υ(\mathbf{r}₁, \mathbf{r}₂,\mathbf{R}_α,\mathbf{R}_β)$.

Nuclei much heavier than electrons and move much more slowly. Assume nuclei are fixed in space (“clamped”) and electrons move in static field of those electrons. Equivalent to assuming that nuclear kinetic energy is decoupled from electron dynamics. Only change is that \[\hat{h} = -\frac{1}{2}∇^2 - ∑_α \frac{Z_α}{|\bm{r}-\bm{R}_α|} \]

Schrödinger equation becomes parameteric in nuclear positions; solutions $E(\mathbf{R}_α,\mathbf{R}_β)$ define a potential energy surface (PES).

\[E_\text{PES}(\mathbf{R}_α,\mathbf{R}_β) = E_\text{Schr} +\frac{1}{2}∑_α,β\frac{Z_α Z_β }{|\mathbf{R}_β - \mathbf{R}_α|}\]

Model chemistry

Essentially always start with \begin{equation} \left \{ \hat{h} +v_{\text{Coulomb}}[ρ] + v_\text{exchange}[ψ_i] + v_\text{correlation}[ψ_i]\right\}ψ_i(\mathbf{r}) =ε_i ψ_i(\mathbf{r}) \end{equation} label:fock

Standard models of today all treat the one-electron and Coulomb pieces exactly and treat the electron-electron interactions at various levels of approximation.

	$v_{\text{exchange}}$	$v_{\text{correlation}}$
Wave function theory (WFT)
Hartree	self-interaction	neglect	historic
Hartree-Fock	exact	neglect	superceded
MPn, CC	exact	perturbative	state-of-the-art
CI	exact	variational	specialized
Density functional theory (DFT)
Hartree-Fock-Slater	$[ρ^{4/3}]$	neglect	historic
Local density approximation	$[ρ^{4/3}]$	$[ρ]$	general purpose solids
(LDA)
Generalized gradient approximation	$[ρ,∇\rho]$	$[ρ,∇\rho]$	general purpose solids/surfaces
(GGA)
“Improved” GGA	$[ρ,∇\rho]$	$[ρ,∇\rho]$	general purpose
(RPBE, BEEF, Mxx)
Hybrid	$≈$ exact	$[ρ,∇\rho]$	general purpose molecules
(B3LYP, PBE0, HSE06)			specialty solids/surfaces
Meta GGA	$[ρ,∇\rho,∇^2ρ]$	$[ρ,∇\rho,∇^2ρ]$	developing

The choice of the electronic structure model is the most fundamental approximation in applying these methods. Determined from experience and need.

Specification in GAMESS (https://www.msg.chem.iastate.edu/GAMESS/GAMESS.html) is a bit arcane. Default is Hartree-Fock. To specify DFT model, use

$CONTRL DFTTYP =   Slater (HFS), SVWN (LDA), PBE (GGA), B3LYP (Hybrid), M06 (Minnesota optimized)

Beyond Hartree-Fock

Many methods available. See manual for full description. Most common is second-order perturbation theory, “MP2”:

$CONTRL MPLEVL=2 $END

If you want a very high quality number, have a big computer and time to wait, try “coupled cluster”:

$CONTRL CCTYP=CCSD(T) $END

Bring back the basis sets

The one-electron equations eq ref:fock give us defining expressions for the energy-optimal orbitals, but they aren’t convenient to solve for anything more complicated than an atom. Expand solutions in a basis set: \[ψ_i(\bm{r}) = ∑_ν C_{μ i}φ_ν(\bm{r}) \] Often atom-centered. You’ll see the term “linear combination of atomic orbitals,” LCAO.

Abbreviate $\hat{f}ψ_i = ε_iψ_i\]. Substitute in \(ψ_i$, mulitple through by a basis function $φ_μ$: \[∑_ν F_{μ ν} C_{ν i} = ε_i ∑_ν S_{μ ν} C_{ν i}, \qquad \bm{FC} = \bm{SC}ε \] where \[ F_{μ ν} = \langle ψ_μ|\hat{f}|φ_ν\rangle\qquad S_{μ ν} = \langle ψ_μ|φ_ν\rangle\] Matrix equation to solve.

Historically interesting, “semi-empirical” methods (MNDO, \ldots) worked by parameterizing the matrix elements against atom properties.

Recall $\hat{f}$ depends on the density, which can be written \[ ρ(\bm{r}) = ∑_{μ ν} P_{μ ν}φ_μ(\bm{r})φ_ν(\bm{r}),\qquad P_{μ ν}=2∑_i C_{μ i} C_{ν i} \] Depending on implementation, pieces of $\hat{f}$ can often be computed just once and reused, eg one-electron integrals $\langle\phi_μ|\hat{h}|φ_ν\rangle$.

Algorithm:

Put your atoms somewhere in space
Select a basis
Pre-compute what you can
Guess some coefficients/density/density matrix
Construct secular matrix elements
Solve secular matrix equation for $C$ and $ε$
Construct and compare new density to old
Update density and repeat, or \ldots
\ldots if less than tolerance, all done!

ALWAYS check to be sure result has converged, to the state you want!

H2O Energy Example

Hartree-Fock calculation on \ce{H2O}, minimal (STO-3G) basis set.

!   File created by the GAMESS Input Deck Generator Plugin for Avogadro
 $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=RHF RUNTYP=ENERGY COORD=CART $END
 $DATA 
Title: H2O energy evaluation
C1
O     8.0    -0.89600     3.13196     0.00000
H     1.0     0.07400     3.13196     0.00000
H     1.0    -1.21933     3.71670     0.70316
 $END

See ./Resources/H2O-STO3G.gamout.

Symmetry

Often problem can be simplified by taking advantage of symmetry of the system.

!   File created by the GAMESS Input Deck Generator Plugin for Avogadro
 $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=RHF RUNTYP=ENERGY COORD=CART $END
 $DATA 
Title: H2O energy evaluation
CNV 2

O     8.0    -0.89600     3.13196     0.00000
H     1.0     0.07400     3.13196     0.00000
H     1.0    -1.21933     3.71670     0.70316
 $END

Results get labeled by symmetry labels.

See ./Resources/H2O-C2V.gamout.

Examples

Dissociating H2+ example

Compute energy vs distance. Should dissociate to H atom and \ce{H+} ion.

Oops, come on, LDA! Illustrates self-interaction problem in LDA. Electron is too eager to be diffuse, spreads out over both atoms when it should localize on one.

Dissociating HHe+ example

Compare Hartree-Fock and LDA for \ce{H-He+} vs distance. (Isoelectronic to \ce{H2}, but avoids any problems with symmetry. Should dissociate to \ce{H+} and He. Does it?

 $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=RHF RUNTYP=ENERGY ICHARG=1 MULT=1 $END
 $DATA 
Title
C1
H     1.0     0.   0.  0.
He    2.0     0.   0.  XXX
 $END

Equilibrium distance? How’s the dissociation state? Bond energy? Truth is about \SI{-0.075}{Hartree}. LDA has advantage of cancellation of errors between exchange and correlation errors. A good thing!

Dissociated He2 example

Open-shell systems

First, some jargon related to unpaired electrons:

# unpaired electrons	$S$	$2S+1$	name
0	0	1	singlet
1	$1/2$	2	doublet
2	1	3	triplet
3	$3/2$	4	quartet

Model has to be generalized somewhat to deal with systems with unpaired electrons. One approach is to construct wavefunctions that are exactly spin-adapted (eigenfunctions of the $\hat{S}$ operator). Possible in the Hartree-Fock world, but messy. More common is to relax that constraint a bit, define different orbital wavefunctions for spin-up and spin-down electrons, called unrestricted or spin-polarized (opposite of non-spin-polarized!). Means that electron density has different spin-up and spin-down parts.

Controlled in GAMESS using the $CONTRL group:

$CONTRL SCFTYP = RHF   non-spin-polarized, default
        SCFTYP = UHF   spin-polarized
        MULT   = 1 (default), 2,...  spin multiplicity = 1 + number of unpaired electrons
        ICHARG = 0 (default), 1,...  net charge
$END

Gaussian basis sets

Gaussian functions ($e^{-ζ|\mathbf{r}|^2}$) are the most popular choice for atom-centered basis sets. They do not efficiently represent molecular wavefunctions, but one- and two-electron integrals in WFT can be solved analytically over Gaussians.

Other choices, like Slater functions ($e^{-ζ|\mathbf{r}|}$) are possible but require numerical quadrature.

Gaussian basis sets have to be created for any given atom and must be used consistently within a set of calculations.

Primitive is a single Gaussian function, possibly multipled by a polynomial to look like s, p, \ldots. Defined by an exponent $ζ$ that determines how extensive (small $ζ$) or compact the function is.
Contraction is a pre-set linear combination of several primitive Gaussians.
Basis set is a predefined set of exponents and contraction coefficients appropriate for some specific atom.

Gaussian basis set nomenclature

Minimal basis contains one contracted function for every atomic orbital. STO-3G is the poster child.
Double zeta contains two contracted functions for every atomic orbital
Split valence is more common, single zeta in core, double zete in valence, typical of Pople basis sets, eg “6-31G”
Triple-split valence would be “6-311G”
Spherical vs Cartesian determines whether a d function have 6 or 5 parts (the sixth of which is an s).
Polarization functions are functions of one angular momentum greater than the highest angular momentum occupied states, eg $p$ function for H, or $d$ function for C. Important to capture the polarization of charge when atoms make molecules. Arcane nomenclature, eg 6-31G(d,p) or 6-31G**.
Diffuse functions are small exponent functions to describe anions or loosely bound electrons. Again argane nomenclature, eg 6-31+G(d,p). Yech.
Correlation consistent and atomic natural orbital are series of basis sets that are constructured to be efficient and to improve systematically.
Complete basis set (CBS) limit is notion of extrapolating energies from a series of systematically improving basis sets. Very common in high accuracy calculations.

Standard basis sets in `GAMESS`

Specified in $BASIS group. Some common choices, in increasing level of sophistication:

Name	Type	Flags
	Pople type	The most venerable and widely used
STO-3G	Minimal	GBASIS=STO NGAUSS = 3
3-21G	Split valence	GBASIS=N21 NGAUSS=3
6-31G(d)	Split valence polarized	GBASIS=N31 NGAUSS =6 NDFUNC=1
6-311+G(d,p)	Triple-split valence	GBASIS=N311 NGAUSS=6 NDFUNC=1 NPFUNC=1
	polarized and augmented	DIFFSP=1

	Polarization-consistent	Good for DFT
PC0	Split valence	GBASIS=PCseg-0 ISPHER=1
PC1	Split valence polarized	GBASIS=PCseg-1 ISPHER=1
PC2	Triple split double polarized	GBASIS=PCseg-22 ISPHER=1

	Correlation-consistent	Good for MP2 and beyond
cc-pVDZ	Split valence polarized	GBASIS=CC2
cc-pVTZ	Triple split double polarized	GBASIS=CC3
aug-cc-pvDZ	augmented with diffuse functions	GBASIS=ACC2

	Effective core potentials	Good for treating heavy atoms
SBKJC	Split valence + core potential	GBASIS=SBKJC
Hay-Wadt	Split valence + core potential	GBASIS=HW

Complicated field, which is why the old standards live on in routine calculations. Optimal approach is to employ a composite model, calibrated by someone else, with well defined set of basis functions and treatments of exchance and correlation. A composite model pieces together results from a number of different calculations to estimate a higher accuracy model.

Electron cores

Low energy “core” electrons typically don’t participate in chemical bonding but can add substantially to computational cost. Generally seek approximations, especially for heavy elements/metals.

Heart of approach is to partition an atom into core and valence parts. Seek ways to express the influence of the core on the valence without actually having to compute the core. Essentially seek to write \[ \hat{v}^\text{ee} ≈ \hat{v}^\text{ee,core} + \hat{v}^\text{ee,val} \] where the core potential is some simpler, composite expression of the influence of the core on the valence. Typically take the electron cores as “frozen” in the pure atomic states, and express influence on valence through angular-momentum-dependent operators parameterized against accurate atomic calculations. Goal is to recover valence wavefunctions with less effort.

Relativistic effects

Relativistic kinetic energy is relativistic total energy minus the rest energy: \[T = \sqrt{p^2c^2 + m_0^2c^4} - m_0c^2\] Taylor expanding about $p^2 = 0$ gives the first-order mass-velocity correction: \[ T ≈ \frac{p^2}{2m_0} - \frac{p^4}{8 m_0^3 c^2}\] Reduces to non-relativistic result when $c→\infty$. Electrons near core move at speeds close to $c$, second term becomes non-negligible and diminishes their energy. Most important for $s$ states that penetrate closest to nucleus; they shield nucleus better and other valence states rise up in energy.

Electron spin and orbital magnetic moments also couple when $l>0$, leads to spin-orbit coupling that splits $p$, $d$, \ldots states into $j = l ± s$ states.

Darwin correction corrects s orbitals for electron and nucleus being at the same point; comes from solution of full Dirac relativistic equation for the atom.

Relativistic effects typically incorporated implicitly, by including in model for core electrons and thus capturing their effect on the valence. Spin-orbit, if necessary, added after the fact.

Implementations

Non-relativistic and relativistic effective core potentials (ECPs) available for many elements. These specify the potential felt by the valence electrons due to the core in terms of radial potential functions and angular projection operators. Typically these have to be combined with basis functions designed to work with them.

Most common are Hay-Wadt (LANL) and Stevens-Basch-Krause (SBK). Other more modern ones also available, like Stuttgart.

Essential to all plane-wave codes, like Vasp, but implemented differently. Will touch on later in class.

Population analysis

The molecular orbitals contain information that can be helpful in understanding structure and bonding:

Charge density most direct representation of electron distribution

\[ ρ(\bm{r}) = ∑_\text{occupied} |ψ_i(\bm{r})|^2 = ∑_μ,νP_μ,νφ_μ(\bm{r})φ_ν(\bm{r}) \]

Moments of charge density (dipole, quadrupole), useful for thinking about molecule-molecule interactions. (Only exactly defined for neutrals!)
Electrostatic potential, or Coulomb potential created by electrons and nuclei. More refined way of thinking about “how spots” on a molecule. Commonly used to parameterize classical forcefields, by seeking set of atom-centered charges that reproduce calculated electrostatic potential, as is done with CHELPG. Not uniquely defined.
Population analyses, which attempt to distribute electrons to individual atoms and possibly bonds based on decomposition of molecular orbitals. Chemically it is intuitively nice to assign charge to individual atoms. There is no single “right” way to do this\ldots the “charge” on an atom in a molecule in not uniquely defined! Consider an occupied molecular orbital ψ made up of two basis functions on two different atoms, α and β:

\begin{eqnarray} ψ &=& c_α φ_α + c_β φ_β
\langle ψ | ψ \rangle &=& c_α^2 + c_β^2 + 2 c_α c_β \langle\phi_α | φ_β\rangle \end{eqnarray} In Mulliken analysis, $c_α^2$ is fraction of ψ assignable to the atom of α, $c_β^2$ fraction assigned to atom of β. Remainder is the “overlap” population, which is split evenly between the two. Summing over all occupied orbitals and subtracting nuclear charges gives gross atomic charges.

In Löwdin analysis, basis functions are pre-orthogonalized, so last term vanishes.

Both approaches very sensitive to choice of basis set. Only sensible to compare within a common model type across molecules.

Localized orbitals is notion of creating linear combinations of ψ that satisfy some constraint for being compact. Leads to orbitals that are more naturally “bonding.”
Natural orbitals a rigorous scheme for orthogonalizing and assigning charge. Based on recognition that there is a set of orthogonal orbitals that optimally describe the density. Localizing these give natural bonding orbitals. See 06Weinhold.pdf.
Bader analysis another method, based on a geometric analysis of the total charge density. Define See Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, 1990.

Molecular orbital (MO) diagrams

Correlate molecular orbitals with their parent fragments. Use to be the thing. Seldom now

Implementation details of SCF methods

Basis is often orthonormalized to eliminate overlap from H-F-R equation; allows equations to be solved by matrix diagonalization.

Initial density matrix $\mathbf{P}$ are obtained by solving an approximate Hamiltonian (like extended Hückel). Always beware! Initial guess can influence final converged state.

Because the number of 2-electron integrals grows as $N^4$, they are sometimes calculated as needed “on-the-fly”, so-called direct SCF.

Hartree-Fock integrals can be computed analytically in a Gaussian basis. Any other choice of basis, or any DFT functional, requires integrals to be computed by quadrature. Used to be a real hang-up. Today, algorithms are very robust to establish grids and do quadrature.

SCF updating

The SCF procedure is an optimization problem: find set of coefficients that minimizes the total energy. As discussed above, success depends on a reasonable initial guess for density matrix and judicious updating. Various strategies can be used to speed and stabilize convergence, like damped mixing of previous cycles.

Second-order SCF is a convergence acceleration method that requires calculation or estimation of the first- and second-derivatives of the energy with respect to the orbital coefficients. See e.g. Chaban et al., Theor. Chem. Accts. 1997, 97, 88-95.

Pulay’s “direct inversion in the iterative subspace,” or “DIIS,” is a popular and powerful acceleration procedure that extrapolates from several previous Fock matrices to predict optimal next Fock to diagonalize. An opportunity for machine learning?

Controlled in GAMESS using the $SCF group.

$SCF DIRSCF=   .T./.F. controls direct scf
     SOSCF=    .T./.F. second-order scf
     DIIS=     .F./.T. direct inversion in the iterative subspace
     DAMP=     .T./.F. damping, on for initial iterations
$END

Can also control initial guess orbitals. Particularly powerful feature is to restart from converged orbitals from a previous calculations (.dat file), controlled with $GUESS group.

$GUESS GUESS = HUCKEL   construct intial guess from a simple Hamiltonian
             = MOREAD   read in orbitals from $VEC group.

\newpage

Potential energy surfaces

The potential energy surface (“PES”) is the sum of the repulsive energy of the nuclei and the kinetic and potential energies of all the electrons:

\begin{equation} E_\text{PES}(\mathbf{R_α},\mathbf{R_β},\ldots) =E_\text{elec} +∑_α=1^N ∑_{β =α +1}^N \frac{Z_α Z_β e^2}{R_{α β}} \end{equation}

Specifying atomic positions

./Images/Internals.pdf

Cartesian

Computationally straightforward but don’t correspond with our physical notion of bonds, bends, etc. Easiest to get out of a piece of software. A molecule has $3 N-6 $ internal degrees of freedom ($3N-5$ if linear), but Cartesians specify $3N$. The extra values correspond to the location of the center of mass and molecular oriendation. Codes will typically center and reorient the Cartesians.

In GAMESS, would specify Cartesian coordinates for \ce{FCH2CH2F} like this:

 $CONTRL COORD=CART $END
 $DATA
FCH2CH2F drag calculation
C1
C     6.0    -3.76764     0.33879     0.03727
C     6.0    -2.35246     0.34495     0.03689
F     9.0    -4.72277     0.58147    -1.18012
F     9.0    -1.59909    -0.68487    -0.83662
H     1.0    -4.04387     1.08375     0.75395
H     1.0    -3.92958    -0.71060     0.16941
H     1.0    -2.03786     0.18875     1.04760
H     1.0    -2.09983     1.28759    -0.40187
 $END

Internal coordinates

These provide a more intuitive representation and can be convenient when building molecules by hand. In codes like GAMESS, most commonly defined using “z-matrix” notation. Specify each atom in terms of its distance, angle, and dihedral angle with three previous atoms.

In Gamess, would specify z-matrix for \ce{FCH2CH2F} like this:

$CONTRL SCFTYP=RHF RUNTYP=ENERGY COORD=ZMT $END
$DATA
FCH2CH2F drag calculation
C1
C
C   1   r1
F   2   r2   1   A1
H   2   r3   1   A2   3   D1
H   2   r4   1   A3   3   D2
F   1   r2   2   A1   3   D3
H   1   r3   2   A2   6   D1
H   1   r4   2   A3   6   D2

r1=1.5386
r2=1.39462
r3=1.11456
r4=1.12
A1=109.54214
A2=111.
A3=110.
D1=120.
D2=-120.5
D3=50.
$END

Particularly convenient when you’d like to “scan” over the value of some coordinate. Variable can be applied to more than one independent coordinate, if the molecule has symmetry.

Features of potential energy surfaces

One-dimensional example

./Images/PES.pdf

Note 3-fold periodicity as expected for rotation about a \ce{C−C} single bond. Note too there are some special points:

Minima are places where energy bottoms out. More formally, first derivative of energy, or slope, or “gradient” $g = 0$, and second derivative, or curvature, or “Hessian” $H > 0$. These are the locally stable conformations of the molecule. Note that lowest energy in this case is not trans, but rather gauche conformations. Are you surprised?
Saddle points are places where energy is maximized. Physicially, corresponds to “transition states” connecting low-energy conformations. Gradient $g = 0$, but curvature $H < 0$.

Many-dimensional PES features

Gradient becomes vector and Hessian a matrix

\[ \bm{g} = \left ( \begin{array}{c} \frac{∂ E}{∂ q_1} \ \vdots \ \frac{∂ E}{∂ q_3N} \end{array} \right )\qquad \bm{H} = \left (\begin{array}{ccc} \frac{∂^2 E}{∂ q_1^2} & \cdots & \frac{∂^2 E}{∂ q_1∂ q_3N} \ \vdots & \ddots & \vdots \ \frac{∂^2 E}{∂ q_1∂ q_3N} & \cdots & \frac{∂^2 E}{∂ q_3N^2} \end{array} \right ) \]

gradient is vector tangent to PES. The force on an object is $\bm{F} = - \bm{g}$, so the gradients are often called the forces. Where gradient (slope) is negative, force is positive, and vice versa. Force always pushes system toward nearest minimum. If the potential is harmonic, then the force constant $k = H$, so the Hessian is also called the “force constant.”
Hessian matrix is real and symmetric. Diagonalization gives eignevalues and eigenvectors. Eigenvectors give “natural” directions along PES (physically, the harmonic vibrational modes), and eigenvalues indicate curvature in that direction.
Minimum on multidimensional PES has gradient vector $\bm{g} = 0$ and all positive Hessian eigenvalues.
First-order saddle point, or transition state, has $\bm{g} = 0$ and one and only one negative Hessian eigenvalues. (Physically, one unique direction that leads downhill in energy.) Must correspond to lowest-energy point connecting two minima.
Minimum energy pathway (MEP) or intrinsic reaction coordinate (IRC) is steepest descent pathway (in mass-weighted coordinates) from saddle point to nearby minima. Path a marble with infinite inertia would follow.
Higher order saddle points have $\bm{g} = 0$ and more than one negative Hessian eigenvalue. Can always lead to lower energy first order saddle point. These generally do not have chemical significance.

In computational chemistry/materials science, it is frequently our job to identify the critical points (minima and transition states on a PES). In liquids, PES much more flat and lightly corrugated. Statistical mechanics becomes mores important.

Each distinct electronic state defines its own PES. Remember that there are multiple PES’s for any atom configuration, corresponding to different electronic states. Sometimes these states can interact, intersect, giving avoided crossings, conical intersections. Lead to more complicated dynamical behavior.

Energy gradients and second derivatives

Gradients

\[ \frac{∂ E_\text{elec}}{∂ q_i} = \frac{∂}{∂ q_i} \langle Ψ | \hat{H} | Ψ \rangle = \langle \frac{∂ Ψ}{∂ q_i} | \hat{H} | Ψ \rangle + \langle Ψ |\frac{∂ \hat{H}}{∂ q_i} | Ψ \rangle + \langle Ψ | \hat{H} | \frac{∂ Ψ}{∂ q_i} Ψ \rangle \]

Hellman-Feynmann theorem says that sum of first and last terms vanish, so in principle only need to compute middle, which involves derivative of Coulomb potential and can be evaluated.
Pulay forces are forces from first and last terms that appear when basis functions are centered on atoms. An advantage of plane-wave basis sets, for which these terms vanish.

Hessian

In some electronic structure models can be computed analytically. More commonly, determined from numerical differentiation of gradients. Implementations typically assume that system is at minimum.

h must be small enough to stay in the harmonic region, but big enough to avoid numerical noise swamping the gradients.

For a molecule with N atoms, to construct complete $ 3N× 3N$ Hessian, have to evaluate gradients $6N$ times for two-sided differencing. Each pair of displacements completes one row of Hessian. Obviously tends to be quite expensive.

To get better precision and accuracy, could calculate more than two displacements, and could fit to a more complicated function than a harmonic potential.

Geometry optimization algorithms

Energy-only

Trudge infers gradient and locates minimum by wandering around

Energy + gradient

Steepest descent, just march down hill. $\bm{R}^′ = \bm{R} - λ \bm{g}$. Safe, very inefficient near minimum. May do line search to adapt λ.
Conjugate gradient is steepest descent plus orthogonalization to previous step. Safe, less very inefficient, common choice when far from minima.

Energy + gradient + higher order

Quasi-Newton Raphson takes advantages of both first and second derivative information: \[ \bm{R}^′ = \bm{R} - \bm{H}^-1(\bm{R})\bm{g}(\bm{R}) \]

Typically do not know Hessian and it is expensive to calculate. Make an initial guess, then update Hessian with gradient information from each geometry step. “Learning” PES as we go. Generally converges very rapidly near minima, where surface is not too anharmonic.

Rational function optimization is similar in spirit, also constructs Hessian, but uses more sophisticated (than quadratic) guess form of PES to update positions.
Direct inversion in the iterative subspace (DIIS) uses sizes of QNR steps as estimates of error and constructs new step from linear combination of previous that minimizes error inferred from previous steps: \[ \bm{err}(\bm{R}) = ∑_i c_i \bm{H}^-1_i \bm{g}_i \]

Generally very efficient in region of minimum. Algorithm can misbehave away from minima, possibly even converging to nearby saddle points, so often started with conjugate gradient steps.

Machine learning

All of these based on some form of assumed model of underlying PES (Taylor expand to second order around minimum). Emerging are methods that “fingerprint” structure and construct and improve energy model with each step. Remains to be seen if a new “standard” emerges.

Global optimizations

Simulated annealling, genetic algorithms, \ldots, more on the exotic side.

Convergence criteria

Typically determine convergence by enforcing maximum on each individual force component

GAMESS offers a limited set of algorithms.

$STATPT METHOD = NR     ! Quasi-Newton Raphson
               = RFO    ! rational function optimization
        OPTTOL = 0.0001 ! convergence criterion in au/bohr.
        HESS = GUESS    ! guess an initial Hessian
             = READ     ! read from $VIB group
$END

Specify `GAMESS` calculation type

Specified in $CONTRL group by RUNTYP flag:

Calculation	RUNTYP=
Single-point energy	Energy
Single-point energy + force	Gradient
Geometry optimization	Optimize
Linear scan over PES	Scan
Energy + second derivative	Hessian
Transition state search	Sadpoint
Intrinisc reaction coordinate	IRC

Note too that specifying EXETYP=CHECK will check your input without actually running the job.

Hessian (force) calculation can be done analytically or by numerical differentiation of forces, depending on electronic structure method:

$FORCE METHOD = ANALYTICAL
              = SEMINUM
       VIBSIZ = 0.01    ! step size (bohr)
$END

Efficient coordinate systems

Optimizations are most efficient in coordinate system that diagonalizes Hessian, so that optimizations along each direction are (nearly) independent. Large off-diagonal Hessian terms imply strong coupling. Cartesian coordinates do not reflect physical forces in system and are generally poor choice/slow convergence for optimizations. Can choose other coordinate systems. Forces/Hessian always computed in cartesians, so any other coordinate system requires transformations back and forth.

Cartesians simplest to implement, consistent performance. Typically only choice in supercell calculations.
Z-matrix are easy to define and use for organic molecules with no rings. For $3N-6$ degrees of freedom, must specify $N-1$ distances, $N-2$ angles, and $N-3$ dihedrals. Typically will converge much faster than Cartesians for small molecules. Implemented in most molecular codes, including GAMESS.

$CONTRL COORD=ZMAT NZVAR = 3N-6 $END

Natural internals are generalizations of z-matrix that are coordinates that approximately diagonalize Hessian. Hard to generate a priori or automatically. Specification in GAMESS is arcane, done in $ZMAT.
Redundant internals are an over-determined set of internal coordinates, like all bond distance, all angles between bonded atoms, all dihedrals. Easily constructed automatically. Mapping between Cartesians and redundant internals becomes more complicated; transform from redundant to Cartesians is overdetermined and has to be solved iteratively. Works very efficiently for molecules though. Implemented (in spirit) in GAMESS.

$CONTRL NZVAR = 3N-6 $END
$ZMAT DLC=.TRUE.  AUTO=.TRUE. $END

	FCH2CH2F	C5H10
Cartesian	13	30
Z-matrix	11	failed
Redundant/DLC	16	failed

Performance of models for geometries

Generally pretty good! Gross geometries of molecules can be computed with good reliability using common approximations (H-F, LDA, GGA). Subtler things (does \ce{FCH2CH2F} really prefer trans or gauche?) can take more care in choice of electronic structure model. https://cccbdb.nist.gov/ is a great place to look for benchmarks.

Vibrational frequencies

Suppose we are at a minimum on a PES. Near that minimum, we can Taylor expand the PES and truncate at second order (“harmonic”, or quadratic, approximation). The Schrödinder equation describing the dynamical motion of the nuclei can be written in terms of displacements from the equilibrium position, $q_i = x_i - x^\text{eq}_i$, $i =1, \ldots, 3N$ and the Hessian $H$: \[\hat{H} = -\frac{\hbar^2}{2}∑_i \frac{1}{m_i}\frac{∂^2}{dq_i^2}+\frac{1}{2}∑_i,j H_ij q_i q_j\] or in mass-weighted coordinates, $ξ_i=\sqrt{m_i}q_i$: \[ \hat{H} = -\frac{\hbar^2}{2}∑_i \frac{∂^2}{dξ_i^2}+\frac{1}{2}∑_i,j ˜{H}_ij ξ_i ξ_j,\qquad ˜{H}_ij=\frac{1}{\sqrt{m_i m_j}}H_ij \] From eigenvalues $κ_i$ and eignevectors $s_i$ (“normal modes”) of mass-weighted Hessian, can transform into $3N$ one-dimensional problems: \[\hat{H}_i=-\frac{\hbar^2}{2}\frac{d^2}{ds_i^2}+κ_is_i^2\] This is one-dimensional harmonic oscillator Hamiltonian, solutions well known.

./Images/HO.pdf

Do this in $3N$ Cartesian space, so 6 (or 5) of the normal modes correspond to translations and rotations of the molecule. If calculation is exact, these will have $κ_i = 0$. Numerical errors may make them somewhat non-zero. If necessary, these can be projected out by transforming Hessian to internal and back to Cartesian coordinates.

Note it is impossible for molecule to just sit at $q = 0$. Nuclei are always vibrating about $x^\text{eq}$. Gives zero point energy \[ZPE = \frac{1}{2}∑_i hν_i\]

For a chemical bond, $κ\approx \SI{500}{N.m^-1}$, $hν ≈ \SI{0.8}{eV} $.

Absorption intensities

Vibrational modes can be probed/observed spectroscopically. Intensity of stimulated absorption from vibrational state i to f given by Einstein coefficient of stimulated absorption: \[B_if=\frac{|μ_if|^2}{6ε_0\hbar^2}\] Arises from coupling of electric field of light with dipole of system. Transition dipole moment, $μ_if)\, given by: \[μ_if = \langle ψ_i|\hat\mu |ψ_f \rangle\] where \(\hat\mu$ is dipole operator, \[\hat\mu = ∑ q_i \bm{r}_i \] where sum runs over all charged particles and $q_i$ is each charge. Dipole moment changes as molecule vibrates. In direction $ξ$, can write \[\hat\mu(ξ(t)) = μ(0) +ξ(t)\left( \frac{dμ}{dξ}\right )_ξ=0 + \ldots\] Use harmonic oscillator wavefunctions for $ψ_i$ and $ψ_f$: \[μ_if = \langle ψ_i|\hat\mu |ψ_f \rangle = \langle ψ_i|μ(0) |ψ_f \rangle + \left( \frac{dμ}{dξ}\right )_ξ=0 \langle ψ_i | ξ|ψ_f\rangle + \ldots\] First integral vanishes for $i≠ f$. Second integral provides gross selection rule that intensity of transition proportional to the dynamic dipole moment along vibrational normal mode $ξ$ and particular selection rule that intensity of transition is zero unless $f = i ± 1$. Latter comes from nature of Hermite polynomials. At normal temperatures, $i = 0$, and the only observable vibrational transitions are $1→ 2$.

Performance of harmonic approximation for vibrational spectrum

Harmonic vibrational frequency systematically overestimate experiment. Convolution of harmonic approximation error: actual PES is not exactly harmonic, and errors intrinsic to electronic structure model. Errors are typically systematic (see eg https://doi.org/10.1021/jp960976r).

	scale factor
HF/3-21G	0.9085
HF/6-31G(d)	0.8929
MP2/6-31G(d)	0.9434
B3LYP/6-31G(d)	0.9613

Relative intensities are generally predicted with good reliability. Model predicts absorption peaks to be delta functions. Peaks always broadened due to a variety of fundamental and instrumental considerations. Common in displaying spectra to arbitrarily broaden using a Lorentzian function.

Transition states

Symmetry

Can exploit symmetry to force calculation to converge to a transition state. Eclipsed form of \ce{FCH2-CH2F} has mirror symmetry ($C_s$), if initialized in that configuration, optimization must preserve symmetry.

 $DATA
  FCH2CH2F eclipsed z-matrix (Cs, or mirror, symmetry)
CS

 C   
 C      1   1.5000000
 F      2   1.4000000  1   109.5421400
 H      2   1.1000000  1   109.5421400  3   120.0000000  0
 H      2   1.1000000  1   109.5421400  3  -120.0000000  0
 F      1   1.5000000  2   109.5421400  3     0.0000000  0
 H      1   1.1000000  2   109.5421400  6   120.0000000  0
 H      1   1.1000000  2   109.5421400  6  -120.0000000  0
 $END

See results in ../Labs/Gamess/FCH2CH2F/SYMMETRY.

Coordinate drag

If reaction state coordinate maps closely onto some internal coordinate, then can do a series of “constrained” optimizations, fixing the coordinate of interest to a series of values and relaxing all other coordinates. Use IFREEZ within GAMESS. For example, define \ce{FCH2CH2F} using z-matrix and freeze \ce{F-C-C-F} dihedral angle at a series of values.

 $CONTRL SCFTYP=RHF DFTTYP=PBE RUNTYP=OPTIMIZE COORD=ZMT NZVAR=18 ISPHER=1 $END
 $BASIS GBASIS=PCseg-1 $END
 $STATPT IFREEZ(1)=12 $END
 $DATA
FCH2CH2F eclipsed TS
C1
C
C   1   r1
F   2   r2   1   A1
H   2   r3   1   A2   3   D1
H   2   r4   1   A3   3   D2
F   1   r5   2   A4   3   D3
H   1   r6   2   A5   6   D4
H   1   r7   2   A6   6   D5

...
D4=120.
...
 $END

See results in ../Labs/Gamess/FCH2CH2F/SCAN.

A quasi-NR optimization started at one of the approximate TS’s will usually converge to the exact TS.

 $CONTRL SCFTYP=RHF DFTTYP=PBE RUNTYP=OPTIMIZE COORD=ZMT NZVAR=18 ISPHER=1 $END
 $STATPT METHOD=NR $END
 $BASIS GBASIS=PCseg-1 $END
 $DATA
FCH2CH2F zmatrix optimization near saddle point, no hessian
C1
 C   
 C      1   1.5256029
 F      2   1.4018052  1   110.3919125
 H      2   1.1066523  1   112.5759655  3   119.5903058  0
 H      2   1.1070172  1   108.5521191  3  -119.1579857  0
 F      1   1.4018888  2   110.3079213  3   120.0000000  0
 H      1   1.1065092  2   112.5372024  6   119.5247291  0
 H      1   1.1069911  2   108.5976341  6  -119.1645584  0
 $END

See results in ../Labs/Gamess/FCH2CH2F/QNR-TS. Always good practice to follow up with a frequency calculation to confirm.

Coordinate dragging can fail when the reaction coordinate is non-linearly related to the multiple internal coordinates. Plus, it is relatively expensive, as it involves a lot of optimizations.

1st order methods (NEB, …)

A scan is a first order method that freezes the gradient along a search direction. A better algorithm would use the information along the search direction to find the TS and the entire MEP. That’s the spirit of a nudged elastic band calculation, for instance. Also of newer methods to generate more sophisticated approximations. Touch on that later....

2nd order methods

Hessian-based methods (like quasi-NR and DIIS) generally work more efficiently, assuming you can find a region reasonably close to the TS and can get a good guess of the Hessian with the appropriate (1!) number of negative eigenvalues. These work like regular optimization methods, but search uphill along the negative eigenvector and downhill along the other directions. The Hessian update scheme has to be appropriately modified to accommodate the negative eigenvalue. The choice of coordinate system can become even more crucial.

Algorithm:

Identify and optimize structures of reactant and product states
Guess a structure for the TS, interpolated either mathematically or empirically between reactant and product
Compute initial TS Hessian (RUNTYP = HESSIAN). Only needs to be approximately correct, so appropriate to use a less-expensive method.
Confirm Hessian has desired number of imaginary modes, or at least an imaginary mode that maps well onto desired reaction coordinate.
Do saddle-point search with guessed Hessian (RUNTYP=SADPOINT, HESS=READ).
Check result with final Hessian calculation

 $CONTRL SCFTYP=RHF DFTTYP=PBE RUNTYP=SADPOINT COORD=CART ISPHER=1 $END
 $BASIS GBASIS=PCseg-1 $END
 $FORCE HESS=READ $END
 $DATA  
H2ON2                                                                           
C1       1
N           7.0      0.9405942384     -0.3241591975     -0.0156283505
N           7.0      0.0107640781      0.6540868381      0.0152029180
O           8.0     -0.9317086136     -0.2406652788     -0.0189127039
H           1.0     -0.2853795803     -0.9336844866      0.3015543918
H           1.0      1.8537652825      0.1690870887      0.0045149723
 $END      
 $HESS
ENERGY IS     -184.6737814637 E(NUC) IS       74.1740447013
 1  1 8.31017204E-01 1.46125961E-01 7.06555414E-04-1.42105871E-01 2.88851646E-02
 1  2-1.06097190E-03-3.64083936E-01 2.22079857E-02-2.60487424E-03-3.19838288E-02
.....
 $END

See Results in ../Labs/Gamess/H2NNO-TS.

Well climbing

Algorithms exist that will climb out of basins and seek transition states....

Intrinsic reaction coordinates

In more complicated systems, it can often be difficult to know exactly what minima a particular transition state corresponds to. The intrinsic reaction coordinate (IRC) or minimum energy path (MEP) is steepestdescent path from TS towards both basins. Starts from TS, steps forward in direction given by gradient, using second order method. Have to use care in selecting algorithm and step sizes to stay on the path. Can be useful for locating variational transition state…configuration where free energy (rather than energy) is maximized.

 $CONTRL SCFTYP=RHF DFTTYP=PBE RUNTYP=IRC COORD=UNIQUE ISPHER=1 $END
 $BASIS GBASIS=PCseg-1 $END
 $IRC PACE=LINEAR STABLZ=.TRUE. NPOINT=10 SADDLE=.T. FORWRD=.T. $END
 $FORCE HESS=READ $END
 $DATA
H2ON2                                                                           
C1       1
 N           7.0   0.9233485039  -0.2399318160   0.0736256617
 N           7.0  -0.0272197420   0.6088503884  -0.0717581187
 O           8.0  -1.0933769256  -0.1081391759   0.0571959479
 H           1.0  -0.1059622027  -1.0519964547   0.2191639975
 H           1.0   1.8912457716   0.1158820221   0.0085037392
 $END      
 $HESS
ENERGY IS     -185.5854194050 E(NUC) IS       73.0033764178
 1  1 8.43287379E-01-4.86363743E-02 1.07610516E-04-3.07227749E-01 1.52208588E-01
.....

Molecular dynamics

Basic idea is to propogate atoms forward in time using some model to compute the forces on atoms and Newton’s laws to describe kinetic energy of atoms. May be done at constant energy (NVE) or, by coupling kinetic energy to an appropriate reservoir, at constant temperature (NVT). Details beyond this course (perhaps beyond this instructor!). Primary point for us is that cost of electronic structure calculations is great enough that typically some adjustments must be made to make force calculations cheap and fast enough to support number of evaluations necessary to do meaningful dynamics.

Cheap parameters

Simplest trick is to back off on the precision and accuracy of calculations to reduce force evaluation cost.

Car-Parrinello dynamics

Treat wavefunction itself as a dynamical variable and propogate forward in time along with the nuclei. Ideally parallels but does not exactly follow Born-Oppenheimer energy surface. Huge advance when originally introduced, has fallen out of favor as more conventional B-O approaches compete more effectively.

Biased sampling

Ask Whitmer. See https://doi.org/10.1021/acs.jctc.8b00192. \newpage

WFT Beyond Hartree-Fock

DFT or WFT calculations of PES’s are expensive. Not uncommon to use an heirarcy of methods:

Characterize PES using some low-cost model (LDA)
Improve structures using some more reliable and expensive method
Improve final energies using some even more reliable and expensive method

WFT has the advantage of offering a series of systematically improving models. Starts with Hartree-Fock ground-state wavefunction and build in “many-body” features/electron correlation by combining with “excited” Hartree-Fock wavefunctions:

./Images/correlation-summary.pdf

static correlation are electron correlation effects that arise from the restrictive form of the H-F wavefunction
dynamic correlation is like it sounds, the “dance” of all electrons about one another
configuration interaction is a variational approach to adding in dynamical correlation by combining H-F determinants. Difficult to apply, rarely used any more.
size consistency is the property that the energy model scales properly with the system size. CI lacks this.
perturbation theory (MPn) is non-variational, size-consistent, and user-friendly approach
coupled-cluster (CCxxx) is a systematic, size-consistent approach to introducing H-F correlation. Improvement on perturbation theory.
quadratic configuration interaction (QCIxxx) is an approximate coupled cluster model.

Heirarcy of models:

HF < MP2 ~ MP3 ~ CCD < CCSD < MP4 < CCSD(T) ~

Reliability also tied to basis set completeness.

“Model chemistry” some linear combination of these, often calibrated against experimental data. “G2” most venerable:

	Hartree-Fock	MP2	MP4	QCISD(T)
6-31G(d)	ZPE	Structure
6-311G(d)		1	2	3
6-311+G(d,p)		4	5
6-311G(2df,p)		6	7
6-311+G(3df,p)		8

\[ \text{QCISD(T)/6-311+G(3df,p)} ≈ 2 + (3-2) + (5-2) + (7-2) + (8-1) - (4-1) - (6-1)\] \[ Δ E(\text{HLC}) = \SI{-2.50}{mHa} * (\text{# electron pairs}) - \SI{0.19}{mHa}* (\text{# unpaired electrons}) \] \[ E(\text{G2}) = \text{QCISD(T)/6-311+G(3df,p)} + 0.8929 * \text{ZPE} + Δ E(\text{HLC})\]

Quality assessed using, eg, mean absolute deviation from some reliable data set. Many elaborations on this same idea in the literature.

\newpage

First-principles thermodynamics

Connection Between QM and Thermodynamics

Internal energy

The internal electronic energy of a single molecule from a WFT/DFT calculation is the energy associated with taking infinitely separated constituent nuclei and electrons at rest and forming a molecule: \begin{equation} 2~\mathrm{H}^+ + 8~\mathrm{O}⁸⁺ + 10~\mathrm{e}^- → \mathrm{H_2O}\qquad E^\mathrm{elec} \end{equation} Calculate $E^\mathrm{elec}$ within the Born-Oppenheimer approximation, so the nuclei are fixed in space at the minimum energy configuration.

Zero-point vibrational energy

\begin{equation} E^0=E^\mathrm{elec} + ZPVE \end{equation} ZPVE can be calculated reliably within the harmonic approximation, according to \begin{equation} \mathrm{ZPVE}=\frac{1}{2}h∑_i=1^3n-6ν_i \end{equation} where $ν_i$ are the harmonic vibrational frequencies, obtained from a vibrational frequency analysis. $E^0$ is the minimum physically meaninfful energy of the molecule.

Finite T energy

Energy can be deposited in a syste as translational and rotational kinetic energy, in excited vibrational modes, in the interaction of a molecule with an external electric or magnetic or gravitational field, or .... If we assume that the energy in these various degrees of freedom are separable, we can write: \begin{equation} E_i=E^0+E^\mathrm{trans}+E^\mathrm{rot}+E^\mathrm{vib} +E^\mathrm{elec*}+E^\mathrm{ext} \end{equation} To fully describe microscopic energetic state of a system, would have to specify all of these.

Typically want collective properties at equilibrium, like the internal energy $U$ or enthalpy $H$ or Gibbs energy $G$, under some external constraints like temperature $T$ or volume $V$. These thermodynamic quantities are averages over the energy states of an ensemble of molecules. The way this averaging is performed is the realm of statistical thermodynamics.

Canonical ensemble

Free variables are the number of molecules $N$, the total volume $V$, and the temperature $T$. Probability for a molecule to be in some energy state $E_i$ above $E^0$ is given by the Boltzmann factor, \begin{equation} P(E_i) \propto e^-E_iβ=e^-E_i/k_BT,\qquad\beta=1/k_BT \end{equation} Defines an exponentially decaying probability function for a state to be occupied at some temperature. Temperature is the characteristic of a system following this distribution.

./Images/boltzmann.pdf

Averages and partition functions

Let’s use this to calculate the internal energy $U$ of a molecule at some temperature. \begin{equation} U(T)=\frac{∑_iE_iP(E_i)}{∑_iP(E_i)} \end{equation} where the denominator ensures that the probability is normalized. \begin{eqnarray} U(T) & =& \frac{∑_iE_i e^-E_iβ}{∑_ie^-E_iβ}
& = & \frac{\frac{∂}{∂\beta}∑_ie^-E_iβ}{∑_ie^{-E_i β}}\ & = & -\frac{∂ ln ∑_i e^-E_iβ}{∂ β} \end{eqnarray} The sum over energy states is evidently a special quantity, called the partition function: \begin{equation} Q=∑_ie^-E_iβ \end{equation} All thermodynamic quantities can be written in terms of the partition function.

Stat Mech applied to stuff

./Images/gls.pdf

Separability

in principle need to sum over all the types of energy states (translational, rotational, vibrational, …) of every molecule. Seemingly impossible task. One simplification is if we can write energy as sum of energies of individual elements (molecules) of system: \begin{align} E_j&=ε_j(1)+ε_j(2) + … + ε_j(N)
Q(N,V,T) &= ∑_j e^-E_jβ \ &=∑_je^{-(ε_j(1)+ε_j(2) + … + ε_j(N))β} \end{align} If molecules/elements of system can be distinguished from each other (like atoms in a fixed lattice), expression can be factored: \begin{align} Q(N,V,T)&=\left ( ∑_j e^-ε_j(1)β\right )\cdots \left ( ∑_j e^-ε_j(N)β\right ) \ &= q(1)\cdots q(N) \ \text{Assuming all the elements are the same:}\ &= q^N \end{align} If not distinguishable (like molecules in a liquid or gas, or electrons in a solid), problem is difficult, because identical arrangements of energy amongst elements should only be counted once. Approximate solution, good almost all the time: \begin{equation} Q(N,V,T)=q^N/N! \end{equation} Sidebar: “Correct” factoring depends on whether individual elements are fermions or bosons, leads to funny things like superconductivity and superfluidity.

This $q(V,T)$ is the molecular partition function, and is calculated by summing over the individual energy states of a single molecule (starting at $E_0$).

Further simplified by factoring into contributions from various ($3N$) molecular degrees of freedom: \begin{eqnarray} q(V,T)&=&\left(∑_\mathrm{trans} e^{-e_\mathrm{trans}β}\right) \left(∑_\mathrm{rot} e^{-e_\mathrm{rot}β}\right) \left( ∑_\mathrm{vib} e^{-e_\mathrm{vib}β} \right) \left( ∑_\mathrm{elec} e^{-e_\mathrm{elec}β}\right)
&=& q_\mathrm{trans}q_\mathrm{rot}q_\mathrm{vib}q_\mathrm{elec} \ U & = & E_0 + U_\mathrm{trans}+U_\mathrm{rot}+U_\mathrm{vib}+U_\mathrm{elec} \end{eqnarray} Similarly for other thermodynamic quantities, for example, \begin{equation} C_v=\left(\frac{∂ U}{∂ T}\right)_V = C_{v,\mathrm{trans}}+C_{v,\mathrm{rot}}+C_{v,\mathrm{vib}}+C_{v,\mathrm{elec}} \end{equation} Thermodynamic quantities are sums of contributions from indvidual degrees of freedom.

Have to somehow model these motions and have to use our quantum mechanical results to parameterize the models.

Ideal gas

Ideal gas of molecules

Assume molecules are indistinguishable and that internal energy is seperable \[ Q_ig(N,V,T) = \frac{(q_\mathrm{trans}q_\mathrm{rot}q_\mathrm{vib}q_\mathrm{elec})^N}{N!} \] \[ F(N,V,T) = F_\mathrm{trans}+ F_\mathrm{rot} + F_\mathrm{vib}+F_\mathrm{elec}\] for any thermodynamic function $F$.

Electronic partition functions $→$ spin multiplicity

Governed by Fermi-Dirac distribution. Electronic degeneracy at normal T.

Vibrational thermodynamics: harmonic oscillator

Energy spectrum above ZPE given \begin{equation} E_v=hvν,\qquad v=0,1,2,… \end{equation} Define characteristic temperature $Θ =h ν/k_B$. \begin{eqnarray} q(T) & = &∑_v=0^∞ e^-v(Θ/T)
& = & \frac{1}{1-e^-Θ/T} \end{eqnarray} (geometric series). Partition function increasing function of T. Look at limits.

\noindent Internal energy: \begin{eqnarray} U(T) &=&-\frac{∂ ln q}{∂ β}
& = & \frac{ε_0}{e^ε_0β-1} \end{eqnarray} Need vibrational spectrum to get these contributions.

Results only as good as H-O model! For low frequency modes, errors can be substantial, esp for entropy. Can apply more sophisticated models.

Rotational thermodynamics: rigid rotor

Characteristic temperature $Θ_\mathrm{rot} = \hbar^2/2 I k_B$. Moment of intertia I only depends on shape or molecule—geometry optimzation.

“High” T $q_\mathrm{rot}(T) ≈ σ Θ_\mathrm{rot}/T$, most often true

Translational states: particle-in-a-box

\[ E_n=\frac{n^2π^2\hbar^2}{2 m L^2} = ε_0n^2,\qquad n=1,2,3,\ldots \] $Θ_\mathrm{trans} = ε_0/k_B$ is typically tiny, allows partition function sum to be approximated by integral: \begin{eqnarray*} q_\mathrm{trans,1D} ≈ ∫_0^∞ e^{-x^2β\epsilon_0}dx = L/Λ
Λ = \frac{h}{\sqrt{2π m k_B T}} \qquad \mathrm{thermal\ wavelength} \ q_\mathrm{trans,3D} = V/Λ^3 \end{eqnarray*}

Thermal wavelength $Λ$ depends only a molecule mass and is of the order the box dimensions at which quantization is evident. Typically a tiny number (eg \SI{1.7e-11}{m} for Ar in a \SI{1}{\liter} volume at \SI{298}{K}. $q_\mathrm{trans}$ is thus enormous: lots of translational freedom. $q$ depends on volume, introduces volume/concentration/pressure dependence into thermo functions. Conventional to define a standard state $V^ˆ$ volume, or corresponding pressure.

	Energy (cm^-1)	Temperature (K)
	Characteristic	Characteristic	States @ RT
translational	$\hbar^2/2 m L^2 ≈ 10^-21$	$10^-21$	$10³⁰$	classical limit

rotational	$≈ 1$	$≈ 1$	100’s	semi-classical

vibrational	$≈ 1000$	$≈ 1000$	1	non-classical

electronic	$≈ 10,000$	$≈ 10,000$	1	non-classical

Chemical reactions and equilibria

Chemical reaction

General chemical reaction $∑_i ν_i A_i = 0$, $ν_i$ stoichiometric coefficients
Thermodynamic change $Δ W^ˆ(T) = ∑_i ν_i W^ˆ_i(T)$, where $W = A, U, S, G, \ldots $
“Standard state” derives from concentration dependence of entropy
“Standard state” corresponds to some standard choice, $(N/V)^ˆ = c^ˆ$, e.g. \SI{1}{mol/l} (T-independent), or $(N/V)^ˆ = P^ˆ/RT$, e.g. \SI{1}{bar} (T-dependent)
Permits functions to be easily computed at other concentrations, e.g. \begin{displaymath} A(T,N/V) = A^ˆ(T) + k T ln\left ( (N/V)/(N/V)^ˆ \right ) =A^ˆ(T) + k T ln \left ( c/c^ˆ \right ) \end{displaymath}
Example: ethane dehydrogenation, \ce{C2H6 -> C2H4 + H2}, 1 bar standard state
Reaction entropy captures contributions of all degrees of freedom
Reaction energy (internal, Helmholtz, …) must also capture difference in \SI{0}{K} electronic energy \begin{displaymath} Δ U^ˆ (T) = U^ˆ_\mathrm{B}(T)-U^ˆ_\mathrm{A}(T)+Δ E(0) \end{displaymath}

Chemical equilibrium

At chemical equilbrium, total free energy minimized with respect to reaction advancement $ξ$ \[G(T,ξ) = ξ (Δ G^ˆ + k T ∑_i ν_i ln P_i/P^ˆ) \]
Equilibrium condition—equate chemical potentials \begin{eqnarray*} μ_A(N,V,T) & = & μ_B(N,V,T)
E_A(0) - k T ln (q_A/N_A) & = & E_B(0) - k T ln (q_B/N_B) \ \frac{N_B}{N_A} = \frac{N_B/V}{N_A/V} & = &\frac{q_B(T,V)/V}{q_A(T,V)/V} e^{-Δ E(0)/kT} \end{eqnarray*}
$q/V = 1/Λ^3$ has units of number/volume, or concentration
Equilibrium constant—convert units to some standard concentration $c^ˆ$ or pressure $P^ˆ$ \begin{eqnarray*} q_A^ˆ(T) & = & (q_A(T,V)/V) (1/c^ˆ)
q_A^ˆ(T) & = & (q_A(T,V)/V)(RT/P^ˆ) \ K_eq(T) & = &\frac{q_B^ˆ(T)}{q_A^ˆ(T)} e^{-Δ E(0)/kT} = e^{-Δ G^ˆ(T)/kT} \end{eqnarray*}

Reaction rates

Rate: number per unit time per unit something (volume, area, …). Reactions hypothesized to occur through sequence of elementary steps.

For example, ozone decomposition, rate second-order at high $P_\ce{O2}$, first-order at low $P_\ce{O2}$

\ce{2 O3 -> 3 O2}
\ce{O3 ->[k_1] O2 + O}
\ce{O2 + O ->[k_-1] O3}
\ce{O + O3 ->[k_2] 2 O2}

See two compute rates of individual steps.

Transition state theory (TST)

Assumptions
1. Existence of reaction coordinate (PES)
2. Existence of dividing surface
3. Equilibrium between reactants and “transition state”
4. Harmonic approximation for transition state
rate proportional to concentration of “activated complex” over reactants times crossing frequency \begin{eqnarray*} r & = & k C_AC_B
& = & k^\ddagger C_AB^\ddagger \ & = & ν^\ddagger K^\ddagger C_A C_ B \ & = & ν^\ddagger \frac{k_BT}{hν^\ddagger}\bar{K}^\ddagger(T) C_A C_B \ & = & \frac{k_B T}{h} \frac{q^\ddagger(T)}{q_A(T) q_B(T)} e^{-{Δ E(0)/k_BT}} C_A C_B \end{eqnarray*}
application to two molecules - vinyl alcohol to acetaldehyde
microscopic reversibility
equilibrium requirement $K_eq(T) = k_f(T)/k_r(T)$

Thermodynamic connection

Relate activated complex equilibrium constant to activation free energy \[ \(\bar{K}^\ddagger(T) = e^{-Δ G^{ˆ \ddagger}(T)/kT} = e^{-Δ H^{ˆ \ddagger}(T)/k_BT}e^{Δ S^{ˆ \ddagger}(T)/k_B} \]
Compare to Arrhenius expression \[E_a = Δ H^{ˆ \ddagger}(T) + kT, A = \frac{k_B T}{h}e^1e^{Δ S^{ˆ \ddagger}(T)/k_B}\]

\newpage

Plane waves and core potentials

Periodic boundary conditions

Free particle moving in one dimension \[φ_G(x) = e^{i G x} \] G can take any value. Eigenfunction of momentum and kinetic energy operators. In particular, \[E_G = \frac{\hbar^2}{2m}G^2\]

Suppose we impose periodic boundary conditions, so that $φ_G(x) = φ_G(x+na), n=± 1, \ldots$. a is a lattice vector. \[e^{i G x} = e^{i G (x + a)} → e^{i G a} = 1\] Places constraints on G: \[ G = \frac{2π}{a} n, n = ± 1, \ldots\] $2π/a$ is a reciprocal lattice vector.

Normalized on domain a, \[φ_G(x) = \frac{1}{\sqrt{a}}e^{i G x}\]

Properties of these plane wave functions:

Periodic, $λ_G = a/n$
Orthonormal, $\langle φ_G | φ_G^′\rangle = δ_G,G^′ $
Definite kinetic energy
Complete \begin{eqnarray} ψ(x) &=& ∑_G φ_G(x) \langle φ_G | ψ \rangle
& = & ∑_G φ_G(x) ψ(G) \end{eqnarray}

$ψ(G)$ is Fourier transform of $ψ(x)$. Both contain exactly the same information. Can take advantage of Fourier transform machinery to transform between two representations. For general function $ψ(x)$, represent on some grid of size N. Size of grid determined by maximum frequency/minimum wavelength included. $Nlog(N)$ operations to transform back and forth.

Integrals can be evaluated in either x or G space: \begin{eqnarray} I & = & ∫_a A^*(x) B(x) d x
& = & ∑_G,G^′ A^*(G) B(G) ∫ e^-iGxe^{iG^′ x} dx \ & = & ∑_G,G^′ A^*(G) B(G) δ_G,G^′ \ & = & a ∑_G A^*(G) B(G) \end{eqnarray} Can leverage to simplify solving any problems that are periodic.

Example

Hydrogen atom in a box

Suppose we place an atom at some point R in each box. In principle, can represent solutions as linear combinations of plane waves. Plane waves become basis functions.

Does where we put the atom matter to the answer? No. Alway related by structure factor.

In principle, would have to let G run to infinity to have complete basis. In practice, would have to cut off at some point, defined by \[ \frac{\hbar^2}{2 m_e} G^2 < \frac{\hbar^2}{2 m_e} G_\text{cut}^2 = E_\text{cutoff} \] Sets minimum wavelength retained in basis, and thus size of basis. For fixed $E_\text{cutoff}$, larger a implies larger basis. Basis size has to be a whole number, so changes discontinuously with $E_\text{cutoff}$.

Recall our basic equation \begin{equation} \left \{ \hat{h} +v_{\text{Coulomb}}[ρ] + v_\text{exchange}[ψ_i] + v_\text{correlation}[ψ_i]\right\}ψ_i(\mathbf{r}) =ε_i ψ_i(\mathbf{r}) \end{equation} Contains kinetic and potential energy terms.

Kinetic energy is diagonal in plane wave basis…easy! \[ \langle φ_G | -\frac{1}{2}\frac{d^2}{dx^2} | φ_G^′ \rangle = \frac{1}{2}G^2 δ_G,G^′ \]

Potential energy terms can take advantage of Fourier transforms. \begin{eqnarray} v(x) & = & ∑_G^{”} v(G”) e^{i G” x}
\langle φ_G | v(x) | φ_G’\rangle & = & ∑_G^{”} v(G”) \langle φ_G | φ_G” | φ_G’ \rangle \ & = & v(G’-G) \end{eqnarray} How many Fourier components to include in the sums? Turns out for a basis of dimension m you need 2m components to specify the potential exactly, but you can generally get away with smaller. The cost and accuracy of the calculation scale with this choice.

Solutions given by $ψ_i(x) = ∑_G c_i(G) φ_G(x) $.

One key thing we would need to evaluate potentials is electron density, $ρ(x)$. Let $f_i$ be occupancy of each orbital. \begin{eqnarray} ρ(x) & = & ∑_i f_i |ψ_i(x)|^2
& = & \frac{1}{a}∑_i f_i ∑_G,G’ c_i^*(G)c_i(G’) e^i(G-G’)x \ & = & ∑_{-2 G_\mathrm{max}}^{2 G_\mathrm{max}} n(G) e^{i G x} \end{eqnarray} Electron density exactly represented in plane wave basis with cutoff four times that of the basis set cutoff.

Three-dimensional periodicity

Same ideas generalize to three dimensional periodicity.

Lattice constants of periodic cell become three vectors, $\bm{a}_1$, $\bm{a}_2$, $\bm{a}_3$. As we’ll discuss a little bit later, these vectors will be from one of the Bravais lattices. Every point \bm{R} is equivalent to any other point $\bm{R}’=\bm{R} + n_1 \bm{a}_1+ n_2 \bm{a}_2+ n_3 \bm{a}_3$.

Periodic box matrix $\bm{h} = [\bm{a}_1, \bm{a}_2,\bm{a}_3]$. Periodic box volume $Ω = \text{det} \bm{h}$.

Reciprocal lattice defined such that $\bm{b}_i⋅\bm{a}_j = 2 π δ_ij$. Given by \[2π(\bm{h}^t)^-1 = [\bm{b}_1,\bm{b}_2,\bm{b}_3]\]

Reciprocal lattice vectors become $\bm{G} = i \bm{b}_1 + j \bm{b}_2 + k \bm{b}_3$. Basis functions become $φ(\bm{r}) = \frac{1}{Ω}e^{i \bm{G}⋅\bm{r}} $.

Locations of atoms $\bm{R}_α$ within this periodic cell can be specified in Cartesians or as fractional coordinates of the lattice vectors. Related by $\bm{R}_α = \bm{h}\bm{f}$.

In general, can use this periodic representation to describe atoms, molecules, whatever. If we want to describe isolated things, box dimensions must be large enough to avoid artificial interactions between periodic images.

In particular, Coulomb potentials decay like $1/r$ and are thus long-ranged. Have to use special tricks (Ewald sums) to evaluate these sums, and have to carefully group electrostatic terms to avoid non-convergent sums.

Key limitations:

Periodic cell must be charge-neutral. The electrostatic energy of an infinite, charged system diverges.
Supercell must not have a net electric field.
The absolute electrostatic potential is not well-defined (there is no “vacuum” to reference the energy of an electron to).

These limitations can be overcome to some extent by applying artificial backgrounds.

GAMESS vs. Vasp

`GAMESS`	`Vasp`
System in infinite vacuum	System periodically replicated to ∞
Any charge	Charge-neutral
Cartesian or internals	atom positions relative to supercell
Gaussian basis set	Plane-wave basis set
Hartree units	eV units
Energy referenced to isolated charge particles	Energy referenced to potential model
1 input file	POSCAR: structure input
	INCAR: program options
	POTCAR: identity of atoms
	KPOINTS: k-point sampling
2 output files	OUTCAR: main output
	OSZICAR: summary output
	CONTCAR: final geometry
	CHGCAR: charge density on 3D grid
	WAVECAR: final wavefunctions
	DOSCAR: summary densities of states

Vasp documentation available at https://cms.mpi.univie.ac.at/wiki/index.php/The_VASP_Manual.

Vasp structure definition

POSCAR file specifies lattice vectors, numbers and types of atoms, and their positions in either Cartesian or fractional coordinates. Vasp uses convention in POSCAR that atoms be specified in groups of like type, and that the order in the POSCAR corresponds to the order of atoms in a composite POTCAR. Here’s a POSCAR for an \ce{N2} in a \ce{10}{\angstrom} cubic box.

10 Ang. cell, diagonal orientation
 10.
 1.0 0.0 0.0
 0.0 1.0 0.0
 0.0 0.0 1.0
2
cartesian
     0. 0. 0.
0.0692820323028 0.0692820323028 0.0692820323028

Vasp generates a CONTCAR at the end of any job of identical format to a POSCAR. XDATCAR contains a trajectory of the same format.

Vasp model specification

Electronic structure

ISPIN = 1, 2: non-spin-polarized or spin-polarized
NUPDOWN = xx: set number of up - down electrons
MAGMOM = xx: specifies initial magnetic moment of each atom
NELECT: specifies total number of electrons (seldom used)
FERWE: specifies occupancy of spin-up levels (seldom used)
FERDO: ditto for spin-down
ISMEAR: specifies how electrons are distributed, or “smeared out,” amongst orbitals/bands near Fermi level
ISMEAR = 0: “Gaussian” smearing, a safe options for molecules and most solids
ISMEAR = 1: “Methfessel and Paxton,” useful for geometry optimizations for metals
ISMEAR = -5: tetrahedron method with Blöchl corrections, useful for accurate single-point energies
SIGMA = 0.05: smearing parameter for ISMEAR. For molecules and insulators, \SI{0.05}{eV} is a sensible value. \SI{0.2}{eV} for metals. Larger values can help convergence but too large gives unphysical results. Aim to keep the electronic “entropy” (difference between free and total electronic energies) \SI{<1}{meV/atom}.

Electron-electron interaction model

GGA: specifies exchange-correlation model. Defaults to POTCAR value, typically PBE GGA. Many options possible. It is only meaningful to compare energies between calculations computed with exactly the same exchange-correlation model.
GGA=RPBE: “revised” PBE
GGA=OR: “optimized” PBE
GGA=BF: “Baysian-error estimated functional”
LHFCALC: turns on exact H-F exchange. Typically quite expensive and used sparingly.
HFSCREEN=0.6: turns on HSE06 “screened” exchange, which interpolates smoothly between exact exchange at short distance and GGA exchange at long range

Dispersion corrections

Generally methods that include some wavefunction-based electron correlation are needed to capture dispersion (non-bonded) interactions well. A pragmatic approach in plane-wave calculations is to add-on a dispersion correction, so that $E_\text{disp-corrected} = E_\text{DFT} + E_\text{disp}$.

IVDW: turns on approximate corrections to introduce a dispersion correction to energy/forces
IVDW=10: “DFT-D2” method of Grimme, adds on Lennard-Jones-type dispersion correction. Function only of atom positions.
IVDW=11: “DFT-D3” method of Grimme, corrects for differences in local coordination.
IVDW=20: Tkatchenko-Scheffler method, uses converged charge density to parameterize D2 model

Precision parameters

Precision controlled by size of basis and grids used to transform wavefunctions and density back and forth between real and reciprocal space. It is only meaningful to compare energies between calculations computed with exactly the same precision parameters.

ENCUT: sets plane-wave cutoff (basis set size), in eV. Defaults to largest value in supplied POTCARs.
PREC: general precision parameter to set ENCUT and corresponding FFT grid and PAW parameters
PREC=NORMAL: for routine calculations
PREC=ACCURATE: for higher precision, eg for forces for a Hessian calculation
LREAL: controls evaluation of parts of PAW potential
LREAL=.FALSE.: (default) evaluates exactly in reciprocal space. More computationally costly, esp for larger systems.
LREAL=Auto: generates optimized real-space representation of PAWs. Typically cheaper, esp for larger systems.
ROPT: accuracy of LREAL real-space method
IDIPOL=1-3: turns on dipole corrections for dipolar supercell

SCF

Vasp uses same essential idea as GAMESS. Guess initial charge density and set of orbital coefficients. Vasp initial guesses are not internally consistent: charge density typically a superposition of atoms, but wavefunctions are random. Density and wavefuctions are optimized iteratively and in parallel until internally consistent. A mixture of algorithms are used, in particular conjugate-gradient, Davidson, and DIIS.

ALGO=Normal: (default) block Davidson optimization, very safe
ALGO=Fast: starts with block Davidson, then switches to DIIS. Good compromise method.
ALGO=VeryFast: only DIIS. More aggressive.
NBANDS: how many orbitals (bands) to calculate. You always want to calculate all the occupied ones and, for numerical reasons, at least some of the unoccupied ones. In general for performance reasons you don’t need or want to calculate all of the empty ones. The default is safe.
NLEM: maximum number of SCF cycles
EDIFF=1e-5: precision criterion to stop SCF, in eV. 1e-4 is generous, 1e-5 normal, 1e-6 more accurate

Core electron treatment

Wavefunctions oscillate too rapidly near core to capture reasonably using plane-waves; cut-off would be too high. Always combine plane waves with some model for core states of the system. Typically assume core electrons are “frozen,” and incorporate relativistic effects. Analogous to core potentials in GAMESS.

Orthogonalized plane-wave (OPW)

If we know the core electron wavefuctions, and recast DFT model to calculate difference between any valence orbital and the core orbitals. Create new OPW basis: \[ \X_G^\text{OPW}(\bm{r}) = φ_G(\bm{r}) + ∑_\text{m} b_mG χ_m(\bm{r}) \] where $χ_m$ are core orbitals that are fixed. Equivalent to modifying the electron-electron interaction by \[ \hat{v}^\text{OPW} = \hat{v}_\text{ee} + ∑_m (ε - ε_m) |χ_m\rangle\langleχ_m| \] The last part is expensive to evaluate, but can be approximated.

Pseudopotential

Replace last part with a potential designed to replicate the influence of the core electrons (and nuclei) on the valence. Loose any description of core, but (hopefully) valence is retained.

Cutoff radius - boundary between core and valence regions
Transferability - ability of PP to be applied in different chemical environment
Local vs non-local - spherically symmetric vs angular-momentum-dependent
Norm-conserving - preserves $|ψ|^2$ of wave-functions
Soft vs hard - location of boundary. Closer in, “harder” the potential, higher the plane-wave cutoff needed. Further out, “softer” the potential, lower the plane-wave cutoff.
Ultra-soft pseudo-potential - (Vanderbilt) relax norm-conserving requirement and fix up with augmentation charges. First potentials that enabled large-scale calculations

Augemented plane-wave (APW)

Partition space into a spherical core and plane-wave valence and patch up at the boundary

Projector-augmented plane-wave (PAW)

Most modern approach, combines advantages of APW and of ultrasoft PPs. Not strictly a PP approach: constructs full wavefunction and density as combination of valence (plane wave) part and precomputed core parts from atomic calculations using same XC potential. Retains full nodal structure of valence wavefunctions. Generation of these potentials is not for the casual user (just as the case for making Gaussian basis sets). The great strength of Vasp is a very complete and reliable library of PAW models for the whole periodic table and for several exchange-correlation functional. Stored in POTCAR.

Wavefunctions and charge densities

Primary outputs of Vasp calculation are wavefunction, density, energy, and forces.

LCHARGE=.TRUE.: turns on creation of the CHGCAR file, which contains the total charge and spin (up-down) densities, evaluated on $NGX × NGY × NGZ $ grid. Not so useful by itself, but can be used to construct charge density differences.
LWAVE=.TRUE.: turns on creation of the WAVECAR file, which contains the final wavefunctions, useful for restarts
LORBIT: controls output of DOSCAR and PROCAR, which contain analysis of total density of states and of each band. Controlled with RWIGS, EMIN, and EMAX

Bader charges can be extracted from CHGCAR.

Exploring potential energy surfaces

A nice feature of plane-wave methods is that the forces on the atoms are gotten essentially for free. No Pulay forces, because ions and basis functions are decoupled.

Generally limited to Cartesian coordinates in PBC calculations. No z-matrix, internals, or redundant coordinates. The cartesian locations of atoms can be fixed. Perhaps ASE offers more options.

Geometry optimizations

Again, same idea as in GAMESS. Use some algorithm to move atoms in response to forces.

EDIFFG=-0.03: convergence criterion, in eV/Å. \SI{-0.05}{eV/\angstrom} for routine, \SI{-0.01}{eV/\angstrom} for more precise.
IBRION=2: conjugate-gradient optimizaiton
IBRION=1: Quasi-Newton-Raphson with DIIS algorithm. Uses a digonal Hessian guess that gets updated. Works well when initial geometry is pretty good, poorly otherwise
NFREE=xx: number of previous steps to include in the DIIS update.
IBRION=3: damped molecular dynamics optimization
NSW: maximum number of geometry steps
POTIM=0.5: tunes scaling of forces to optimize optimization algorithm.
ISW=1: only move the atoms in the optimization

Harmonic frequency calculation

IBRION=5: numerical Hessian and frequency evaluation, moving all atoms
IBRION=6: take advantage of symmetry to reduce number of unique geometries to evaluate forces
NFREE=2: two-sided differentiation
POTIM=0.01: specifies step-size for finite difference calculation
IBRION=7,8: alternative, analytical perturbation-theory-based evaluation of Hessian matrix

Transition states

Commonly found using either a “nudged elastic band” (NEB) method or “dimer” method. NEB requires definition of a “string” of structures, will address later. “Dimer” needs an initial guess of transition state and a (good) guess of the reaction coordinate.

IBRION=44: specifies dimer calculation. See sample input below. Get trial direction from prior Hessian calculation.

ammonia flipping
1.
6. 0. 0.
0. 7. 0.
0. 0. 8.
H N
3 1
cart
       -0.872954        0.000000       -0.504000        ! coordinates for atom 1
        0.000000        0.000000        1.008000
        0.872954        0.000000       -0.504000
        0.000000        0.000000        0.000000        ! coordinates for atom N
       ! here we define trial unstable direction:
        0.000001    0.522103   -0.000009        ! components for atom 1
       -0.000006    0.530068    0.000000
       -0.000005    0.522067   -0.000007
        0.000001   -0.111442    0.000001        ! components for atom N

IMAGES: turns on nudged elastic band

Molecular dynamics

IBRION=0: molecular dynamics, propogating atomic positions forward in time according to forces. Typically combine with PREC=LOW.
NSW: number of MD steps
POTIM=0.5: step size, in fs
SMASS=-3: NVE dynamics, which (should) conserve energy
SMASS=-1: scaled velocity MD, eg for simulated annealling
SMASS>0: NVT dynamics, SMASS specifies the thermostat mass
TEBEG: initial temperature, in K
TEEND: final temperature, if temperature is to be scaled

\newpage

Periodic electronic structure

Isolated vs. periodic systems

Power of periodic boundary condition model is ability to describe things that are extensive: solids, surfaces, even liquids.

./Images/ClusterPeriodic.pdf

Key difference is that atoms communicate across the boundary in the periodic model. Results in discrete states $→$ continuum of electronic states.

Bloch’s theorem and qualitative band structure

Images borrowed from https://doi.org/10.1002/anie.198708461.

Cyclic chain of H atoms of increasing length (Born-van Karman boundary condition):

./Images/CyclicH.pdf Energy levels increasing with chain length. Wavefunctions grow in a regular way, consistent with (and identifiable by) underlying symmetry. Nodes increase systematically. In limit, get an infinite band of levels.

Infinite chain of H atoms separated by a lattice constant a. Electrons move in that periodic potential. Expand wavefunctions in atom-centered 1s orbitals.

./Images/Hchain.pdf

Bloch’s theorem states that the wavefunctions can be written as the product of a cell-invariant part (the H 1s functions here) and a cell-periodic part. The periodic part is indexed by the wavevector k. k takes as many values as there are periodic units N: \[ k = \frac{π}{a}\frac{n}{N},\qquad n=-\frac{N-1}{2}, \ldots, 0, \ldots, \frac{N-1}{2} \] If N is infinite, then k is a continuous and real variable. The space of unique values of k is called the first Brillouin zone. The periodic phases of the basis functions correspond to an underlying wavelength associated with k: \[λ_k = 2π/k\] By the de Broglie relationship, then, k relates to the momentum of an electron in that energy level. Each k is degenerate such that $E(k) = E(-k)$ (moving with the same momentum to the left and right gives the same energy).

The continuum of energy levels in k is called the band structure. Conventionally plotted vs $|k|$. The width of the band is called the dispersion:

./Images/banddispersion.pdf

The highest filled energy state, in this case (k=π/2a\), is the /Fermi energy. Band structure on the right would correspond to a metal, as the energy difference between different momentum states is zero.

Band dispersion depends on cell overlap, dispersion depends on type and orientation of orbitals:

Band folding

Multi-dimensional periodicity

Density of statues

Bravais lattices

Quantitative supercell calculations

Brillouin zone integration

k-point sampling

Fermi smearing

\newpage

Practical supercell calculations

\newpage

Surfaces

Surface planes

Slab models

Surface energy

Surface potentials and Fermi energies

Surface adsorption

Coverage-dependent adsorption

Reaction barriers

\newpage

Implicit solvation

Density functional theory

Electron density ρ as fundamental quantity

Thomas-Fermi-Dirac model

Hartree-Fock-Slater model

Hohenberg-Kohn theorems

Kohn-Sham construction

Exchange-correlation functionals

So Kohn et al. showed that the DFT approach is theoretically well-grounded and provided one way to practically apply it. Promise is that if we can find an approximation to the (unknown) true $v_\text{xc}$ with the right balance of simplicity and accuracy, we will have one sweet theory. Has to incorporate both exchange, like Slater tried to do, and correlation.

How to proceed? Lots of approaches, and jargon here is at least as bad as in wavefunction-based methods. Perdew 2006 describes the “Jacob’s ladder” of approximations:

LDA

One well-defined limit is the homogeneous electron gas, and this is the usual starting point for modern approximate DFT methods. Assume exchange and correlation potentials at any given point depend only on the value of ρ there (or spin-up and spin-down ρ, if spin-polarized). We know from Slater and Dirac’s work what the exchange potential is for this system.

It is possible to determine numerically the correlation energy for a given density from quantum Monte Carlo calculations. Ceperley and Alder (PRL 1980, 45, 566) did this to very high accuracy, and others (Vosko, Wilk, and Nusair, “VWN”, and Perdew and Wang, “PW”) fit these numerical results to analytical models in ρ. This combination of local exchange and correlation defines the LDA model.

LDA offers modest improvement over HFS for molecules. “Homogeneous” approximation pretty severe for an atom or molecule. Nonetheless, works surprisingly well for structures and charge distributions, but has problems in calculating accurate bond energies, typically overbinding. Also tends to underestimate the HOMO-LUMO gap in molecules and analogous band gap in solids.

FilesExpand file tree

comp-notes.org

Latest commit

History

comp-notes.org

File metadata and controls

Introduction

What do we care about?

How are we going to figure these out? With only a computer?

Our goals

Reading resources

Software tools

Notebooks

Molecular methods

Supercell methods

Great for getting started

Refresher on Quantum Mechanics

Why quantum mechanics?

Postulates of non-relativistic quantum mechanics

Notes on constants and units

Example: Energy states of a particle in a box

Hydrogen atom: simplest chemical “thing”

Schrödinger equation

Analytical solutions

Variational principle

Basis functions

Secular equations

Spin

(Two is too) many electrons

Schrödinger equation for He

The Hartree atom

The Pauli principle

Slater determinants and Hartree-Fock

Hartree-Fock equation

Basis of wavefunction theory (WFT)

Hartree-Fock-Slater

Notes

Basis of density functional theory (DFT)

Implementations

GAMESS

FDA

Performance

Correlation

Practical electronic structure

Born-Oppenheimer approximation

Model chemistry

Beyond Hartree-Fock

Bring back the basis sets

H2O Energy Example

Symmetry

Examples

Dissociating H2+ example

Dissociating HHe+ example

Dissociated He2 example

Open-shell systems

Gaussian basis sets

Gaussian basis set nomenclature

Standard basis sets in GAMESS

Electron cores

Relativistic effects

Implementations

Population analysis

Molecular orbital (MO) diagrams

Implementation details of SCF methods

SCF updating

Potential energy surfaces

Specifying atomic positions

Cartesian

Internal coordinates

Features of potential energy surfaces

One-dimensional example

Many-dimensional PES features

Energy gradients and second derivatives

Gradients

Hessian

Geometry optimization algorithms

Energy-only

Energy + gradient

Energy + gradient + higher order

Machine learning

`GAMESS`

`FDA`

Standard basis sets in `GAMESS`

Specify `GAMESS` calculation type