Algorithm for solving the eigenvalue reponse equation to obtain excitation energies

(1)

Master's thesis, 30 hp | Educational Program: Physics, Chemistry and Biology Spring term 2016 | LITH-IFM-A-EX-16/3145-SE

Algorithm for solving the

eigenvalue response equation to

obtain excitation energies

Daria Burdakova

Examinator, Patrick Norman Tutors, Joanna Kauczor, Dirk Rehn

(2)

Datum

2016-04-04

Avdelning, institution

Division, Department

Department of Physics, Chemistry and Biology

Linköping University

URL för elektronisk version

ISBN

ISRN: LITH-IFM-A-EX--16/3145--SE

_________________________________________________________________

Serietitel och serienummer ISSN

Title of series, numbering ______________________________

Språk Language Svenska/Swedish Engelska/English ________________ Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport _____________ Titel Title

Algorithm for solving the eigenvalue response equation to obtain excitation energies

Författare

Author Daria Burdakova

Nyckelord

Computational physics, numerical spectroscopy, response theory, numerical algorithm, excitation energies

Sammanfattning Abstract

Light-matter interactions lead to a variety of interesting phenomena, for example photosynthesis which is a process fundamental to life on earth. There exists many different spectroscopic methods to measure light-matter interactions, for example UV/Vis spectroscopy, that can provide information about electronically excited states. However, numerical methods and theory are important to model and gain understanding of these experiments. Quantum chemistry provides that understanding, giving the possibility to numerically calculate molecular properties like excitation energies.

The aim of this thesis was to implement a reduced-space algorithm in Dalton, to solve an eigenvalue equation obtained by response theory, for the calculation of excitation energies of molecular systems.

There already was a similar algorithm in Dalton, that was able to perform these calculations. However, in a different module of Dalton used mainly for complex response theory, an algorithm to obtain eigenvalues was missing. The new implementation was similar to the existing one, except for the division of the reduced space into even and odd parts used in the complex response module.The thesis starts with a quick introduction of light-matter interactions and proceeds with a description of many-body theory, including numerical methods used in that field. In the end of the theoretical part, the eigenvalue equation, used to calculate excitation energies, is derived. In the following section, the reduced-space algorithm is described. In the end of the thesis, numerical results obtained with the algorithm are presented, including a small basis set and method study. The comparison with the existing implementation of the similar algorithm verified the successful implementation of the algorithm presented in this thesis.

(3)

(4)

Abstract

Light-matter interactions lead to a variety of interesting phenomena, for exam-ple photosynthesis which is a process fundamental to life on earth. There exists many different spectroscopic methods to measure light-matter interactions, for example UV/Vis spectroscopy, that can provide information about electronically excited states. However, numerical methods and theory are important to model and gain understanding of these experiments. Quantum chemistry provides that understanding, giving the possibility to numerically calculate molecular properties like excitation energies. The aim of this thesis was to implement a reduced-space algorithm in Dalton, to solve an eigenvalue equation obtained by response theory, for the calculation of excitation energies of molecular systems.

There already was a similar algorithm in Dalton, that was able to perform

these calculations. However, in a different module of Dalton used mainly for

complex response theory, an algorithm to obtain eigenvalues was missing. The new implementation was similar to the existing one, except for the division of the reduced space into even and odd parts used in the complex response module.

The thesis starts with a quick introduction of light-matter interactions and proceeds with a description of many-body theory, including numerical methods used in that field. In the end of the theoretical part, the eigenvalue equation, used to calculate excitation energies, is derived. In the following section, the reduced-space algorithm is described. In the end of the thesis, numerical results obtained with the algorithm are presented, including a small basis set and method study. The comparison with the existing implementation of the similar algorithm verified the successful implementation of the algorithm presented in this thesis.

(5)

(6)

Acknowledgments

First I would like to thank my family for their support! I am also very grateful to my examiner Patrick Norman for taking time to sit down with me and go through what I have learned to see If I had understood everything correctly, for answering my questions by trying to guide me towards the answer. Though it sometimes took me some time to get on the right track it was very giving and I learned a lot! A big thanks to Thomas Fransson for helping me with LaTex and explaining a lot of things for me! A really big thanks to my supervisor Joanna Kauczor for taking a lot of time going through the algorithm with me on a white board, step by step. For helping me with the code, always making me feel welcome to stop by and ask questions! And last but not least I would like to thank my second supervisor Dirk who helped debugging the code always offered his help and guidance.

(7)

(8)

Chapter 1 Introduction

Light-matter interaction leads to a variety of interesting phenomena which have been observed and studied in a wide range of scientific fields like biology, chemistry and physics. An example that is usually brought up in biology classes at school is photosynthesis. It is a fundamental process to life on Earth in which light hits a green leaf and the energy of the light is transformed into chemical energy in the form of glucose [1]. Another interesting example is laser cooling. This is a technique used to cool molecules to a temperature close to absolute zero

by using one or several lasers. A simple description of this process is that a

molecule absorbs a photon and one or several phonons which lead to excitation of an electron. Afterward the electron returns to the ground state by emitting one or several photons with a total energy that is higher than the incident photon. This way a molecule can lose its thermal energy [2].

Various methods can be used to measure or calculate how light interacts with matter. One can calculate a variety of properties like the polarizability (ability for a molecule to be polarized) and excitation energies. To obtain excitation energies experimentally, one can either consider the emission or absorption spectrum of the molecule. The emission spectrum is created when electrons in a molecule jump from a higher to a lower state which results in an emission of photons. The light that is emitted from the molecule consists only of certain wavelengths and from these wavelengths it is possible to obtain the excitation energies of the molecule. On the other hand the absorption spectrum shows which frequencies a molecule can absorb. Similar to the case above, the specific frequencies that are absorbed provide information about the excitation energies. One of the methods used for measuring emission spectra of a molecule is the fluorescence spectroscopy. In flu-orescence spectroscopy, a beam of light is used to excite electrons in a molecule which leads to an emission of light from the molecule. To obtain an absorption spectrum, a microscope in transmission mode can be used. The process is the fol-lowing: light passes through the sample and hits a detector that gives information about which wavelengths have been absorbed. Although there are numerous ways to obtain molecular properties experimentally, in some cases it can be beneficial to calculate the desired properties numerically. For example if the system that is

(11)

measured is very unstable there may not be enough time to do the measurements. Therefore, it is important to have accurate mathematical models and numerical methods that can be used to calculate various properties of light-matter interaction in a molecule.

The excitation energies can be obtained by solving the Schr¨odinger equation.

The problem quickly becomes more complex when the number of particles in-creases (many-body system). Therefore, it is important to have efficient numerical methods to deal with this issue.

A method that can be used to approximate the wave function and find the corresponding energy is the Hartree–Fock method (HF). In HF, a wave function is expressed by a Slater determinant whose elements are atomic orbitals [3]. To obtain the energies, the HF equation is to be solved. It is an approximation to the

Schr¨odinger equation but due to the simplifications it is easier to solve.

Another method used to model a many-body system is the density functional theory (DFT). The main idea of DFT is that an electronic density can be associated with a specific ground state energy. This is very useful for large systems since the complexity of the electron density does not depend on the number of electrons.

Response theory is used to see how the electronic structure and other properties of the molecule are affected by an external field. The first step is to introduce a perturbed Hamiltonian and then use expansions of energy or other observables to obtain the response equations. From the response equations different properties are obtained. This thesis is about solving the linear response equation for a special case that gives the excitation energies. This allows for obtaining molecular excitation spectrum.

(12)

Chapter 2 Electronic structure theory

2.1 The Pauli exclusion principle

The Pauli exclusion principle states that no more than one electron can occupy the same quantum state. Each quantum state is characterized by four quantum

numbers n, l, ml and ms. Each orbital is characterized by a set consisting of the

first three quantum numbers. This means that two electrons in the same orbital

have to have different spin – ms. This requirement is satisfied when the wave

functions are antisymmetric. By the antisymmetry it is meant that the total wave function of a system has to be antisymmetric with respect to interchanging two electrons. When two electrons are interchanged, the sign of the wave function changes as follows [4]:

Ψ(¯n1, ¯n2, . . .) = −Ψ(¯n2, ¯n1, . . .), (2.1)

2.2 The Born–Oppenheimer approximation

The basic idea of the Born–Oppenheimer (BO) approximation is to take advantage of the fact that the electrons in a system are moving much faster than the nuclei, because the nuclei is massive compared to the electrons. The electrons can be viewed as moving in a potential field created by fixed nuclei. This greatly simplifies

solving the non-relativistic time independent Schr¨odinger equation (SE):

ˆ

H |Ψi = E |Ψi , (2.2)

where ˆH is the Hamiltonian operator. For a static system, when the

Hamilto-nian operator acts on an arbitrary wave function, the energy is obtained. The Hamiltonian for a system of N electrons and M nuclei is given as:

ˆ H = − N X i=1 1 2∇ 2 i− M X A=1 1 2MA ∇2 A− N X i=1 M X A=1 ZA riA + N X i=1 N X j>i 1 rij + M X A=1 M X B>A ZAZB RAB , (2.3) 3

(13)

where RA and ri are the position vectors of nuclei and electrons, respectively and

RABis the distance between the A−th and the B−th nucleus. Similarly, riAis the

distance between the i−th electron and the A−th nucleus and Rij is the distance

between two electrons, the i−th and the j−th. ZA is the atomic number of the

A−th nucleus and MAis the ratio of the mass of the A−th nucleus to the electron

mass. ∇2

i and ∇2A are Laplacian operators that define differentiation with respect

to the coordinates of the i−th electron and the A−th nuclei.

The first and second terms in the equation are the operators for the kinetic energy of the electron and the nuclei, respectively. The third term stands for the Coulomb attraction between the nuclei and the electrons. The last two terms represent the electron-electron and the nuclei-nuclei repulsions.

It is possible to neglect the second term and assume that the last term is constant. The constant term can be neglected, because addition of a constant term to an operator does not affect its eigenfunctions, but it correspondingly shifts the operator eigenvalues. The remaining term is the electronic Hamiltonian:

ˆ Helec= − N X i=1 1 2∇ 2 i − N X i=1 M X A=1 ZA riA + N X i=1 N X j>i 1 rij . (2.4)

The electronic Hamiltonian describes the motion of N electrons in a field created by M nuclei, that are viewed as point charges.

Rewriting the SE for the electronic Hamiltonian: ˆ

Helec|Ψeleci = Eelec|Ψeleci . (2.5)

The electronic wave function depends explicitly on the electronic coordinates and parametrically on the nuclear coordinates. The electronic energy depends parametrically on the nuclear coordinates. Each eigenstate gives rise to a so-called Born-Oppenheimer surface [3]. In the figure below one can see an example of two BO surfaces. They correspond to the ground state and the first excited state plotted for a water molecule under symmetrical stretching. The energies are plotted against the distance between the atoms (R) in the water molecule.

One can see that in a region around the equilibrium separation (ca 0.95 ˚A)

there is a large separation between the two states and the Born-Oppenheimer approximation works very well in this region [5].

(14)

2.3 Slater determinants 5

0.5

1.0

1.5

2.0 R (Angstrom)

0

2

4

6

8

10

12

14 Energy (eV)

First excited state

Ground state

Figure 2.1. The ground (blue) and the first excited state (green) energy vs the distance between the atoms in a water molecule.

2.3 Slater determinants

The Slater determinants represent the wave function of a multi-electron system [3]. For an N particle system, it is defined as:

|Ψi = Ψ(¯x1, ¯x2, . . . , ¯xN) = 1 √ N ! χ1(¯x1) χ2(¯x1) · · · χN(¯x1) χ1(¯x2) χ2(¯x2) · · · χN(¯x2) .. . ... . .. ... χ1(¯xN) χ2(¯xN) · · · χN(¯xN) (2.6)

Switching the coordinates of two particles corresponds to switching the two rows in the determinant that corresponds to those particles. Interchange of any two rows of the determinant results in the change of sign. This ensures that the Pauli principle is satisfied.

2.4 Variational principle

The variational principle is used to calculate the approximate energy for the ground and excited states with approximate wave functions [6]. The chosen wave

(15)

func-tions depend on one or several parameters. The parameters are fixed so that the expectation value of the energy is minimized.

Assuming that the SE is to be solved and ψ is a normalized trial solution, the following can be stated:

E0≤

D

ψ| ˆH|ψE, (2.7)

where E0is the ground state energy. The wave function is varied until the

expec-tation value of the Hamiltonian is minimized and an approximation to the ground

state is obtained [7]. When the proper wave function for the ground state ψ0 has

been found, it can be used to find the first excited state. This can be done by

choosing a new trial function ψ orthogonal to ψ0 (hψ|ψ0i = 0). Now the new

equation looks as follows:

E1≤

D

ψ| ˆH|ψE. (2.8)

As previously, the wave function is chosen so that it minimizes the expectation value of the Hamiltonian. This process can be repeated to calculate subsequent excited states. It should be mentioned that not all systems satisfy the variational principle. It is advantageous to be able to use a variational method because then the approximate energy is approached from above, and the lowest energy is always the best approximation.

(16)

2.5 The Hartree–Fock method 7

2.5 The Hartree–Fock method

Hartree–Fock is a method used in computational physics to determine the energy and the wave function of a system. This method is built on three important approximations, namely, the Born–Oppenheimer approximation, the Hartree–Fock approximation and the Linear Combination of Atomic Orbitals approximation (LCAO) [3]. In the Hartree–Fock approximation the ground state of an N -electron system is described by the Slater determinant:

|Ψeleci = |χ1χ2...χNi . (2.9)

From the variational principle it follows that the optimal wave function of the

functional in equation (2.7) is the one that yields the lowest energy (E0). In

equation (2.10) ˆH is the full electronic Hamiltonian.

EHF =DΨ0 ˆ H Ψ0 E . (2.10)

After choosing spin orbitals that will be included in the Slater determinant in equation (2.6), the energy in equation (2.7) is minimized with respect to these spin orbitals. This yields the Hartree–Fock equation:

f (i)χ(i) = χ(i), (2.11)

where f (i) is the Fock operator,

f (i) = − N X i=1 1 2∇ 2 i − N X i=1 M X A=1 ZA riA + vHF(i), (2.12)

vHF_{(i) is the average potential that the i th electron is experiencing in the presence}

of other electrons and therefore depends on the spin orbitals of those electrons. Since the Fock operator is a one electron operator, the many-electron problem reduces to a single-electron problem (The latter is much easier to solve.) [3]. The HF equation is non-linear and must therefore be solved iteratively.

2.5.1 Self consistent field procedure

The Self Consistent Field (SCF) procedure is iterative and is composed of five main steps, which are presented in Figure 2.2. In the first step an initial guess of the spin orbitals is made. The next step is to calculate the average potential

vHF_{(i) within which the eigenvalue equation is solved. When solving the}

eigen-value equation, a new set of spin orbitals is obtained. This set is used to calculate a new average potential. This is repeated until the difference in the potential between two successive iterations is smaller than a chosen threshold value. Then the orbitals that are used to calculate the Fock operator will be the same as its eigenfunctions [3].

The main advantage with SCF is that it is cheap. SCF is a good starting point for more advanced methods. The downsides are that it is accurate only when used on local interactions and does not include electron correlation. This leads to a

(17)

Figure 2.2. SCF procedure.

difference between the true total energy and the HF energy. This difference is called correlation energy.

(18)

2.6 Post-Hartree–Fock 9

2.5.2 Excited determinants

Excited Slater determinants are used in Post-Hartree–Fock methods which are improvements to the HF method. For calculating the excited determinants, the HF ground state is used as a reference [3]:

|Ψi = Ψ( ¯χ1, ¯χ2, . . . , ¯χa, ¯χb, . . . , ¯χN). (2.13)

When an electron is excited from the occupied orbital χato a virtual orbital χr

the occupied orbital is replaced by the virtual in the determinant and we obtain a singly excited determinant:

|Ψi = Ψ( ¯χ1, ¯χ2, . . . , ¯χr, ¯χb, . . . , ¯χN). (2.14)

If two electrons transfers from the occupied orbitals χa and χb to the virtual

orbitals χrand χswe will get a doubly excited determinant:

|Ψi = Ψ( ¯χ1, ¯χ2, . . . , ¯χr, ¯χs, . . . , ¯χN). (2.15)

Excited Slater determinants are used in Post-HF which will be described in the next section.

2.6 Post-Hartree–Fock

Post-HF methods are based on HF. They introduce some improvements to the

original method. In the original HF only one Slater determinant is used and

there is no correlation interaction. To obtain correlation interaction, more Slater determinants are needed. This is implemented in Configuration Interaction (CI) [8] and Coupled Cluster (CC) [9].

The CI wave function is a linear combination of Slater determinants:

ΨCI= a0ΨHF + X S aSΨS+ X D aDΨD+ ... = N X i=0 aiΨi. (2.16)

The trial wave function ΨCI is a combination of Slater determinants. The

subscripts S and D stand for second and double excitation. The expansion coef-ficients corresponding to the determinants have to be chosen so that the energy is minimized. To save computation time the CI is truncated to only allow certain excitations. The different CI methods are characterized by the allowed excitations. For example, in CISD only single and double excitations are allowed. With the CI procedure, an eigenvalue matrix equation is obtained:

Ha = ESa, (2.17)

where a is the coefficient vector and E is the eigenvalue matrix. H is the hamil-tonian matrix and S is the overlap matrix. Their elements are:

(19)

Si,j= hΨi|Ψji = δi,j. (2.19)

By solving the eigenvalue equation, we obtain the eigenvalues that correspond to the energies of the ground and excited states. One of the downsides is that the method is very time consuming and is only suitable for small systems. Another drawback of CI is that it does not exhibit the property called size extensivity. Size extensivity means that the energy of two infinitely separated atoms is twice as big as the energy of one single atom of the same sort. To keep size extensivity, an alternative approach is to use a non-linear expansion of Slater determinants. A method of doing this is CC, where the wave function is:

|ΨCCi = e

ˆ T_|Ψ

0i (2.20)

there Ψ0 is the reference wave function and is usually chosen to be the HF wave

function. ˆT is defined as:

ˆ

T = ˆT1+ ˆT2+ ... + ˆTN, (2.21)

ˆ

T1 operating on the wave function gives the singly excited determinant and ˆT2

gives the doubly excited determinant. A Taylor expansion of the operator eTˆ

gives mixed terms:

The CC wavefunction is also truncated, and different methods are characterized by the allowed excitations. In CCSD, only single and double excitations are involved. The CC eigenvalue equation is:

HeT|ΨHFi = EeT|ΨHFi, (2.23)

In practice, using this equation is unlivable and hΨHF|e−

1

T projected on the

equa-tion. The energies of the ground and excited states are obtained as eigenvalues of this equation. Comparing CC and CI, CCSD performs better than CISD. CC also has size extensivity but like CI, it is very time consuming and can not be used on large systems.

2.7 Density functional theory

The Density Functional theory is based on the theorems presented by Hohenberg and Kohn. The first theorem states that electronic density determines the ground state of the electronic energy [10]. The second theorem states that, given a func-tional, the correct electron density of the ground state minimizes the energy of the functional. The use of the electronic density to determine the energy of the ground state is very helpful for calculations made on larger systems. Since a wave function depends on the coordinates and the spin of the electrons in the system, the wave function quickly becomes very complex with the increasing number of

(20)

2.7 Density functional theory 11 electrons in the system. This problem does not occur for the electron density of the system which does not depend on the number of electrons. So, instead of a many-body system with interacting particles in an external field, one can consider non-interacting particles in an effective field. The effective field consists of the ex-ternal field and exchange and correlation interactions between the electrons. The exchange and correlation interaction have been proven to be difficult to model, which is solved to some extent by certain functionals.

Some of the frequently used functionals are Local Density Approximation (LDA), Generalized Gradient Approximation (GGA) and Local Spin Density Ap-proximation (LSDA) [11]. LDA is the simpler functional, which contains

approx-imations to the exchange-correlation energy functional. These approximations

depend only on the local electronic density:

E_XCLDA(n) =

Z

XC(n)n(r)d3r. (2.24)

An improvement to this approximation can be made by taking into account the electron spin, which is done for the LSDA functional:

EXCLSDA(n↑, n↓) =

Z

XC(n↑, n↓)n(r)d3r. (2.25)

The extension of the parameters by adding the electronic density gradient gives the GGA functional [12]:

E_XCGGA(n↑, n↓) =

Z

XC(n↑, n↓, ∇n↑, ∇n↓)n(r)d3r. (2.26)

LDA and GGA work well in most cases except for weakly bound systems. Further improvement can be made by implementing hybrid functionals which combine correlation exchange from HF theory with other sources. Some commonly used hybrid functionals are B3LYP, CAMB3LYP, PBE0 and PW91. The most popular of the three is B3LYP [13]:

EXCB3LY P = a0LDA + 0.2(ExHF − E LDA x ) + ax(ExGGA− E LDA x ) + E LDA c + ac(EcGGA− E LDA c ). (2.27)

where a0, ax and ac are usually chosen to be 0.2, 0.72 and 0.81, to get a better

agreement with experiments. CAMB3LYP introduces an improvement to the long range exchange function in B3LYP [14]. PBE0 combines HF and PBE with a ratio 1:3 [15]: E_XCP BE0= 1 4E HF x + 3 4E P BE x + E P BE c . (2.28)

In general DFT has a good balance between price and accuracy. The disadvantage is that the use of this method requires experience and there is no systematic way of improving it.

(21)

(22)

Chapter 3 Response theory

3.1 Molecular properties and perturbation theory

In response theory, the effect of a perturbation on a molecule is analyzed. In general, the perturbation can be something else than an electric field but for the purpose of this thesis we focus on external electric perturbations. When a molecule is subjected to an electric field, the Hamiltonian for the system becomes:

ˆ

H = ˆH0− ˆµF (t), (3.1)

where ˆH0 is the Hamiltonian of an unperturbed system, ˆµ is the electric dipole

operator and F (t) is the amplitude of the electric field.

There are two cases that should be considered. In the first case, the electric field is constant and in the second case the electric field is oscillating. In the first case, called the static case, the energy levels are well defined. In the second case, however, the electric field is oscillating and the energy levels are ill-defined, which makes the problem more complicated.

3.1.1 Static external field

For the static case, the external field F (t) is constant. When the electric field does not depend on time, one can obtain the energy levels for the molecule by solving the time-independent SE. The energy of the system can also be obtained by a Taylor expansion around zero perturbation strength, which works for small perturbations [16]: E = E0+ ∂E ∂Fα _F α=0 Fα+ 1 2 ∂2_E ∂Fα∂Fβ _F α=Fβ=0 FαFβ+ . . . . (3.2)

Various molecular properties can be derived from this expansion. The first deriva-tive of the energy with respect to the electric field corresponds to the electric dipole moment. The polarisability of the molecule can be obtained from the sec-ond derivative. The first- and the secsec-ond-order hyperpolarisability can be obtained by considering the third and the forth derivatives, respectively.

(23)

3.1.2 Time-dependent external field

Since most of experiments are performed with an oscillating electric field, the second case becomes more interesting. We assume that the perturbation is slowly turned on at t = 0:

F (t) = Fωsin(ωt) × erf(at), (3.3)

where Fω_{is the electronic field amplitude, ω is the frequency of the oscillation, a}

is a constant and erf(at) is an error function that ensures that the field is slowly turned on at t = 0.

The perturbation operator has the following form: ˆ Vα(t) = X ω ˆ V_αωF_αωe−iωtet, (3.4)

where α is the molecular axis, is a positive infinitesimal and ˆVω

α is an operator.

Since the perturbation has to be real, the conditions F_α−ω = [F_α−ω]∗ and ˆV_α−ω =

[ ˆV_α−ω]† = ˆVω

α should be satisfied. When the electric field becomes time dependent

the molecular properties will also vary with time. Like, for example, in the case of polarization:

µ(t) = µ0+ αF (t), (3.5)

where µ0 is the permanent polarizability and the second term comes from the

oscillations of the electric field. The Taylor expansion of the polarizability gives:

µ(t) = µ0+ αF (t) + 1

2βF

2_{(t) +}1

6γF

3_{(t) + ...,} _(3.6)

where β and γ are the first- and second-order hyperpolarizability, respectively. In general a molecular property can be obtained as an expectation value of an

operator ˆΩ:

hψ(t)| ˆΩ|ψ(t)i =hψ(0)| ˆΩ|ψ(0)i

+ hψ(1)| ˆΩ|ψ(0)i + hψ(0)| ˆΩ|ψ(1)i (3.7)

+ hψ(2)| ˆΩ|ψ(0)i + hψ(1)_{| ˆ}_Ω|ψ(1)_{i + hψ}(0)_{| ˆ}_Ω|ψ(2)_i

+ . . . ,

By inserting the perturbation operator in this equation, we obtain:

hψ(t)| ˆΩ|ψ(t)i =h0| ˆΩ|0i +X ω1 hh ˆΩ; ˆVω1 β iiF ω1 β e −iω1t_et _(3.8) +1 2 X ω1ω2 hh ˆΩ; ˆVω1 β , ˆV ω2 γ iiF ω1 β F ω2 γ e−i(ω1 +ω2)t_e2t + . . . ,

(24)

3.1 Molecular properties and perturbation theory 15 where hh ˆΩ; ˆVω1 β ii, hh ˆΩ; ˆV ω1 β , ˆV ω2

γ ii are the linear and the quadratic response

func-tions, respectively.

3.1.3 Reference parametrization

The form of the linear response function can be deduced by considering the ref-erence parametrization. When solving the SE, |ψi can be parametrized, i.e. ex-panded in a set of known functions. A common parametrization is:

|ψ(t)i =X

n

cn(t)|ψni, (3.9)

where

cn= hn|ψ(t)i. (3.10)

For the time evolution of the parameters cn, the parametrization can be rewritten:

|ψ(t)i =X

n

dn(t)e

−iEnt

~ |ni, (3.11)

The new parametrization is motivated by the fact that, in absence of a perturbation

field, cn(t) = cn(0)e

−iEnt ~ .

Substituting the new parametrization into the perturbed SE and projecting hm| on it, gives: i~_δtδdm(t) = X n Vmn(t)dn(t)e −iEnt ~ (3.12)

where Vmn are elements of the perturbation operator ˆV , given by:

Vmn= hm|

X

ω

ˆ

V_αωF_αωe−iωtet|ni. (3.13)

The parameters in equation (3.11) can be expanded in the order of perturbation:

dn(t) = d(0)n + d

(1)

n (t) + d

(2)

n (t) + . . . , (3.14)

The N th order solution to the differential equation is:

d(N )_m (t) = 1 i~ t Z −∞ X n Vmn(t0)d(N −1)n (t 0_)e(Em−En)t0/~_dt0_, _(3.15)

and the zeroth order solution can be deduced as:

d(1)_m(t) = 1 i~ t Z −∞ X n Vmn(t0)d0ne(Em−En)t 0 /~_dt0 = 1 i~ t Z −∞ X n hm|X ω1 ˆ Vω1 α F ω1 α e−iω1 t0_et0_|niδ 0neiωmnt 0 /~_dt0 = −1 ~ X ω1 hm| ˆVω1 α |0iFαω1 ωm0− ω1− i ei(ωm0−ω1)t_et , (3.16)

(25)

where d(0)n = δ0n is used as the initial condition and the transition angular

fre-quency (Em− En)t0/~ is denoted by ωmn.

An expansion in the order of perturbation is now performed for the wavefunc-tion:

|ψ(N )_{(t)i = |ψ}(0)_{(t)i + |ψ}(1)_{(t)i + |ψ}(2)_{(t)i + . . . ,} _(3.17)

where the N th order wavefunction is obtained by combining equations (3.11) and (3.14):

|ψ(N )_{(t)i =}X

n

d(N )_n (t)e−iEnt~ |ni, (3.18)

Equations (3.16) and (3.18) together with (3.7) will now be used to deduce the linear response equation:

h ˆΩi(1)= hψ(1)| ˆΩ|ψ(0)i + hψ(0)_{| ˆ}_Ω|ψ(1)_i = − h0|eiE0t/~_ΩˆX n 1 ~ X ω1 hn| ˆVω1 α |0iFαω1 ωn0− ω1− i

ei(ωn0−ω1)t_et_e−iEnt/~_|ni

−X n 1 ~ X ω1 hn| ˆVω1 α |0i[Fαω1]∗ ωn0− ω1+ i

e−i(ωn0−ω1)t_et_eiEnt/~_{hn| ˆ}_Ωe−iE0t/~_|0i

where h ˆΩi(1) _{is the first order correction to the expectation value in equation}

(3.7). Comparing this to equation (3.8), we can conclude that the sum over state expression of the linear response function is:

where the term in the perturbation operator is chosen so that i = 0.

3.1.4 Physical properties

As mentioned above it is possible to retrieve physical properties from the response functions. For example, combining the expression for the linear response function together with equation (3.8) one can obtain the expression for the polarizability tensor: ααβ= −hhµα; µbiiωb = − X n (µ 0m α µm0m ωb− ωm − µ 0m β µ m0 α ωb+ ωm ) (3.21)

(26)

3.2 Response equation 17 By considering the residues of the poles in the this equation, one can obtain the transition moments from ground state to an excited state |mi:

lim ωb→ωm (ωb− ωm)hhµα; µβiiωb= µ 0m α µ m0 β (3.22) lim ωb→−ωm (ωb+ ωm)hhµα; µβiiωb= −µ m0 α µ 0m β (3.23)

Apart from the transition moments, in general, various properties can be calculated from non-linear response functions expressed in order of perturbation [17]:

Hyperpolarizability Physical effect

α(0;0) static polarizability

α (−ω;ω) frequency-dependent polarizability

β (0;0,0) static first hyperpolarizability

β (-2ω;ω,ω) second harmonic generation

β (-ω,ω,0) dc-Pockels effect

β(0;−ω,ω) optical rectification

γ (0;0,0,0) static scond hyper polarizability

γ(-3ω;ω,ω,ω) third harmonic generation

γ(-2ω;ω,ω,0) dc-second harmonic generation; electric field induced

SHG

γ(-ω;ω,−ω,ω) intensity dependent refractive index; degenerate four

wave mixing

γ(−ω1;ω1,-ω2,ω2) ac Kerr effect; optical Kerr effect

γ(-ω;ω,0,0) dc Kerr effect; electro-optical Kerr effect

γ(0,ω;-ω,0) dc optical rectification; electric field induced optical

rectification

3.2 Response equation

The sum over state expression of the linear response function can be written as a matrix equation: hh Â; ˆBii = −X n (h0| Â|nihn| ˆB|0i En0− ω +h0| ˆB|nihn| Â|0i En0+ ω ) = A[1](E[2]− ωS[2])−1B[1] (3.24) where E[2]=           E1 . ._.

₀

EN E1

0

. ._. EN           (3.25) and S[2]= 1 0 0 −1

(27)

A and B in equation (3.24) are known, but the matrix inverse is hard to compute

in practice due to large dimensions. Therefore, it is beneficial to set x = (E[2]−

ωS[2])−1B[1] to obtain an equation of the form Ax = b. Motivated by this, the

response equation is rewritten in the matrix form:

(E[2]− ωbS[2])xb = gb, (3.26)

where gb _{is the generalized property gradient and x}b _{is the response vector. In}

the non-exact perturbation theory, E[2]is no longer diagonal. In SCF theory and

DFT E[2] has a block structure [18]:

E[2]=

A B

B∗ A∗

E[2] can be expressed as a sum of E[0] and the exchange correlation E[1], where

E[0] is a diagonal matrix containing the orbital energy differences. To obtain the

excitation energies one can solve the eigenvalue equation:

E[2]xb= ωbS[2]xb, (3.27)

where ωb is the excitation energy, and xb is the corresponding transition moment.

A detailed description of how this equation is solved in practice will be given in the next chapter.

(28)

Chapter 4 Algorithm

In this section, we describe the algorithm used to solve the response equation:

(E[2]₀ − ωbS[2])xb= gb. (4.1)

4.1 General case

The response equation can be written in a more general form:

Ax = b. (4.2)

While this equation may look reasonably easy to solve, this is not the case. Due

to the size of A it is not possible to calculate A−1. Instead this equation has to

be solved iteratively in a reduced space, as will be described below. As the result we get an approximate solution.

After n iterations the subspaces xn _{and σ}n _{consist of:}

x(n)= {x1, x2, ..., xn}. (4.3)

σ(n)= {σ1, σ2, ..., σn}. (4.4)

where σi= Axi. x1is given in a way specific to the problem and xnis generated by

the algorithm recursively. One iteration of the algorithm consists of the following steps:

i) Equation (4.2) is solved in a reduced space for the coefficients (α):

Aredαred= bred, (4.5)

where

[Ared]i,j= xTiσj, (4.6)

and

[bred]i= xTib. (4.7)

ii) The optimal vector is constructed: 19

(29)

˜ Xopt = n X i=1 αixi. (4.8)

The optimal vector is the solution vector that is the best in a certain iteration. iii) The residual is constructed

r = b − A ˜Xopt= b −

n

X

i=1

αiσi. (4.9)

The residual is the difference between the right-hand side and the the left-hand side where the exact solution has been replaced by the optimal vector. This residual becomes smaller and smaller in each iteration when the optimal vector comes closer and closer to the correct solution. The last step in equation (4.9) can be justified

by keeping in mind that σi= Axi:

A ˜Xopt= A(

n

X

i=1

αixi) = Aα1x1+ Aα2x2+ ... + Aαnxn=

n

X

i=1

αiσi. (4.10)

The problem is considered to be solved to a desired accuracy, when the norm of the obtained residual becomes smaller than a given threshold value. If this holds the algorithm terminates. Otherwise another iteration has to be performed and a new vector has to be added to the subspace in equation (4.3). There are two ways of choosing the new vector. One can simply use the latest residual as the new vector. The downside with this method is that it has a low convergence rate. The other method is to precondition the residual before it is added to the subspace. Here the latter method will always be used.

iv) The vector xn+1is calculated:

bn+1= P(rn). (4.11)

where

P(rn) = (A0)−1rn. (4.12)

A0 is an approximation of A which is easier to invert. For example A0 could

be the diagonal of A. Before adding the new vector to the subspace it has to be orthonormalized against the vectors in the subspace in equation (4.3). This is done by the Gram–Schmidt process:

xn+1= bn+1− n X i=1 hbn+1, xji hxj, xji xj. (4.13)

When the new vector xn+1 is obtained, σn+1 is calculated and they are both

added to the subspaces xn and σn. All these steps are repeated until convergence

(30)

4.1 General case 21

Figure 4.1. Flowchart of the different steps in the algorithm.

To gain a better understanding of the steps of the algorithm, the first two iterations are presented in the subsections below.

First iteration

The start vector x1is added to the subspace {x}. The corresponding σ is obtained

by the linear transformation :

σ1= A x1. (4.14)

It is necessary to clarify that for getting σ1, A is not multiplied by x1, but is

obtained by calculating the perturbed Fock matrix from the density matrix [18]. Equation (4.2) is solved in the reduced space (i):

(31)

Aredxred= bred. (4.15)

where in this iteration

Ared= x1σ1. (4.16)

and

bred= x1b. (4.17)

In step (ii) the optimal vector is calculated: ˜ Xopt,1= n X i=1 αixi= α1x1. (4.18)

The residual can now be obtained (iii) in accordance with:

r1= b − A ˜Xopt,1= b − α1σ1. (4.19)

After this step the residual is usually still too large compared to the threshold

value. This means that in step (iv) another vector is added to xn_:

x2= P(r1) = A−10 r1. (4.20)

Before adding x2to the subspace, it has to be orthogonalized against x1and then

normalized. After that the second σ is constructed:

σ2= Ax2. (4.21)

Now x2 and σ2 can be added to the corresponding subspaces.

Second iteration

In the second iteration our subspaces look the following:

x(2)= {x1, x2}. (4.22)

σ(2)= {σ1, σ2}. (4.23)

In step (i) the reduced equation:

xT 1σ1 xT1σ2 xT₂σ1 xT2σ2 α1 α2 =x T 1b xT₂b (4.24) is solved. The optimal vector is generated in step (ii) by the formula:

˜ Xopt,2= n X i=1 αixi= α1x1+ α2x2. (4.25)

The residual is calculated as step (iii):

r2= b − A ˜Xopt,2= b − α1σ1− α2σ2. (4.26)

The residual is compared to the threshold value and if it is small enough the algorithm terminates here, otherwise another iteration will be needed.

(32)

4.2 Solving the response equation 23

4.1.1 Eigenvalue equation in the general case

In this section we consider the eigenvalue equation:

Ax = λx. (4.27)

As in the general case, after n iterations we have the subspaces:

xn= {x1, x2, ..., xn}. (4.28)

σn= {σ1, σ2, ..., σn}. (4.29)

The eigenvalue equation is solved in the reduced space both for αredand λ:

Aredαred= λnαred. (4.30)

where Ared is given by equation (4.6). The optimal vector is found in the same

way as previously: ˜ Xopt,n= n X i=1 αixi. (4.31)

Since we have a different type of equation, the residual is constructed in a different way, namely, as:

r = A ˜Xopt,n− λnX˜opt,n= n X i=1 αiσi− λn n X i=1 αixi. (4.32)

As previously, if the residual is larger than a preset threshold value, a new is added to the basis.

4.2 Solving the response equation

After considering the general case, we focus on a more specific case of response equation. In this case, equation (4.2) turns into equation (3.26):

(E[2]− ωS[2])x = G. (4.33)

The algorithm has the same structure as in the previous case but with minor changes. One of the changes is that the trial vectors that were previously denoted ny x are now replaced by b as in the literature. The iteration begins by adding a

vector bi to the subspace {b}. When a vector is added it has to be orthonormal

(normalized and orthogonal) to the subspace. The next step is to calculate the

corresponding σ and ρ denoted σi and ρi [18]:

σi= E[2]bi. (4.34)

ρi= S[2]bi. (4.35)

(33)

(Ered− ωSred)xred= Gred, (4.36) where [Ered]i,j = bTiE [2]_b j= bTi σj, (4.37) [Sred]i,j= bTiS[2]bj= bTi ρj, (4.38) [Gred]i = bTiG. (4.39)

The following step in the iteration is to calculate the optimal vector ˜Xopt:

˜ Xopt= n X i=1 xredbi. (4.40)

The residual is constructed as: r = G − ( n X i=1 αiσi− ω n X i=1 αiρi). (4.41)

If the residual is smaller than a given threshold value the algorithm terminates.

Otherwise, a new vector bi+1 is added to {b}:

bi+1= P(r), (4.42)

where P is a preconditioning matrix:

P = (E0− ωS[2])−1. (4.43)

Here E0is a diagonal matrix containing the orbital energy differences. These steps

are repeated until convergence is obtained.

4.3 Solving the response equation for arbitrary

number of frequencies

It is beneficial to solve the response equation for several frequencies at the same time since for close-lying frequencies the solution can be found in the same reduced

space. In the beginning of the iteration, new vectors bi are added to the vector

space {b}. For each vector bi that is added to {b}, the corresponding vectors σi

and ρiare calculated [18]. The σiand the ρivectors are obtained in the same way

as in equations (4.34) and (4.35). Now ω is a vector containing N frequencies. The reduced equation is solved for each of the N frequencies to obtain the coefficient vectors:

(Ered− ω1Sred)αred,1= G

(Ered− ω2Sred)αred,2= G

.. . (Ered− ωNSred)αred,N= G

(34)

4.3 Solving the response equation for arbitrary number of frequencies 25 Where the coefficient vectors are of the same dimension as the number of vectors

in {b}, N. After obtaining the coefficient vectors the optimal vectors ˜Xopt(w) are

calculated for each frequency:

˜ Xopt(ω1) = n X i=1 αred,1(i)bi ˜ Xopt(ω2) = n X i=1 αred,2(i)bi .. . ˜ Xopt(ωN) = n X i=1 αred,N(i)bi (4.45)

The next step is to calculate the residuals r(w):

rω1= G − ( n X i=1 αred,1(i)σi− ω1 n X i=1 αred,1(i)ρi) rω2= G − ( n X i=1 αred,2(i)σi− ω2 n X i=1 αred,2(i)ρi) .. . rωN = G − ( n X i=1 αred,N(i)σi− ωN n X i=1 αred,N(i)ρi). (4.46)

The norms of the residuals are then calculated and compared to a given thresh-old value. If a residual corresponding to a certain frequency is larger than this threshold value, the iterative process will continue for this frequency. If, however, a residual is small enough then in the subsequent iteration the optimal vector corresponding to this frequency will not be changed.

In the next iteration, N new vectors may be added to the reduced space. The potential vectors are obtained by preconditioning the residuals :

b1= Pω1(r(w1)) b2= Pω2(r(w2)) .. . bN = PωN(r(wN)) (4.47)

(35)

Pω1 = (E0− ω1S [2]₎−1 Pω2 = (E0− ω2S [2]₎−1 .. . PωN = (E0− ωNS [2] )−1 (4.48)

The new vectors are orthogonalized against the basis vectors. If the part of a vector that is orthogonal to the basis vectors is big enough it will be added to the basis. This is done subsequently until all the residuals are below the threshold value.

4.3.1 Eigenvalue equation in the special case

We consider the eigenvalue equation:

E[2]X = ωS[2]X. (4.49)

It is solved in the reduced space for the N frequencies:

Eredαred,1= ωSredαred,1. (4.50)

The optimal vectors are obtained in exactly the same way as in the previous case, namely, from equation (4.45):

˜ Xopt(ω) = n X i=1 αred,ibi. (4.51)

The next step is to calculate the residuals rwwhich will be different since the right

hand side is now equal to zero:

rω= E[2]X˜opt− ωS[2]X˜opt = ( n X i=1 αred,iσi− ω n X i=1 αred,iρi). (4.52)

The norms of the residuals are then calculated. If the algorithm did not converge, new vectors are added. The new vectors are calculated in the same way as in the previous case, see equation (4.48). A problem with this way of solving the eigenvalue equation is that complex eigenvalues may be obtained. This can be solved by adding symmetric vectors as described in the next section.

4.3.2 Splitting the trial vectors into symmetric and

anti-symmetric parts

In order to avoid complex eigenvalues, the vectors are split into a symmetric and an antisymmetric part. The steps in the iteration are the same as before but the difference is that the solution vector is now divided into a symmetric (gerade, denoted g) and an antisymmetric (ungerade, denoted u) part [19]:

(36)

4.3 Solving the response equation for arbitrary number of frequencies 27 The response equation is separated into a symmetric and an antisymmetric part,

keeping in mind that E[2]does not affect the symmetry while S[2]changes it. Thus,

we have:

E[2]bi,g= σg,i. (4.54)

E[2]bi,u= σu,i. (4.55)

S[2]bi,g= ρu,i. (4.56)

S[2]bi,u= ρg,i. (4.57)

So if S[2] acts on a symmetric vector, it becomes antisymmetric and vise versa.

The response equation is split into a symmetric and an antisymmetric equation:

E[2]Xg− ωS[2]Xu= Gg. (4.58)

−ωS[2]_X

g+ E[2]Xu= Gu. (4.59)

By combining these equations, we obtain the following matrix equation: E[2] −ωS[2] −ωS[2] _E[2] Xg Xu =Gg Gu . (4.60)

In the first step of the iteration, the trial vectors are split into a symmetric and an antisymmetric part:

b = bg+ bu (4.61)

If the symmetric part of a vector is too small to be added to the reduced space, the antisymmetric part is also discarded. The next step is to solve the reduced equation: Egg_red −ωSgu_red −ωSug_red Euu_red Xred,g Xred,u =Gred,g Gred,u (4.62) where

[Egg_red]i,j = bTi,gσj,g (4.63)

[Euu_red]i,j = bTi,uσj,u (4.64)

[Sgu_red]i,j = bTi,gρj,u (4.65)

[Sgu_red]i,j= [Sug_red]j,i (4.66)

Now the optimal vectors can be obtained as: ˜ Xg= X i [Xred] g ibi,g (4.67) ˜ Xu= X i [Xred]uibi,u (4.68)

Using these vectors the residuals are calculated as:

rg= Gg− ( X αi,gσi,g− ω X αi,gρi,g) (4.69) ru= Gu− ( X αi,uσi,u− ω X αi,uρi,u) (4.70)

As in the previous cases if the residual is small enough the algorithm terminates, otherwise, vectors are added.

(37)

4.3.3 Initial trial vectors

Initial trial vectors for solving the general response equation

The iterative algorithm for solving the response equation starts with producing proper initial vectors [18]. In the general case, the trial vectors are obtained in the following way: b1(ω1) = (E [2] 0 − ω1S[2])G b2(ω2) = (E [2] 0 − ω2S[2])G .. . bN(ωN) = (E [2] 0 − ωNS[2])G

This is a system of N equations for the N different frequencies. The first of the

generated vectors is b1, and after being normalized it generates the space of the

trial vectors {b}.

Initial trial vectors for solving the eigenvalue equation

It is a bit more complicated to find initial trial vectors for the eigenvalue equation,

so it is done in several steps. At the beginning, frequencies ωiare calculated:

ωi=

E[2]_0,i

abs(S[2]_i )

. (4.71)

A diagonal matrix consisting of the calculated frequencies is then created:    ω1 0 . ._. 0 ωN   . (4.72)

After this step the frequencies are sorted in decreasing order and the columns of the matrix are permuted accordingly, so that the largest frequency is in the first column and the smallest in the last. The frequencies are then replaced by ones and this final matrix is used to calculate the trial vectors. For simplicity this process is illustrated below for a four by four matrix:

    ω1 0 0 0 0 ω2 0 0 0 0 ω3 0 0 0 0 ω4     . (4.73)

Let us assume that the following inequality hold for the frequencies ω2 > ω4 >

ω1> ω3. The frequencies in the matrix will then be shifted by shifting the columns:

    0 0 ω1 0 ω2 0 0 0 0 0 0 ω3 0 ω4 0 0     . (4.74)

(38)

4.4 Eigenvalue equation with vector division 29 This matrix is then used to create a matrix with only ones and zeros. The ones being in the same places as the frequencies in the previous matrix:

    0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0     . (4.75)

This matrix is now used to create the trial vectors. The symmetric and the anti-symmetric trial vectors are:

bi,g= 1 2(vi+ v p i), (4.76) bi,u= 1 2(vi− v p i). (4.77)

Where vi is the i–th column in the matrix.

4.4 Eigenvalue equation with vector division

In the final stage, the eigenvalue equation is solved. In this section, it is presented in the same way as it is done numerically in this thesis. Similar to the section before, the solution vector will be divided into a symmetric and an antisymmetric part:

E[2](Xg+ Xu) = ωS[2](Xg+ Xu). (4.78)

The response equation can now be accordingly separated into a symmetric and an antisymmetric equation:

E[2]Xg= ωS[2]Xu, (4.79)

E[2]Xu= ωS[2]Xg. (4.80)

This gives the following matrix equation: E[2] −ωS[2] −ωS[2] E[2] Xg Xu =0 0 . (4.81)

In the first step of the iteration, all the trial vectors are split into symmetric and antisymmetric parts:

b = bg+ bu. (4.82)

As previously, if the symmetric or the antisymmetric component of a vector pair is too small, the vector pair will be neglected. The next step is to solve the reduced

equation for both Xi (i = g, u) and ω:

Egg_red −ωSgu_red −ωSug_red Euured Xred,g Xred,u =0 0 , (4.83)

(39)

where

[Egg_red] = bT_i,gσi,g, (4.84)

[Euured] = b

T

i,uσi,u, (4.85)

[Sgu_red] = bT_i,gρi,u, (4.86)

[Sgu_red]i,j= [S

ug

red]j,i. (4.87)

Then the optimal vectors can be obtained: ˜ Xg= X i [Xred]gibi,g, (4.88) ˜ Xu= X i [Xred]uibi,u. (4.89)

Using the optimal vectors, the residuals are calculated:

rg= E[2]X˜g− ωS[2]X˜u, (4.90)

ru= E[2]X˜u− ωS[2]X˜g. (4.91)

As before if the residuals are small enough the algorithm terminates.

4.4.1 First iteration

The first iteration will be shown for three frequencies. Since we have three fre-quencies, in each iteration the maximum number of vector pairs added is three:

bg₁, bu1, b g 2, b u 2, b g 3 and b u

3. Let us assume that b

g

1 and b

u

3 are too small. This

means that only bg₂ and bu2 will be added to the vector spaces bg and bu:

bg= {bg₂} (4.92)

bu= {bu2} (4.93)

Whenever new vectors are added the corresponding σ and ρ are calculated. The reduced space is now defined by:

Egg_red −ωSgu_red

−ωSug_red Euu_red

. (4.94)

The next step is to calculate the optimal vectors: ˜ Xg= X i [Xred] g ibi,g= σ g 2b g 2, (4.95) ˜ Xu= X i [Xred]uibi,u= σ2ub u 2. (4.96)

In the next step, the residuals will be calculated:

rg= E[2]X˜g− ω S[2]X˜u (4.97)

ru= E[2]X˜u− ωS[2]X˜g (4.98)

(40)

4.4 Eigenvalue equation with vector division 31 Second iteration

The vectors that could be added now are bg₄, bu₄, bg₅, bu₅, bg₆ and bu₆. Assuming

that bg₆ is too small, only bg₄, bu₄, bg₅ and bu₅ will be added:

bg= {bg₁, bg₂, bg₃}, (4.99) bu= {bu1, b u 2, b u 3}. (4.100)

The corresponding σ and ρ are calculated. The reduced space is now defined by:

Egg_red −ωSgu_red

−ωSug_red Euu_red

. (4.101)

The optimal vectors are obtained in the same way as previously: ˜ Xg= X i [Xred]gibi,g= σg1b g 1+ σ g 2b g 2+ σ g 3b g 3, (4.102) ˜ Xu= X i [Xred]uibi,u= σ1ub u 1+ σ u 2b u 2+ σ u 3b u 3. (4.103)

The next step is to calculate the residuals:

rg= E[2]X˜g− ωS[2]X˜u, (4.104)

ru= E[2]X˜u− ω S[2]X˜g. (4.105)

If the residuals are small enough, another iteration is not needed. Otherwise new vectors are to be added and another iteration will be performed.

(41)

(42)

Chapter 5 Numerical results

5.1 Evaluation of the algorithm

An algorithm for calculating excitation energies was already available in Dalton. In this thesis the same algorithm was implemented in Dalton using symmetric trial vectors (see section 4.3.2). The new implementation of the algorithm is motivated by the necessity to substantially improve the modular structure of the old one.

Calculations have been performed on a variety of molecules to make sure that the same excitation energies were obtained using both paired and symmetric vec-tors. Besides comparing the results, it was also important to compare the efficiency of the implementations. The convergence and the size of the reduced space of the two implementations were compared iteration by iteration, and it was concluded that they were also identical. This is due to the fact that the symmetric and the paired trial vectors span the same space. The vectors are basically the same, only rotated with respect to each other. In Table 5.1 the convergence of the algorithm is illustrated. The calculations were performed on ethene with CAMB3LYP and

10−5 a.u. as the threshold value. The table shows how the first three excitation

energies and their corresponding residuals evolve in the different iterations. In the

Figure 5.1. The structures of a) ethene b) benzene

(43)

Table 5.1. The eigenvalues, the residuals and the size of the reduced space are shown iteration by iteration. The calculations were performed on ethene using aug-cc-pVDZ and CAMB3LYP with 10−5a.u. as threshold. The energies are given in a.u.

Iter. E |Ru+ Rg| 1 0.26112177 0.4227155·10−1 0.28564442 0.4276224·10−1 0.28573595 0.4065476·10−1 2 0.25462927 0.14102·10−1 0.27713470 0.987249·10−2 0.27961106 0.211188·10−1 3 0.25444844 6.67672·10−3 0.27698169 3.66254·10−3 0.27944302 3.03620·10−3 4 0.25441237 2.72279·10−3 0.27694201 4.12168·10−3 0.27943546 6.93229·10−4 5 0.25440839 4.74864·10−4 0.27691170 3.02254·10−3 0.27943466 3.17348·10−4 6 0.25440826 9.29062·10−5 0.27690359 4.30645·10−4 0.27943453 9.14210·10−5 7 0.25440825 2.55844·10−5 0.27690350 3.88370·10−5 0.27943453 1.06030·10−5 8 0.25440825 3.21639·10−6 0.27690350 5.04346·10−6 0.27943453 2.88791·10−6

first iterations, the excitation energies change significantly, while in the following iterations, only the last decimals are changed. Another thing that is interesting to look at, is the size of the reduced space and the number of iterations it takes to reach convergence which is presented in Table 5.2. Table 5.2 shows that the number of iterations does not significantly change upon increasing the size of the molecule. This is an attractive result since it is important to have an algorithm that is effective for bigger systems.

(44)

5.1 Evaluation of the algorithm 35

Table 5.2. The size of the reduced space and the number of iterations required to converge is presented for different molecules. The calculations were performed with CAMB3LYP and the basis set aug-cc-pVDZ. The threshold was set to 10−5 a.u. and three excitation energies were requested.

Ng+ Nu Nit Ethene 48 8 Butadiene 72 17 Benzene 60 12 Cyclopentadiene 60 10 Furan 56 11 Formaldehyde 42 8

(45)

5.2 Choice of basis set study

In this section we consider testing the algorithm with different basis sets. The functions that are commonly used are Gaussian functions:

ΦGT O_abc = N xaybzce−ζr2, (5.1)

where N is a normalization constant, a, b and c specify the angular momentum, ζ controls the orbital width and ’GTO’ stands for ’Gaussian type orbitals’. Ideally, basis sets should be infinite but that is computationally impractical. There are different families of the basis sets. One of the commonly used is the Dunning family consisting of the basis sets cc-pVXZ where X stands for D,T, Q and so on (double, triple, quadruple-zeta) and ’cc-p’ stands for ’correlation-consistent polarized’. Larger values of X include larger shells of polarized functions. To further improve the cc-pVXZ basis sets, one can augment them (aug-cc-pVXZ) by introducing diffuse functions. The diffuse functions play an important role when dealing with excited states since the wave functions of excited states are diffuse. For diffuse functions, the exponent ζ is small which allows the electrons to be far away from the core [20]. Figures 5.2 and 5.3 show the first six excitation energies of ethene and benzene, respectively, calculated with CAMB3LYP using different basis sets. The basis sets cc-pVQZ and aug-cc-pVQZ are omitted for benzene because of convergence failure. It can be seen that the excitation energies change significantly for the cc-pVXZ basis sets as X goes from D to Q. Due to the diffuse nature of excitation energies, the augmented basis sets give a significantly better result than the non-augmented basis sets. Since the difference between aug-cc-pVTZ and aug-cc-pVQZ is small, it is not necessary to go beyond aug-cc-aug-cc-pVTZ.

5.3 Choice of method

Testing the algorithm with different methods, HF and DFT with different func-tionals will be used. The DFT funcfunc-tionals that will be considered are LDA, BLYP, CAMB3LYP, PW91 and PBE0 (For a further description of the functionals see Section 2.7). In Table 5.4 the six first excitation energies of benzene are listed for HF and the various DFT functionals. The calculations are performed with the

basis set aug-cc-pVTZ and 10−5 a.u. as threshold value. The results are

com-pared with experimental data [21]. Comparing the performance of the different functionals, it can be seen from the table that HF roughly overestimates the ex-citation energies and gives the biggest error. The DFT functional which gave the lowest error is PBE0, followed by LDA and B3LYP. LDA, being one of the best performing functionals in these calculations, is merely by chance since LDA is the simplest functional and is usually inferior to the other functionals. Figures 5.2, 5.3 and Table 5.3 shows that the choice of method and functional greatly affects the results. Figure 5.4 shows the first excitation energies of ethene calculated with HF, B3LYP, CAMB3LYP and PBE0 using different basis sets. The experimental value of the excitation energy is taken from [22]. It is clear that to reach the reference value, both an accurate method and a good basis set are needed. To illustrate this,

(46)

5.3 Choice of method 37

cc-pVDZ

cc-pVTZ

cc-pVQZ

aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVQZ

0.26

0.28

0.30

0.32

0.34

0.36 Excitation energies (au)

1

_B

3u 1

_B

1g 1

_B

1u 1

_B

2g 1

_A

g 1

_B

3u

Figure 5.2. The first six energy levels of ethene calculated for different basis sets with CAMB3LYP.

consider HF combined with a big basis set, aug-cc-pVDZ. It does not come close to the reference value, because despite a good basis set, the unreliable method results in an error of approximately 0.01 eV. The same refers to PBE0 using cc-pVDZ, where the functional performs well but, because of the small basis set, the result is very far from the desired value. When comparing results from simulations with data acquired during experiments, it is necessary to take into account that they are not performed under the same conditions. For example, the simulations are performed in the gas phase at 0 K and do not include zero-point vibrations.

(47)

cc-pVDZ

cc-pVTZ

aug-cc-pVDZ

aug-cc-pVTZ

0.20

0.22

0.24

0.26

0.28

0.30

0.32 Excitation energies (au)

1

_B

2u 1

_B

1u 1

_E

1g 1

_A

2u 1

_E

1u 1

_E

2u

Figure 5.3. The first six energy levels of benzene calculated for different basis sets with CAMB3LYP.

Table 5.3. The excitation energies of benzene calculated using different DFT functionals and HF. The excitation energies are given in eV.

HF LDA B3LYP CAMB3LYP PW91 PBE0 Expt.

11B2u 6.582 5.268 5.406 6.454 5.255 5.494 4.90 11E1g 6.582 6.021 6.005 6.454 5.805 6.152 6.20 11B1u 6.998 6.060 6.005 6.946 5.805 6.308 6.33 11A2u 7.136 6.060 6.049 7.048 5.972 6.308 6.93 11_E 2u 7.345 6.849 6.930 7.048 6.816 7.037 6.94 11_A 1u 7.398 6.849 6.930 7.170 6.816 7.037 6.95 Mean diff. 0.63 0.3 0.3 0.47 0.4 0.25

(48)

5.3 Choice of method 39

cc-pVDZ

cc-pVTZ

cc-pVQZ

aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVQZ

0.270

0.275

0.280

0.285

0.290

0.295

0.300 Excitation energies (a.u.)

HF

B3LYP

CAMB3LYP

PBE0

Figure 5.4. The first excitation energy of ethene calculated with different basis sets and methods.

(49)

Figure 5.5. The structure of bacteriochlorin

Table 5.4. The excitation energies of bacteriochlorin, given in eV. The experimental data is taken from reference [14]

6-31G* aug-cc-pVDZ Expt. 1.9226 1.83 1.6 2.5249 2.47 2.3 3.9159 3.78 4.2588 4.14 4.3638 4.16 4.5078 4.28 4.5122 4.35 4.5431 4.45 4.7228 4.46 5.0589 4.48

5.4 Bacteriochlorin

In this section, we consider calculating excitation energies of bacteriochlorin. Bac-teriochlorin is used in photodynamic therapy (PDT) to treat cancer [23]. In PDT a photosensitizer is administered to the patient. A photosensitizer is a molecule that can alter another molecule through a photochemical process. Inside the body the photosensitizer accumulates in tumors and can then be illuminated with light. The light excites the photosensitizer molecules leading to a transfer of energy or electrons to oxygen. This leads to a creation of singlet oxygen or hydroxyl radicals that can reduce tumours. The first ten excitation energies of bacteriochlorin have been calculated and are presented in Table 5.4. The experimental data is taken from [14].

The basis set 6−31G∗, being a small basis set, gives excitation energies that are slightly further away from the reference values than in case of the bigger basis set,

(50)

5.5 Conclusions 41 aug-cc-pVDZ. Since the difference between excitation energies given by 6 − 31G and aug-cc-pVDZ is not that big, in this case, using 6 − 31G would be preferable since it requires less computational time.

5.5 Conclusions

The numerical results allows us to conclude, that the implementations using paired and symmetric vectors give identical results. The performance of the two imple-mentations is almost the same, with a small difference explained by a difference in convergence criteria. In the implementation, where symmetric trial vectors are used, the residual norm is directly compared to the threshold value, while for paired vectors the residual is first divided by the length of the corresponding response vector.

The basis set study (Figures 5.2 and 5.3) shows that increasing the value of X in the cc-pVXZ basis sets greatly improves the results, but the biggest improvement comes from augmenting the basis sets with diffuse functions. Considering the cc-pVXZ basis sets, it can be seen that it is not necessary to go beyond aug-cc-pVDZ.

From the method study it is possible to conclude that having a big basis set is not enough, an accurate method is also needed.

The calculations performed on bacteriochlorin, illustrate that the implemented algorithm performs well on bigger systems.

(51)

(52)

Bibliography

[1] A. D. Bryant and N. Frigaard, “Prokaryotic photosynthesis and phototrophy illuminated,” Trends Microbiol., vol. 14, no. 11, pp. 488–496, 2006.

[2] X. L. Ruan and M. Kaviany, “Advances in laser cooling of solids,” Heat Transfer, vol. 129, pp. 3–10, 2006.

[3] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry. Dover Publica-tions, 1996.

[4] P. A. Tipler and R. A. Llewellyn, Modern Physics. W. H. Freeman and

company, 2007.

[5] T. Jecko, “On the mathematical treatment of the Born-Oppenheimer approx-imation,” J. Math. Phys., vol. 55, no. 053504, 2014.

[6] D. J. Griffiths, Introduction to Quantum Mechanics. Pearson Education In-ternational, 2005.

[7] T. Sommerfeld, “Lorentz trial function for the hydrogen atom: A simple, elegant exercise,” J. Chem. Educ., vol. 88, pp. 1521–1524, 2011.

[8] D. Maurice and M. Head-Gordon, “Analytical second derivatives for excited electronic states using the single excitation configuration interaction meth-ods,” Mol. Phys., vol. 96, pp. 1533–1541, 1999.

[9] J. D. Purvis and R. J. Bartlett, “A full coupled cluster singles and dou-bles model: The inclusion of disconnected triples,” J. Chem. Phys., vol. 76, pp. 1910–1918, 1982.

[10] P. Hohenberg and W. Kohn, “Inhomogeneous gas,” Phys. Rev., vol. 136, no. 3B, 1964.

[11] K. Burke and L. O. Wagner, “DFT in a nutshell,” Int. J. Quant.Chem., vol. 113, pp. 96–101, 2012.

[12] P. J. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, “Atoms, molecules, solids, and surfaces: Applica-tions of the generalized gradient approximation for exchange and correlation,” Phys. Rev., vol. 46, pp. 6671–6687, 1992.

Algorithm for solving the eigenvalue reponse equation to obtain excitation energies