Quantum Chemistry of Solids / 21-Basis Sets and Pseudopotentials in Periodic LCAO Calculations

8.3 Relativistic E ective Core Potentials and Valence Basis Sets

321

a scheme is proposed in which most of the electrons are treated in the standard way (Dirac–Hartree–Fock equations with only Coulomb interactions and calculation of all the other terms with the aid of first-order perturbation theory), while for the core electrons of large-Z atoms or ions, the generalized relativistic HF equations are used. Therefore, one applies the solutions of the generalized relativistic HF equations for the construction of the relativistic Slater determinant in the case of core electrons, while for the rest of the electrons the one-electron functions are obtained from the standard Dirac–Hartree–Fock equations. Unfortunately, the application of this fully relativistic theory (involving four-component spinors built in terms of a large and small component) is practically too di cult not only for solids but even for molecular systems (this application is limited to small-size molecules).

As was noted in the preceding sections the relativistic calculations of molecular systems are made in quasirelativistic (two-component) approximations based on Hamiltonians operating on the large component only. Between di erent alternatives to decouple large and small components for molecules and solids the Douglas–Kroll–Hess (DKS) approach [548, 549] is the most popular as this approach is variationally stable (and therefore can be incorporated into SCF calculations) and allows its accuracy to be improved [550]. In the DKS approach the decoupling of large and small components of the Dirac Hamiltonian is made through a sequence of unitary transformations

Hm = UmUm−1 · · · U1U0HDU0†U1† · · · Um† −1Um†

(8.62)

Each transformation Un is chosen in such a way that the o diagonal blocks of the Dirac Hamiltonian are zero to a given order in the potential. The decoupled, blockdiagonal transformed Hamiltonian

HL	O(V n)
DKHn
Hn = O(V n)	HDKHnS		(8.63)

The large-component block, HDKHL n is the two-component e ective Hamiltonian to be included in the electronic-structure calculation. The order n of transformation can be systematically increased to improve the accuracy of transformation. Two classes of operators are involved in the construction of the DKHn Hamiltonian: a function of the square of the momentum f (P 2), p = −i and those involving the potential V, (σp)V (σp), where σ are the Pauli spin-matrices [550]. The explicit form of HDKHn up to fifth order was derived in [551, 552]. The DKH approach leads to a two-component formalism, which is separable in a spin-averaged (spin-free) and spin-dependent part. In the scalar-relativistic (SR) approximation the spin-dependent part is ignored. The DKH transformation for the one-electron operators in scalarrelativistic HF LCAO calculations with periodic boundary conditions is used in [553] and implemented in CRYSTAL88 code. The fact that scalar-relativistic e ects are short range was made use of so that the electrostatic potential calculated by the Ewald method was not relativistically corrected. The results of SR HF LCAO calculations [553] of silver halides AgX (X=F, Cl, Br) with the fcc NaCl structure are presented in Table 8.6.

The reoptimized uncontracted correlation corrected basis sets were used. The diffuse exponents (< 0.04) in the basis set for Ag were discarded. The results given in Table 8.6 demonstrate the influence of scalar-relativistic e ects in the AgX crystals

322 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

Table 8.6. Scalar-relativistic e ects in AgX (X=F, Cl, Br) crystals, [553]

Property	AgF	AgCl	AgBr	Ag
˚	5.175	5.937	6.207	4.456
Lattice parameter (A)
	(5.191)	(5.990)	(6.260)	(4.546)
Binding energy (kJmol−1)	–371.2	–385.9	–333.4	–53.0
	(–407.5)	(–417.7)	(–362.8)	(–39.4)
Bulk modulus (kbar)	459	306	273	298
	(460)	(297)	(259)	(322)

properties (the results of nonrelativistic HF LCAO calculations are given in brackets, the binding energy is given with reference to the HF ground-state energies of the free atoms). It is seen that lattice parameters show small relativistic e ects – they decrease (up to 2% for Ag) as compared to the nonrelativistic values. The binding energies demonstrate a remarkable bond destabilization (up to 8.9% for AgF and 34.5% for Ag as compared to the nonrelativistic values). Bulk moduli change up to 7.5% (for Ag). Though scalar-relativistic e ects of silver compounds are expected to be relatively small (compared with the compounds of heavier elements), remarkable changes of the HF binding energy and of the bulk modulus are observed. We note that in the relativistic HF calculations the correlation e ects are not included explicitly in the Hamiltonian (the correlation correction is included only in AO basis sets). In Table 8.7 properties of AgCl crystal are compared at various levels of theory. Table 8.7 demonstrates good agreement for the lattice parameter and bulk modulus between HF calculations with RECP (Hay–Wadt RECP is used) and the scalar (one-component) relativistic HF DKH approach. Slight di erences are seen for the binding energy.

Table 8.7. Crystal properties of fcc AgCl at various levels of theory [553]

Method	Total energy	Lattice	Binding energy	Bulk modulus
	(a.u.)	const. (A)˚	(kJ mol−1)	(kbar)
Nonrel. HF	–5656.601	5.990	–417.7	297
Rel. HF	–5771.344	5.937	–385.9	306
HF RECP	–159.709	5.943	–376.5	307
DFT ECP	–160.870	5.610	–514.6	560
Exper.	–	5.510	–519.8	535

We noted that it is necessary to simultaneously include the relativistic and correlation e ects in the Hamiltonian. Such inclusion for periodic systems is usually made in the relativistic DFT approaches.

The fully relativistic LCAO method for solids, based on the DKS scheme in the LDA approximation was represented in [541]. The basis set consists of the numericaltype orbitals constructed by solving the DKS equations for atoms. This choice of basis set allows the spurious mixing of negative-energy states known as variational collapse to be overcome. Furthermore, the basis functions transform smoothly to those in the nonrelativistic limit if one increases the speed of light gradually in a hypothetical way.

8.3 Relativistic E ective Core Potentials and Valence Basis Sets

323

This is also an important feature to avoid the variational collapse. The proposed approach was applied to Au and InSb crystals. The basis functions were chosen so that they have enough variational flexibility: not only were used AOs of neutral atoms but also those of positive ions Au+, Au2+, Au3+, In2+, Sb2+. In the comparison for Au the results of scalarand two-component relativistic calculations (lattice constant and bulk modulus) indicates that the spin-orbit coupling plays a minor role in the structural properties of Au. The comparison of results of the relativistic calculations with those of the nonrelativistic calculations shows that the lattice constant is overestimated by 5% and the bulk modulus is underestimated by 35% in nonrelativistic calculations. This strongly shows the importance of the inclusion of the relativistic e ects in the study of the structural properties of Au. In contrast with the case of Au, the results of nonrelativistic calculations of InSb are not so poor; the error in the lattice constant is 1% and the error in the bulk modulus is 10%. This should be due to the fact that In and Sb are not so heavy that the relativistic e ects do not play an important role in studying the structural properties of InSb.

The implementation of scalar relativity in an all-electron LDA linear combinations of the Gaussian-type orbitals (LC-GTO) method for solids was reported by Boettger [542] and applied for crystalline Au with a fcc lattice. This initial implementation was later extended to include spin-orbit coupling terms, produced in the second-order DKH transformation. The GTO Au atomic basis set of 19s14p10d5f primitive GTOs was derived from scalar-relativistic (SR) and nonrelativistic (NR) calculations of paramagnetic atom. These calculations show that bulk and one-electron properties obtained by the LCGTO method, are very close to those obtained with allelectron, scalar-relativistic DFT techniques using other basis sets (LAPW or LMTO SR methods). The e ciency of the LCGTO approximation in relativistic calculations of solids was demonstrated also in the calculations of equilibrium volumes and bulk moduli for the light actinides Th through Pu [555]. It is concluded that two independent (both scalar and quasirelativistic) electronic-structure DFT methods with the di erent choice of basis set (LCGTO and LAPW) are in good agreement with each other. Table 8.8 shows this agreement for fcc Th and Pu metals giving numerical values of the equilibrium volume V (a.u) and the bulk modulus B (GPa), for comparison also results of FPLMTO calculations are given. The agreement between LCGTO and FPLAPW results is important, keeping in mind the primary advantage of LCAO methods over other existing DFT approaches (FLAPW, FLMTO)- the possibility of LCAO methods to treat both molecular and periodic systems and in this way to bridge the gap between quantum chemistry and solid-state physics.

Comparison of the nonrelativistic and scalar-relativistic results for fcc Au reveals the large impact that relativity has on the lattice constant (6%) and bulk modulus (57%) [542]. The most important qualitative change in the band structure of fcc Au is the more than 2-eV lowering of the s-band relative to the bottom of d- bands. In addition, the overall width of the d-bands is increased by more than 15% due to a relativistic delocalization of the d- states. The spin-orbit coupling included LCGTO DFT-GGA calculations were made for fluorite structure actinide oxides MO2 (M=Th,U,Pu) and their clean and hydroxylated surfaces, [556], magnetic ordering in fcc Pu [557] and bulk properties of fcc Pb [558].

In the calculations using nonrelativistic LCAO codes for solids the relativistic effects are included implicitly in RECPs and basis sets. In the relativistic calculations

324 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

Table 8.8. Equilibrium volume V (a.u) and bulk modulus B (GPa) for Th and Pu with FP-LMTO, LCGTO and FPLAPW calculations, using DFT LDA and GGA Hamiltonians, with and without spin-orbit (SO) e ects included [555]. Experimental data for fcc Th:

V =221.7 a.u., B=58 GPa

			fcc Th			fcc Pu
		LDA	+SO	GGA	+SO	LDA	+SO	GGA	+SO
FPLMTO	V	-	182.1	-	199.9	109.2	-	119.2	-
	B	82.6	61.5	-	-	214	-	143	-
LCGTO	V	200.5	199.3	216.9	214.8	111.5	125.7	121.2	141.4
	B	64.8	71.7	58.8	63.4	218	124	170	97
FPLAPW	V	199.7	204.1	219.3	218.1	111.9	120.2	122.3	133.4
	B	61.0	78.6	56.7	73.1	194	143	153	121

of periodic systems RECPs of the Stuttgart-K¨oln group (energy-consistent RECPs) are often applied. These RECPs (and the corresponding basis sets) are obtained in the relativistic DHF atomic calculations and are tabulated as the linear combination of GTFs, see www.theochem.uni-stuttgart.de. As was noted in Sect. 8.1, the basis sets generated for molecular calculations have to be adapted for the periodic systems calculations. As an example we mention the valence basis sets for lanthanide 4f-in-core relativistic pseudopotentials, adapted for crystal orbital ab-initio calculations [559]. The crystalline calibration calculations were made for Ln2O3 (Ln=La–Pm) crystals

with hexagonal space group P 3m1 and one formula unit per cell. The calibration within the HF and DFT schemes was made using the CRYSTAL03 code [23]. To verify these basis sets the calculated geometries and cohesive energies were compared with experimental data. The valence basis sets adapted for crystal calculations restore the calculated cohesive energy of Ln2O3 (Ln = La–Nd) to more than 88% of the experimental data within the a-posteriori HF correlation scheme in combination with gradient-corrected functionals. Good agreement has also been found between the conventional DFT results and the experimental cohesive energy with a deviation of only a few per cent.

The implementation of the LCGTO relativistic two-component DKH approximation in a fully self-consistent all-electron DFT approach for molecules (this implementation requires changes to be made in the Hamiltonian [550]), allows the RECP results (when the nonrelativistic calculations are made with the use of relativistic pseudopotentals) to be compared with those obtained in scalar-relativistic and two-component relativistic calculations. Such a comparison was made in [560] for the bond dissociation energy in the UF6 molecule, obtained from the relative energies of the fragments (UF6 = UF5 + F), corrected for zero-point energy and spin-orbit interaction. The small-core RECP PBE0 and B3LYP calculations give the bond-dissociation energy as 68.34 and 69.59 kcal/mole, respectively, all-electron two-component relativistic PBE0 and B3LYP calculations – the values 68.94 and 70.39 kcal/mole, while the experimental values are given between 69.5 and 73 kcal/mole. Note that the large-core RECP result is o by more than 50% so that to obtain a good agreement with experiment one needs to work with a hybrid density functional and small-core RECP.

8.3 Relativistic E ective Core Potentials and Valence Basis Sets

325

This conclusion was taken into account in RECP calculations of crystals. RECP hybrid DFT calculations were made for PtO crystal [421]. For Pt, 60 core electrons [Kr]4d104f 14 were replaced with the small-core relativistic Stuttgart–Dresden e ective core potentials RECPs. A valence basis set was optimized for Pt considering 5s, 5p, 5d, 6s, and 6p orbitals. This basis set circumvents numerical linear dependencies in periodic calculations and it is of similar quality to the Stuttgart basis set for most practical applications. Benchmark calculations in PtO using the 6-31G(d), the 6-311G(d), and 8-411G(d) basis sets for oxygen showed that these three bases yielded similar results. RECPs of the Stuttgart group were also applied in periodic hybrid DFT calculations of PuO2, Pu2O3 [422], PtN [420], UO2 [543].

The molecular calculations [550, 560] allow us to conclude that the DKH approximation in combination with hybrid DFT functionals provides a reliable tool for the prediction of structural and thermochemical properties of molecules.

The all-electron DFT LCAO approach with Gaussian basis sets was extended to scalar-relativistic calculations of periodic systems [561]. The approach is based on a third-order DKH approximation, and similar to the molecular case, requires only a modification of the one-electron Hamiltonian. The e ective core Hamiltonian is obtained by applying the DKH transformation to the nuclear–electron potential VN . Considering that relativistic e ects are dominated by the short-range part of the Coulomb interaction, it is proposed to replace the nuclear–electron Coulomb operator used to build the DKH Hamiltonian by a short-range Coulomb operator

Z
V ζ (r) = − r erf c(ζr)	(8.64)

where the complimentary error function screens the long-range tail of the potential. The use of the complimentary error function simplifies the calculation of the matrix elements between GTOs, but in general, any screening function could be employed. The correct nonrelativistic long-range behavior of the Coulomb operator is retained if the core Hamiltonian is defined as

Hcore = HDKHnζ

− VNζ + VN

(8.65)

where the supraindex stands for quantities involving the short-range potential given by (8.64). The parameter ζ switches from a Hamiltonian where only the free-particle (averaged) kinetic energy is relativistically corrected (ζ → ∞), to the full-range DKHn

Hamiltonian (ζ = 0). In this way, the relativistic term HDKHnζ − VNζ for periodic systems can be evaluated using VNζ and pVNζ p obtained by adding contributions from a finite number of cells in direct space. It should be pointed out that the summation of VNζ and pVNζ p for ζ = 0 diverges. The divergent term can be removed from the summation by adding a uniform charge-canceling background, as is usually done in

Ewald-type summations. For the relativistic term HDKHnζ −VNζ the contribution from the uniform background is exactly zero, since the DKH transformation does not change

a uniform potential. The choice of ζ only a ects the relativistic correction to the core Hamiltonian, while the nonrelativistic electrostatic contributions are entirely done using FMM, see Chap. 7. The tests on bulk systems by changing the value of ζ from 0.05 a.u to 0.5 a.u., showed that the lattice constants, bulk moduli, and bandgaps are rather insensitive to ζ in this range [561], (ζ =0.1 has been chosen in the discussed

326 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

calculations). In this approach, relativistic e ects are included through the nucleus– electron interaction, while the electron–electron, Coulomb, the portion of Hartree– Fock exchange, and DFT exchange-correlation potential remain nonrelativistic. Such a strategy has been successfully employed for molecules and solids earlier.

The benchmark scalar-relativistic calculations [561] were made for the bulk metals (Pd, Ag, Pt and Au) and the large bandgap semiconductors AgF and AgCl. It was shown that scalar-relativistic e ects reduce the lattice constant by 0.06–0.10 ˚A for the 4d metals (Pd and Ag), and by 0.14–0.22 ˚A for the 5d metals (Pt and Au). For the 4d metals, scalar-relativistic e ects increase the calculated bulk moduli by 20–40 GPa, while for the 5d metals this increase is between 60 GPa and 100 GPa. For both AgF and AgCl crystals scalar-relativistic e ects decrease the energy gap – by 1.0 eV for AgF and 0.9 eV for AgCl.

The DKH approach allows relativistic all-electron DFT calculations to be performed using the traditional LDA and GGA approximations, as well as meta-GGA and hybrid density functionals. Slater-type basis sets were used in the relativistic DFT calculations, based on the ZORA Hamiltonian and applied for heavy metals and their surfaces [562, 563].

Concluding this chapter we note that the relativistic correlated calculations of solids are mainly made in the DFT scheme. The post-DHF methods for solids wait their further development and implementation in computer codes for periodic systems.