Distorted Space and Multipoles in Electronic Structure Calculations

Full text

(1)Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 601. Distorted Space and Multipoles in Electronic Structure Calculations FREDRIK BULTMARK. ACTA UNIVERSITATIS UPSALIENSIS UPPSALA 2009. ISSN 1651-6214 ISBN 978-91-554-7407-2 urn:nbn:se:uu:diva-9553.

(2)

(3)

(4)

(5)

(6)

(7)

(8) ! " # $ %& '( ) ) ) * ! +

(9)

(10) ,

(11)

(12) -

(13) ! . / "! %&! 0

(14) #

(15) -

(16) 0 1

(17) ! 2.

(18)

(19) !

(20)

(21)

(22)

(23) $! $3 ! ! 40.5 &67&7((376367%! +

(24)

(25) )

(26)

(27)

(28)

(29) ) ! +

(30)

(31) , 82*9:

(32) ; 2*9 8

(33) 2*9:

(34) 2*9< 8 :

(35) , )

(36)

(37) ! =. )

(38) ) 2*9

(39)

(40) )

(41) ) )

(42)

(43) ) .

(44)

(45) ! 2

(46) ), , )

(47)

(48) ) ! +

(49)

(50) 82:

(51)

(52) )

(53)

(54) ! 2

(55) )

(56) /

(57) ,

(58)

(59) 7

(60)

(61) )

(62)

(63)

(64) > ) ,

(65) ) ! 2 )

(66) ,

(67)

(68) 7 ))

(69) ! .

(70)

(71)

(72)

(73) 7 ))

(74) .

(75)

(76)

(77)

(78) )

(79)

(80) ) 2< )

(81)

(82)

(83) ! 4

(84)

(85) )

(86) ) 2<

(87)

(88)

(89) , , )

(90)

(91) )

(92) )

(93)

(94)

(95) )

(96) )

(97)

(98)

(99) )

(100)

(101) ! 2

(102) .

(103) )

(104) ) =

(105) ; )

(106) ,

(107)

(108)

(109)

(110)

(111)

(112) ) ? ;

(113)

(114)

(115)

(116) ! *

(117)

(118)

(119)

(120) )

(121)

(122) #

(123) !

(124) "

(125) # $ " %&

(126) ' " ()* +" ! , -./" " 012-+3+ " @ "/ . / %& 4005 $(7$%3 40.5 &67&7((376367%

(127) '

(128)

(129) ''' 7&((A 8 'BB

(130) !/!B C

(131) D

(132) '

(133)

(134) ''' 7&((A:.

(135) if you ﬁnd any grammatical errors you may keep them.

(136) List of Papers. This thesis is based on the following papers, which are referred to in the text by their Roman numerals. I F Bultmark, K Dewhurst, DJ Singh, L Nordström Tests of the efficiency of an augmented distorted planewave basis in electronic structure calculations J. Phys.: Condens. Matter 20, 235241 (2008) II Fredrik Bultmark, Francesco Cricchio, Oscar Grånäs, Lars Nordström Multipole decomposition of LDA+U energy and its application to actinides compounds in manuscript III Francesco Cricchio, Fredrik Bultmark, and Lars Nordström Exchange energy dominated by large orbital spin-currents in δ -Pu Phys. Rev. B 78, 100404 (2008) IV Francesco Cricchio, Fredrik Bultmark, Oscar Grånäs, Lars Nordström Itinerant non-collinear magnetic multipole moments of rank five, triakontadipoles, as the hidden order in URu2 Si2 in manuscript V O. Grånäs, F. Bultmark, F. Cricchio, L. Nordström Analysis of dynamical exchange and correlation in terms of coupled multipoles in manuscript VI Y. Krockenberger, K. Mogare, M. Reehuis, M. Tovar, M. Jansen, G. Vaitheeswaran, V. Kanchana, F. Bultmark, A. Delin, F. Wilhelm, A. Rogalev, A. Winkler, L. Alff Sr2 CrOsO6 : End point of a spin-polarized metal-insulator transition by 5d band filling Phys. Rev. B 75, 020404 (2007) Erratum: Sr2 CrOsO6 : End point of a spin-polarized metalinsulator transition by 5d band filling Reprints were made with permission from the publishers..

(137) Related work. I C. Caleman , C. Ortiz, E. Marklund, F. Bultmark,M. Gabrysch, F. G. Parak, J. Hajdu, M. Klintenberg and N. Tîmneanu Radiation damage in biological material: Electronic properties and electron impact ionization in urea EPL 85, 18005 (2009) II F. Bultmark, J. Lööf, L. Hermansson, H. Engqvist Continuous hydration of low w/c calcium-aluminate cement submitted to Journal of Materials Science III Stina Uppsala, 25-11-1994. IV Tora Danderyd, 31-01-2000. V Vaste Danderyd, 21-10-2005. VI Vide Danderyd 17-10-2008.

(138)

(139) Contents. 1 2. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Electronic structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 One particle equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Density functional theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Local density approximation . . . . . . . . . . . . . . . . . . . . . . 2.3 Augmented planewaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Muffin tin spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Core electrons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Linearised basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Semicore states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.1 Local orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 APW+lo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Augmented distorted planewaves . . . . . . . . . . . . . . . . . . . . . . . 2.10 Adaptive coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Distorting function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Distorted augmented plane waves . . . . . . . . . . . . . . . . . . . . . . 2.13 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Magnetism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.0.1 Magnetic order . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Intra-atomic non collinear magnetism . . . . . . . . . . . . . . . . . . . 4 Multi pole LDA+U . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Multi pole momentum tensor formulation . . . . . . . . . . . . . . . . 4.2 Interpretation of the tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Saturation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Polarisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Hund’s and Katt’s rules . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 9 11 13 14 15 16 16 18 18 20 20 21 21 22 23 24 26 31 33 33 36 39 40 44 45 48 50 59.

(140)

(141) 1. Introduction. In 1897 Thomson [1] proved that electric current consists of negatively charged particles moving under the influence of the electric field - the electrons, as proposed by Stony in 1891. The electrons where later found to come from the studied materials themselves. Matter had small negatively charged particles in it. Since matter in general is not electrically charged there had to be something with positive charge too. In 1911 Rutherford interpreted the experiment of Geiger and Marsden [2] and suggested the existence of small, positively charged nuclei around which the electrons are moving. Shortly after Bohr [3] published his model of quantised energy levels for the electrons, followed by the works of Schrödinger and Heisenberg which led to the theories we today call quantum mechanics. In quantum mechanics the particles are represented by wavefunctions. The wavefunction is a mathematical construction that holds all information that can be known about a system. To extract information operators are applied to the wavefunction. This abstract mathematical model has proven to correctly predict and explain many experimental results and is considered a good model of the interactions in the small size and/or large energy regime. The most important quantum mechanical operator is the total energy operator E. The total energy can also be described as the sum of the kinetic energy and the potential energy, the Hamiltonian, H = T + K. Through this relation the Schrödinger equation is formulated, HΨ = EΨ , (1.1) where Ψ is the total wavefunction of the system. The total energy of a system is a function of the configuration of all particles in the system. Different configurations has different energies. One (or sometimes more) configuration has the smallest total energy. This configuration is of special interest since it is the state in which you would find the system if it was isolated from the rest of the universe for some time, the rest configuration or the ground state. A real system in the real world is never left without interaction with the rest of the universe, and thus a system in its ground state can rarely be observed, but many of the properties of a system in its ground state will be preserved when the system is put in contact with the surroundings. When searching for the ground state configuration, any constant in the Hamiltonian may be omitted since it will have no influence on which configuration has the lowest energy. In theoretical materials physics the properties of macroscopic materials are predicted by calculations or simulations of materials. It is both impossible 9.

(142) and unimportant to keep track of single particles in a huge collection of interacting particles. By approximating matter as organised in infinite perfect crystals, symmetry arguments can transform the problem into a problem concerning only a few tens or hundreds particles obeying periodic boundary conditions. The many particle equations may under certain circumstances be transformed into one particle equations which may be solved numerically. Usually this is done by expressing the wavefunctions in some set of basis functions, ψ = ∑N cN φN , and varying the coefficients cN looking for the minimum value of HΨ. In this thesis augmented planewave basis sets will be discussed and a new basis set is suggested. Implementation of the basis set in planewave codes and numerical results are presented. Throughout the thesis Hartree atomic units are used, for which h¯ = 1, the electron mass me = 1, the Bohr radius a0 = 1, the electron charge e = 1. This means that the speed of light c = 1/α ≈ 137 and the dielectric constant 4πε0 = 1. Energy is given in Hartree (≈ 27.2eV) and distance in Bohr radii (a0 ≈ 0.529Å). Hopefully I manage to use the following conventions for wavefunctions Ψ φ ψ. total many electron wave f unction sinlge electron wave f unction basis f unction. (1.2) (1.3) (1.4). Mathematic notation used in the thesis [a] = (2a + 1) [ab] = (2a + 1)(2b + 1) 1 2 3 [1 2 3 ] 1 2 3 = Y1 m1 Y2 m2 Y3 m3 dΩ 4π 0 0 0 m1 m2 m3 a b c : Wigner 6 j symbol d e f. 10.

(143) 2. Electronic structure. We wish to solve the many particle Schrödinger equation HΨ(r1 , . . . , rN , r1 , . . . , rM ) = EΨ(r1 , . . . , rN , r1 , . . . , rM ) ,. (2.1). where r1 , . . . , rN describe the position of the N electrons and r1 , . . . , rM denote the position of the M nuclei. The total wavefunction Ψ holds all information of the electronic structure (and thus of the properties) of the material. The Hamiltonian H describes the particles interaction with each other and E is the total energy of the system. Consider the many particle Hamiltonian describing any system of interacting electrons and nuclei, H = Tn + Te +Vnn +Vne +Vee ,. (2.2). where Tn is the kinetic energy of the nuclei, Te is the kinetic energy of the electrons, Vnn is the energy of the nuclei-nuclei interaction, Vne is the energy of the nuclei-electron interaction and Vee is the energy of the electron-electron interaction. The most significant difference between the nuclei and the electrons is that the nuclei are heavy and move slowly in an approximately harmonic potential causing crystal vibrations. The electrons are light and move fast around approximately stationary nuclei. In the Born-Oppenheimer approximation we let the nuclei be stationary, that is Tn = 0 and Vnn = constant. The kinetic. Figure 2.1: Periodic array of identical cells.. 11.

(144) energy operator for N electrons in quantum mechanics is N 1 Te = ∑ − ∇2i , 2 i. (2.3). where ∇i is acting on ri , the position of electron i. The Coulomb attraction of the N electrons by the M nuclei is N M. Vne = ∑ ∑ − i. j. Zj , |ri − r j |. (2.4). where Z j is the charge of nucleus j. The Coulomb repulsion of the electrons by the other electrons is Vee =. 1 1 , ∑ ∑ 2 i j=i |ri − r j |. (2.5). where the factor 12 takes care of the fact that in the sum we count each pair of electrons twice. The many particle Schrödinger equation can now be written. N Zj 1 1 2 M 1 ∑ − 2 ∇i − ∑ |ri − r j | + 2 ∑ |ri − r j | Ψ(r1 , r2 . . . , rN ) = i j j=i = EΨ(r1 , . . . , rN ). (2.6). In general we can not solve this equation, because the number of particles is too large. We need more approximations. If we have an infinite periodic crystal, it is easy to see that the potential from the nuclei Vne must be periodic. If r is a vector between a point in one of the cells in Fig. 2.1 to the same place in one of the other cells, a translation of the lattice by r will reproduce the same potential. The one electron solutions to the Schrödinger equations for such periodic potential will also be periodic, apart for a phase shift, φ (ri + ri ) = φ (ri )eik·ri .. (2.7). If the total wavefunction is a product of single electron wavefunctions, the total wavefunction will also be periodic with a phase that is a product of all the phases of the single electron wavefunctions. If we identify the smallest possible cell that with periodic boundary conditions build up the whole crystal we have found the primitive unit cell. If we solve the problem for the primitive unit cell, with periodic boundary conditions satisfied, we have the solution everywhere.. 12.

(145) 2.1. One particle equations. It is still computationally demanding to solve the many particle problem from equation 2.6. It would be much simpler to solve N (coupled) single particle equations. There are several ways of mapping the many particle problem in Eqn. 2.6 to single particle equations. In the Hartree approximation [4] we assume that the total wavefunction can be written as a product of individual electron wavefunctions Ψ(r1 , r2 , . . . , rN ) = φ1 (r1 )φ2 (r2 ) . . . φN (rN ). The total energy of the system is the expectation value of the Hamiltonian, N. Ψ| H |Ψ 1 ∗ ∗ ∗ HΨ = = φ1 (r1 )φ2 (r2 ) . . . φN (rN ) ∑ − ∇2i − Ψ| Ψ 2 i=1. M Zj 1 1 ++ ∑ φ1 (r1 )φ2 (r2 ) . . . φN (rN )dr1 dr2 . . . drN . −∑ 2 j=i |ri − r j | j=1 |ri − r j | (2.8) Assuming normalised one electron wavefunctions φi | φi = 1.To find the ground state wavefunctions we should minimise the expectation value of the Hamiltonian with respect to the one electron wavefunctions φi , d HΨ = ∑ i. ∂ HΨ dφi = 0 . ∂ φi. (2.9). The total number of electrons should be constant |φi |2 dr = 1, ∀i, and thus ∂ d1 = ∂ φi. . |φi |2 drdφi = 0 .. (2.10). Summing these two differentials, multiplying the last with Lagrange multipliers εi (undetermined constants), we get

(146). ∂ ∗ ∑ ∂ φi HΨ + εi φi (ri )φi (ri )dri dφi = 0 i M Zj 1 ∂ (2.11) ∑ ∂ φi φi∗ (ri ) − 2 ∇2ri − ∑ |ri − r j | + i j=1. 1 1 ∗ + ∑ φ j (r j ) φ j (r j )dr j − εi φi (ri )dri = 0 . 2 j=i |ri − r j | Since each of the i terms in the sum above regards an independent variable φi , each of the terms must be zero,. φ ∗ (r )φ (r ) Zj 1 2 M 1 j j j j − ∇ −∑ + ∑ dr j φi (r) = εi φi (r) . (2.12) 2 2 j=i |r − r j | j=1 |r − r j |. 13.

(147) This is the Hartree equation for occupied one electron states φi (r) with eigenvalues εi . If the spin degree of freedom is included for each electron, the Pauli principle states that the total wavefunction should be anti symmetric with respect to interchanging two particles and reversing their spins. A variable separated trial function that fulfils this is a Slater determinant φ (r , s ) φ (r , s ) · · · φ (r , s ) 1 2 2 1 N N 1 1 1 φ2 (r1 , s1 ) φ2 (r2 , s2 ) · · · φ2 (rN , sN ) . Ψ(r1 , s1 , r2 , s2 , . . . , rN , sN ) = .. .. .. .. . . . . φN (r1 , s1 ) φN (r2 , s2 ) · · · φN (rN , sN ) (2.13) If the same arguments as before are used on the new trial function we end up with the Hartree-Fock equations [6]. 1 |φi (r)|2 Zr + (2.14) εi = ∑ φi∗ (r) − ∇2 φi (r)dr − ∑ 2 |r − r| r i. ∑ ij. . |φi (r)|2 |φi (r )|2 drdr − ∑ δsi ,s j |r − r | ij. φ ∗ (r)φ (r )φ ∗ (r )φ (r) i j i j. |r − r |. drdr .. The first term is the kinetic energy and electrostatic interaction with the nuclei of the single electron wavefunction φi (r, si ), the second term is the (repulsive) Coulomb interaction of the single electron wavefunction φi (r, si ) with all the other electrons. The last term is the exchange term that reduces the repulsive interaction between electrons of the same spin. This reduction is needed due to the Pauli principle which states that electrons with the same spin cannot be arbitrarily close to each other. The error in energy in the Hartree-Fock approximation is called the correlation energy, and if we knew the correlation potential, it could be added to equation 4.1 which then would give the exact solutions.. 2.2. Density functional theory. Another way to map the many particle problem to a single particle problem is density functional theory (DFT). In 1964 Hohenberg and Kohn [9] showed that the electron density uniquely defines the potential of a system. The potential in turn defines the Hamiltonian which can be solved for eigenvalues and eigenstates. The occupied eigenstates determines the density. The sum of the eigenvalues of the occupied states gives the total energy. This means that the total energy can be written as a functional of the electron density only. Later Kohn and Sham [10] showed that minimising the total energy leads to an independent set of single particle equations, (T +V )Φn = En Φn .. 14. (2.15).

(148) The eigenstates of the Kohn-Sham equation (equation 2.15) are called pseudo wavefunctions since they are not solutions for any specific electron (there is no assumption of separation of variables), just a solution that together with the other occupied eigenstates gives the correct charge density. The Kohn-Sham eigenvalues are in the same way not real energies for real electrons in a material, but the total energy of the pseudo particle system is minimised for the same density as the real electron density. In principle we can not expect the spectral function (density of states) to correspond to the real one, but surprisingly often the spectral function of the pseudo particles looks similar to the real spectral function obtained from experiments. This might seem strange at first but it is actually in many ways a better approach than the Hartree-Fock approximation, since in general DFT is exact - if we only knew the exact form of the total energy functional. The unknown part of the total energy functional is called the exchange-correlation energy.. 2.2.1. Local density approximation. The most common exchange-correlation functional in DFT calculations is the functional obtained by the local density approximation. This crude approximation suggested by Slater [8], that assumes that the electron density is so slowly varying that it can be treated as a uniform electron gas with density ρ = ρ(r). This approximation leads to the following Hamiltonian for the Kohn-Sham equations (eq 2.15),.

(149) 3ρ(r) 1/3 ρ(r ) dr − 3 . (2.16) H = T +Vne (r) + |r − r | 8π The results agree surprisingly well with both experiments and more exact calculations and to this day it is still a widely used approximation. Whether we choose to follow the Hartree-Fock route to wave-functions for real particles and with exact exchange potential or the DFT route to pseudowavefunctions using approximated exchange correlation potentials we end up with Schrödinger-like equations, HΦn = En Φn ρ=. (2.17). ∑ |Φn |. 2. (2.18). occ.n. (2.19) The total energy is minimised by the ground state density (ground state wavefunctions). If we expand the wavefunctions in a complete basis, with linearly independent basis functions Φn = ∑ cin ψi .. (2.20). i. 15.

(150) The problem transforms to finding the expansion coefficients cin that minimises the total energy.. 2.3. Augmented planewaves. There are many possible basis sets to choose from. The perhaps most obvious basis set is to use the atomic orbitals for all atoms in the sample. This is called LCAO (linear combination of atomic orbiltals) basis. Another basis set that might come into mind is the complete set of planewaves. The LCAO basis has the advantage that the basis functions resemble the solution, which speeds up the process of finding the right coefficients. Unfortunately the LCAO basis is not complete and can not fully describe the electronic structure of a solid. The planewaves are in principle a complete basis set - we have to truncate the series somewhere but we can get as exact solutions as we like by just adding more basis functions. Unfortunately the only system that has a wavefunction that resembles planewaves is the electron gas in a constant potential, all other systems need a lot of planewaves to be described properly. The best basis set would be a set of functions that looks like the eigenstates of the atomic problem close to the atom and is complete in the interstitial.. 2.4. Muffin tin spheres. The Coulomb potential from the nuclei in a solid is almost spherical close to the atoms and slowly varying (almost constant) between the atoms. In a simplified picture the potential would look like a muffin tin - cups of spherical potentials embedded in a constant potential (Fig. 2.2). In the muffin tin (MT) approximation we let each atom in the unit cell be surrounded by a sphere. The radii of the spheres should be chosen so that they do not overlap. The part of space not included in any sphere is called the interstitial region. Inside the spheres the potential is almost spherical. In the interstitial region the potential is almost flat. We know that solutions to the Schrödinger equation in a spherical potential are radial functions times spherical harmonics, ψ(r) = ∑m am f (r)Ym (ˆr). In a constant potential the solutions to the Schrödinger equation are planewaves ψ(r) = Ω−1/2 ∑G cG eiG·r , with Ω as the unit cell volume. This made Slater [7] propose augmented planewaves defined as planewaves in the interstitial and spherical functions inside the MT spheres , r ∈ interstitial Ω−1/2 eiG·r φG (r) = . (2.21) G ∑m am f (r)Ym (ˆr) , r ∈ MT The radial functions f (r) are solutions to the radial Schrödinger equation . l(l + 1) 1 d2 + +V (r) − E r f (r) = 0 , (2.22) − 2 dr2 2r2 16.

(151) Figure 2.2: Potential in the muffin tin approximation. where E are parameters that determine the radial functions for the solutions of each l character. The expansion coefficients cG and am should be chosen so that the wavefunction is continuous at the MT sphere boundary. Consider an atom at position Rα with a MT radius RMT . A general point on the sphere boundary will be rMT = Rα + RMT rˆ α . If we expand a planewave in spherical harmonics at the MT sphere boundary we will see that 4π ∗ ˆ (G)Ym (ˆrα ) = φG (RMT ) = √ eiG·Rα ∑ i j (GRMT )Ym Ω m. 4πi iG·Rα ∗ ˆ ˆ MT ) , =∑ √ e j (GRMT )Ym (G) Ym (R Ω m and the matching coefficients aG m should be chosen as. 4πi iG·Rα 1 G ∗ ˆ √ e j (GRMT )Ym (G) . am = f (RMT ) Ω. (2.23). (2.24). Note that the matching coefficient is undefined if the radial function disappears at the sphere boundary, fl (RMT ) = 0. If there are bands in the region of the energies that produce radial functions that are zero at the MT sphere, good solutions might be hard to find. The APW basis set works well if we have only one band for each -character. If there are several bands of the same -character, like the 2s and 3s, or Eg and T2g d-bands, the energy dependence of the basis set will become a problem. 17.

(152) 2.5. Core electrons. For large atoms with many electrons, the lowest lying states for the atom will be well localised inside the MT spheres. Electrons in these states will experience a strong Coulomb attraction to the nucleus and the presence of other nuclei will be screened by the outer electrons. The high kinetic energy of these core electrons makes relativistic effects important and new terms appear in the Hamiltonian that will only make a difference for the core states. Instead of solving the Schrödinger equation for these states, we solve the Dirac equation, −iα · ∇ + c2 β +V (r) Ψ = EΨ α=. 0 σ1 σ1 0. ,. 0 σ2 σ2 0. ,. 0 σ3 σ3 0. . ,β =. I 0 0 −I. . (2.25). The solutions to the Dirac equation can be written. gκ (r)χκ μ (ˆr) , Ψ(r) = −i fκ (r)σ χκ μ (ˆr). (2.26). where κ is the relativistic quantum number, μ its projection on the z-axis, χκ μ (ˆr) are spin-angular functions, that is spherical harmonics times two component spinors. The potential used is the first component of a spherical harmonics expansion of the true potential, V (r) = V (r)| Y00 . Some distance from the nuclei, the wavefunctions of the core electrons will be small, and at the MT sphere the core wavefunctions and their derivatives will be assumed to be zero. In this way we do not need to match them to any interstitial function or core wavefunction of any other atom. We will also assume that the core wavefunctions are orthogonal to any valence state. This is not entirely true, since the core wavefunctions will be small but finite outside the MT spheres. The overlap between core and valence states will in general be small. The charge density of the core states outside the MT spheres, the core leakage, may be calculated. The core leakage should be small, otherwise electrons has to be moved from the core to the valence. The use of core states will solve many of the problems of the energy dependent basis set for the valence states. All states that are (almost) confined to the MT spheres are eliminated from the valence state problem.. 2.6. Linearised basis functions. The APW basis set require that the energies E in equation 2.22 equal the band energy. In general there are several valence bands of the same -character due to the non spherical potential in the interstitial region (the crystal field splitting) and relativistic effects in the MT spheres. A method of solving the prob18.

(153) lem of the energy dependence of the basis set was invented by Andersen [11] and the basic idea is to approximate the radial functions of equation 2.22 with functions that are linear in the energy parameter, g (r, E ) = f (r, E0 ) + f (r, E0 )(E − E0 ) Here we need the energy derivatives of the radial (r, E) = ∂ f (r, E), evaluated at the energy E 0 satisfying fm ∂E . 2 d ( + 1) − 2+ +V (r) − E r f (r, E ) = r f (r, E ) . dr 2r2 The solutions are now written , r ∈ interstitial Ω−1/2 eiG·r G φG r) = . G ∑m am f (r) + bm f (r) Ym (ˆr) , r ∈ MT. (2.27) functions. (2.28). (2.29). G Here we have two parameters, aG m and bm , to match the basis function at the MT sphere boundary. To determine these two matching coefficients both the basis function and its first derivative are made continuous at the sphere boundary,. 4π ∗ ˆ ˆ MT ) = φG (rMT ) = √ eiG·rα ∑ i j (GRMT )Ym (G)Ym (R Ω m . 4πi iG·rα ∗ ˆ ˆ MT ) = =∑ √ e j (GRMT )Ym (G) Ym (R Ω m G G ˆ MT ) = ∑ am f (RMT ) + bm f (RMT ) Ym (R m 4π iG·rα ∂ φG (r) d j (Gr) ∗ ˆ ˆ MT ) = √ e = i Ym (G)Ym (R ∑ ∂ r r=rMT dr Ω r=R m MT . 4πi iG·rα d j (Gr) ∗ ˆ ˆ MT ) = =∑ √ e Ym (G) Ym (R dr Ω r=R m MT . (r) d f d f (r) Ym (ˆrM T ) . = ∑ aG + bG m m dr dr r=RMT r=RMT m Since the spherical harmonics are orthogonal to each other, these relations generates the system of equations ⎧ ˆ √ j (GRMT )Y ∗ (G) ⎪ f (RMT ) + bG f (RMT ) = 4πi aG ⎪ m m m Ω ⎨ (2.30) ⎪ ∂ j (Gr) ⎪ 4πi ∗ ˆ ⎩ aG ∂ f (r) + bG ∂ f (r) = √Ω ∂ r Ym (G) , m ∂ r m ∂ r r=RM T. r=RMT. G from which the matching coefficients aG m and bm may be determined. This generates smooth wavefunctions but the radial functions inside the MT. 19.

(154) spheres loose some freedom which causes them to differ from the physical solutions. In general a larger basis set is required to represent the solution in the LAPW method. In return we get a basis set that does not depend on energy.. 2.7. Semicore states. All states confined to the spheres are put in the core. The LAPW basis set can handle all states of the same -character, if the energy gap between the bands is not too large. In some systems there are states of the same character, but different n quantum number character that has to be treated as valence states. The states in the lower band are called the semicore states, the upper band contains valence states. The energy difference between the semicore states and the valence states is often too large to treat them with the same LAPW basis set. One solution to this problem is to diagonalise two Hamiltonians, one made with a basis set constructed with the lower bands energy as linearisation energy and another with the upper band as linearisation energy. This is of course an unattractive solution to the problem since it forces us to make twice the computational effort. If we could add another basis function of the same to the basis set we could get away with half the work.. 2.7.1. Local orbitals. A possible choice of new basis functions are the local orbitals suggested by Singh [13]. That is localised basis functions in each MT sphere that resemble the semicore states 0 r ∈ interstitial lo χm = am f (r, E1 ) + bm f (r, E1 ) + cm f (r, E2 ) Ym (ˆr) r ∈ MT (2.31) The second energy E2 is chosen to match the band energy of the semicore state. This LAPW+lo basis set will be able to produce correct solutions for both the semicore states and the valence states. Two of the matching coefficients are used to satisfy the condition that the function and its first derivative is zero at the MT sphere boundary. The third coefficient can be chosen to unity or determined by matching each local orbital to a set of fictive planewaves. The latter method is often preferred since it produces local orbitals with the same inversion symmetry as the planewaves. The local orbitals appear as an extended basis set beside the augmented planewaves and matrix elements of the Hamiltonian and overlap matrices has to be calculated for the local orbitals as well as for the other basis functions. This makes the size of the secular matrix larger and the computation slower, but much faster than solving a secular problem for each energy window.. 20.

(155) 2.8. APW+lo. With a sufficient number of planewaves, augmentation functions and local orbitals the LAPW+lo basis will produce correct solutions. Due to the constraint on the first radial derivative of the augmentation functions, the sufficient number of planewaves is usually quite large compared to the number of basis functions needed in an APW calculation. The APW basis functions describe the solutions better, but demand the right energy parameter. In 2000 Sjöstedt et. al. [17] showed that it is possible to combine the better physical description of the valence states of the APW basis functions with the relative independence of the energy parameter of the LAPW basis functions. This is done by adding another set of local orbitals, 0 r ∈ interstitial lo = . (2.32) φm am f (r, E ) + bm f (r, E ) Ym (ˆr) r ∈ MT One of the coefficients am and bm are chosen to satisfy the condition lo (R T ) = 0, the other can either be put to unity or matched to a set of φm M fictive planewaves as in the previous section. These local orbitals resemble the LAPW radial functions in the sense that they have a component of the energy derivative of the radial functions. This will decrease the sensibility of lo local orbitals the basis functions to the choice of energy parameters. The χm of the previous section are still needed to treat the semicore states. In total 2 there are three different basis functions inside the spheres, (apw max + 1) APW apw 2 1 basis functions, φm (r) = aapw r), (lo to max + 1) local orbitals m f (r, E )Ym (ˆ . lo1 lo1 1 reduce the energy dependence, φm = alo r), m f (r, E ) + bm f (r, E ) Ym (ˆ and Ns.c. = ∑m δm(m)sc local orbitals to treat the semicore states, lo2 lo2 lo2 = sc ) Y (ˆ 2 alo f (r, E ) + b f (r, E ) + c f (r, E The energy χm m r). m m m parameters E may be the same for all radial functions with the same -character, but there is always a possibility to choose different energy parameters to avoid linear dependence of the basis functions. The energy parameters can be adjusted to the center of the band in each iteration as in the APW method, but the error in energy is only of the order (E − E )4 and in general it is better to keep the basis set unchanged to speed up convergence. Some caution should be used when choosing the energy parameter for the semicore local orbitals though, since the semicore states form very narrow bands and the description of the semicore states will be poor if the energy parameter fall outside the band.. 2.9. Augmented distorted planewaves. All improvements of the APW basis set discussed in the previous section has concerned the augmentation part of the basis functions. The solution inside the MT spheres often determine the important properties of materials and in 21.

(156) many cases the interstitial region is not important at all. With the tight binding approximation, where the electrons of one atom only interacts with electrons of the nearest neighbours, it is possible to explain many properties of several materials. Sometimes the approximation fails due to important interactions in the interstitial region. Then we need a basis set with high resolution in the interstitial region, such as the APW basis. Since the planewaves are a complete basis set, the basis have any resolution you want. It is just a question of using a large enough number of basis functions. The number of basis functions needed can become quite large. It would be appreciated to find a smaller basis set that would be able to describe the charge density in the interstitial region with the same precision as the corresponding set of planewaves.. 2.10. Adaptive coordinates. In 1993 Gygi [14] proposed another way of changing the planewave basis set. By using planewaves on a distorted space mesh r , the planewaves can be made to be more similar to atomic centered basis functions in the vicinity of atoms, 1 ∂ r 1/2 iG·r e , (2.33) ηG (r) = √ Ω ∂r where Ω is the unit cell volume and ∂∂rr is the determinant of the Jacobian matrix of the coordinate transform. Let each atom distort space around it through a distorting function F and element specific parameters, A and R0 ,

(157) r − ri (r − ri ) F i . r = r + ∑ Ai (2.34) Ri R i. 0. 0. A uniform grid of space points in the undistorted flat space will by this distortion transform to a grid with a higher density of space points in regions where the derivative of the distortion is large. A constant distortion will just shift the grid. The radial density of space points around a single atom, located at the origin will be.

(158)

(159) A r r r dr = 1+ F + F . (2.35) dr R0 R0 R0 R0 The largest density of space points will be found where the second radial derivative is zero.

(160)

(161) d 2 r A r r r = 2 2F + F = 0, (2.36) dr2 R0 R0 R0 R0. 22.

(162) Figure 2.3: 2F + F (solid) and 3F + F (dashed) for the distorting functions suggested by Gygi. Panel I refers to equation 2.38 and panel II refer to equation 2.39. In both cases the maximum density of space point is found at the origin only. R0 = 1. and the third derivative is negative.

(163)

(164) A d 3 r r r r = 3 3F + F < 0. dr3 R0 R0 R0 R0. (2.37). Note that in order to find a maxima at the atomic position, equations 2.36 and 2.37 reduces to the conditions F (0) = 0 and F (0) < 0.. 2.11. Distorting function. In his implementation, Gygi [15] did not use MT spheres and had to make sure his basis functions had many nodes in the vicinity of the atoms. He uses the distorting functions F(x) = sech(x) , (2.38) and. 2 1 F(x) = tanh(x)e−x , x. (2.39). where x = Rr0 is used for simplicity. Both distorting functions fulfil the criteria F (0) = 0 and F (0) < 0, and give maximum density of space points at the atomic positions. There are no other maxima in either of the two functions as shown if Fig. 2.3. When the MT spheres are introduced, the planewaves does 23.

(165) Figure 2.4: 2F + F (solid) and 3F + F (dashed) for the distorting function from equation 2.4. The maximum density of spacepoints is found at r ≈ 0.557RMT if R0 = 1.. not exist inside the spheres and there is no reason to gather a high density of space points inside the spheres. Instead the maximum density should be around the MT sphere boundary. We suggest the distorting function F(x) = sech(x)tanh(x) .. (2.40). The criteria for space point density maximum is fulfilled at a point r = 0, as shown in Fig. 2.4. The higher density of space points in the regions of interest has to be compensated for by a lower density of space points in other regions. Other functions may be considered as good distorting functions, but here we will investigate the following two, F(x) = sech(x) F(x) = sech(x)tanh(x) .. (2.41) (2.42). Outside the MT spheres they are similar, and we expect them to generate similar basis functions in the interstitial region.. 2.12. Distorted augmented plane waves. To increase the efficiency of the distorted planewaves, MT spheres may be introduced in which spherical basis functions are used. This would generate the distorted augmented planewave (DAPW) basis set. The augmentation parts 24.

(166) of the basis functions are constructed just as before and consists of spherical harmonics times radial functions. As before the augmentation wavefunctions are matched to the interstitial wavefunctions. Here the matching coefficients will not be sums over Bessel functions as in equations 2.24 and 2.30, but over the spherical harmonics expansion of the distorted planewaves at the sphere boundary, RG m (RMT ) =. ∗ ηG (r)|r=RMT Ym (ˆr)dSMT .. (2.43). This is the first complication of the distorted augmented planewave method. Other complications are the surface contribution to the kinetic energy ?? and the solutions to Poissons’s equation ?? in the interstitial region. The first complication may be solved by either distorting the MT spheres, so that the terms in the expansion are Bessel functions again, or by calculating the spherical harminics expansion in equation 2.43. The spherical harmonics transform of a general function on a sphere is numerically expensive, but this approach may still be better than to work with distorted spheres, since the spherical potential inside the MT spheres would not be spherical anymore. On the other hand the non-spherical part of the potential inside the MT spheres has the symmetry of the lattice, so by doing all calculations in the distorted space the basis functions will have the symmetry of the crystal. By working in distorted space the second problem (the surface contribution to the kinetic energy) will be solved in the usual way and likewise the solutions to Poisson’s equation. In the last case it is useful to approximate the distorted plane waves by their Fourier transform on the large Fourier mesh used for the charge density and potential, ηG (r) ≈. G max. ∑. . cGG eiG ·r .. (2.44). G. The coefficients cGG will form a N by M matrix U, where N is the number of G vectors in the Fourier mesh, and M is the number of distorted plane waves. This matrix will have the property that it can transform a Hamiltonian for N (L)APW basis functions to a Hamiltonian for M distorted (L)APW basis functions, (2.45) HDAPW = U † HAPW U . Further the eigenstates of this Hamiltonian may be transformed back into the planewave basis by a matrix multiplication, Ψ pw = U † Ψd pw . The use of the transformation matrix U is a fast way of implementing the distorted augmented planewaves, although setting up the Hamiltonian will take as long time as before and the method will not be as efficient as it could be.. 25.

(167) Figure 2.5: Density of nodes for the atom centred function (equation 2.38) and the MT sphere boundary centred function (equation 2.40). 2.13. Results. When using the DAPW basis set, there are two parameters for the interstitial basis functions, R0 and A. First we will use R0 to make it easier to compare the two distorting functions from equation 2.38 and 2.40. If we look at the density of space point generated by the two distortions in Fig. 2.5, we see that for a certain choice of RMT , the distance where the density cuts the dotted line (that is where we have the planewaves back), can be made equal for the two distorting functions. Let R0 = c( f )d0 , c ≈ 0.8336 for function 2.38 and 0 c ≈ 0.6835 for function 2.40 and the dimensionless quantity RdMT will be the number of RMT from the atom the density of space points equal one (as in flat space). Closer to the atom the density will be higher, further away it will be lower than in flat space. The distorting strength A has to be optimised by testing different distortions and comparing the convergence. We start by looking at atoms, that is a calculation of a large simple cubic structure (lattice constant a = 20a0 ) with only one atom. We compare the self-consistent total energy calculated with different basis set in the interstitial region in Fig. 2.6. The curves represent basis sets constructed using different values of d0 . On the x-axis we have the reduced distorting strength AR = A/d0 , a more suitable parameter than A by itself since it gives a better measure of the strength of the distortion around d0 . From Fig. 2.6 we see that we get the same result using a considerably smaller basis set when we use the DAPW basis sets as if we 26.

(168) Figure 2.6: Atomic calculations of total energy versus reduced distorting strength AR , for a reduced DAPW+lo basis set with |G + k| < 2.5. E0 is the total energy of the reference state with full APW+lo basis set with |G + k| < 3. The panels with roman numeral I refer to distortion of type Eq. (2.38) while II refer to Eq. (2.40). The 0 : 0.5 (◦), calculations have been performed for several values of the parameter. RdMT 1.0 (), 1.5 (

(169) ), 2.0 (×), and 3.0 ( ).. 27.

(170) Figure 2.7: Bulk calculations of total energy versus reduced distorting strength AR , for a reduced DAPW+lo basis set with |G + k| < 3. E0 is the total energy of the reference state with full APW+lo basis set with |G + k| < 3.5. The panels with roman numeral I refer to distortion of type Eq. (2.38) while II refer to Eq. (2.40). The calculations 0 : 0.8 (◦), 1.0 (), 1.2 have been performed for several values of the parameter. RdMT (

(171) ), 1.5 (×), and 2.0 ( ).. 28.

(172) Figure 2.8: Calculations for the compound LiNbO3 with the same setup as in Fig. 2.7. 0 with values 0.8 (◦), 1.0 (), 1.2 (

(173) ), and 1.5 (×) are Calculations with parameter RdMT displayed.. 29.

(174) would use a quite large APW basis set. With this large lattice constant and small R( MT ) |G + k| < 3 corresponds to 2800 planewaves, compared to the 1600 planewaves corresponding to |G + k| < 2.5. It is commonly estimated that the computer time needed to solve the secular problem scales as the number of basis functions cubed. That would mean that the diagonalisation of the smaller matrix is up to five times faster than the diagonalisation of the larger matrix. The range of parameters that generates basis sets with optimal convergence are 1 d0 2 and 34 AR for both distorting functions and all atoms. In terms of space point density this means that the density should be higher than for flat space in a region not further from the atom than twice the MT radius. Beyond this point the density of space points is lower than in flat space. There is no detectable difference between the atom centered distorting function and the MT sphere centered distorting function. This is expected since the part of the distorted plane wave inside the MT sphere is not used. In sparse structures such as single atoms we have a large interstitial region and the efficiency of the basis set depend heavily on the possibilities for the basis functions to combine in such a way that the wavefunctions are zero in a large part of space, but at the same time make a good representation of the physical wavefunction close to the atom. That is why we need so many planewaves in flat space. The DAPW basis set is expected to perform best for open structures. In Fig. 2.7 calculations of the same atomic species as in Fig. 2.6 in their solid state configuration (with 20% tetragonal distortion to enlarge the interstitial region) are shown. Here we see that the DAPW basis set is not as easily optimised as before - in general we do not hit the lower dashed line anymore. On the other hand almost the same parameters still holds. The best DAPW basis functions are generated using 1 d0 1.5 and 12 AR 34 for both distorting functions and all atoms. In the solid state configuration the distortion is not allowed to reach as far from the atoms as in the case of free atoms. This is of course due to the fact that the space point density in the vicinity of an atom is disturbed by the distortion from the neighbouring atoms. For the same reason the distorting strength must be kept slightly lower than for the free atoms. In Fig. 2.8 we show calculations of a compound, LiNbO3 in its ground state solid structure. We used the same distorting function and parameters A and d0 for all three species in each calculation, although any combination may be used. The optimal parameters are about the same as for the previous solid calculations, except that the range may be allowed to increase slightly. The perovskite structure of LiNbO3 is rather open and the fraction of interstitial region is larger than for the denser fcc and bcc structures. This is likely to be the reason for the slightly improved results compared to the single species solid calculations.. 30.

(175) 2.14. Conclusions and outlook. In this work a new set of basis functions for electronic structure calculations is presented. The distorted augmented planewaves, composed of the distorted planewaves in the interstitial region and the augmentation part of the well known APW basis set inside the MT spheres. The construction of the new basis set by distorting space around atoms and define the planewaves on the distorted mesh is discussed and the difficulties of implementing the method are mentioned. The simplified implementation by transformation matrices is explained and two distorting functions are selected to demonstrate the expected increase of efficiency of the new basis set. Test calculations show that for sparse systems even the simplified implementation of the new basis set gives optimal basis functions using the same parameters for the distortion regardless of atomic species. For solid structures the efficiency of the new basis set decrease for close packed structures, but universal distortion parameters can be chosen that generates the best possible basis set regardless of structure and chemical composition. There are many possible distorting functions that have not yet been tried, and several should be investigated to explore the possibility of finding even better basis sets. A full implementation of the method may also show increased efficiency, although there are a few difficulties left to take care of.. 31.

(176)

(177) 3. Magnetism. According to Hund’s first rule, almost all atoms are magnetic. Any open shell with n electrons will have a spin magnetic moment mspin with magnitude n ≤ 2 + 1 n mspin = (3.1) 2 + 1 − n n ≥ 2 + 1 Already for diatomic molecules the spin moment usually vanish, since the energy levels of a molecule are spin degenerate molecular orbitals that to a good approximation can be written as a linear combination of atomic orbitals (LCAO). Each molecular level will be populated by two electrons, one of each spin. If the number of electrons is even the molecule will have equally many spin up and spin down electrons. When more atoms are added to the system and we go towards a crystalline solid the energy levels broaden and form bands. If a band crosses the Fermi level the total energy of the system may be lowered by breaking some symmetry related to this band and split the band. Any symmetry that splits the band will do if the energy gain by the splitting is larger than the energy loss in kinetic energy. It could be a broken crystal symmetry as in the case of Jahn-Teller like distortions, where the t2g and eg states for d-electrons may split in case of odd d-occupation. Or it could be a broken time reversal symmetry, allowing the spin up and spin down bands to split. This is the Stoner model for ferromagnetism ??, which states that if the exchange energy gained by spin splitting is larger than the kinetic energy loss (by increasing the Fermi level), the band will split and form a ferromagnetic metal, see Fig 3.1. 3.0.1. Magnetic order. Another picture of magnetism is the (classical) Heisenberg model, where each atomic site is home to a localised atomic moment with direction ei . The magnetic sites interact through the Heisenberg Hamiltonian H = − ∑ Ji j ei · e j .. (3.2). i= j. If the Heisenberg parameter Ji j > 0 the interaction is ferromagnetic, if Ji j < 0 the interaction is antiferromagnetic. The minimum energy arrangement of the magnetic moments gives us the ground state magnetic configuration. In case of negative Ji j ’s this configuration might be complicated, for example in the case of a triangular lattice the ground state configuration is non-collinear, see 33.

(178) Figure 3.1: Stoner model for ferromagnetism. The total energy of the system is lowered by the spin splitting if the gain in exchange is larger than the loss in kinetic energy.. Fig. 3.2 The x- and y-components will show up as spin mixing terms in a two component spinor picture, and the Hamiltonian will have the form ↑↑ H ↑↓ H = Bxc · σ , (3.3) H mag = H ↓↑ H ↓↓ where σ is a vector containing the three Pauli spin matrices 0 1 0 i 1 0 σx = , σy = , σz = 1 0 −i 0 0 −1. (3.4). Another special non-collinear arrangement of localised magnetic moments is the spin spiral structure. Here magnetic moments rotate in a symmetric way in some direction of the crystal, see Fig. 3.3. In case the spin orbit coupling is sufficiently weak to be neglected, the symmetry of the spin spiral can be used by generalising Bloch’s theorem and combine a spatial translation by vector T with a rotation in spin space by an angle φ around the spiral axis q . ?? (3.5) T˜T,φ = TT Rφ . 34.

(179) Figure 3.2: Heisenberg model for moments on a triangular lattice. The antiferromagnetic interaction leads to a non-collinear arrangement of the magnetic moments.. Since spin space is decoupled from real space we may consider the spiral axis q to define the z-direction in spin space, regardless of the direction in real space. Then the magnetisation density m(r) and field b(r) can be replaced by the translationally invariant quantities u(r) = e−iq·r [mx (r) + imy (r)] h(r) = e. −iq·r. [bx (r) + iby (r)]. The magnetic part of the Hamiltonian may then be written as Hmag = eiq·r h(r)σ− + h.c. + bz (r)σ z ,. (3.6) (3.7). (3.8). where σ− = 12 (σx − iσy ). The eigenstates of this Hamiltonian will be generalised Bloch spinors of the shape αn,k (r)eik·r (3.9) ψn,k (r, q) = βn,k (r)ei(k+q)·r The magnetic interactions are fairly weak compared to the terms in the Hartree Hamiltonian. It is possible to treat the magnetic interactions as a small correction to the eigenstates of the non magnetic Hamiltonian by adding them in a second variational stage using the eigenstates of the spin independent 35.

(180) Figure 3.3: Spin spiral. The direction of the magnetic moment varies as mˆ = eiR·q as we move through the crystal.. Hamiltonian as basis functions for the magnetic calculation, in a similar way as the spin orbit interaction is commonly implemented. In case the spiral qvector is commensurate with the k-point mesh used for the first variational step, we simply use the eigenstates at k as spin up and k + q as spin down when we set up the second variational Hamiltonian. We just have to remember that the magnetisation density calculated is not the regular m(r), but rather the translationally invariant u(r). If the spiral vector q is incommensurate to the k-point set all points k + q must be added in the first variational step.. 3.1. Intra-atomic non collinear magnetism. In the Heisenberg model the magnetic moments are vectors sitting on points in a lattice, This picture is known as the atomic moment approximation (AMA). Many magnetic ions really have a smooth magnetisation density that can be modelled by AMA. But some systems are different. The magnetisation density close to the atom can be rapidly varying, not often in magnitude but in direction. In Fig. 3.4 the magnetisation density of US is shown. US is a ferromagnetic compound discussed in papers II and IV, with an atomic moment in the [111]-direction of the crystal. The magnetisation density is rotating on the surface of a sphere around the U atom. The strong spin orbit coupling ties the spin density to certain orbital quantum numbers. The spin density in each direction will have the shape of the selected f-orbitals. In total each 36.

(181) Figure 3.4: The magnetisation density in US. The x-,y- and z-components are shown in different colors. The lobes represent isosurfaces of the magnetisation density.. component form a triakontadipole (that is a ’32-pole’). The shape of a single triakontadipole is shown in Fig. 3.5. The polarisation of such multi poles of quite high order will be discussed further in the chapter Multi pole LDA+U.. 37.

(182) Figure 3.5: The magnetisation density in the z-direction minus the average in Californium. The dominating part of the magnetisation density is uniform since Californium has a large magnetic moment, but the magnetisation density is not uniform.. 38.

(183) 4. Multi pole LDA+U. The commonly used local density approximation (LDA), although being the simplest possible approximation to the exchange correlation energy, works surprisingly well for many systems [25]. Since LDA automatically fulfils the sum rule of the exchange correlation hole and that both the exchange energy and the correlation energy are underestimated and of opposite sign and thus partly cancel each other, it is still the most commonly used functional for electronic structure calculations. LDA can by construction only take local effects into account. In some systems the non-local effects are strong, such as in transition metal oxides or systems containing Lathanide or Actinide atoms. To include the non-local effects it is common to add an orbital dependent operator to the Hamiltonian that acts on a few selected states only (usually atomic orbitals with a certain angular momentum quantum number ). The energy of this Hartree-Fock like contribution is given by EHF =. 1 ∑ [ρac ρbd − ρad ρbc ] ab |g| cd , 2 abcd. (4.1). where ρab is the density matrix acting as an occupation matrix, a, b, c and d runs through all the selected states and g is the interaction. If the interaction is the bare Coulomb interaction the energy becomes the Hatree-Fock energy for the selected states. Hartree-Fock is by definition correlation free and known to overestimate the exchange energy. By screening the Coulomb potential a more physical Yukawa potential can be used [28, 29] V (r) =. e−λ r . r. (4.2). The screening parameter λ may be interpreted as the Thomas-Fermi wave vector of an electron in the selected state. The screening may be set to a reasonable value by hand, used as a free parameter or calculated by constrained DFT methods.. 39.

(184) 4.1. Multi pole momentum tensor formulation. The interaction term ab |g| cd of Eqn. 4.1 may be expanded in spherical harmonics g(r12 ) =. ∞. ∑ gk (r1 , r2 ). k=0. k 4π ∑ Ykq∗ (ˆr1 )Ykq (ˆr2 ) . 2k + 1 q=−k. (4.3). Now the interaction term consists of a spin term S(sa , sb , sc , sd ) = δsa ,sc δsb ,sd , radial terms F. (k). =. . (4.4). gk (r1 , r2 )R2 (r1 )R2 (r2 )r12 r22 dr1 dr 2 ,. (4.5). and angular terms Ω(k) =. ×. k. 4π ∑ 2k + 1 q=−k . . ∗ Ym (ˆr )Y ∗ (ˆr )Y (ˆr )d rˆ1 a 1 kq 1 mc 1. (4.6) ∗ Ym (ˆr )Y (ˆr )Y (ˆr )d rˆ2 . b 2 kq 2 md 2. The angular integrals can be written as Wigner 3j symbols [20]

(185) 1/2 4π Ym1 (ˆr)Ykq (ˆr)Ym2 (ˆr)d rˆ = 2k + 1 = (2 + 1). k 0 0 0. . k m1 q m2. . (4.7). Using Eqn. 4.4 to 4.7 the interaction term may be expressed as 2 2 k ab |g| cd = δsa ,sc δsb ,sd []2 ∑ F (k) 0 0 0 k=0. ×. k. ∑ (−1)ma +mb +q. q=−k. . k −ma −q mc. . k −mb q md. (4.8) .. Instead of using the density matrix as occupation matrix we may choose to use tensor occupation numbers w defined as the expectation value of a multi. 40.

(186) pole momentum tensor operator [21, 22, 23, 24] wkx = Tr ωxk σ0 ρ . (4.9) . . ωxk ≡ mb |ωxk |ma = (−1)−mb. k −mb x ma. . ω k ,. . where σ0 is the unit matrix in spin variables. We will call ω k = n−1 k .. (4.10). (4.11). There are several choices for nk , and each choice gives a different type of irreducible spherical tensor. With the choice nk = . (2)! . (2 − k)!(2 + k + 1)!. (4.12). wk corresponds to a kind of multi-pole momentum tensor.[26, 27] While the choice nk = 1 .. (4.13). will give Racah’s unit tensor. The formalism presented here is valid for any choice of irreducible tensor. It is worth noting here that the spherical tensor, C k might seem to be a good choice, but it does not work since then n−1 k would vanish for odd k. By rewriting σ0 as s 0 s , (4.14) σ0 = δs1 ,s2 = (−1)s−s1 −s2 0 sa we may construct a double tensor k p k p wkp xy = wx wy = Tr ωx χy ρ. χyp. = (−). s−sb. s −sb. p s y sa. (4.15) . sχ k s ,. . where it is convenient to choose −1/2 . sχ k s = n−1 sp = (p + 2)!. (4.16). (4.17). 41.

(187) Then wk = wk0 . The transform from density matrix to multi pole momentum occupation numbers can be inverted k mc − (4.18) ρac = ∑(2k + 1)nlk (−) −mc x ma kx × ∑(2p + 1)nsp (−)sc −s py. s −sc. p s y sa. wkp xy .. (4.19). The transform is thus 1-to-1, which can also be seen through the fact that the density matrix is a Hermitian 2(2 + 1) × 2(2 + 1) matrix and thus have (2 + 1)(4 + 3) real and (2 + 1)(4 + 1) imaginary independent components. 2s 2 The multi pole momentum tensors has ∑k=0 (2k + 1) × ∑ p=0 (2p + 1) = (4 + 1 + 1)(2 + 1)/2 × (4s + 1 + 1)(2s + 1)/2 = (4 + 2)2 independent components, with s = 1/2. The multi pole momentum tensors contains the same information as the density matrix, but has a more intuitive meaning. With the introduction of the double tensors wkp , both the direct and the exchange term can be put in simple forms EH =. . ∑. F (k) I(, k, k) wk0 · wk0. (4.20). t=0 k=2t. EX = −. . ∑ t=0. F (k). 2. ∑. k1 =0. 1. J(, k, k1 ) ∑ wk1 p · wk1 p ,. (4.21). p=0. k=2t. with (2 + 1)2 2 nk1 I(, k, k1 ) = 2. 42. . k 0 0 0. 2 (4.22).

(188) and (2 + 1)2 (2k1 + 1) (−)k1 n2k1 J(, k, k1 ) = 4. . k 0 0 0. 2 . k1 k. . (4.23) 2k1 + 1 (−)k1 = I(, k, k1 ) 2. . k1 k. = I(, k, k1 )B(, k, k1 ) . (4.24). Then the full HF energy expression can be written EHF =. . ∑. F (k). t=0. 2. 1. k1 =0. p=0. ∑ I(, k, k1 ) ∑ [δ (k, k1 )δ (p, 0) − B(, k, k1 )] wk1 p · wk1 p .. k=2t. (4.25) When the spin-orbit coupling is strong (which is the case for many of the strongly correlated systems we are addressing with this method) it is useful to couple the two indices to a third. This brings us a three index tensor k r p 1 k1 pr (−1)k1 −x+p−y Nk−1 wt = ∑ wkxy1 p 1 pr −x t −y xy =∑. k1 r p −x t −y. ab. (−1)k1 −x+p−y (4.26). Nk1 pr nsp nk1. xy. × ∑(−)−mb. . k1 −mb x ma. . (−1)s−sb. s −sb. p s y sa. ρab. Here the symbol Nabc is a generalisation of the bare Wigner-3j symbols a b c , but which is non-zero also for g = a + b + c odd, 0 0 0 . (g − 2a)!(g − 2b)!(g − 2c)! 1/2 Nabc = ig (g + 1)! (4.27) g!! × . (g − 2a)!!(g − 2b)!!(g − 2c)!! 43.

(189) It reduces to the ordinary bare 3j-symbol for even g. The Hartree and exchange contributions to the Hartree-Fock energy written in terms of the three index tensors is then EH =. 4. ∑ F (k) I(, k, k)wk0k · wk0k. 2k=0. (4.28) EX = −. ×. 1. 4. 2. 2k=0. k1 =0. ∑ F (k) ∑ J(, k, k1 ) k1 +p. ∑ ∑. (2r + 1)(−1)k1 +p+r Nk21 pr wk1 pr · wk1 pr ,. (4.29). p=0 r=|k1 −p|. since 2 =1 (2k + 1)Nk0k. . ∀k .. k s sb −s ρab = ∑[k]nlk (−1) [p]nsp (−1) ∑ −mb x ma −sc py kx k r p × ∑[r]Nkpr (−1)x+y−k−p wtkpr . −x t −y rt mb −. 4.2. (4.30). p s y sa (4.31). Interpretation of the tensors. So far it may seem as the expressions only gets more complicated when working with multi pole momentum tensors rather than the density matrix. It may be true that the mathematical expressions are longer, but the tensor occupation numbers tells us a lot about the system. Tensors with even k represent multi poles of the charge (with p = 0) or magnetisation (with p = 1) densities. Tensors with odd k has odd time reversal symmetry and thus represent charge (p = 0) and spin (p = 1) currents. Note that it is only the angular part of the mentioned quantities that is considered, since the radial dependence has been integrated out and is contained in the Slater parameters F ( j) . The larger k the larger is the order of the multi pole. w000 is proportional to the charge in the selected -shell, w011 is proportional to the total spin moment of the -shell and w101 is proportional to the total orbital moment of the -shell. w110 is proportional to L · S, an orbital spin current associated with the spin orbit coupling. w111 is proportional to L × S and will not show up that often. The five components of w112 contain other linear combinations of Li S j not included in the scalar and vector products. In total the w11x x=[0,1,2] tensors has nine components which will all be linearly independent combinations of the 44. .

(190) components in the product ⎞ ⎞ ⎛ ⎛ Lx # Lx Sx Lx Sy Lx Sz $ ⎟ ⎟ ⎜ ⎜ LT S = ⎝ Ly ⎠ Sx Sy Sz = ⎝ Ly Sx Ly Sy Ly Sz ⎠ Lz Lz Sx Lz Sy Lz Sz. (4.32). In a similar way w202 contains five charge quandrupoles that looks like 3dorbitals and w21 the 15 linear combinations of spin quadrupoles, and so on.. 4.2.1. Saturation. We may see the tensor occupation numbers as generators of density matrices. It follows from Eqn.4.31 that by occupying the tensor component wtkpr , we generate a density matrix which should be added to the density matrix of the system before the polarisation. An unpolarised system (here a less than half filled s-shell) has a density matrix ⎞ ⎛ ↑ ↓ ns ⎜ ⎟ ρ (0) = ⎝ ↑ 1 0 ⎠ ns ≤ 1 , (4.33) 2 ↓ 0 1 that corresponds to occupying w000 = ns . The only possible polarisation of an s-electron system id spin polarisation. Putting w011 0 = ns generates the density matrix 1 0 n s ρ (1) = , (4.34) 2 0 −1 and the total density matrix is then ns ρ (0) + ρ (1) = 2. . 2 0 0 0. .. The magnetic moment is given by 1 0 2 0 ns m = Tr[σ · ρ] = Tr eˆz = ns eˆz . 2 0 −1 0 0. (4.35). (4.36). Spin polarisation in the x- and y-directions are obtained by occupying w011 1 and . It is necessary to occupy both since to generate a Hermitian (physical) w011 −1 density matrix kpr ∗ t wkpr −t = (−1) (wt ) .. (4.37). It may seem as this reduces the degrees of freedom from (2[])2 to 2[] × 2( + 1) − 1, but the tensor occupation number may have both a real and an 45.

(191) imaginary component (except for wkpr which is always √ real), and then we are 0 √ 2 011 011 back at (2[]) again. If we put w1 = 2 and w−1 = − 2 we get the density matrix 0 1 n s , (4.38) ρ (1) = 2 1 0 and magnetic moment ns m = Tr[σ · ρ] = Tr 2. . 0 1 1 0. . 1 1 1 1. eˆx = ns eˆx .. (4.39). 011 ∗ √ If instead we put w011 = 2i the generated density matrix is −1 = − w1 0 −i n s ρ (1) = , (4.40) 2 i 0 and the magnetic moment is 0 −i 1 −i ns eˆy = ns eˆy . m = Tr[σ · ρ] = Tr 2 i 0 i 1. (4.41). If we diagonalise the density matrix, the eigenvalues should be in the interval [0, 1]. If we start from an unpolarised density matrix n I, (4.42) ρ (0) = 2[] # $∗ and add the density given by setting wtkpr = (−1)t wkpr = α, where α starts −t from zero and grows in magnitude to finite positive or negative numbers. Some eigenvalues of the total density matrix will become larger, others smaller. As soon as one eigenvalue hits zero or one, the density matrix is saturated by this polarisation and we can not polarise it more. With only one electron in the shell it is sufficient to require that the generated total density matrix is positivedefinite. With more electrons we have to check that all eigenvalues are less than one too. To find the number max(|α|), that is the saturation value of wtkpr we start from a density matrix with one electron. With only one electron in the system we know that at saturation n kpr I + αXt = 0 , detρ = det (4.43) 2[]. 46.

No results found