The Kohn-Sham potential in the
strong-interaction limit: comparison of
exact and approximate approaches
by
Tim Daas
A thesis submitted in partial fulfillment for the degree of Bachelor of Science
in the
Faculteit der Exacte Wetenschappen Theoretical Chemistry
Bachelor Thesis Scheikunde
The Kohn-Sham potential in the strong-interaction limit:
comparison of exact and approximate approaches
Door
Tim Daas
29 Juni 2017
Studentnummer
10796746
Onderzoeksinstituut
Vrije Universiteit
Onderzoeksgroep
Theoretische Chemie
Verantwoordelijk docent
Prof. Dr. Gori Giorgi
Begeleider
Abstract
Faculteit der Exacte Wetenschappen Theoretical Chemistry
Bachelor of Science by Tim Daas
The strong interacting limit, called the SCE limit, is a well-studied asymptotic case in DFT that can help to increase the number of compounds and reactions that, when coupled to the Kohn-Sham system, can be studied by computational chemistry. A good approximation for the Exc and the energy density in the strong interacting limit
is a Gradient Expansion, called the PC Model, which costs less computational time than the exact SCE functionals. However, the PC Model has not been studied as an approximation for the exchange correlation potentials, vxc, which is essential for
computational simulations, and this is the object of the thesis. From the work done here, it can be concluded that the PC Model approximation works decently in the strong interacting limit for vxc, but has its limitations. In the Helium isoelectronic series, the
PC model worked well at intermediate radii, but it lost its accuracy at low and high radii. The latter one is caused by the wrong asymptotic behaviour of the PC Model. For the neutral atoms the same result is obtained, but for Hooke’s series the PC Model approximation does not hold, due to the fast decay of the density. In order to partially solve the asymptotic behaviour of the PC model and have approximation for the physical vxc, the SPL approximation was done, where the Kohn-Sham system was coupled to the
SCE-limit. Important values and functionals in the SCE limit needed in SPL were obtained from the PC model. For H– and He, SPL does not approximate the exact potentials well, but for Be and, especially Ne, SPL worked well as an approximation. The errors of the SPL approximation come from the correlation potentials, which are more important for smaller atoms. Overall, the PC model works well in the SCE limit as an approximation for the exact vxc, but it has its limitations. The same is true for
I want to thank Msc. Sara Giarrusso for being my daily supervisor, helping me get started, when I got stuck and proofreading my thesis. I also wanted to thank prof. dr. Paola Gori Giorgi for giving me the chance to do a project in her group, for proofreading and for assisting me in my endeavours too. Thirdly, I wanted to thank Mike Daas for reading my paper as an outside source and checking my English. Lastly I wanted to thank the rest of the research group of theoretical chemistry, especially dr. Mehdi Farzanehpour, dr. Klaas Giesbretz, Msc. Juri Grossi, dr. Ivan Infante, Bsc. Derk Kooi, Prof. dr. Michael Seidl, Wouter Schreuder and drs. Stefan Vuckovic.
In de scheikunde zijn computersimulaties nodig omdat sommige reacties te gevaarlijk zijn, te duur zijn, te snel gaan of een combinatie van deze drie redenen. Computersimu-laties zijn vooral gemaakt op basis van de kwantummechanica, dit om een fundamenteler beeld te krijgen van hoe reacties verlopen. Alleen kunnen deze simulaties niet gemaakt worden voor systemen groter dan He of H2, zonder dat benaderingen nodig zijn. Bij deze systemen zijn namelijk meer dan drie elektronen die elkaars lading voelen en dus willen afstoten. Om deze bewegingen te beschrijven, is veel tijd en computerrekenkracht nodig. Een van de beste manieren, binnen de scheikunde, om dit soort systemen te beschrijven, is Density Functional Theory. Dit betekent dat alles beschreven kan worden als functies van de elektronendichtheid. Wanneer je deze dichtheid weet, kan je alles simuleren, maar ook hier zorgen meer-elektronen-systemen voor complicaties.
De meest gebruikte manier, om toch deze systemen snel te kunnen berekenen, is om te postuleren dat elektronen geen lading hebben, waardoor je de kinetische energien exact kan berekenen, wat het Kohn-Sham-Systeem wordt genoemd. In vergelijking met de echte wereld is er nu geen interactie, maar dit wordt opgelost door een moeilijk exact berekenbare term, de exchange correlation energy, Exc. Later worden benaderingen
van de Exc toegevoegd die het afstotende effect beschrijven, bijvoorbeeld verbieden dat
er elektronen mogen zijn in een bepaalde ruimte om een ander elektron. Dit is niet hoe de elektronen zich fysisch gedragen, maar zorgt er wel voor dat simulatie snel en goed gedaan kunnen worden. Alleen systemen met lage elektronendichtheden, zoals H-, kunnen niet beschreven worden omdat hier de interactie belangrijker is.
Er is veel onderzoek gedaan naar het verbeteren van het Kohn-Sham-Systeem, maar o.a. de VU is bezig met een ander systeem. In plaats van te veronderstellen dat de elektronen geen interactie hebben, hebben ze nu juist een oneindig sterke interactie, SCE-limiet, en geen kinetische energie. Je zou denken dat de elektronen nu oneindig ver uit elkaar zouden gaan, maar beperkingen zijn toegevoegd zodat de elektronen dezelfde dichtheid moeten geven als in de werkelijkheid. Je kan je dit voorstellen als elektronen die stil staan en op een bepaalde plaats worden gehouden door ijzeren staven, die de oneindige
interactie voorstellen. Dit systeem is net als het geen-interactie-system, makkelijker te berekenen. alleen zijn benaderingen ook hier nodig. Om de Exc makkelijker te
berekenen, is het gesplitst in lokale potentialen die bij elkaar opgeteld de Exc geven.
Een bepaalde benadering, genaamd het PC Model, werkt heel goed voor energiedichthe-den in het SCE-limiet Deze benadering zegt dat je de elektronen kan plaatsen in een denkbeeldige vloeistof met een positieve lading gelijk aan het atoomnummer. Omdat de elektronen niet bewegen, kan je de ruimte indelen in ronde neutrale cellen met elk een elektron erin. Deze beschrijving kost weinig computerkracht, omdat het alleen afhangt van de elektronen dichtheid en de afgeleide daarvan. Het zou handig zijn als deze be-nadering ook werkt voor de lokale potentialen.
Ik heb de kwaliteit van dit model onderzocht aan de hand van drie series. Voor de eerste serie, met atomen/ionen met 2 elektronen, werkte het goed alleen minder goed dichtbij de atoomkern er op grote afstanden ervan. Voor de neutrale atomen werkte het beter en voor een niet bestaand systeem, waar de elektronen worden beschreven als veren, werkte het amper. Ik heb daarna het SCE-limiet en Kohn-Sham-Systeem gekoppeld aan elkaar om de echte situatie te beschrijven en deze koppeling werkte ook goed voor de potentialen. Door mijn resultaat zouden simulaties sneller gedaan worden als men het PC Model koppelt met het geen-interactie-systeem, al zijn er beperkingen en deze moeten worden verbeterd.
Abstract ii Acknowledgements iii Samenvatting iv 1 Introduction 1 2 Theory 4 2.1 DFT . . . 4 2.1.1 Introduction . . . 4 2.1.2 Kohn Sham . . . 6 2.2 SCE . . . 9 2.2.1 SCE Exact . . . 10 2.2.2 PC Model . . . 12 2.2.3 SPL . . . 15 2.3 Studied Systems . . . 17
2.3.1 Helium Isoelectronic Series . . . 17
2.3.2 Neutral Atoms . . . 17
2.3.3 Hooke’s Atom. . . 18
3 Computational Details 19 4 Results and Discussion 21 4.1 PC Model . . . 21
4.1.1 Helium Isoelectronic Series . . . 21
4.1.2 Neutral Atoms . . . 23
4.1.3 Hooke’s Atom. . . 24
4.2 SPL . . . 25
5 Conclusions 28
A Results for the Helium Isoelectronic series 30
B Results for the Neutral atoms 38 vi
C Results for the Hooke’s series 43
D Results for the SPL potentials 50
Introduction
Many important reactions done in chemistry are either too fast to study, too danger-ous, too costly or a combination of the three, which means that computer simulations are needed in those cases. Furthermore, computer calculations can also be used to understand the mechanism of certain reactions. Due to important breakthroughs by the likes of Planck, Einstein, Schr¨odinger and Bohr, quantum mechanics was born and thus a fundamental theory of particles, including the molecules used in chemistry. Due to high computational cost needed to calculate chemical systems, approximations are needed [1]. Luckily, chemists and physicists came up with important approximations that could relieve stress on the computer and still get reasonable theoretical data that are consistent with the experimental ones. Nonetheless, still many chemically interest-ing systems cannot be explained or computed accurately and many errors arise in these kinds of calculations.
An important way of calculating chemical reactions is DFT (Density Functional Theory), which uses electron density dependent functionals instead of wave functions. DFT is in principle exact, according to the inventors Hohenburg and Kohn [1, 2]. However, there are still unknown functionals within DFT due to the limitations that the many body calculation presents us with. Those unknown functionals are studied extensively and all the unknown terms are collected into one, called the exchange correlation functional. If this functional were known, all DFT calculations would be exact. Many studies within theoretical chemistry are either trying to find new approximations or improving existing functionals to have higher accuracy. Most of these functionals work well for specific systems, for example predicting spectroscopic spectra or catalytic reactions, but fail in other chemical systems.
In DFT, most of the calculations are based on the Kohn-Sham system, which is exactly solvable [1, 3]. In this system electron-electron interactions are mapped into single-particle potentials and the electrons are postulated to only have kinetic energy. Although this potential can be calculated exactly in principle, in practical calculations must be approximated. To increase the speed of the calculation on the Kohn Sham system, easy calculable functionals or other approximations are made for the exchange correlations energies to fix the errors in this system.
Even though the current approximations predicted many energies and reactions with high accuracy, they have some flaws, such as ignoring the existence of the hydride-ion and homogeneous bond breaking of small molecules [4,5]. Some of these problems can be solved by introducing the other limit, where there is no kinetic energy and only strong electron-electron correlation, also known as the strictly correlated electron limit (SCE) [6]. The SCE framework coupled to the Kohn-Sham formalism, works better on systems with low electron densities, which are systems such as the hydride-ion or homogeneous bond breaking, thus increasing the amount of usable systems.
Much research has been done into the exact SCE functionals and the approximations thereof, most notably the PC Model, to check how good the approximations work for the energies and energy densities of some small or symmetrical systems [7]. The PC (point charge plus continuum) Model is rather interesting, because it is only dependent on the density and its gradient and can easily be calculated [6]. Although the PC Model is good in approximating the exact SCE energies, it has never been used to study the exchange correlation potentials of the exact SCE, which are more useful for Kohn-Sham DFT. If the PC model is a good approximation for exact potentials, then it can be used for more complex systems and, in turn, increase the speed of other calculations. Furthermore, if the PC model works well, it can be implemented via the SPL method, where the SCE is coupled to the Kohn Sham system, to approximate the physical system.
The main research question of this thesis is: Does the PC Model work as an approxi-mation for the exchange correlation potentials in the strong interaction limit and how well does it work? To test this, the PC model will be studied against the exact SCE potentials for three different systems. The first one is the Helium isoelectronic series where every atom has two electrons only, the second series is the neutral atoms consist-ing of Hydrogen, Helium, Lithium, Beryllium, Boron, Carbon and Neon and the third the Hooke’s atom. This last system does not exist physically, but gives a system that is exactly solvable and thus interesting to test the PC Model against. If the PC Model does not work, other approximations, such as SPL, will be tested to see if they work better than the PC Model. This approximation will not be tested versus the SCE ex-change correlation potentials, but versus the physical ones. The thesis starts with an
introduction to DFT, SCE, SPL and the used systems, which is needed to understand the theoretical background of the study, followed by a discussion of the computational procedure. The fourth chapter will show and discuss the results and the thesis will end with a definitive conclusion and outlook on the strong interaction limit, especially on the PC Model.
Theory
2.1
DFT
2.1.1 Introduction
After the introduction of Quantum Mechanics and the Schr¨odinger equation, where the latter one can be found in equation 2.1,
Hψ = Eψ, (2.1) most quantum systems are in principle exactly solvable [1]. However, just like in gravi-tational calculations on the solar systems, many body problems are still unsolvable due to high computational time needed to calculate exact values. The many body problem in quantum chemistry is found in the ˆEee, shown in the atomic Schrodinger equation
below, due to the n!/2 different interactions between n electrons,
ˆ Hψ = ( ˆT + ˆVN e+ ˆVee)ψ = −~2 2m X j ∇2j−X j,l Zle2 |~rj− Rl|+ 1 2 X j6=j0 e2 |~rj − ~rj0| ψ, (2.2)
where ˆVee is the electron-electron interaction due to the coulomb interaction, ˆVNe the
nucleus-electron attraction between a negative electron and a positive nucleus and ˆT the kinetic energy of the electrons. The kinetic energy of the nucleus is already removed by the Born Oppenheimer approximation. If systems are used with more than one atom, an extra term is added called ˆVN N, which is the nucleus-nucleus interaction. As a further
complication, electrons are fermions with spin 1/2, which means that the wavefunction Ψ, form equation 2.2, must be anti-symmetric. All of this means that any attempt to
calculate energy of either an atom or molecule exactly, is too computational demanding to be done routinely.
Another way to calculate the energy and other important properties of matter is Density Functional Theory (DFT) where the electron density, ρ(~r), takes the role of the wave function, which was proven to be exact by Hohenburg and Kohn in 1964 [2]. In other words, they needed to proof their first theorem: ”the external potential, Vext(~r), is
(within a constant) a unique functional of ρ(~r); since, in turn Vext(~r) fixes ˆH we see
that the full many particle ground state is a unique functional of ρ(~r).” [1,2] The big advantage of DFT is that it is easier to use but still in principle exact, whereas Hartree-Fock, a less computationally demanding wave function method, is not [1]. The electron density is defined as
ρ(~r) = N X
s1,...,sN
Z
d3r2. . . d3rN|Ψ(~r, ~r2, . . . , ~rn; s1, . . . , sn)|2, (2.3)
which means that it is the multiple integral over all but one spatial coordinates and the sum of the spin coordinates, so that ρ(~r) describes the probability to find one electron in the volume element with any spin coordinate. However, constraints must be added to the density to give it a physical interpretation, which are
Z
ρ(~r)d~r1= N, (2.4)
where equation2.3shows that the integration over the whole volume must give the total number of electrons.
The expectation value of the local external potential can, by using the first theorem from Kohn and Hohenburg, [1,2] be rewritten as
hΨ| ˆVext|Ψi =
Z
d~rρ(~r)v(~r). (2.5) By now adding the electronic density into the Hamiltonian of the wave function, the DFT Energy Functional can be derived as
E[ρ(~r)] = F [ρ] + Vext[ρ] = F [ρ] +
Z
d~rρ(~r)v(~r), (2.6) where Vext comes from the nucleus-electron interaction, which depends on the number
of electrons (N), the atomic number (ZA) and the atomic position (RA) and where
F , sometimes called the Hohenberg-Kohn or universal functional [1], is independent of those variables and contains the kinetic energy, T , and the electron-electron energy
contribution, Vee, visualized in equation2.7below,
F [ρ] = T [ρ] + Vee[ρ] = min
Ψ→ρhΨ|Vee+ T |Ψi . (2.7)
The universal functional can also be written as the expectation value of Vee+T , visualized
in the last part of equation 2.7, where Ψ is the wave function that minimizes T + Vee
and yield the density ρ(~r). This density can be found by using the second theorem of Hohenburg and Kohn, the Variational Principle [2], and the Euler Lagrange equation on equation 2.6, so that,
δE[ρ] δρ~(r) =
δF [ρ]
δρ(~r) + v(~r) = µ, (2.8) is obtained, where µ is the Lagrange multiplier that ensures that equation 2.4 still holds [5].
Because the second term on the right-hand side of equation2.5can be calculated, know-ing the F [ρ] exactly, will yield exact solution for a certain system. [1] However, F [ρ] is not completely known, which hinders the exactness of DFT. F [ρ] can be split into three parts,
F [ρ] = T [ρ] + Uh[ρ] + W1[ρ], (2.9)
where T [ρ] is the kinetic energy of the electrons, Uh[ρ] the Hartree energy term which
covers the classic coulomb repulsion between electrons, and the indirect energy, W1[ρ],
which describes the interactions between electrons. Although knowing the indirect en-ergy would yield exact solution, this is not the case for most systems. Even in the small amount of systems where it can be calculated exactly, it still is computationally demanding. For even bigger systems, reasonable approximations are needed to get the wanted properties.
2.1.2 Kohn Sham
An interaction strength scaled functional, which uses the adiabatic connection framework proposed by [8–10], can be introduced as
Fα[ρ] = Tα[ρ] + αVeeα[ρ] = min
Ψ→ρhΨ|α ˆVee+ ˆT |Ψi , (2.10)
to obtain more systems to study. These new systems can, in certain cases, be used as an approximation for the physical systems. From equation2.10 it can be seen, that the electron-electron interaction term is scaled by α so that at α = 1 it would give back the universal functional of equation2.7 [5]. When α = 0 in equation2.10, the most widely used approximation of DFT is obtained, called Kohn Sham DFT. This is a system with
the same density as the physical one but with non-interacting electrons, which means that only the kinetic energy of the electrons must be minimized [1,3]. If α = 0 is inserted into equation2.10, the following formula can be obtained,
F0[ρ] = min
Ψ→ρhΨ| ˆT |Ψi = Ts[ρ], (2.11)
which shows that only the kinetic energy must be calculated exactly. However, the system does not contain any contributions from the electron interactions, what the physical universal functional, F [ρ], obviously does have.
If Ts[ρ] is subtracted from the physical universal functional, the leftover functional is
then called the Hartree Exchange Correlation functional, which consists of a Hartree term, U0
h[ρ], and an exchange correlation functional, Exc[ρ]. Due to the Hartree term
being universal for every system, the Hartree functional is also equal to Uh[ρ] [11]. All
of this will give an equation similar to equation2.8,
F [ρ] = Ts[ρ] + Uh[ρ] + Exc[ρ], (2.12)
where a part of the kinetic energy of T , due to the lack of interaction, has been placed into indirect energy. This means that Exc[ρ] is an collection of all unknown terms
of Kohn-Sham DFT, which are the exchange and correlation energies. The two main properties that the exchange correlation functionals must inhibit are the two errors that the Hartree term makes. These errors are the self-interaction of an electron due to the smeared out negative background and the quantum mechanical behaviour of electrons i.e. Pauli exclusion and Coulomb repulsion. Although the exchange correlation functional is thus unknown, the adiabatic connection gives us a way to calculate Exc[ρ] and in turn
vxc. Because all Ψα[ρ] are made to give back the physical density for all values of α, we
can now define the coupling constant integral, ˆWα, as
Wα[ρ] = hΨα[ρ]| ˆVee|Ψα[ρ]i − Uh[ρ], (2.13)
which integrated over α from 0 to 1, gives us a way of calculating Exc[ρ] with
Exc[ρ] =
Z 1
0
dαWα[ρ]. (2.14)
If equation2.12 is substituted into the Euler Lagrange equation of2.8, it results into, δE[ρ]
δρ(~r) = δTs[ρ]
where vKS is the unique effective Kohn Sham potential, which is a Lagrange multiplier
in equation 2.15 [5]. It must consist of the local external potential, vext, and the local
Hartree and exchange correlation potentials, vH and vxc respectively. This is in formula
form,
vKS(~r) = vext(~r) + vH[ρ](~r) + vxc[ρ](~r). (2.16)
By definition, the vxc can be obtained from ˆExc[ρ] via equation2.15,
δExc[ρ]
δρ(~r) = vxc[ρ](~r). (2.17) Now approximations can be made for the exchange correlation potential so that simula-tions and calculasimula-tions can be done on chemically and physically relevant systems using the Kohn Sham Hamiltonian,
ˆ
H = ˆT + ˆVKS, (2.18)
where ˆVKS is the Kohn-Sham Potential operator, which represents the one-body
poten-tial, which ensures that the ground state density of the KS Hamiltonian is the same as the physical density [3].
In order to get an expression for the exchange correlation Energy at any α, Excα[ρ(~r)], let’s go back to the Kohn Sham system [1, 11]. There we mentioned that Ts ignored
some kinetic energy of the real T due to having no electron interaction and this was placed into the Exchange Correlation functional. A part of the exchange correlation functional can be visualized as the ignored part of the kinetic energy at a certain α, Tα,
due to only having the exact non-interacting kinetic energy of Ts. If we now neglect
the kinetic energy, then the exchange correlation energy at a certain α consists only of the ignored part of the Eeea[ρ(~r)] at a certain α due to the aforementioned errors of the Hartree term, Uhα[ρ(~r)]. It was already mentioned that the Hartree term is universal so that the Hartree term at any α is equal to the one of the physical system. If we multiply the Uh[ρ(~r)] and Eeea[ρ(~r)] by α, then the exchange correlation functional, at α = 0,
consists only of the ignored kinetic part, due to the absence of any form of electron interaction. Furthermore, at α > 0, the ignored kinetic energy term is still present. However, the exchange and correlation terms must also be accounted for. This is done by a scaling parameter, α, to get physically logical systems. This can be summarized by the following formula,
Excα[ρ] = (Tα[ρ] − Ts) + α(Eαee[ρ] − Uh[ρ]). (2.19)
Approximations are needed to be able to calculate accurate exchange correlation poten-tials so that DFT can be used to simulate atoms, molecules or even reactions [1,5]. The
simplest approximation is Local Density Approximation (LDA), which has the following formula,
Exc[ρ] =
Z
d~rρ(~r)LDAxc (ρ(~r)), (2.20) where LDAxc (ρ(~r)) is the local exchange correlation energy per particle, which depends on the local density, ρ(~r). LDA works fine for condensed matter and solid-state physics but for chemistry Generalized Gradient Approximations (GGA) were invented. They use the density and its gradient but also allow the gradient to expand in classically forbidden regions e.g. regions only accessible by quantum tunnelling. This increases the accuracy and opens the world of computational chemistry to a wider range of possible systems to study. Further improvement on functionals, such as hybrids, is out of the scope of this thesis, but the reader can gather more information about this in Kohn et al. [1].
2.2
SCE
Although the Kohn Sham limit can give accurate results for many systems, it has some problems due to the approximations that have to be made in order to obtain exchange correlation functionals and potentials [4]. One of those problems are the errors that arise in the calculation of the fundamental conduction gap, which is important for condensed matter physicists working on semiconductors and chemists working on spectroscopy [5,
12]. Homogeneous bond stretching in small molecules, also important in chemistry, does not work in the Kohn Sham formalism either. The same goes for the existence of the hydride-ion and fluoride.
An interesting limit to study is the case where α → ∞, because it has been proven to produce exact solutions [13]. The Hamiltonian of this so called SCE (Strictly Correlated Electron) limit can now be written for every α as
ˆ
H = ˆTs+ α ˆVee+ ˆVextα , (2.21)
where ˆVextα can be understood as the potential that ensures that the density at a certain α equals the physical density. Just like in equation 2.5, ˆVextα can be written as the sum of multiple vextα terms [14]. In order to get the physical density back, the external potential must compensate the infinitely strong electron-electron interaction. This goes proportionally to α, because higher interaction means more compensation needed, which gives
lim
α→∞
vextα ([ρ], ~r)
α = vSCE([ρ], ~r), (2.22) where vSCE can be defined as the SCE external potential.
If we now go back to equation 2.10 and assume that hΨα| ˆVee|Ψαi grows faster than
hΨα| ˆT |Ψαi [13], equation 2.10can be simplified to
Fα[ρ] = α min
Ψ→ρ hΨα| ˆVee|Ψαi (2.23)
and using 2.11 gives,
Fα[ρ]
α = Wα[ρ] + Uh[ρ] = V
SCE
ee (2.24)
when α → ∞ [4,5,13]. This means that we will only consider Wα[ρ] or VeeSCE, because
we will also consider Fα[ρ] at the same time. VeeSCE can be understood as the SCE version of the Ts in the Kohn Sham system due to both being exactly solvable. It can
also be understood as the sum of all the vextα terms in the α → ∞ limit, thus the sum of all vSCE.
This limit can be interpreted physically as the situation where the electrons, placed in a certain space, will minimize their interaction energy by finding the configurations with the lowest energy where their kinetic energy is zero [4]. So, in this limit there is infinite electronic correlation between the electrons but zero kinetic energy. This means that the movement of one electron affects the position of all the other electrons drastically in a certain space.
2.2.1 SCE Exact
Most of the mathematical details will be skipped, but the reader can find the details in references [5, 13]. Following from equation 2.24 and 2.25, we can see that VeeSCE is the minimum expectation value of Vee with the constraint that the density is smooth [4].
However, the square of the limiting wave function, which is the probability of finding the electrons at positions ~r1, . . . , ~rn, is not smooth and has the form of Dirac delta functions,
because the electrons can only be found at certain specific positions due to the strong correlation between them. This can be described by the so-called co-motion functions, because the placement of one electron affects the position of the rest of the electrons. These N co-motion functions can be mathematically understood as
~
ri= fi[ρ](~r) ∀ i ∈ {2, 3, . . . , N }, (2.25)
where fi[ρ](~r) are the co-motion functions and f1(~r) = ~r [5,13]. The co-motion
func-tions must also obey cyclic group properties in order to give physically relevant results. Because the wave function is associated with the density by Ψ → ρ, the co-motion
functions must be as well. This relation is given in equation 2.25 below,
ρ(~r)d~r = ρ(fi(~r))dfi(~r) ∀ i ∈ {2, 3, . . . , N }, (2.26)
which means that the probability of finding the reference electron at ~r in the volume element d~r is the same as finding a second electron at position fi[ρ](~r) in the volume
dfi[ρ](~r).
To calculate the important properties of the SCE limit in DFT exactly, a more usable formula than 2.26 must be found [5, 13]. This can be done for spherically symmetric densities, where equation2.26 can be rewritten as
4πx2ρ(x)dx = 4πfi2(x)ρ(fi(x)) dfi(x) dx dx. (2.27) This can be simplified by introducing a new variable, called the cumulant, as
Ne(r) =
Z r
0
ds4πs2ρ(s), (2.28) which gives the following formula for the co-motion function for odd numbers of particles,
f2k+1(r) = Ne−1[2k + Ne(r)] for r ≤ aN −2k Ne−1[2N − 2k − Ne(r)] for r > aN −2k, (2.29)
where k is an integer, which runs from 1 to N −12 for odd N [13]. For even numbers of electrons, k runs from 1 to N −22 and in this case, the comotion functions can be obtained by f2k(r) = Ne−1[2k − Ne(r)] for r ≤ a2k Ne−1[Ne(r) − 2k] for r > a2k (2.30) and an extra third formula,
fN(r) = Ne−1[N − Ne(r)]. (2.31)
Equation 2.31can now be simplified for a 2 electron system, which yields,
f (r) = Ne−1[2 − Ne(r)]. (2.32)
According to the article written by Seidl et al. [13] the SCE potential, for N = 2, is now, vSCE[ρ](~r) =
Z r
dx
where the mathematical derivation can be found in reference [13]. This potential, vSCE,
compensates the net force acting on an electron at position, ~r, from the N − 1 other electrons at their positions, fi. Using a shortcut found in Malet et al. [15], it is possible
to couple this to the VeeSCE[ρ], which gives δVeeSCE[ρ]
δρ(~r) = −vSCE[ρ](~r). (2.34) This is equivalent to equation 2.17, because vSCE− vh is the SCE version of vxc, and
VeeSCE[ρ] that of the Exc. The Hartree potential, vh, can be calculated by the following
shortcut for spherical symmetric systems as vh(r) = Ne(r) r + 4π Z ∞ r dr0r0ρ(~r0). (2.35)
Now, by using the cumulant and the co-motion functions, it is possible to obtain exact values for the unknown exchange correlation potential in the SCE limit for N = 2. This has never been studied before, thus giving an interesting starting point for this study. In the SCE formalism, other systems with more electrons can be calculated exactly too. However, they are out of the scope of this thesis theoretically. If interested, the article written by Gori Giorgi et al. [4] would provide an excellent supplementary account to the exact mathematical derivations and other exactly solvable systems within the SCE formalism. Even though, the SCE limit is exactly solvable by using the co-motion functions, it is not known how to do it in general for non spherically symmetric systems. There are numerical algorithms, but they are computationally too expensive.
2.2.2 PC Model
The Point Charge Plus Continuum Model (PC Model) has been introduced by Micheal Seidl in 1999 [6], before the exact SCE was known, in order to calculate properties in the strong interaction limit. Instead of using the density at every point in space, it uses either local or semi-local functionals which are easier to compute. For a more complete derivation of the PC model readers are redirected to reference [6]. In this model the W∞[ρ] is rewritten as the electrostatic energy, ˆEes, of N electrons in a certain state,
smeared smoothly over a background of positive charge with a density equal to the physical density [5]. This is visualized in the equation,
ˆ
Ees= ˆEee+ ˆEeb+ ˆEbb= hΨ[ρ] |Vee| Ψ[ρ]i − Uh[ρ] = W∞[ρ], (2.36)
where ˆEee = hΨ[ρ] |Vee| Ψ[ρ]i is the electronic repulsion energy per definition, ˆEeb =
interaction. Due to how the extreme interaction limit works, it is expected that the elec-trons move into relative positions, which can be approximated as neutral cells containing one electron. Now the Ees can be approximated by calculating the electronic energy of
each cell, Ecell, and summing these up. Using this, the W∞[ρ] can be approximated by
the following formula,
W∞[ρ] =
Z
d~rρ(~r)Ecell([ρ, ~r], (2.37)
which can be used to calculate F∞[ρ].
In the same article [6] values for Ecell are proposed in order to have a LDA and GEA
(Gradient Expansion Approximation) functional for the strong interacting limit. The GEA functional is
EcellGEA = Aρ(~r)1/3+ B|∇ρ(~r)|
2
ρ(~r)7/3 (2.38)
with A = −109(4π3 )1/3 and B = 3503 (4π3 )1/3. The LDA only uses the first term by its definition of being only dependent on the density. Both of these are tested versus the exact SCE energy densities and they work reasonably well [7]. However, their ability to approximate the exchange correlation potential needed in Kohn-Sham DFT calculations are yet to be examined and are thus interesting to study. This can be examined by rewriting equation2.14, where instead of integrating from 0 to 1, the integration is done from 0 to α in the α → ∞ limit,
Excα[ρ(~r)] = Z α
0
dα0W∞[ρ] = αW∞[ρ]. (2.39)
If the limit of Excα[ρ(~r)]/α is taken for α → ∞, then ExcSCE is obtained as ExcSCE= lim
α→∞
Excα[ρ(~r)]
α = W∞[ρ]. (2.40) By applying equation2.17 to equation 2.40, equation 2.40results in
δW∞[ρ]
δρ(~r) = v
SCE
xc [ρ](~r), (2.41)
which can be understood as the functional derivative of the PC model, which gives a way to calculate the vxc of the PC Model. This can be found in the formula below,
δW∞[ρ] δρ(~r) = 4A 3 ρ(~r) 1/3− 2B∇2ρ(~r) ρ(~r)4/3 + 4B 3 |∇ρ(~r)|2 ρ(~r)7/3 . (2.42)
The results of equation2.42 can be compared to the exact SCE results in order to find out if the PC model is a good approximation for vSCE in DFT.
instead of going to ±∞ at larger r. For the PC model this can be checked fairly easily by inserting the following equations in the PC model:
ρ(r → ∞) = rβe−αr (2.43) , α = 2√2I (2.44) , β = Z − N + 1√ 2I − 1 (2.45) with I the ionization energy of an atom, Z the atom number and N the number of electrons. For the neutral atoms β can be rewritten as: β = α2 − 1, which means that β can be eliminated from equation 2.43 and this is also true for Hooke’s atom, where Z can be replaced by n, the principal quantum number. For the helium isoelectronic series, however, β can be rewritten as β = 2(Z−N +1)α − 1. By now inserting equation2.43
into 2.42 and taking the first and second spatial derivative of ρ(~r), the asymptotic behaviour can be examined. The graph below contains the result for the neutral neon atom, which is also true for every other atom/ion:
2 4 6 8 10 r a -60 -40 -20 vXC
Figure 2.1: The asymptotic behaviour of the PC Model exchange correlation poten-tials for Neon vs r/a.
The graph shows that PC Model has the incorrect asymptotic behaviour, because the exact vSCE(~r) graphs follow the behaviour of the coulomb potential at large values of r,
and thus goes asymptotically to zero, whereas the PC model goes to −∞. If this defect happens after the density has decayed to insignificantly low values, it will not affect the accuracy of the PC Model approximations. Therefore, this must be researched before the PC Model potentials can be used as an approximation of the exact SCE potentials.
2.2.3 SPL
The reason why the SCE limit is interesting, apart from having the ability to simulate properties that the Kohn Sham system cannot, is that it can be coupled to the Kohn-Sham systems using the adiabatic connection described earlier [4, 5, 11, 16]. In the Kohn Sham system the only information that is known, is W0[ρ] and its derivative
W00[ρ]. From these two values an extrapolation can be done in order to get the coupling integrand value at α = 1 and, by using equation2.14, the exchange correlation energy. However, the accuracy is very low due to only having one point and slope, which gives a straight line. By also implementing a second system, for example α = ∞, which can be calculated exactly, an interpolation can be made with an higher accuracy than only using α = 0. If the strong limit is implemented, thus W∞[ρ] and the next leading term
W∞0 [ρ], then an interpolation can be done. In the adiabatic connection, the coupling
constant, α, scales the electron-electron interaction, which means that the interpolation is also coupled to α. This interpolation is named SPL after its inventors Seidl, Perdew and Levy. This is visualized by the graph below:
Figure 2.2: The SPL graph obtained by using the adiabatic connection framework for the helium atom [16]. SPL was called ISI in earlier papers, which is also visible in this figure. The two dotted lines with a formula attached to them are extrapolations of one α value and its derivative. The other dotted line is the asymptotic value of W∞
and the bolded line is obtained by SPL. The integral obtained by SPL differs from the not exact Kohn Sham extrapolation.
In the article of Seidl et al. [6] a formula was proposed for SPL, which can be seen in the equation below,
WαSP L[ρ] = W∞[ρ] + W0[ρ] − W∞[ρ] p1 + 2X[ρ]α , X[ρ] = W00[ρ] W∞[ρ] − W0[ρ] . (2.46)
The derivation and the requirements needed to obtain the SPL formula can be found in reference [6]. In that article, the accuracy of the SPL functional was also tested for He, Ne8+, Be, Ne and two other electron densities. Equation 2.17 can be applied to WαSP L[ρ], just like to the W∞[ρ] of the PC Model, so that the vxc[ρ](~r) can be obtained.
Although this has never been tested, in this thesis another form of SPL formula will be obtained, because W00[ρ] is computationally more demanding than W∞0 [ρ]. This can be
done by altering equation 2.46 so that instead of using the W00[ρ], W∞0 [ρ] can be used, which gives WαSP L[ρ] = W∞[ρ] + W0[ρ] − W∞[ρ] p1 + Y [ρ]α , Y [ρ] = (W0[ρ] − W∞[ρ])2 W0 ∞[ρ]2 . (2.47)
It was already proven that the PC Model potentials did not have the right asymptotic behaviour, which could cause problems in the accuracy of the model. The SPL potentials could have the right asymptotic behaviour, because it has the functional derivative of W∞0 , see equation 2.49, which goes to ∞ and can counteract the PC model part.
However, SPL does not approximate the SCE limit, but the physical system, which is much more useful in physics and chemistry. Therefore, the second form will be tested against the exact physical potentials and its accuracy will be checked. Before WαSP L[ρ] can be used to study the potentials, equation 2.47must be integrated via equation2.14
so that it gives the exchange correlation energies, ExcSP L[ρ] as
ExcSP L[ρ] = (W0[ρ] − W∞[ρ])W∞[ρ] + 2 −1 +q1 +(W0[ρ]−W∞[ρ])2 (W0 ∞[ρ])2 (W∞0 [ρ])2 W0[ρ] − W∞[ρ] . (2.48) The potential vSP L[ρ] can be found by letting equation 2.17 work on ExcSP L[ρ], which
gives the following formula, vSP L[ρ] = δExcSP L[ρ] δρ[~r] = ∂ExcSP L[ρ] ∂W0[ρ] δW0[ρ] δρ[~r] + ∂ExcSP L[ρ] ∂W∞[ρ] δW∞[ρ] δρ[~r] + ∂ExcSP L[ρ] ∂W0 ∞[ρ] δW∞0 [ρ] δρ[~r] . (2.49) The partial differential terms can be calculated easily by applying them to equation2.48. The first functional derivative is just vx from equation2.35, the second can be found in
equation 2.41and the third will be explained later, but can be found in equation 2.51. An equation for W∞0P C[ρ] has been proposed in the article written by Seidl et al. [17],
which is W∞0 [ρ] = Z d~rCρ(~r)3/2+ D|∇ρ(~r)| 2 ρ(~r)7/6 , (2.50)
with C = 1.535 and D = −0.02558. The derivation of the exact W∞0 [ρ] is out of the scope of this thesis, but a good, extensive mathematical derivation for the exact W∞0 [ρ]
functional is δW∞0 [ρ] δρ(~r) = 3C 2 ρ(~r) 1/2− 2D∇2ρ(~r) ρ(~r)7/6 + 7D 6 |∇ρ(~r)|2 ρ(~r)13/6. (2.51)
The interpretation of W∞0P C[ρ] is the Zero Point Oscillations (ZPO) of the SCE electrons, due to the Heisenberg uncertainty principle, which says that a particle must have a non-zero velocity if the particles position is known.
2.3
Studied Systems
2.3.1 Helium Isoelectronic Series
The Helium isoelectronic series is the series of atoms (or ions) containing N = 2 electrons but with a different atomic number. In this thesis only the atoms and ions from H– to Ne8+ and Zcrit were considered. Zcrit is the minimal atomic number needed to bind two
electrons and thus has the lowest electron density possible for a two-electron system. Even though it is not a physical system, it is still an interesting system to study the PC Model for, due to the high electron correlation of Zcrit = 0.9110289. The helium
isoelectronic series is used, because a two-electrons-systems can be calculated exactly using the SCE formalism, which means that exact results can be compared to the PC Model results. The Hamiltonian of this series, in atomic units, is
ˆ H = −1 2∇ 2 r1− 1 2∇ 2 r2− Z r1 − Z r2 + 1 |r1− r2| , (2.52) where the first two terms are the kinetic energy of the two electrons, the second two terms the nucleus-electron interaction with the nucleus at the origin and the last term the electron-electron interaction. It is expected that the density will be closer to the nucleus and narrower spread at higher Z-values, due to the higher positive pull felt by the electrons.
2.3.2 Neutral Atoms
The second series is the neutral atoms consisting of H to C and Ne. This series is chosen, because neutral atoms are chemically and physically more important than the helium isoelectronic series. Furthermore, the calculations cost a lot more computational power than the helium isoelectronic series so that having a good approximation is much more important here. The Hamiltonian of this series can be found in section 2.1 as equation 2.2.
2.3.3 Hooke’s Atom
The last series that will be examined is the Hooke’s atom, which is a two-electron series where the external potential of the nucleus-electron interaction is described by Hooke’s law or the harmonic oscillator potential [18]. Now the Hamiltonian of equation2.2can be written as ˆ H = −1 2∇ 2 1− 1 2∇ 2 2+ ω2r2 1 2 + ω2r2 2 2 + 1 |r1− r2| , (2.53) where again the first two terms are the kinetic energy of the two electrons, the second two terms the nucleus-electron interaction where ˆVext= −ω
2r2
2 is used and the last term the
electron-electron interaction. This can easily be rewritten in other coordinates, which are the center mass, R = r1+r2
2 , and the difference vector, r = r2− r1. The Hamiltonian
is now ˆ H = −∇2r−1 4∇ 2 R+ ω2r2 4 + ω 2R2+1 r = ˆHr+ ˆHR, (2.54) which can be split due to having no cross terms between the two coordinates. According to Taut [18], only certain values of ω can be used to calculate exact ground state densities, which are coupled to the principal quantum number, n. In this thesis the studied n values are N = 2, 3, 5, 6, 7, 9, 10, 11, 12, 13.
Computational Details
For the helium isoelectronic series, the calculations are done via the outline proposed in the SCE exact and PC model subsections, 2.2.1 and 2.2.2 respectively. An accurate density, placed on a grid, was used and the rest was calculated exactly using the afore-mentioned formulas. The results of the neutral atoms were obtained by using a fit of accurate densities instead of real ones, due to high computational time it costs to calcu-late important potentials and properties from the exact one. Moreover, the vSCE(~r) was
taken from earlier research done on the isoelectronic series for the same reasons [13]. For the Hooke’s atom, the densities were taken from earlier research by Mirschrink et al. [7] and the following calculations were done via the recipe proposed in the SCE Exact and PC Model subsections. For the Helium isoelectronic series, a new density with a finer grid (so more points) was also studied to remove numerical errors that could possibly arrive with less points on the grid. This was done for the first two elements in the series, H and He. Zcritwas also studied using the finer grid, but was not studied for the original
one.
Then the results of both the PC Model and the exact one were plotted and, if needed, a constant value was added to the PC Model to verify if a constant shift, which does not affect the final energies, makes the results more accurate. Then the relative and absolute errors for all values of r/a, with a being the Bohr radius, were calculated so that the accuracy of the PC Model could be measured more easily. The absolute error will be calculated using the following formula,
AbsErr(r) = |vpcxc(r) − vscexc(r)| , (3.1)
where the difference between the PC Model and the exact one is calculated for every r. For the relative error the equation below is used, which is
RelErr(r) = |v pc xc(r) − vxcsce(r)| −vsce xc(r) , (3.2)
where the difference between the two is divided by vxcsce. The minus sign is added due to vsce
xc (r) being negative, which would give negative relative errors.
For the SPL the functional derivatives of W0, W∞ and W∞0 are needed to be able to
calculate the functional derivative of WαSP L. The first functional is equal to exchange energy or, for a two-electron system, equal to minus one half of the Hartree Energy. The formulas for the second and third term can be found in equation 2.35 and 2.48
respectively. The WαSP L will be tested for the Helium isoelectronic series versus the exact physical exchange correlation potentials.
Results and Discussion
4.1
PC Model
The results will consist mostly of graphs and plots in order to show difference between the PC model potentials and the SCE exact potential at all values of r/a within the density of the atom. The relative and absolute errors are also reported for every atom. Not all graphs are reported in the results, due to the sheer number of graphs made. The rest of the graphs can thus be found in appendices A, B, C, and D.
4.1.1 Helium Isoelectronic Series
For the helium isoelectronic series most of the graphs can be found in appendix A. For most members of the series, the results in the appendix indicated that the PC model works well within the region where the density is significantly high. At small r-values, the PC Model has unexplainable peak in the graph. This is due to the LDA part of the PC model which lays significantly lower at small r-values than the exact graph. At higher values, within the density, the PC model does not converge to 0 but goes slowly to −∞, which causes the PC model to underestimate the exact potentials. This is true for all members, because the second and third term of equation 2.43 do not converge to 0. In short, the PC model works well for r-values in the middle of the density, but the LDA term at small r-values and the two other terms at high r-values decrease the overall accuracy of the PC model as approximation of the SCE exact potentials. The same can be extracted from the absolute and relative error graphs, where absolute errors at intermediate r-values only deviate between the 0.05 and 0.2 and for relative errors between the 0.05 and 0.1, whereas at higher r-values errors can be found of 1.0 and higher. For smaller r-values the errors are less significant but still much higher
than the errors at intermediate r-values. This shows that the PC Model works well for the intermediate r–values but lacks accuracy in the extreme parts. The LDA of the PC Model is much worse because it either underestimates the Exact SCE potentials at small values of r- or overestimates at intermediate and higher r-values. This problem can be seen in the graphs of every member of the helium isoelectronic series, which points toward it being a universal problem of the LDA.
Using a finer grid improves the PC model by removing the peak at small r-values, which means that it is a numerical error of the old grid. This has been studied for H– and He only, but it should also remove the numerical error peaks for other Z values. Even though the density graphs are identical for He, the PC Model graphs differ at higher r-values. This effect is relatively small and happens due to numerical anomalies of the old grid. However, it does not affect the results of the PC Model due to it existing only close to the end of the density for both He and H–.
The most important member of the Helium isoelectronic series, from SCE limit point of view, is H– due to having a low density and thus being strongly correlated. The PC model decently approximates the exact SCE potentials in the hydride-ion and thus can be used in other systems with the same low density and strong correlation. The relative and absolute errors are low for most of the densities and the vxc plots are also similar
at intermediate r-values. Another important low density, high correlation system was the above-mentioned Zcrit. The graphs of the PC Model and the exact SCE of Zcrit are
almost identical, which means that PC Model potentials work better for systems with low densities as the Zcrit and H–.
In all graphs of the PC Model a gap is visible at the nucleus between the LDA/GEA and the exact one. This gap grows for higher Z values and can thus scale with Z. This has been researched and the results can be found in the figure below. From the R2-value,
2 4 6 8 10 Z 0.5 1.0 1.5 2.0 2.5 Vxc TrendLine Zinf(x)
Figure 4.1: The energy difference between the exact one and the PC Model LDA/GEA at the nucleus vs the atomic number Z.
0.9992, it can be seen that the exchange correlation potential differences scale linearly with Z. This can be used to improve the PC Model potentials, because when a correction
is known for one atom, the correction can be scaled for the other members. The graph for Z = ∞ can also be seen in the figure above. This is the high density, low interaction limit due to fact that equation 2.53 can be rewritten by a coordinate change as
ˆ H = Z2 −1 2∇ 2 ~ s1− 1 2∇ 2 ~ s1 − 1 ~ s1 − 1 ~ s2 + 1 Z 1 | ~s1− ~s2| , (4.1) where ~r = ~s/Z. Here the last term goes to zero, when Z → ∞, which proves the low interaction. In this limit, the atomic number scales linearly with the exchange correlation potentials, according to reference [19], which gives us a way to check the scaling of the gap. Although the Z = ∞ grows faster than the linear trend line, they have the same order and from this it can be concluded that the energy difference scales linearly with Z.
It is known that the SCE potentials scale linearly with Z in the helium isoelectronic series. This means that the scaling of the PC Model with Z must also be linear, in order to be a good approximation in the SCE formalism. This can be done mathematically by inserting ρZ(~r) = Z3ρ(Z~r) into equation 2.43, which shows that the PC Model scales
linearly with Z. In the two graphs in the appendix, the scaling of the PC Model and exact exchange correlation potentials can be seen at the nucleus. These values are divided by Z, because when they scale linearly with Z, this should result into a constant. However, the wrong scaling is found, which is caused by the fact that the used density is the real one, which is not scaleable and thus shows deviation from the proposed scaling. Nonetheless, for bigger members the effect is less and both the exact and PC Model graphs should go asymptotically to the vlim(0)-values of -1.36357 and -1.66290 respectively. Because
both scale linearly with Z, their difference should be too and this is already found in the above paragraph.
4.1.2 Neutral Atoms
In Appendix B, the results and graphs of the neutral atoms can be found. For He, Li and Be the graphs of the PC model look similar to the graphs of the exact one, but lose accuracy at low and high values of r. The graphs of the other atoms need to be shifted by a constant before they align with the exact one. This means that the PC model is less accurate for the neutral atoms then for the Helium Isoelectronic series due to a needed shift for some atoms. Furthermore, all the graphs of the PC model for neutral atoms contain a strange offset at intermediate r-values caused by the terms containing the gradient and Laplacian of the density. However, this could be a numerical error, but it is most likely a shell effect.
It can still be found that the PC Model approximation does not work for small or high r-values due to overestimating the exact SCE potentials at lower r-values and underestimating at high r-values. The deviation at higher values is caused by the wrong asymptotic behaviour that the GEA’s have. From the absolute and relative errors the same conclusions as above can be reached, because at intermediate r-values both errors are significantly lower than anywhere else. Only for the heavier atoms, the relative and absolute errors are significantly big, which can be drastically improved by adding a constant. The error graphs can also be found in the appendix.
The LDA does not work as well as the GEA functionals, because it still overestimates the absolute SCE exact potential at higher values of r. However, when shifted it works much better than the GEA, due to having the right asymptotic behaviour. It is not researched in this thesis but the constant added could be either a function of the atomic number, a function of the density or something else. It is therefore interesting to find the shifting constant for every atom and research their scaling for Z or ρ.
4.1.3 Hooke’s Atom
The results of the Hooke’s atoms can be found in the graphs of appendix C. Since the Hooke’s atom densities decay faster than normal exponential decays, which means that the gradient terms are now larger, the PC Model potentials are now more sensitive to changes in the density. This effect can be seen in the figures in appendix C, where the PC Model does not work for Hooke’s atom potentials. The figures show that the GEA and LDA overestimate the exact potentials for all n- and r-values, except at higher r-values where the GEA has the wrong asymptotic behaviour. Furthermore, both the LDA and the GEA do not have the same shape as the exact potentials, so that no constants can be added for a better alignment. Where the GEA and LDA show ups and downs in their graphs for n = 3 and higher, the exact one has the same S-shape as in the other two series. This is again a failure of the PC model approximation, which causes its usefulness to decline.
The relative errors also show that the PC model does not work for Hooke’s atom, because the relative errors are mostly found between the 0.2 and 0.4, whereas in the other series it was around 0.1 or lower. Even though the Hooke’s atom is not a physical system and thus not of interest when applied to chemistry, it shows that the PC model has its limits and cannot calculate accurate potentials for every system. This means that another approximation for the exchange correlation potentials is needed to approximate Hooke’s atom or other similar systems.
4.2
SPL
The graphs showing the results of SPL can be found in appendix D. First the asymptotic behaviour of SPL must be examined, which can be found in the graph below: The
1 2 3 4 5 6 r a -10 -8 -6 -4 -2 vXCSPL
Figure 4.2: The asymptotic behaviour of SPL exchange correlation potentials for Helium vs r/a.
SPL exchange correlation potentials do not have the right asymptotic behaviour either, because they go to −∞ just as the PC model potentials. This comes from the fact that the W∞0 goes as e
1
6αr, whereas W∞ goes as −e 1
3αr, which means that the latter one goes
much faster to its limit than the former one. This means that the SPL does not fully fix the problem that the PC Model has, but it can delay the inevitable so that values inside the density are not affected. Moreover, the SPL approximation can still be used as an approximation for the physical system, whereas the PC model can only be used for the SCE limit.
Secondly, it is important to find out if WSP Lworks as an approximation for the exchange
correlation energies of the physical system. To test this, the exact exchange energies were subtracted from the SPL functionals, which yields EcSP L. The SPL correlation energies can then be compared to the exact ones, so that the accuracy of WSP L can be found. The comparison between ESP Lc and the exact Ec can be found in table
4.1 below. It can be seen that both the SPL-SCE and SPL-PC correlation energies differ significantly from the exact ones, especially for bigger atoms. Furthermore, the SPL-PC approximation works better for He and Be, where SPL-SCE one works better for Neon and the strongly correlated systems H– and Zcrit. From the table it can be
concluded that this version of SPL does not approximate the correlation energies well, due to the lack of W00. However, SPL can still work as an approximation for the exchange correlation energies and potentials, because the absolute values of the correlation energies are relatively low in comparison to exchange energies, which means that the effect of the wrong correlation energy is limited.
Table 4.1: Comparison between the SPL correlation energies obtained via either the exact procedure or the PC Model scheme and the exact correlation energies for Zcrit,
H–, He, Be and Ne. The exact Ex and Ec were obtained from reference [11] and the
W∞SCE and W∞0SCE from references [11,14].
Ex W∞P C W∞SCE W∞0P C W∞0Exact EcSP L−P C EcSP L−SCE EExactc
Zcrit -0.298 -0.448 -0.451 0.0943 0.0985 -0.044 -0.046 -0.049
H− -0.380 -0.556 -0.565 0.139 0.146 -0.041 -0.0434 -0.045 He -1.025 -1.46 -1.50 0.619 0.621 -0.044 -0.055 -0.042 Be -2.67 -3.96 -3.94 2.58 2.59 -0.071 -0.068 -0.096 Ne -12.08 -19.98 -20.04 22.94 22 -0.22 -0.244 -0.39 We can interpret W∞as the exchange correlation energy in the SCE limit, which means
that removing the exchange energy would yield the correlation energy in the SCE limit. It can be seen that even for the most correlated systems, e.g. Zcrit, both the PC and
the exact SCE correlation energies overestimate the absolute correlation energy, -0.15, -0.153 and -0.049 respectively. This means that the PC model and the exact SCE do not work as approximation for the physical correlation energies of atoms, even if the system is strongly correlated like Zcrit.
The results of the exchange correlation potentials for H–, He, Be and Ne have been obtained and the graphs can be seen in appendix D. But before those can be interpreted, W (λ) must be investigated first. When W (λ) is not a smooth function of λ, it cannot approximate the exchange correlation energies and in turn not the exchange correlation potentials, because they both scale smoothly with λ. Fortunately, W (λ) is a smooth function, which is shown in figure 4.3 for the Neon atom:
20 40 60 80 100 λ -17 -16 -15 -14 -13 -12 Wλ
Figure 4.3: The plot of W (λ) versus λ for Neon, which gives a smooth graph. The values for W0, W∞and W∞0 are obtained from reference [11]
For the two-electron-systems, H– and He, the SPL potentials overestimate the physical potentials at all r-values. This effect is the most significant for H– at small r-values, but can also be found at other r-values for both the atoms. It can be seen in the graph of Helium that, again, the wrong asymptotic behaviour affects the results of the
helium isoelectronic series. However, this effect is much weaker due to the functional derivative of the W∞0 counteracting the W∞ term. The two bigger neutral atoms, Be
and Ne, show something different. Here the SPL is almost completely similar to the exact exchange correlation potentials, especially for Neon. For Beryllium, the SPL approximation holds for almost any r-value, only close to the nucleus and at the peak of the exchange correlation potentials does it differ from the exact one. It can also be seen that the wrong asymptotic behaviour of the SPL does not affect the neutral atoms inside the density. This means that the SPL approximation for the exchange correlation potentials is a decent approximation for the helium isoelectronic series and a better one for the neutral atom, especially for high atomic numbers.
From the figures in the appendix for Helium, Beryllium and Neon it can be extracted that the graphs of the exact vx, the exact vxc and vxcSP L are very similar to each other.
This means that for the neutral atoms, the exchange part has the highest contribution, especially for bigger atoms, hence it is the main contributor to the exchange correlation potential. SPL is thus a good approximation due to its ability to approximate the exchange part of the exchange correlation potential. Therefore, the errors of the SPL must lay in the correlation potential, which are also calculated and the resulting graphs can be found in appendix D. This has been done by subtracting the exact exchange potentials from both the vSP L
xc and the exact vxc. It can be seen that for all three
atoms, the correlation potentials cannot be accurately approximated by SPL. For Helium and Beryllium, the SPL approximation overestimates the exact correlation energy, but has approximately the same shape as the exact one, apart from the wrong asymptotic behaviour. For Neon, the errors are much bigger and the shape of the SPL graph does not correlate with the exact one either. This means that the SPL approximation works decently for the exchange correlation potential, but errors in the exchange correlation potential come from overestimating the correlation potential.
Conclusions
It can be concluded that the PC Model can work as an approximation for the SCE exchange correlation potentials, but it has its limitations. It cannot approximate the Hooke’s atom, but it worked well for both the helium isoelectronic series and the neutral atoms. However, for the latter a constant needs to be added for atoms with higher atomic numbers. Especially for systems with a low density, and thus high correlation, such as the negative ions with both N = 2 and Z = 1 and Z = Zcrit respectively, the PC
Model works extremely well and it is thus a good approximation for strongly correlated systems. The biggest problem of the PC Model is its asymptotic behaviour, which goes to −∞ instead of going to 0. This means that it incorrectly shows the behaviour at high r-values, which causes problems especially for atoms with high atomic numbers in the helium isoelectronic series. The small differences at small radii are less important due to their relatively short range.
Where the PC Model was a good approximation for exchange correlation potentials in the strong interaction limit, the SPL approximation works decently for the physical system. It still does not show the right asymptotic behaviour in comparison to the exact exchange correlation energy, but its effect is weaker within the density than the PC Model one. It works really well for the neutral atoms, especially with high atomic numbers, but it overestimates the exact potentials for members of the helium isoelectronic series. However, SPL does not approximate the correlation potentials good enough, meaning that the errors in the exchange correlation potentials are caused by the correlation part. Future research can be done into the exchange correlation energies in the SCE limit for small molecules, such as H2 and LiH, which are computationally more demanding but also chemically more interesting. Furthermore, one can increase the accuracy of the PC model by inventing a GGA functional form of the PC Model and testing this for the helium isoelectronic series and the neutral atoms. This is important because GGA’s have
the right asymptotic behaviour and can thus solve the problems the PC model has for these potentials. Moreover, better approximations are needed for the Hooke’s atom in order to approximate the SCE exchange correlation energy due to the failure of the PC Model. The SPL approximation must also be examined more carefully for more neutral atoms and members of the helium isoelectronic series to give a more general conclusion, but also to gain more knowledge of its accuracy. Furthermore, more research into better approximations for the physical correlation potentials is also needed, due to the failure of the SPL approximation for these potentials.
The PC model can thus be used to calculate the exchange correlation potentials for systems in the SCE limit, such as H–. Although the exact SCE potentials use non-local densities, the PC Model uses semi-local densities, which are much easier to calculate and thus take less time. This will speed up many computational chemistry calculations on or similar to systems in SCE limit. However, before the PC Model can be used in com-putational chemistry, a GGA must be invented and tested first. These functionals are promising due to the already high accuracy of the PC Model apart from the asymptotic behaviour. The GGA will give computational chemistry a broader range of reactions and compounds to simulate, without causing calculations to take days. For good ap-proximations on the physical exchange correlation potentials, the SPL approximation can be used, which is also easily calculable. Its GGA, again, can speed up the reactions on physical systems without losing much accuracy on the exact systems. This can hope-fully explain some systems that could, before the introduction of the PC Model or SPL, not be coupled via the Kohn Sham formalism with the standard approximations.
Results for the Helium
Isoelectronic series
2 4 6 8 10 12 14 r a 0.1 0.2 0.3 0.4 0.5 0.6 DensityFigure A.1: The density vs r/a graph of H– 2 4 6 8 10 r a -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.2: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the H–
2 4 6 8 r a 0.05 0.10 0.15 0.20 0.25 RelativeError
Figure A.3: The relative error graph of H– versus r/a. 2 4 6 8 10 r a 0.05 0.10 0.15 0.20 0.25 AbsoluteError
Figure A.4: The absolute error graph of H– versus r/a.
1 2 3 4 5 r a 0.5 1.0 1.5 Density
Figure A.5: The density vs r/a graph of He 1 2 3 4 5 r a -5 -4 -3 -2 -1 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.6: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the He
0.5 1.0 1.5 2.0 2.5 3.0 r a 0.1 0.2 0.3 0.4 0.5 0.6 RelativeError
Figure A.7: The relative error graph of He versus r/a. 1 2 3 4 r a 0.2 0.4 0.6 0.8 1.0 AbsoluteError
Figure A.8: The absolute error graph of He versus r/a. 1 2 3 4 5 r a 0.5 1.0 1.5 2.0 2.5 Density
Figure A.9: The density vs r/a graph of Li+ 0.5 1.0 1.5 2.0 2.5 r a -4 -3 -2 -1 vXC
Figure A.10: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the Li+
0.2 0.4 0.6 0.8 1.0 1.2 1.4 r a 0.05 0.10 0.15 0.20 RelativeError
Figure A.11: The relative error graph of Li+ versus r/a.
0.5 1.0 1.5 2.0 r a 0.2 0.4 0.6 0.8 AbsoluteError
Figure A.12: The absolute error graph of Li+ versus r/a.
1 2 3 4 5 r a 1 2 3 4 Density
Figure A.13: The density vs r/a graph of Be2+ 0.2 0.4 0.6 0.8 1.0 1.2 1.4 r a -5 -4 -3 -2 -1 0 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.14: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the Be2+
0.2 0.4 0.6 0.8 1.0 r a 0.05 0.10 0.15 RelativeError
Figure A.15: The relative error graph of Be2+ versus r/a.
0.2 0.4 0.6 0.8 1.0 1.2 r a 0.2 0.4 0.6 0.8 AbsoluteError
Figure A.16: The absolute error graph of Be2+ versus r/a.
1 2 3 4 5 r a 1 2 3 4 5 Density
Figure A.17: The density vs r/a graph of B3+ 0.2 0.4 0.6 0.8 1.0 r a -6 -4 -2 0 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.18: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the B3+
0.1 0.2 0.3 0.4 0.5 0.6 0.7 r a 0.02 0.04 0.06 0.08 0.10 0.12 0.14 RelativeError
Figure A.19: The relative error graph of B3+versus r/a.
0.2 0.4 0.6 0.8 1.0 r a 0.2 0.4 0.6 0.8 AbsoluteError
Figure A.20: The absolute error graph of BB+ versus r/a.
1 2 3 4 5 r a 1 2 3 4 5 6 Density
Figure A.21: The density vs r/a graph of C4+ 0.2 0.4 0.6 0.8 1.0 r a -10 -8 -6 -4 -2 0 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.22: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the C4+
0.2 0.4 0.6 0.8 r a 0.1 0.2 0.3 0.4 0.5 0.6 0.7 RelativeError
Figure A.23: The relative error graph of C4+ versus r/a.
0.2 0.4 0.6 0.8 1.0 r a 0.5 1.0 1.5 2.0 2.5 3.0 AbsoluteError
Figure A.24: The absolute error graph of C4+versus r/a.
1 2 3 4 5 r a 1 2 3 4 5 6 7 Density
Figure A.25: The density vs r/a graph of N5+ 0.2 0.4 0.6 0.8 r a -10 -8 -6 -4 -2 0 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.26: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the N5+
0.1 0.2 0.3 0.4 0.5 r a 0.02 0.04 0.06 0.08 0.10 0.12 RelativeError
Figure A.27: The relative error graph of N5+ versus r/a.
0.2 0.4 0.6 0.8 r a 0.5 1.0 1.5 2.0 2.5 AbsoluteError
Figure A.28: The absolute error graph of N5+ versus r/a.
1 2 3 4 5 r a 2 4 6 8 Density
Figure A.29: The density vs r/a graph of O6+ 0.2 0.4 0.6 0.8 r a -12 -10 -8 -6 -4 -2 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.30: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the O6+
0.1 0.2 0.3 0.4 0.5 r a 0.05 0.10 0.15 0.20 0.25 0.30 RelativeError
Figure A.31: The relative error graph of O6+ versus r/a.
0.2 0.4 0.6 0.8 r a 1 2 3 4 5 6 7 AbsoluteError
Figure A.32: The absolute error graph of O6+ versus r/a.
1 2 3 4 5 r a 2 4 6 8 Density
Figure A.33: The density vs r/a graph of F7+ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 r a -15 -10 -5 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.34: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the F7+
0.1 0.2 0.3 0.4 r a 0.05 0.10 0.15 RelativeError
Figure A.35: The relative error graph of F7+ versus r/a.
0.1 0.2 0.3 0.4 0.5 0.6 r a 0.5 1.0 1.5 2.0 2.5 3.0 AbsoluteError
Figure A.36: The absolute error graph of F7+ versus r/a.
1 2 3 4 5 r a 2 4 6 8 10 Density
Figure A.37: The density vs r/a graph of Ne8+ 0.1 0.2 0.3 0.4 0.5 0.6 r a -15 -10 -5 0 vXC vXCSCE(r) vXCPCG(r) vXCPCLDA(r)
Figure A.38: The exact vXCSCE, vXCPC and vXCLDA versus r/a graph of the Ne8+
0.1 0.2 0.3 0.4 r a 0.1 0.2 0.3 RelativeError
Figure A.39: The relative error graph of Ne8+ versus r/a.
0.1 0.2 0.3 0.4 0.5 0.6 r a 1 2 3 4 5 6 7 AbsoluteError
Figure A.40: The absolute error graph of Ne8+ versus r/a.
1 2 3 4 5 r a 0.5 1.0 1.5 Density Finer Grid Grid Old
Figure A.41: The density vs r/a graph of He for the old grid
vs a finer grid. 1 2 3 4 5 vXC -4 -3 -2 -1 r a OldGrid NewGrid
Figure A.42: The vXCPC graphs for the old and a finer grid versus r/a graph of He.
