Quantum Chemistry of Solids / 14-Hartree-Fock LCAO Method for Periodic Systems

4.3 Density Matrix of Crystals in the Hartree–Fock Method		135
Pr,r (k + bm) = Pr,r (k)		(4.106)
Using the hermiticity of the DM
ρ(R, R ) = ρ (R , R)		(4.107)
we find that
ρr,r (Rn) = ρr ,r(−Rn);	Pr,r (k) = Pr ,r(−k)	(4.108)

From (4.96), the normalization condition for the DM in k-space can easily be found to be

Pr,r(k) d3r = n (4.109)

Va

Here, the integration is performed over a primitive cell of volume Va.

In the case where the many-electron wavefunction is calculated in the singledeterminant approximation, the spinless DM is duodempotent:

	ρ(R, R )ρ(R , R )d3R = 2ρ(R, R )		(4.110)
VN
Similar relations also hold for the matrices ρr,r (Rn) and Pr,r (k)
Va
Rm	d3r ρr,r (Rm)ρr ,r (Rn − Rm) = 2ρr,r (Rn)		(4.111)
		d3r Pr,r (k) Pr ,r (k) = 2Pr,r (k)	(4.112)
	Va

The one-electron density matrix has the same point-symmetry group F as the crystal structure.

Let operator gˆ = tvR correspond to the coset representative tvR in the coset decomposition of the crystal structure space group G over translation subgroup T (see (2.15)). In this notation the density matrix point symmetry can be written in the form

gρˆ (R, R ) = ρ(ˆg−1R, gˆ−1R ) = ρ(R, R )						(4.113)
For the density matrix ρrr (Rn) in the direct space, we can write
gρˆ rr (Rn) = ρgˆ	1r+Rn,gˆ	1r	+R (R−1Rn + Rm	−	Rn) = ρrr (Rn)	(4.114)
−	−		n

where Rm is the translation vector that returns gˆ−1r back to the reference unit cell. In a similar manner, the R n vector returns R−1r back to the reference cell. Using relation (4.105) for the DM in the reciprocal space, we obtain

gˆPrr (k) = exp iR−1k(r − Rn) PR−1r+Rn ,R−1r +R n (R−1k) = Prr (k) (4.115)

Hence, we have a useful relation between matrix elements of the DM for the star of the wavevector k:

136	4 Hartree–Fock LCAO Method for Periodic Systems
	Prr (Rk) = exp ik(Rn − R)) PR−1r+Rn ,R−1r +Rn	(4.116)
	In practical calculations of the band structure of crystals by the Hartree–Fock

method, the one-electron density matrix Prr (k) in the reciprocal space can be calculated only for a rather small finite set of special points (see Sect. 4.2). Let us consider such a set of special points {kj }, j = 1, 2, . . . , N0. Then, in the calculation of the density matrix of the infinite crystal the integration over the Brillouin-zone volume VB is changed by the sum over the special points chosen

1

VB

1

VB

d3kPrr

(k) =

N0

kj

Prr (kj )

(4.117)

Thus the the approximate density matrix ρ˜(R, R ) is introduced

ρ˜(R, R ) =

1

j

N0

k

Prr (kj )

j

ρ˜r,r (Rn) =

1

exp (−ikRn) Pr,r (kj )

(4.118)

N0

k

The approximate density matrix ρ˜(R, R ) satisfies the translation-symmetry conditions. This matrix also satisfies the point-symmetry condition, provided that the whole star of each vector kj is included in the set of special points {kj }. However, the approximate density matrix ρ˜(R, R ) does not satisfy other properties of the exact DM. It is easy to see that an arbitrary finite sum over vectors kj , of delocalized Bloch functions ϕikj (r) does not decrease with |r| → ∞. Therefore, provided that the R vector in expression (4.97) for the DM is fixed, the approximate density matrix ρ˜ does not decrease and does not approach zero as |R −R | → ∞. This incorrect asymptotics

of the approximate DM is the reason for the divergence in the calculation of ρ2 and for the violation of the idempotency conditions. Indeed, using the orthonormality of the set of Bloch functions ϕik, it is easy to verify [95] that

ˆ2		N	ˆ	→ ∞ at N → ∞
ρ¯	= 2 N0 ρ¯				(4.119)

The incorrect asymptotic behavior of the approximate density matrix ρ˜ gives rise to divergences in calculations of the average values of some physical quantities. In particular, for the exchange energy per unit cell Kex, we obtain the divergence

Kex = N		d3R	d3R |ρ˜R		R	|		→ ∞ at N → ∞	(4.120)
1				(R, R )			2
	VN	VN		| −		|

It was empirically found that, in a calculation of the electron structure of crystals with nonlocal exchange, the integration of the Hartree–Fock exchange potential over the direct lattice requires a so-called exchange-interaction radius to be introduced. The exchange-interaction radius cannot be arbitrarily large but must correspond to the used number of points k in the Brillouin zone.

This fact was explained in [95]. It was shown that the use of a finite number of special points and simple cubature formulas for integration over the Brillouin zone

4.3 Density Matrix of Crystals in the Hartree–Fock Method

137

leads to an incorrect behavior of the o diagonal matrix elements of the DM at large distances. This leads to nonphysical divergences in the calculation of the Hartree– Fock nonlocal exchange potential and the Hartree–Fock exchange energy of a crystal. Similar divergences take place in the calculation of the square of the DM. This naturally violates the idempotency of the exact one-electron matrix required of the Hartree–Fock one-determinant method. In the theory of special points, the divergence problem arises only in the cases where the total exchange energy and the exchange potential of the crystal depend on o diagonal elements of the DM. Thus, in contrast to the Hartree–Fock method, this problem is nonexistent in methods based on densityfunctional theory; it is precisely the approximation where the theory of special points was first applied. In calculations of the electron structure with nonlocal exchange, this theory must be modified to avoid the above-mentioned divergences (see Sect. 4.3.3). The considered properties of DM (translation and point symmetry, duodempotency and normalization properties) do not depend on the basis choice (LCAO or plane waves). In the next section we consider the density matrix of a crystal in the LCAO approximation.

4.3.2 The One-electron Density Matrix of the Crystal in the LCAO Approximation

In Sect. 4.1.5 the Hartree–Fock LCAO approximation for periodic systems was considered. The main di erence of the CO LCAO method (crystalline orbitals as linear combination of atomic orbitals) from that used in molecular quantum chemistry, the MO LCAO (molecular orbitals as linear combination of atomic orbitals) method was explained. In the CO LCAO approximation the one-electron wavefunction of a crystal (CO - ϕik(R)) is expanded in Bloch sums χµk(R) of AOs:

	ϕik(R) =				Ciµ(k)χµk(R)	(4.121)
				µ
where	1
					exp(ikRn)χµ(R − Rn)	(4.122)
	χµk(R) =	√
			N	Rn

In (4.121) and (4.122), the index µ labels all AOs in the reference primitive cell (µ = 1, 2, . . . , M ) and the index i numbers the energy bands (i = 1, 2, . . . , M ). The appeaance of energy bands for periodic systems is the result of translation repeating of the primitive unit cells AO over the infinite crystal of the cyclic system modeling the infinite crystal. In the calculation of the molecular electronic structure the MOs are filled by electrons, taking into account the degeneracy of levels due to the pointsymmetry group. In crystals, COs are filled by electrons and the degeneracy of oneelectron states is defined by irreps of the space-symmetry group of a crystal. This means that the one-electron energy bands are filled by electrons and for nonconducting crystals the valence bands are filled and the conduction bands are empty.

The Bloch sums (4.121) of AOs, as well as the AOs, do not constitute an orthonormal basis; that is, the overlap integrals

138 4 Hartree–Fock LCAO Method for Periodic Systems
Sµν (k) = d3Rχµk(R)χνk(R)
sµν (Rn) =	d3Rχµ(R)χν (R − Rn)	(4.123)

are not equal to δµν .

The crystalline orbitals ϕik(r) compose an orthonormal set from which the following orthonormality relations can be derived for the elements of the matrix C(k) involved in (4.57):

{C(k)S(k)C(k)}ij = δij

(4.124)

In terms of the Bloch sums of AOs, the DM elements Pµν (k) can be expressed as

Pµν (k) = ni(k)Ciµ(k)Ciν (k) (4.125)

i

Within the LCAO approximation, the DM elements in the coordinate space are given by an expression similar to (4.101),

ρµν (Rn) =	1		exp(−ikRn)Pµν (k)	(4.126)
	N k

and the DM in the reciprocal space is related to the DM in the direct (coordinate) space through a relation similar to (4.105)

Pµν (k) = exp(ikRn)ρµν (Rn)

(4.127)

Rn

In the reciprocal space, the analog of the normalization condition (4.109) for the DM in the AO representation is the relation

Sp (P(k)S(k)) =	Pµν (k)Sµν (k) = n	(4.128)
	µν

The normalization condition for the DM in the coordinate space (analog of (4.95)) is

Sp [ρ(Rn)s(−Rn)] = n

(4.129)

Rn

The idempotency relation for the density matrix P(k) in the reciprocal space (with allowance for the nonorthogonality of the basis) has the form (ni(k) = 0, 2)

P(k)S(k)P(k) = 2P(k)	(4.130)
In the coordinate space, the idempotency relation for the DM is written as
m
ρ(Rm)S(Rm − Rm)ρ(Rn − Rm) = 2ρ(Rn)	(4.131)

Rm ,R

In various semiempirical versions of the Hartree–Fock approximation (see Chap. 6), the orthonormal set of L¨owdin atomic orbitals (LAOs) rather than the nonorthogonal AO basis is used; the LAOs are defined as

4.3 Density Matrix of Crystals in the Hartree–Fock Method		139
χµLk(R) =	Sµν−1/2(k)χνk(R)	(4.132)

ν

In this basis, the normalization condition for the density matrix PL(k) is simplified; instead of (4.129) we have

Sp PL(k) = n, Sp ρL(0) = n

(4.133)

In the LCAO basis, the idempotency relations for the DM in the coordinate and reciprocal spaces are similar to (4.128) and (4.129) for the DM in the coordinate representation; in the reciprocal space, we have

PL(k)PL(k) = 2PL(k)	(4.134)
and in the coordinate space, the relation has the form
ρL(Rm)ρL(Rn − Rm) = 2ρL(Rn)	(4.135)

Rm

In particular, for the reference primitive cell (Rn = 0), with allowance for the hermiticity of the DM, we have

ρL(Rn)ρL (Rn) = 2ρL(0)		(4.136)
Rn
The o diagonal elements of the DM in the AO basis determine the quantities
WAB (Rn) =
	|ρµνL (Rn)|2	(4.137)

µ A,ν B

that can be considered as the extension to the crystals of the Wiberg indices introduced for molecules [96]. The Wiberg indices WAB (Rn) can be interpreted as chemical-bond indices (orders) between atoms A and B, belonging to the reference and Rnth primitive cells, respectively [97, 98]. These indices are subject to a relation that is a consequence of the idempotency of the DM. To derive this relation, we consider the diagonal matrix elements of (4.136) and carry out summation over all AO

indices of atom A. The result is
	WAB (Rn) = 2ρAL −		ρµµL (0)	(4.138)
	B =A Rn		µ A
Here, ρL	is the total electron population (in Lowdin’s sense) of atom A,
A
	ρAL =	ρµµL (0)		(4.139)

µ

Let us define the covalence CA of atom A as the sum of the chemical-bond orders (Wiberg indices) between atom A and all other atoms of the crystal. Using the

idempotency relation (4.138), we have
CA =	WAB (Rn) = 2ρAL −	ρµµL (0)	(4.140)
	B,Rn =A,0	µ A

140 4 Hartree–Fock LCAO Method for Periodic Systems

It follows from (4.140) that the covalence of atom A in the crystal can be calculated either by summing the bond indices between atom A and all other atoms of the crystal or by using only single-center DM elements related to atom A. This property of the covalence is a consequence of the idempotency of the DM.

The Wiberg indices (4.137) and atomic covalence (4.140) are called the local properties of the electronic structure of periodic systems. These properties also include AO populations, atomic charges (electrovalencies) and total atomic valences [97]. The analysis of the local properties of the electronic structure in molecular quantum chemistry is very popular as it gives useful information about the chemical bonding. The local properties of the electronic structure of crystals are considered in more detail in Chap. 9. The above consideration holds for the density matrix of the basic domain of the crystal; that is, it is assumed that the number N of primitive cells in this domain is so large that the introduction of cyclic boundary conditions virtually does not a ect the density matrix of the infinite crystal.

In actual practice, the SCF calculations of periodic systems are made with the approximate DM ρ˜(R, R ), which is calculated by summation in (4.102) over a finite (relatively small) number of k-points. As a result, the divergence of the summation in the direct lattice can appear in the nonlocal exchange part of the Fock matrix. This problem is considered in the next section.

4.3.3 Interpolation Procedure for Constructing an Approximate Density Matrix for Periodic Systems

Let the electronic structure of a crystal be calculated using the one-electron wavefunctions found at a finite number L of points {kj } in the BZ (j = 1, 2, . . . , L). This raises the question of how the sum over k points in the BZ should be approximately calculated in (4.102) for the one-electron DM of the crystal.

We consider a set of points {kj } generated by the largeunitcell (LUC) – small brillouin zone (SBZ) method, see Sect. 4.2.2. In this method, the primitive lattice vectors ai(i = 1, 2, 3) are transformed with the aid of a matrix l whose elements are integers (see (4.77)):

3
i	(4.141)
AjL = ljiai, L = | det l|
=1

The basis vectors ALj determine an LUC and a new Bravais superlattice. The LUC thus constructed has volume VL = LVa and consists of L primitive cells. The superlattice vectors A are linear combinations (with integral-valued coe cients) of the basis vectors ALj . The matrix l in (4.141) is chosen such that the point symmetry of the new superlattice is identical to that of the original lattice (the corresponding transformation (4.141) is called a symmetric transformation, see Sect. 4.2.1). The type of direct lattice can be changed if there are several types of lattice with the given point symmetry. The LUC is conveniently chosen in the form of a Wigner–Seitz (WS) cell, which possesses the point symmetry of the lattice.

We introduce the periodic boundary conditions for the crystal domain coinciding with the LUC; that is, we assume that all translations through the superlattice vectors A are equivalent to the identity translation. Thus, we have a system of finite size, i.e. a cyclic cluster belonging to the symmetry group GL = TLF, see Sect. 2.2.3

4.3 Density Matrix of Crystals in the Hartree–Fock Method

141

(we consider only symmetric transformations). Here, the subgroup TL includes L translations through the vectors R0n of inner translations of the original direct lattice that lie within the LUC or fall on its boundary. The lattice sites lying on the boundary of the LUC are connected by superlattice vectors A. These lattice sites should be counted only once, because they belong simultaneously to several LUCs.

For the cyclic cluster thus constructed, the following orthogonality relations hold:

	1
			exp(ikj Rn0 ) = δk,b	(4.142)
	L		R0
			n
1	j
L	k		exp(−ikj Rn0 ) = δRn0 ,A	(4.143)

These relations are a generalization of the analogous equations for the basic domain of the crystal because the cyclic cluster is obtained with the help of transformation (4.141), in which the matrix l can be nondiagonal. The vectors kj in (4.142) and(4.143) label L di erent irreducible representations of the group TL and can be found from the relation

exp(ikj A) = 1, (j = 1, 2, . . . , L)

(4.144)

Equation (4.142) is a consequence of the orthogonality of the characters of irreducible representations of the translation group to the character of the unit representation (k = 0), while (4.143) means that the characters of a regular representation of the group are equal to zero for all elements of the group except for the identity element (i.e. except for the identity translation and equivalent translations through the superlattice vectors A).

Let the density matrix P(k) be known at a finite set of points determined by the LUC-SBZ method and, therefore, satisfying (4.144). Our aim is to approximate the DM at an arbitrary point k in the BZ. The interpolation procedure suggested in [70] and discussed in this subsection is appropriate for calculations in both the coordinate representation and the AO basis. For this reason, we drop the indices on the DM, keeping in mind that these indices are r and r in the coordinate representation and µ and ν in the AO or orthogonalized-AO representation.

The expansions of the density matrix P(k) given by (4.105) and (4.127) can be rewritten in the form

P(k) = exp(ikRn0 )ρ(Rn0 ) +	exp(ik(Rn0 + A))ρ(Rn0 + A) (4.145)
Rn0	A =0Rn0

where the translations (R0n + A) lie in the basic domain of the crystal. As mentioned in Sect. 4.3.2, the o diagonal elements of the density matrix ρ(R0n + A) fall o with distance as Wannier functions (exponentially in the case of insulators). Therefore, as the LUC grows in size and the values of |A| become su ciently large, the second term in (4.145) will be small in magnitude. With this in mind we will approximate the

˜

density matrix P(k) of the crystal as follows [70]. In the expansion given by (4.145), we drop the sum over the superlattice sites with A = 0 and take the remaining expression as an interpolation formula for determining the DM at any k point in the BZ; we rewrite this expression in the form

142 4 Hartree–Fock LCAO Method for Periodic Systems

˜		0	(4.146)
	0 0
P(k) =	exp(ikRn)ρ	(Rn)
	R0
	n

This formula was proposed in [70] for the interpolation of an arbitrary periodic function f (k). The interpolation coe cients ρ0(R0n) (the number of which is L) can be found from the condition

˜		0	(4.147)
	0 0
P(kj ) =	exp(ikj Rn)ρ	(Rn) = P(kj )
	R0
	n

Using the orthogonality relation (4.143), the interpolation coe cients can be found
to be	1			j
ρ0(Rn0 ) =				exp(−ikj Rn0 )P(kj )	(4.148)
	L
				k

The coe cients ρ0(R0n) can also be represented as a sum of the DM elements ρ(R0n)

over the superlattice sites. Indeed, substituting (4.146) for P(kj ) into (4.148) and
using (4.143), we have
ρ0(Rn0 ) =	ρ(Rn0 + A)	(4.149)

A

It should be noted that the matrix ρ0(R0n) can be defined for all vectors Rn of the Bravais lattice by using the appropriate extensions of (4.148) and (4.149). It is easy to see that ρ0(R0n) is a periodic function of period A. Substituting (4.148) for the coe cients ρ0(R0n) into (4.149), we obtain an interpolation formula for the DM in the reciprocal space,

˜	j
P(k) =				P(kj )Ωj (k)	(4.150)
		k
	1
Ωj (k) =				exp(i(k − kj )Rn0 )	(4.151)
	L
				R0
				n

Here, Ωj (k) are interpolation weights, the sum of which is equal to unity (the normalization natural for weighting factors). Indeed, using (4.143), we find

			1			(4.152)
	Ωj (k) =	exp(ikRn0 )	L		exp(−ikj Rn0 ) =	exp(ikRn0 )δRn0 ,0 = 1
kj	R0			kj	R0
	n				n

For the appropriate DM in the coordinate space, one can write equations similar to (4.118) and (4.133):

	1					1
ρ˜(Rn) =			exp(−ikRn)P(˜ k) =		ρ0	(Rn0 ) "			exp(ik(Rn0 − Rn))#	(4.153)
	N k			R0			N	k
				n

According to (4.143), the expression in the square brackets on the right-hand side of (4.153) is equal to unity if the vector Rn belongs to the set of vectors {R0n} (i.e. if this vector lies within the LUC or on its boundary) and vanishes otherwise. Therefore, the appropriate DM can be represented in the form

4.3 Density Matrix of Crystals in the Hartree–Fock Method

				ρ˜(Rn) = ω(Rn)ρ0(Rn)
where ω(Rn) is the so-called weighting function,
							{	}
ω(R	) =	1		exp(ik(R0		R	)) = 1, if Rn {Rn0	}
n		N		n	−	n	0, if Rn / Rn0
			k

143

(4.154)

(4.155)

We note that the procedure described above for interpolating the DM in the BZ is not uniquely determined, because the LUC (i.e. the set of vectors {R0n}) can be variously chosen for the same superlattice. Furthermore, the LUC can be chosen di erently for the pairs of DM indices r, r and µ, ν. In this book, the LUC is taken to be the Wigner–Seitz cell, because only this cell has a symmetry identical to the point symmetry of the superlattice in all cases. In order to correlate the LUC with a cyclic cluster, we choose the LUC to be dependent on the pair of DM indices as follows. In the coordinate representation, the LUC (VA-region) is centered at the point (r − r ); therefore, we have

ω

rr

(R

) = ω(R

n

+ r

−

r) =

1,

if

Rn + r

r VA

(4.156)

n

0,

if

Rn + r

− r / VA

−

In the AO representation, the LUC is centered at the point dµ − dν where dµ and dν are the position vectors of the two atoms to which the AOs with indices µ and ν belong, respectively. Thus, we have

ω

µν

(R

) = ω(R

n

+ d

ν −

d

) =

1,

if

Rn + dν

dµ VA

(4.157)

n

µ

0,

if

Rn + dν

− dµ / VA

−

The weighting function ωµν (Rn) of (4.157) introduced into expressions for the DM specifies the cyclic boundary conditions and the cyclic cluster. Indeed, let an arbitrarily chosen LUC be fixed and let us consider the orbitals of atoms A and B in this LUC (µ A, ν B ). Out of all matrix elements ρ˜µν (Rn + A) with indices µ and ν kept fixed and the vector A running over the superlattice, only one matrix element is nonzero. For this matrix element, the vector (dν + Rn + A) (the position vector of atom B) falls into the Wigner–Seitz cell centered at atom A site. This matrix element exactly equals the matrix element ρ0µν (Rn).

According to (4.154), the approximate density matrix ρ˜(Rn) found by interpolation in the BZ contains the weighting function of (4.155)–(4.157) as a factor. This function ensures the proper behavior of the o diagonal elements of the approximate DM as |Rn| → ∞. As already mentioned, the matrix without a weighting factor is a periodic (not vanishing at infinity) function

ρ0(Rn + A) = ρ0(Rn)

(4.158)

However, this DM is frequently used in many calculations based on the Hartree–Fock approximation or its semiempirical (with nonlocal exchange) versions for crystals (the CNDO and INDO methods, see Chap. 6). In those calculations, all summations over the lattice sites are usually truncated by introducing artificially interaction ranges. The non decaying density matrix ρ0(Rn) gives rise to a divergent exchange term in

144 4 Hartree–Fock LCAO Method for Periodic Systems

the Fock matrix. In other words, as the corresponding interaction range increases in size at a fixed number of involved k points, the results do not converge to a certain limit and the total energy of the system sharply decreases.

In order to avoid these divergences,the exchange interaction range should be chosen such that the corresponding sphere di ers only slightly from the Wigner–Seitz cell of the superlattice, which generates precisely the set of k points used in the calculations. In this case, the size of the summation domain in the coordinate space is in accord with the number of k points used in the calculations and the exchange term in the HF operator does not diverge. The approximate density matrix ρ˜0(Rn) is not subject to these drawbacks, and the balance between the size of the summation domain in the coordinate space and the number of {kj } points involved occurs automatically.

It should be noted that in the versions of the density-functional theory in which the exchange-correlation term depends on the electron density alone, both approaches are equivalent. Indeed, the electron density ρ(R, R) = ρrr(0) depends only on the diagonal elements of the DM; therefore, the weighting function of (4.155) and (4.156) is equal to unity in this case.

In general, the approximate DM does not satisfy all the conditions to which the exact DM is subject. Let us elucidate which of the relations presented in Sect. 4.2.2 holds for the approximate DM and which do not.

The normalization conditions (4.95), (4.128), and (4.129) are very important. The approximate DM in the coordinate representation and in the LAO basis meets these conditions because the weighting function for the diagonal elements of the DM is equal to unity. In the nonorthogonal basis, a modified normalization condition is satisfied,

Sp [ρ˜(Rn)˜s(−Rn)] = n

(4.159)

Rn

where s˜(Rn) is an approximate overlap-integral matrix, which is obtained by interpolating in the BZ in much the same way as the approximate DM was obtained and has the form

s˜µν (Rn) = ω(Rn + dν − dµ)sµν0 (Rn)						(4.160)
		1		j
sµν0 (Rn) =	L			exp(−kj Rn)S(kj ) =	s(Rn + A)	(4.161)
				k	A

It is easy to verify that in all cases the approximate DM is Hermitian, i.e. it obeys relations identical to (4.107) and (4.108).

In general, the approximate DM is not idempotent, because (4.112) holds only at points k = kj (j = 1, 2, . . . , L) and is not satisfied at other points of the BZ. For this reason, (4.111) and (4.136) in the coordinate space do not generally hold. However, in the important particular case where the vector Rn in these equations is zero, the idempotency relation is satisfied. In the coordinate space and in the LAO

representation, we have
	ρ˜(Rm)ρ˜(−Rm) = 2ρ˜(0)	(4.162)

Rm

Therefore, the important relation for the Wiberg indices (4.140) is also satisfied. We note that the matrix ρ0(Rn) obeys the relation