Introduction to Electronic Structure Theory Dage Sundholm Department of Chemistry University of Helsinki Finland

CSC Spring School, 20.3.2013

What is Electronic Structure Theory • Electronic structure theory (EST) describes the motion of

electrons in atoms and molecules. • The motion of electrons and nuclei of molecules can be

separated into 1. Electronic motion 2. Vibrational motion of the nuclei 3. Rotational motion of molecules • In most electronic structure calculations, one is mainly

interested in (1) • In this brief introduction, (1) will be discussed, (2) is

mentioned, and (3) is omitted.

Born-Oppenheimer Approximation • Nuclei are much heavier than electrons • The lightest nucleus i.e., the proton is 1836 times heavier

than the electron. • The heavy nuclei move much slower than the fast electrons

• The electronic motion can always adjust to the nuclear

positions • The electronic Schrodinger ¨ equation can be solved for

given fixed positions of the nuclei • The decoupling of electronic and nuclear motions is called

Born-Oppenheimer approximation (BOA)

Born-Oppenheimer Approximation • The decoupling of electronic and nuclear motions is called

Born-Oppenheimer approximation (BOA) • BOA is the birth of molecular structure • Calculating the electronic energy for a number of nuclear

positions yields the potential energy surface of the molecule • At the non-relativistic Schrodinger ¨ level, the Hamiltonian is 2 X ~2 X 2 X ZA e2 ˆ =−~ H ∇2i − ∇A − 2m 2MA 4π0 riA i

+

A

i,A

X ZA ZB e2 X e2 + 4π0 RAB 4π0 rij

A>B

i>j

Hamiltonian 2 X ~2 X 2 X ZA e2 ˆ =− ~ H ∇2i − ∇A − 2me 2MA 4π0 riA i

+

A

i,A

X e2 X ZA ZB e2 + 4π0 RAB 4π0 rij

A>B

i>j

where • me is the electron mass, • e is the charge of the electron, • MA is the mass of nucleus A, • ZA is charge of nucleus A, • riA is the distance between electron i and nucleus A, • RAB is internuclear distances, • rij is interelectronic distances.

Hamiltonian • In EST calculations, atomic units (a.u.) are used. • • • •

ˆ =− H

The electron mass me = 1 The charge of the electron e = 1 Dirac’s constant ~ = 1 1 Coulomb’s constant 4π =1 0

X1 i

2

∇2i −

X 1 X ZA X ZA ZB X 1 ∇2A − + + 2MA riA RAB rij A

i,A

A>B

i>j

• The 1st term is the kinetic energy of the electrons • The 2nd term is the kinetic energy of the nuclei • The 3rd term is the electron-nucleus Coulomb attraction • The 4th term is the nucleus-nucleus Coulomb repulsion • The 5th term is the electron-electron Coulomb repulsion

Hamiltonian ˆ =− H

X1 i

2

∇2i −

X 1 X ZA X ZA ZB X 1 ∇2A − + + 2MA riA RAB rij A

i,A

A>B

i>j

¨ The Schrodinger equation reads ˆ HΨ(r , R) = EΨ(r , R) The decoupling of electronic and nuclear motions (BOA) yields Ψ(r , R) = Ψ(r ; R)Φ(R) • Ψ(r , R) is a function of electron and nuclear coordinates. • Ψ(r ; R) is a function of electron coordinates for given

parameters R. • Φ(R) is a function of the nuclear coordinates.

Electronic and Nuclear Schr¨odinger Equations ¨ The electronic Schrodinger equation in BOA   X1 X ZA X ZA ZB X 1 −  Ψ(r ; R) ∇2 − + + 2 i riA RAB rij i

i,A

A>B

i>j

= Ee (R)Ψ(r ; R) ¨ The nuclear Schrodinger equation in BOA ! X 1 − ∇2 + Ee (R) Φ(R) = En Φ(R) 2MA A A

Computational Procedure • The electronic Schrodinger ¨ equation is solved for given R • The nuclear positions are adjusted to obtain the lowest

energy. • The molecular structure with the lowest energy is called

equilibrium structure or equilibrium geometry. • The equilibrium structure is THE molecular structure. • Check whether it is a minimum or a saddle point! • For a minimum, the second derivative of the energy with

respect to nuclear displacements should have a positive curvature in all dimensions. • A first estimate of the nuclear motion (vibrational

frequencies) of the molecule is the by-product.

Example of potential energy surfaces

Wave packages on potential energy surfaces Chemical reactions and other events involving nuclear motion can be modeled as wave packets moving on the potential energy surface.

Computer Exercise In the Computer exercise this afternoon you optimize the molecular structure of some small molecules. • Optimize the molecular for your favourite small molecule • Investigate how the molecular structure depends on the

level • Check whether you have found a minimum or a transition

state • Search for the transition state of a chemical reaction • Simulate the vibrational spectrum

Mathematics: Dirac Notation Electronic structure theory needs appropriate mathematical tools. A wave function ψ(r) can be represented as a vector |Ψi =

n X

ai |Φi i

i=1

where |Ψi is a state vector, ai are expansion coefficients, and |Φi i are fixed basis vectors of the Hilbert space.   a1  a2 
  |Ψi =  .  hΨ| = a1∗ a2∗ · · · an∗  ..  an

Dirac Notation The scalar (inner) product is defined as

hΨa |Ψb i =

n X

ai∗ bi = ha|bi

i=1

Z

∞

Ψ∗a (r )Ψb (r )dr

hΨa |Ψb i = ∞

ˆ is characterized by its effect on the basis vectors An operator A n X

ˆ ji = A|Ψ

|Ψi iAij

i

ˆ = Aij = hi|A|ji

Z

∞

∞

ˆ j (r )dr Ψ∗i (r )AΨ

Slater determinants The Slater determinant is a useful but clumsy notation χ1 (x1 ) χ2 (x1 ) · · · χN (x1 ) 1 χ1 (x2 ) χ2 (x2 ) · · · χN (x2 ) Ψ= .. .. .. .. N! . . . . χ1 (xN ) χ2 (xN ) · · · χN (xN ) The equivalent Dirac notation is simpler and more practical Ψ = |χ1 (x1 )χ2 (x2 ) · · · χN (xN )i The notation using occupation number (ni ) vectors is Ψ = |n1 , n2 , · · · nN i where ni is 0 or 1. Appropriate mathematics is developed to simplify expressions involving occupation number vectors.

Occupation number vectors

• The abstract occupation number vectors can be used to

represent Slater determinants. • The main advantage is that an efficient algebra can be

developed for them. |Ψi = |k1 , k2 , · · · , kN i ki = 1 for occupied orbital and ki = 0 for unoccupied orbital Inner product: hk |mi = ΠN i=1 δki mi |ci =

P

k

ck |k i,

hc|di =

P

k

ck∗ dk ,

|kihk | = 1

Second Quantization The whole machinery of the second quantization is obtain from the algebra of the occupation number vectors. Creation operator: aP† |k1 , k2 , · · · , 0P , · · · , kN i = δkP ,0 ΓkP |k1 , k2 , · · · , 1P , · · · , kN i kQ with ΓkP = ΠP−1 Q=1 (−1)

Annihilation operator: aP |k1 , k2 , · · · , 1P , · · · , kN i = δkP ,1 ΓkP |k1 , k2 , · · · , 0P , · · · , kN i Anticommutation relations: † ]+ = 0 [aP , aQ ]+ = 0 [aP† , aQ

[aP† , aQ ]+ = δPQ

Operators in Second Quantization ˆ Hamiltonian and Energy (E = hHi): X X ˆ = H hpq ap† aq + gpqrs ap† ar† as aq pq

pqrs

One-electron interaction (e.g., dipole moment hˆ µi with µ ˆ = {x, y, z}): X µpq ap† aq µ ˆ= pq

Two-electron interaction (e.g., Coulomb and exchange energy 1 h r12 i): X 1 1 = ( )pqrs ap† ar† as aq r12 pqrs r12

The N-body Wave Function The most simple ansatz for the N-body wave function is the Hartree Product Ψ(x1 , x2 , · · · , xN ) = χ1 (x1 )χ2 (x2 ) · · · χN (xN ) where xi are 4N coordinates, i.e., the three Cartesian coordinats plus the spin. χi (xi ) are spin orbitals. χi (xi ) = φi (ri )|αi or χi (xi ) = φi (ri )|βi |αi and |βi denote spin up and spin down, respectively. The Hartree product does not fulfil the antisymmetry condition of the electrons (fermions).

Accurate N-body Wave Functions John C. Slater introduced determinants as N-body wave function. The Slater determinants fulfil the antisymmetry condition of fermions. How can the wave function be improved? The use of large one-particle basis sets makes it possible to approach the basis-set limit. By writing the wave function as a linear combination of many Slater determinants improves the quality of the wave function. By considering all possible Slater determinants in a given ¨ one-particle corresponds to the solution of the Schrodinger equation in that basis.

The Hartree-Fock Model The Hartree-Fock or the mean-field model is obtained by assuming that the wave function consists of one single Slater determinant.

In the Hartree-Fock model, the electrons are assumed to move in the average field of all other electrons and nuclei. The electronic motion is uncorrelated.

The Hartree-Fock Model The equations of the Hartree-Fock model is obtained by expressing the energy as expectation value of the Hamiltonian for a one-Slater-determinant wave function ˆ HF i EHF = hΨHF |H|Ψ Differentiating the energy expression with respect to linear parameters of the one-particle functions (orbitals) φi Ensuring that the orbitals are orthonormal hφi |φj i = δij The Hartree-Fock Lagrangian reads L(φi ) = EHF (φi ) −

X ij

ij hφi |φj i − δij




The Fock Equation The Hartree-Fock energy for closed-shell systems EHF = 2

X

hii +

i

X

2giijj − gijji




ij

ˆ φ (1)dr1 hij = φ∗i (1)h R R ∗ 1 j∗ −1 gijkl = φi (1)φk (2)r12 φj (1)φl (2)dr1 dr2 R

The molecular orbitals (MO), φi , are expanded in a set of non-orthonormal atomic orbitals (AO) χµ φi =

X

χµ Cµp

µ

Inserting it into the energy expression and differentiation with respect to Cµp yields the Fock equation

The Fock Equation The Fock equation for non-orthonormal AOs reads X

AO Fµν Cνi = i

ν

X

Sµν Cνi

ν

where the Fock matrix has been introduced X
AO 2gµνii − gµiiν = hµν + Fµν i

The AO Fock matrix can be calculated in the AO basis   X 1 AO AO Fµν = hµν + Dστ gµνστ − gµστ ν 2 στ

The Hartree-Fock Iterative Procedure

1. Guess C 2. Calculate DAO = CDMO CT 3. Calculate FAO 4. Diagonalize FAO C = SC 5. Determine an error vector e.g., e = FDS − SDF 6. Check the convergence 7. Introduce damping using e.g., DIIS 8. Make a new guess for C and go to 2 9. End when converged

Hellmann-Feynman Theorem Assume that the Hartree-Fock wave function (Ψ) is perturbed by an external perturbation V : The perturbed energy can be written as: ˆ + λV |Ψ + λ∆Ψi E(λ) = hΨ + λ∆Ψ|H The first-order correction to the energy is obtained by differentiating E(λ) with respect to λ dE(λ) ˆ + λV |Ψ + λ∆Ψi + hΨ + λ∆Ψ|V |Ψ + λ∆Ψi = 2Reh∆Ψ|H dλ

dE(λ) dλ λ=0

ˆ = 2Reh∆Ψ|H|Ψi +hΨ|V |Ψi | {z } =0

Open-Shell Systems For molecules with unpaired electrons or for general open-shell systems, different orbitals are used for describing the motion of the α (spin-up) and the β (spin-down) electrons. The unrestricted Hartree-Fock (UHF) model is a generalization of the restricted Hartree-Fock (RHF) model for closed-shell molecules. The matrix sizes of the UHF are twice the size of the RHF ones. UHF is not unproblematic as it might introduce spin contamination effects, convergence problems, negative gap between the highest occupied molecular orbital and the lowest unoccupied molecular orbital (HOMO-LUMO gap). Correlation effects are often larger for open-shell systems than for the closed-shell ones.

Density Functional Theory Hartree-Fock does not consider any electron correlation effects by definition ECorrelation = EExact non−relativistic − EHartree−Fock Density functional theory (DFT) is the most simple model that consider electron correlation effects. DFT is structurally quite similar to Hartree-Fock (HF) The effctive HF potential ˆ vˆeff = vnuc + vJ + K The effctive DFT potential vˆeff = vnuc + vJ + vxc

Density Functional Theory DFT is formally exact. However vxc is not known. The art of DFT is to find an accurate approximation to vxc DFT is based on the electron density, ρ(r ), which is the probability to find an electron in dr . The electron density is related to the wave function Z Z Z ρ(r1 ) = n dr2 dr3 · · · drn |Ψ(r1 , r2 , . . . , rn )|2 The density uniquely determines the Hamiltonian, because R • ρ(r )dr = n • The cusps at the nuclei yield the nuclear positions • The slope of the nuclear cusps yields the nuclear charge

Density functional theory Simple DFT approach: • Take the exact ρ(r )

ˆ • Use ρ(r ) to obtain H • Construct and solve the Schrodinger ¨ equation • The exact ground-state energy is obtained

The proceduce can be formulated mathematically as E0 = E[ρ(r )] which is an energy functional that is variational E[˜ ρ] ≥ E[ρ]

Hohenberg-Kohn The energy can be divided into something that explicitly depend on the electron density and the universal part (the rest) Z E[ρ] = dr ρ(r )(vnuc + velectrostatic ) + FHK The energy functional E[ρ] = T [ρ] + J[ρ] +

0 Exc [ρ]

Z +

dr ρ(r )vnuc

The kinetic contribution is: T [ρ] = Ts + Tc [ρ]

with

Ts = −

1 X hψi |∇2 |ψi i 2 i∈occ

The Kohn-Sham energy is: Z EKS = Ts + J[ρ] + Exc [ρ] +

dr ρ(r )vnuc

Kohn-Sham Equations In

Z EKS = Ts + J[ρ] + Exc [ρ] +

dr ρ(r )vnuc

everything is known except Exc [ρ], which is small compared to the rest. Different DFT functionals have their own approximation to Exc [ρ] Minimization of EKS with respect to the orbitals yields the Kohn-Sham equations   1 − ∇2i + vnuc + vJ + vxc ψi = i ψi 2

with

vxc =

δExc [ρ] δρ

Density functionals

1

LDA:

vxc ∝ ρ(r ) 3

GGA:

~ vxc = f (ρ, ∇ρ)

Hybrid

~ + EHF,x vxc = f (ρ, ∇ρ)

Meta GGA

~ τ ) where τ = P 1 |∇ψi (r )|2 vxc = f (ρ, ∇ρ, i 2

B3LYP Exc [ρ] = (1−a)ExS [ρ]+aExHF +bExB88 [ρ]+cEcLYP +(1−c)EcVWN [ρ]

a = 0.20

b = 0.72

c = 0.81

Stairway to heaven • The exact functional • ··· • Stairway to heaven • Orbital dependent functionals • Exact exchange functionals • Parametrized functionals • Meta GGA • Hybrid functionals • GGA functionals molecules • LDA functionals formal • Hartree-Fock Slater, the first • Thomas Fermi

Configuration Interaction Configuration interaction (CI) is the most simple ab initio method to consider electron correlation. The wave function is expanded as a linear combination of Slater determinants (configurations). |ΨCI i =

X

CI |Ii

I

The expansion coefficients are determined by minimizing the CI energy, with the constraint that the states are orthonormal. ˆ CI i ECI = hΨCI |H|Ψ

with

hΨCI |ΨCI i = 1

LCI = ECI −  (hΨCI |ΨCI i − 1)

Configuration Interaction The CI wave function considers electron correlation

The CI wave function can be written as |ΨCI i = CHF |ΨHF i + CS |ΨS i + CD |ΨD i + CT |ΨT i + · · · + CN |ΨN i ˆ CI = ECI ΨCI The CI eigenvalue equation: HΨ

Configuration Interaction |ΨS i All configurations that are obtained by single excitations from the |ΨHF i reference. |ΨD i All configurations that are obtained by double excitations from the |ΨHF i reference. |ΨT i All configurations that are obtained by triple excitations from the |ΨHF i reference. |ΨN i All configurations that are obtained by all N excitations from the |ΨHF i reference. |ΨCI i including all terms up to |ΨN i is the full configuration interaction (FCI) or exact diagonalization wave function i.e., ¨ it is the solution of the Schrodinger equation in the given basis set.

Excited Slater Determinants

The Slater determinants of the CI model are obtained by single (S), double (D), triple (T), . . . alternations of the occupation of the Fock space.

Matrix Elements of the CI Hamiltonian Matrix

Brillouin’s theorem says that single excited determinants do not directly couple with the Hartree-Fock reference. The Hamiltonian matrix is banded because only Slater determinants which are conncted via singles and double excitations yield non-vanishing matrix elements.

Configuration Interaction CI is a straightforward method to consider electron correlation effects. Truncated CI wave functions are not size-extensive i.e., the energy of two non-interacting fragments is not twice the energy of one of the fragments. The only exception is the FCI wave function. FCI calculations are computationally expensive. The largest FCI calculation as far as I know had 2500 billion (2.5 1012 ) Slater determinants. FCI calculations are indispensable as benchmark, whereas truncated CI has lost its importance as production tool due to the size problem.

Slater-Condon rules

Most matrix element of the CI Hamiltonian vanish, because configurations whose occupation number vectors differ in more than two position pairs do not directly couple. The reason is that the Hamiltonian contains only one- and two-body interaction terms The non-vanishing matrix elements of the CI Hamiltonian can be identified using second quantization. The obtained expressions are called Slater-Condon rules

One-electron Slater-Condon rules |ki is a occupation number vector containing the occupation of the orbitals for a given Slater determinant. kP = 1 for occupied orbitals kP = 0 for unoccupied orbitals fPQ are the matrix elements of the one-electron operator ˆf Identical bra and ket : hk|ˆf |ki =

X

kP fPP

P

bra and ket differ in one pair : hk|ˆf |li = ΓkI ΓlJ fIJ bra and ket differ in two or more pairs : hk|ˆf |li = 0 with kQ ΓkP = ΠP−1 Q=1 (−1)

Two-electron Slater-Condon rules ˆ gPQRS are the matrix elements of the two-electron operator g 1 2

Identical bra and ket

ˆ |ki = hk|g

differ in one pair

ˆ |li = ΓkI ΓlJ hk|g

differ in two pairs

ˆ |li = ΓkI ΓkJ ΓlK ΓlL (gIKJL − gILJK ) hk|g

P

ˆ |li = 0 differ in more than two pairs hk|g with P−1 ΓkP = ΠQ=1 (−1)kQ

PR

P

kP kR (gPPRR − gPRRP ) R

kR (gIJRR − gIRRJ )

CAS and RAS Models The Complete Active Space and Restricted Active Space

Multireference Configuration Interaction Multireference configuration interaction (MRCI) belong to the category of ”As Full CI as possible” models.

RAS can be considered as a MRCI with orbital optimization. The reference can also consist of a few selected important configurations with single and double replacements to the complementary space.

Multireference Perturbation Theory The CASPT2 method is a multireference (MR) perturbation theory method with a CASSCF wave function as the reference. MRPT should be formulated such that when |Φ0 i → |HF i MRPT → Møller-Plesset PT. This is not trivial because the MR Fock operator is not unique. In CASPT2, the spin-averaged first-order density matrix of the CAS function is used to define the Fock operator. The open-shell and closed-shell configurations are then treated unequally, which is fixed with semi-empirical level shifts The NEVPT2 model avoids the problem by adding two-electron terms to H0 , but it yields orbital labelling dependent results.

Multireference Perturbation Theory In the single reference case, H0 is the sum of one-electron Fock operators, which is diagonal in canonical orbital basis. The CAS reference is not an eigenfunction of the n-electron Fock operator, that is, FCAS is not a proper H0 . This is fixed by using projection operators H0 = PFP + QFQ where P projects to the CAS reference. The Fock operator of the MR case is not diagonal. Unitary transformations in the individual spaces remove off-diagonal elements in each orbital block, but not between the inactive, active and secondary orbital blocks.

Multireference Perturbation Theory The first-order wave function is obtained as in the single reference case as (H0 − E0 ) Ψ1 = (H0 − H) Ψ0 which has to be solved iteratively. The first order wave function is expanded in the first-order interaction space Ψ1 =

X

qs ˆ ˆ Tpr Epq Ers Ψ0

pqrs

ˆ pq is the spin-averaged excitation operator that moves where E an electron from the q orbital of the active or inactive spaces to an empty orbital p of the secondary space.

Multireference Perturbation Theory The configurations of Ψ1 comprise those who have non-vanishing Hamiltonian matrix elements with Ψ0 This is called internal contraction that reduces significantly the number of Tqs pr amplitudes Internal contraction leads to complicated expressions involving higher-order density matrices and non-orthogonal configurations including linear dependencies. CASPT2 is able to handle bond breaking, reactions, transition metals, and excited states. CASPT2 yields pure spin states.

CASPT2 Problems The internal contraction scheme leads to intruder states as for MR coupled-cluster, MRCI and other similar approaches. CASPT2 is not a black box method. It requires experience and a deep understanding of the electronic structure of the investigated system to obtain reliable results. All orbitals that should be in the active space cannot be included due to the huge computational costs. The borderline between active and secondary orbitals is not obvious. CASPT2 with the minimal reference based on chemical intuition often yield unreliable results.

Coupled-Cluster Models The coupled-cluster ansatz for the wave function is ΨCC = exp(Tˆ )Φ0

where

Φ0 = ΦHF

in most cases

Tˆ 2 Tˆ 3 Tˆ n + + ··· 2! 3! n! The cluster operator Tˆ consists of single Tˆ1 , double Tˆ2 , triple Tˆ3 replacement operators. exp(Tˆ ) = 1 + Tˆ +

Tˆ = Tˆ1 + Tˆ2 + Tˆ3 + · · · + Tˆn The maximum is n = number of electrons, which is identical to full CI When the series is cut at Tˆ2 , it is the CCSD model and a cut at Tˆ3 yields the CCSDT model, etc.

Coupled-Cluster Models The cluster operators Tˆk operating on the HF reference generates all Slater determinants for single (S), double (D), triple (T), ... replacements. Tˆ1 =

occ X vir X

tia aa† ai

a

i occ vir

1 X X ab † † tij aa ab ai aj Tˆ2 = 4 ij

ab

tia and tijab are the singles and doubles amplitudes, respectively ˆ + exp(Tˆ ) = 1

     1 ˆ2 1 ˆ3 ˆ ˆ ˆ ˆ ˆ T1 + T2 + T1 + T3 + T2 T1 + T1 · 2 6 | {z } | {z } | {z } single double triple 

Coupled-Cluster Models To keep the expressions as general as possible the coupled¨ cluster Schrodinger equation is multiplied from left by exp(−Tˆ ) ˆ H|Ψi = E|Ψi Ψ = exp(Tˆ )Φ0 ˆ exp(Tˆ )|Φ0 i = E exp(Tˆ )|Φ0 i H ˆ exp(Tˆ )|Φ0 i = E exp(−Tˆ ) exp(Tˆ )|Φ0 i exp(−Tˆ )H ˆ exp(Tˆ ) exp(−Tˆ )H | {z } similarity transformed Hamiltonian

|Φ0 i = E|Φ0 i

Coupled-Cluster Models The coupled cluster equations for the amplitudes are obtained ¨ by projecting the coupled-cluster Schrodinger equation against the excited Slater determinants ˆ exp(Tˆ )|Φ0 i = 0 hµ| exp(−Tˆ )H hµ| are the Slater determinant that enter the coupled-cluster state with connected amplitudes Tˆ1 , Tˆ2 , Tˆ3 , · · · hµ| = hΦ0 |τµ† ;

abc hµ| = hΦai |, hΦab ij |, hΦijk |, · · ·

where τµ are second-quantization operators that creates all possible single, double, triple, ... replacements from the Hartree-Fock reference (Φ0 )

Coupled-cluster energy The coupled-cluster energy is obtained by projecting against the Hartree-Fock reference ˆ exp(Tˆ )|Φ0 i E = hΦ0 | exp(−Tˆ )H The coupled-cluster energy can be expressed as

ECC = EHF +

X

fia tia +

ia

1X 1X hij||abitijab + hij||abitia tjb 4 2 iajb

iajb

hpq||rsi = hpq|rsi − hpq|sr i Z Z hpq|rsi =

−1 ∗ φ∗p (1)φr (1)r12 φq (2)φs (2)dr1 dr2

fia are elements of the Fock matrix in the MO basis

Approximate Coupled-Cluster SD and Related Methods Coupled-cluster energy:

Second-order Møllet-Plesset perturbation theory

Coupled-cluster equations:

Similarity-transformed Hamiltonian:

Approximate Coupled-Cluster and Related Methods Jacobi matrix:

Jacobi matrix for CC2:

Jacobi matrix for CIS(D∞ ):

Jacobi matrix for ADC(2):

Approximate Coupled-Cluster and Related Methods

Time-dependent density functional theory TDDFT response equations:

Van der Waals Interaction The van der Waals (vdW) interaction originates from fluctuating dipoles due to the correlated motion of the electrons. Hartree-Fock and the commonly used DFT functionals are not able to take vdW interaction into account. Second-order Møller-Plesset perturbation theory (MP2) is the cheapest ab initio method that include vdW interaction. MP2 overestimate vdW interaction. The spin-component-scaled (SCS) MP2 yields often similar vdW energies as the computationally more expensive CCSD(T) method. A simple semi-empirical patch to consider vdW interactions is Grimme’s D3 method, which consider vdW interaction at molecular mechanics (MM) level. B3LYP+D3 is good and pragmatic level when optimizing the molecular structure of large molecules such as biomolecules.

Basis sets Basis sets can be divided into the • Gaussian-type (GTO) basis functions (in most quantum

chemistry programs) • Slater-type (STO) basis functions (in ADF), efficient for

DFT • Plane-wave basis sets, the favourite in the physics

community (VASP) • Numerical basis functions (GPAW uses grid and

atom-centered basis functions) • Fully numerical (Octopus)

GTOs are better than STOs in ab initio correlation calculations, because GTOs have the wrong asymptotic form at the nuclei ¨ and at infinity. (Almlof)

Classes of Gaussian-Type Basis Sets • Cartesian basis sets χi (r ) = x l y m z n exp(−αi r 2 ), (6d,10f) • Spherical basis sets χi (r ) = r l Ylm (θ, ϕ) exp(−αi r 2 ), (5d, 7f) • Even-tempered basis sets • Uncontracted primitive basis sets with constant ratio between the exponents αi−1 /αi = αi /αi+1 • General contracted basis sets • All basis functions comprise all primitive basis functions • Segmented contracted basis functions • As general contracted, but for small blocks of primitive basis functions • Completeness-optimized uncontracted basis sets • Optimize the exponents to obtain as complete basis set as possible in a given range of exponents

Basis sets Basis sets can be downloaded from https://bse.pnl.gov/bse/portal A few personal principles concerning the choice of basis set. • Use as large basis set as possible • Use more than one basis set • The series of Dunning basis sets is very useful when

aiming at benchmarking. • The Dunning basis sets are large, which put limitations on

the size of the molecule and the level of theory. • Diffuse basis functions might be needed when the

molecule is not very large or when it is negatively charged • Special basis sets might be need for some properties.

Basis sets A few personal principles concerning the choice of basis set. • The MOLCAS ANO basis sets should be used only in

combination with MOLCAS and other programs that can take the advantage of general contracted basis sets. • I prefer the Karlsruhe basis sets • the def2-TZVP basis sets is usually enough for my

purposes, def2-SVP in molecular structure optimization of large molecules. • When diffuse basis functions are needed, the def2-TZVPD

is much more cost efficient than Dunning’s aug- basis sets • The Pople basis sets are not very accurate • Use Stuttgart pseudopotentials instead of LANL2DZ ECPs

Direct Diagonalization Large eigenvalue problems cannot be solved by constructing the matrix and diagonalizing it. The largest matrix that fit into the memory of a modern computer is say 64 GB = 2 · 109 double precsion words. The size of the largest matrix is then less than 105 The largest CI problem that has been solved consisted of 2.5 · 1012 Slater determinants The CI matrix has to be diagonalized using direct methods Direct methods means that the linear transformation H C = σ is performed without construction of H.

Davidson-Olsen diagonalization The Hamiltonian matrix can be written as H = H0 + H1

(1)

The eigenvector and the eigenvalue is analogously written as C = C (0) + C (1)

E = E (0) + E (1)

and

(2)

The approximate eigenvalue E (0) is obtained as T

E (0) =

C (0) HC (0) C (0)T C (0)

(3)

inserting Eq. (1) and Eq. (2) into HC = EC and neglecting quadratic terms yields (H0 − E (0) )C (1) = −(H − E (0) )C (0) + E (1) C (0)

(4)

Davidson-Olsen diagonalization

The correction vector C (1) is the obtained as C (1) = −(H0 − E (0) )−1 [(H − E (0) )C (0) − E (1) C (0) ]

(5)

The energy correction E (1) is obtained by multiplying Eq. (5) T from the left by C (0) and require orthonormality of the C (0) and C (1) vectors T

E

(1)

C (0) (H0 − E (0) )−1 (H − E (0) )C (0) = C (0)T (H0 − E (0) )−1 C (0)

which can be inserted in Eq. (5)

(6)

Davidson-Olsen diagonalization The correction vector C (1) obtained using C (1) = −(H0 − E (0) )−1 × "

T

(H − E

(0)

)C

(0)

C (0) (H0 − E (0) )−1 (H − E (0) )C (0) (0) − C C (0)T (H0 − E (0) )−1 C (0)

# (7)

is orthogonal against C (0) by construction regardless of how the partitioning of H into H0 + H1 is done. • The algorithm requires that the inverse (H0 − E (0) )−1 times

a vector must be calculated. • In the Davidson diagonalization algorithm H0 is the

diagonal of the H matrix and the E (1) correction term is omitted.

Davidson-Olsen diagonalization The Davidson-Olsen formulation, which is correct to first order, has the following advantages as compared to the slightly simpler Davidson diagonalization procedure • The update vector C (1) is orthogonal against C (0) • A better H0 does not introduce any new formal problems

not even for the case H0 → H. • Block diagonal H0 or other more accurate approximations

to H can be used. • An algorithm with only two vectors (C (0) and C (1) ) leads to

convergence in only a few iterations. • Very large CI problems can therefore be solved

The zeroth-order Hamiltonian The better H0 the better is the convergence.   Hblock LT H0 = L D The matrix times vector operation is the H0 x = y  x=

x1 x2



 y=

Hblock x1 + LT x2 = y1 Lx1 + Dx2 = y2

y1 y2



The zeroth-order Hamiltonian x2 = D−1 (y2 − Lx1 ) Hblock x1 + LT D−1 (y2 − Lx1 ) = y1 

T

−1

 L x1 = y1 − LT D−1 y2

T

−1

Hblock − L D

x1 = Hblock − L D

 −1  T −1 L y1 − L D y2







x=

x1 x2

 y=

y1 y2



Direct Inversion of the Iterative Space One of the most important developments in electronic structure calculations is the Direct Inversion of the Iterative Space (DIIS) method for damping the iterative procedure of the solution of the Fock and Kohn-Sham equations. In the update of the new Fock matrix, the Fock matrices are considered FnAO

=

n X

ωi FiAO

i

where FiAO are the Fock matrices of the previous iterations and the DIIS weights ωi are obtained by solving the DIIS equation.

Direct Inversion of the Iterative Space The DIIS equation is obtained by constructing the DIIS Lagrangian ! m X L = ω † Bω − λ 1 − ωi , i

with Bij = hei |ej i. constructed from the error vector ei = FiAO Di S − SDi FiAO Di is the density matrix of iteration i and S the overlap matrix

Direct Inversion of the Iterative Space The DIIS equation is a matrix equation of the size of the iteration space or of the m last iterations.     0 B11 B12 · · · B1m −1 ω1  B21 B22 · · · B2m −1   ω2   0      .. .. .. ..   ..  =  .. ..  .  .   . . . . .      Bm1 Bm2 · · · Bmm −1   ωm   0 −1 λ −1 −1 · · · −1 0

      

Solution of the DIIS equation yields ωi which combined with FkAO , k = 1, m yields the new extrapoleted (damped) FiAO

The End Now this is not the end. It is not even the beginning of the end. But it is, perhaps, the end of the beginning.

Sir Winston Churchill, Speech in November 1942 British politician (1874 - 1965)