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

4.2 Special Points of Brillouin Zone

125

√

a1 = (a, 0, c) a2,3 = (−a/2, ±a 3/2, c)

The basic translation vectors of the new lattice composed of supercells for symmetrical transformation (4.77) have the same form (the parameters a and c of the initial lattice are replaced with the parameters s1a, s2c)

√

a1 = (s1a, 0, s2c) a2,3 = (−s1a/2, ±s1a 3/2, s2c)

Inserting them in (4.77) one obtains nine equations for nine elements of the matrix l. The solution of this system is

l

ij

=

s2 + 2s1

δ

ij

+

s2 − s1

(1

−

δ

ij

)

3

3

where i,j=1,2,3. As the matrix elements lij must be integers let us assign (s2 −s1)/3 = n2 and (s2 + 2s1)/3 = n1 + n2. The matrix of the symmetrical transformation for a rhombohedral lattice with the corresponding value of L may be found in Appendix A and [86] where the matrices of the symmetrical transformations for all three-dimensional crystal lattices are given. Let us explain the peculiarities of the consideration made for each of the seven possible crystal systems.

In the triclinic crystal system an arbitrary matrix with integer elements defines a symmetrical transformation (any transformation seems to be symmetrical because of the low point symmetry of the lattice).

In the monoclinic crystal system there are two lattices: simple (P) and basecentered (A) each of which is defined by five parameters. Therefore the matrices of symmetrical transformations are determined by five integers.

In the hexagonal crystal system there is only one lattice type (P), but the basic translation vectors may be oriented in two di erent ways relative to the basic translation vectors of the initial lattice: either parallel to them or rotated through an angle of π/6 about the z-axis. Therefore, two types of symmetrical transformation are possible in this case (with two parameters for each).

In the orthorhombic crystal system the base-centered lattice merits special attention because of di erent possible settings. Let the initial base-centered lattice have the setting C. The transition to base-centered lattices with settings C and A (or B) gives di erent results (see Appendix A). The change of setting for the transition to other types of lattice does not give new supercells.

In the tetragonal crystal system there are two types of Bravais lattice (P and I). All their symmetrical transformations may be obtained from the symmetrical transformations for orthorhombic lattices if one sets n1 = n2 and takes into account that base-centered and face-centered orthorhombic lattices become simple and bodycentered tetragonal ones, respectively.

The matrices of symmetrical transformations for all the types (P, F, I) of cubic lattices may also be obtained from the matrices of symmetrical transformations for orthorhombic lattices if one sets n1 = n2 = n3 = n.

4.2.2 Special Points of Brillouin-zone Generating

In self-consistent calculations of the electronic structure of crystals both in the basis of plane waves and in the basis of localized atomic-type functions, one needs in every

126 4 Hartree–Fock LCAO Method for Periodic Systems

stage of the self-consistent procedure to evaluate an approximate electron-density matrix by integration over the Brillouin zone (BZ):

ρ(r, r ) = BZ Qr,r (k)dk,	Qr,r (k) =		ψv (k, r )ψv (k, r)	(4.78)
		v

where the sum is over occupied electronic states of the crystal. In practice, ψ(k, r) (and ϕ(k) ≡ Qr,r (k)) are calculated in some finite meshes of k-points. Therefore, one needs for the evaluation of the integral (4.78) to construct the interpolation formula of the numerical integration based on an interpolation procedure by means of some set of interpolation functions.

The functions ϕ(k) ≡ Qr,r (k)) are periodic in the reciprocal space with periods determined by the basic translation vectors bi of the reciprocal lattice and having the full point symmetry of the crystal (F is a point group of order nF of the crystal).

				3
ϕ(k + bm) = ϕ(k) = ϕ(f k),		f F,	bm =	i		(4.79)
					mibi
				=1
The plane waves exp(ik · an) seem to be the most			convenient as interpolation
					3
functions for the integrand in the BZ integration, where an =					i=1 niai	are direct
lattice translation vectors and a	i are primitive translations. It is			easy with plane waves

to take into account the translational and point symmetry of the crystal.

In the problem of BZ integration, nodes for the interpolation (k-points of the BZ) were called special points (SP).

Many procedures for special points generation have been proposed [87–90]. The supercell (large unit cell–small Brillouin zone) method appears to be the most general and fruitful in practical applications [86,87] as it gives as a rule the most e cient sets of special points.

Let Bi(Γ1) (i = 1, 2, 3) and bj (Γ2) (j = 1, 2, 3) be basic translation vectors of the reciprocal lattices corresponding to direct ones determined by basic translation

vectors a	i(Γ1) and aj (Γ2),	respectively. The transformation (4.77) of the direct lattices
			)
is accompanied by the following transformation of reciprocal lattices:
		)	i	)
	bj (Γ2) =			(l−1(Γ2Γ1))ij Bi(Γ1)	(4.80)

For the symmetrical transformation (4.77) the transformation (4.80) is also symmetrical as it does not change the point symmetry of the reciprocal lattice. The symmetrical transformation is compatible with the change of the reciprocal-lattice type in the limits of the same crystal system too. The vectors bj defining the small Brillouin zone are very important in the theory of special points [86]. Let f (K) be the function with a point symmetry F to be integrated over the initial Brillouin zone where the wavevector K varies. Usually, the point-symmetry group F either coincides with the crystal

Pm(K))	)	)	× Ci (otherwise) [86].
class F of the crystal (if F contains the inversion I) or F = F

The function f (K) may be expanded in Fourier series over symmetrized plane waves

f (K) = fmPm(K) (4.81)

m

		4.2 Special Points of Brillouin Zone	127
	1
Pm(K) =	nG	g G exp(iK · gam)	(4.82)

3

where am = i=1 miai is some translation vector of the direct lattice, the integer m = 0, 1, 2, . . . numerates the stars of vectors gam (g F ) in order of increasing length. The functions Pm(K) have the properties [88]

1	BZ Pm(K)dK = δm0, m = 0, 1, 2 · · ·	(4.83)
VBZ

Let m be a subset of m corresponding to the stars of vectors gam (g F ). As is proved in [86], for symmetrical transformation (4.77) L points in the initial Brillouin

zone (related to the initial basic translation vectors ai)
Kt(k) = k + j	qtj bj	(4.84)

where qtj are integers and k is an arbitrary vector in the small Brillouin zone (related to the basic translation vectors aj ), satisfy the relation:

L
Pm(Kt(k)) = L					P	m		(k)δmm		(4.85)
t=1
		m
or
N
wsPm(Ks(k)) =					P		m		(k)δmm	(4.86)
s=1
			m

In the latter relation s numbers the di erent irreducible wavevector K stars that contain the points (4.84) and ws = Ls/L (Ls is the number of points (4.84) belonging

to the sth star,		N	(k)	are usually chosen in the irreducible
		s=1 Ls = L). The points Ks
part of the	Brillouin zone.

Relations (4.81)–(4.83) give the following formula of the appoximative numerical integration:

	N
1
VBZ	BZ f (K)dK ≈ s=1 wsf (Ks(k))	(4.87)

Let M = so be the number denoting the set of vectors gam with the smallest (nonzero) lengths. If f (K) is some linear combination of Pm(K) with m < M , then the formula (4.87) appears to be exact. The number M characterizes the accuracy of the numerical integration formula. The special choice of k can either increase the accuracy M , or change the number N of points K(sk) (or both), [89].

To generate the set of points K(sk) for any of the 14 Bravais lattices it is su cient to find the inverse of the corresponding matrix from Appendix A, to pick out according to (4.84) L points in the Brillouin zone related to basic translation vectors ai and to distribute them over stars. The distribution of these points over stars depends on the symmetry group F of the function f (K) and can not be made in general form.

The sets of special points for numerical integration over the Brillouin zone of cubic crystals are given in Tables 4.1–4.3. They are obtained by symmetrical increasing of

128 4 Hartree–Fock LCAO Method for Periodic Systems

the unit cells of cubic lattices for L ≤ 32. The tables contain the sets (4.84) with k = 0 and some e cient sets with k = 0. Some sets were found earlier in [87,89]. The types of lattices (TL) composed of supercells are indicated in the second columns. The symbols of the matrices Q = l−1 of the transformation of basic vectors of reciprocal lattices are given in the third column. The matrices Qi themselves are given in Appendix B. The system of symmetrized plane waves Pm(K) depends on the crystal class. Therefore, the number M that characterizes the e ciency of the set of special points may be di erent for di erent crystal classes [86]. In Tables 4.1–4.3 the numbers M

are given for crystal classes Td, O, Oh and T, Th (in parentheses).

In addition, in Table 4.1 for L = 4 two special points (1/2, 0, 1/4) and (1/2, 1/4, 0) are related to di erent stars in crystal classes T, Th and each of them is a special point with the weight 1/4; in crystal classes Td, O, Oh they are related to the same star and this three-special-point set contains only one of these points with the weight 1/2.

Table 4.1. Special points of the Brillouin zone for the simple cubic lattice generated by

the symmetrical transformation; B1 =									2π	(1, 0, 0), B2 =				2π		(0, 1, 0), B3 =				2π				(0, 0, 1); Ks =
									a					a						a
α1B1 + α2B2 + α3B3 ≡ (α1, α2, α3)
L	T L	l(−1)	M	Special points Ks and weights ws in {}
2	F	21 Q2	2	(0, 0, 0) {21 } ( 21 , 21 , 21 ) {21 }
4	I	21 Q3	3	(0, 0, 0) {41 } ( 21 , 21 , 0) {43 }
			4	( 41 , 0, 0) {41 } [( 21 , 0, 41 ) {41 } ( 21 , 41 , 0) {41 }] ( 21 , 21 , 41 ) {41 }
8	P	21 Q1	4	(0, 0, 0) {81 } ( 21 , 0, 0) {83 } ( 21 , 21 , 0) {83 } 21 , 21 , 21 ) {81 }
			15(19)	( 81 , 81 , 81 ) {81 } ( 83 , 81 , 81 ) {83 }									( 83 , 83 , 81 ) {83 } ( 83 , 83 , 83 ) {81 }
16	F	41 Q2	7	(0, 0, 0) {	1		} ( 41 , 41 , 41 ) {21 } ( 21 , 0, 0) {									3	} ( 21 , 21 , 0) {	3			}
					16											16		16
				(1, 1, 1) {	1		}
					16
27	P	31 Q1	9	(0, 0, 0) {	1		} ( 31 , 0, 0) {92 } ( 31 , 31 , 0) {94 } ( 31 , 31 , 31 ) {												8			}
					27														27
32	I	41 Q3	12	(0, 0, 0) {	1		} ( 41 , 41 , 0) {83 } ( 21 , 0, 0) {									3	} ( 21 , 41 , 41 ) {83 }
					32											32
				( 21 , 21 , 0) {		3		} ( 21 , 21 , 21 ) {			1	}
						32					32

4.2.3 Modification of the Monkhorst–Pack Special-points Meshes

The most popular special-points (SP) sets are meshes obtained by a very simple algorithm proposed by Monkhorst and Pack (MP) in [90]:

3
i	2pi − n − 1
kp(n) = up(ni )bi, n = 1, 2, ..., pi = 1, 2, . . . , n, up(ni ) =	2n	(4.88)
=1

4.2 Special Points of Brillouin Zone

129

Table 4.2. Special points of the Brillouin zone for the face-centered cubic lattice ) generated by the symmetrical transformation

L	T L	l(−1)	M	Special points Ks and weights ws in {}
4	P	21 Q3	2	(0, 0, 0) {41 } (0, 21 , 21 ) {43 }
			2	( 41 , 41 , 21 ) {1}
			2	( 21 , 21 , 21 ) {1}
8	F	21 Q1	4	(0, 0, 0) {81 } ( 21 , 21 , 21 ) {21 } (0, 21 , 21 ) {83 }
				(0, 41 , 41 ) {41 }
			4						( 41 , 41 , 21 ) {21 } ( 41 , 41 , 21 ) {41 }
			8(10)	( 41 , 41 , 41 ) {41 }					( 21 , 21 , 41 ) {43 }
16	I	41 Q4	6(7)	(0, 0, 0) {	1		}		( 41 , 41 , 21 ) {43 } (0, 21 , 21 ) {						3	}
					16										16
				(0, 41 , 41 ) {83 }
			6(7)						( 21 , 21 , 21 ) {41 } ( 41 , 41 , 21 ) {83 }
27	F	31 Q1	9	(0, 0, 0) {	1		} ( 31 , 31 , 31 ) {			8		} (0, 31 , 31 ) {92 } ( 31 , 31 , 32 ) {94 }
					27					27
32	P	41 Q3	8	(0, 0, 0) {	1		} (0, 41 , 41 ) {				3		} ( 41 , 41 , 21 ) {83 } ( 21 , 21 , 21 ) {81 }
					32						16
				(0, 21 , 21 ) {		3		} ( 41 , 21 , 43 ) {		3		}
						32				16
			8	( 41 , 41 , 41 ) {41 } ( 21 , 21 , 41 ) {43 }

)B1 = 2aπ (−1, 1, 1), B2 = 2aπ (1, −1, 1), B3 = 2aπ (1, 1, −1)

K s = α1B1 + α2B2 + α3B3 ≡ (α1, α2, α3)

These meshes are widespread in modern Hartree–Fock and density-functional theory calculations of crystals and are automatically generated in the corresponding computer codes. But MP meshes of SP have at least one shortcoming, which may be understood if one oversees the problem from the point of view of the more general large unit cell–small Brillouin zone (LUC–SBZ) method of SP generation.The MP method is a particular case of the LUC–SBZ method. In [91] the modification of the MP meshes was suggested that makes faster and more regular convergence of the results of the self-consistent calculations of the electronic structure of crystals.

MP meshes (4.88) are a particular case of the meshes (4.84), which corresponds to the transformation (4.77) with the diagonal matrices

lji(n) = nδji, L = n3,

and a special choice of k in (4.84):

1 (b1 + b2 + b3) k = 2n

0

	n = 1, 2, ...,				i, j = 1, 2, 3	(4.89)
			1	(1, 1, 1) for even n
=			2n			(4.90)
			0		for odd n

The k-vector choice (4.90) is not the best one. As an example, we demonstrate this for cubic crystals. In this case, transformation (4.89) assures at least the cuto

130 4 Hartree–Fock LCAO Method for Periodic Systems

Table 4.3. Special points of the Brillouin zone for the body-centered cubic lattice ) generated by the symmetrical transformation

L	T L	l(−1)	M				Special points K s and weights ws in {}
2	P	21 Q2	2	(0, 0, 0) {21 } ( 21 , 21 , 21 ) {21 }
			2	( 41 , 41 , 41 ) {1}
			2	(0, 0, 21 ) {1}
			2	( 41 , 41 , 41 ) {1}
			3	(		1	, 41 , 41 ) {1}
						12
4	F	41 Q5	3	(0, 0, 0) {41 } ( 41 , 41 , 41 ) {21 } ( 21 , 21 , 21 ) {41 }
			3	( 41 , 41 , 41 ) {21 } (0, 0, 21 ) {21 }
			6	( 81 , 81 , 81 ) {21 } ( 83 , 83 , 83 ) {21 }
8	I	21 Q1	5	(0, 0, 0) {81 } (0, 0, 21 ) {43 } ( 21 , 21 , 21 ) {81 }
			5	( 41 , 41 , 41 ) {43 } ( 41 , 41 , 41 ) {41 }
			6	(0, 0, 41 ) {41 } (0, 41 , 41 ) {21 } ( 41 , 41 , 21 ) {41 }
16	P	41 Q2	6	(0, 0, 0) {				1		} ( 41 , 41 , 41 ) {83 } (0, 0, 21 ) {83 } ( 41 , 41 , 41 ) {81 }
								16
				( 21 , 21 , 21 ) {					1		}
									16
			6	( 81 , 81 , 81 ) {81 } ( 81 , 81 , 83 ) {21 } ( 85 , 83 , 83 ) {81 } ( 83 , 81 , 81 ) {41 }
27	I	31 Q1	10(11)	(0, 0, 0) {				1		} (0, 0, 31 ) {94 } ( 31 , 31 , 31 ) {92 } (0, 31 , 31 ) {																	8		}
								27																			27
32	F	81 Q5	12	(0, 0, 0) {				1		} ( 81 , 81 , 81 ) {41 } ( 41 , 41 , 41 ) {														3	} (0, 0, 21 ) {	3		}
								32																16		16
				( 81 , 83 , 83 ) {41 } ( 41 , 41 , 41 ) {												1	} ( 21 , 21 , 21 ) {							1	}
																16								32

)B1 = 2aπ (0, 1, 1), B2 = 2aπ (1, 0, 1), B3 = 2aπ (1, 1, 0)

K s = α1B1 + α2B2 + α3B3 ≡ (α1, α2, α3)

length Rcut = n ·|ai| whatever the k. The k-vector choice (4.90) increases slightly the cuto length Rcut for even n in an fcc lattice (see Table 4.4).

In all the other cases of cubic crystals (in bcc and sc lattices and for odd n in fcc lattice) the meshes MP (4.88) do not increase the standard cut-o length Rcut = n·|ai| and are equivalent to the meshes with k=0 that contain the point Γ and consist of the whole k-stars.

In [91] another choice of the k-point in (4.88) was proposed and the corresponding SP sets were called modified MP (MMP) meshes. They depend on the type of cubic crystal lattice.

1. ) = 1 (1, 1, 1) for sc lattice (MMP1=MMP2); k

4n

4.2 Special Points of Brillouin Zone

131

Table 4.4. Characteristics of MP, MMP1 and MMP2 k-points meshes for sc and bcc cubic lattices (Rcut and M are the cuto length and the precision of the meshes, N is the number of k-points in the IBZ; Rcut is in a units, a is a cubic lattice constant)

						4Rcut2	(in	a2), M, N
				sc							bcc
n	Γ M=MP			MMP1=MMP2			Γ M = MP			MMP1			MMP2
1	4	1	1	16	4	1	3	1	1	4	2	1	8	3	1
2	16	4	4	64	15	4	12	5	3	16	6	2	32	12	6
3	36	9	4	144	40	10	27	10	4	36	14	5	72	30	18
4	64	15 10		256	79	20	48	19	8	64	26	8	128	59	40
5	100 25 10			400	141	35	75	32	10	100	45	14	200 104		75
6	144 40 20			576	224	56	108 51		16	144	70	20	288 163 126
7	196 56 20			784	341	84	147 72		20	196 102 30			392 221 196
8	256 79 35			1024 488		120	192 99		29	256 143 40			512 265 288

2. k =		1			(1, 1, 1) (MMP1) and k =			1		( 1/3, 1, 1) (MMP2) for bcc lattice;
			n
)						4n
	4 1				)					1−
3. k) =						(1, 1, 1) (MMP1) and k) =					(0.4588, 0.3112, 0.1477) (MMP2) for fcc
			2n						2n

lattice.

The main characteristics of the proposed meshes for sc and bcc lattices are given in Table 4.4 and for fcc lattice in Table 4.5.

Table 4.5. Characteristics of Γ M, MP, MMP1 and MMP2 k-points meshes for fcc cubic lattice and convergence of PW DFT binding energy E per unit cell (in ev) for SiC crystal for Γ M, MP and MMP2 k-points meshes, [91](Rcut and M are the cuto length and the precision of the meshes, N is the number of k-points in the IBZ; Rcut is in a units, a is a cubic lattice constant).

		Γ M				MP				MMP1			MMP2
n	4Rcut2	M N		E	4Rcut2	M N		E	4Rcut2	M N		E	4Rcut2	M	N
1	2	1	1	−2.931	2	1	1	−2.931	4	2	1	12.843	6	3	1
2	8	4	3	12.236	16	8	2	14.990	16	8	2	14.990	24	13	6
3	18	9	4	14.615	18	9	4	14.615	36	20	6	15.100	54	33	18
4	32	17	8	14.987	64	40	10	15.110	64	40	10 15.110		96	67	40
5	50	29 10 15.076			50	29	10	15.076	100	71	19 15.111		150	121	75
6	72	47 16 15.100			144	113 28		15.111	144	113 28 15.111			216	180 126
7	98	68 20 15.108			98	68	20	15.108	196	163 44 15.111			294	227 196
8	128	98 29 15.110			256	207 60		15.111	256	207 60 15.111			384	253 288

132 4 Hartree–Fock LCAO Method for Periodic Systems

MMP2 meshes are the best ones as they give the maximal value of the cuto length Rcut and the precision M for the given transformation (4.89). MMP1 meshes for bcc and fcc lattices give the smaller value of the cuto length and the precision, but at the same time consist of the smaller number of points in the IBZ. The one-SP meshes MMP2 (for n = 1) are mean value points of Baldereschi, [92].

The idea of the proposed modification of the MP meshes can be also used to generate MMP meshes for all noncubic lattices.

The e ciency of MMP k-points meshes is demonstrated for an fcc lattice in Table 4.5, where PW DFT binding energies E per unit cell (in eV) for SiC crystal (fcc lattice) are given for di erent choices of k-meshes.

DFT-PW calculations [91] were performed on a SiC crystal for di erent SP meshes to investigate which one corresponds to the best convergence of the results for the total energy per unit cell. The experimental value a = 4.35 ˚A was taken for the fcc lattice constant of SiC, the electron–ion interaction was described by ultrasoft, Vandebildt-type pseudopotentials [93] , the cuto energy for the plane-wave basis set was taken to be 1000 eV. The generalized gradient corrections (GGA) DFT method was used (see Chap. 7).

In Table 4.5 are given the calculated binding energy values (defined as the difference between the sum of Si and C atomic energies with Vanderbildt-type pseudopotentials and the total energy per primitive unit cell) for three types of meshes mentioned above and used in calculations.

As is seen, the faster convergence takes place for the modified Monkhorst–Pack schemes. For comparison, the experimental 13.0 eV and Hartree–Fock 9.0 eV may be found in [94].

The suggested modification of Monkhorst–Pack special-points meshes for Brillouinzone integration is essentially useful for crystals with many atoms in a primitive unit cell or for point-defect calculations in a supercell model. In these cases each step of the self-consistent procedure in Hartree–Fock or DFT calculations is time consuming, so that the higher e ciency of k-point meshes shortens the computing time. It is also important for the lattice parameters or atom-position optimization in crystals, when the self-consistent procedure is repeated for di erent atomic structures. The special point set generation considered can be applied both in LCAO and plane-wave calculations to approximate the DM of a crystal.

4.3 Density Matrix of Crystals in the Hartree–Fock Method

4.3.1 Properites of the One-electron Density Matrix of a Crystal

As is well known, the energy of a system within the single-determinant Hartree–Fock approximation can be expressed in terms of the one-electron density matrix (DM). The one-electron spinless DM ρ(R, R ) is defined as

ρ(R, R ) = ψ(R, R2, ...RM ) ψ (R , R2, ...RM ) d3R2, d3R3...d3RM (4.91)

VN

where the electron position vectors R and R vary within the basic domain of a crystal of volume VN , i.e. a cyclic cluster consisting of N primitive unit cells. The DM of

ni(k) ϕik(R) ϕik(R )

4.3 Density Matrix of Crystals in the Hartree–Fock Method

133

the infinite crystal is obtained from the DM of the basic domain by taking the limit N → ∞. In its translational symmetry, the DM is periodic on the direct lattice:

ρ(R, R ) = ρ(R + Rn, R + Rn)

(4.92)

where Rn is an arbitrary translation vector of the Bravais lattice.

We will represent the electron position vector R in the form (r, Rn), where Rn specifies the primitive cell within which the end of the vector R lies and r is the position vector of an electron within this primitive cell. Therefore, we have R = r + Rn. Using (4.92), the density matrix can be written in the form

ρ(R, R ) = ρ(r + Rn, r + Rn ) = ρ(r, r + Rn − Rn) = ρr,r (Rn − Rn) (4.93)

The notation ρr,r (Rn) for the one-electron DM in the coordinate representation implies that the indices r and r of the matrix vary continuously only within the primitive cell. Therefore, there is an analogy between the properties of the DM in the coordinate representation and the properties of the DM represented in terms of a set of basis functions, for example, in terms of Bloch sums of atomic orbitals (AOs) or plane waves.

As is known, the diagonal elements of the one-electron DM in the coordinate representation are equal to the electron density:

ρ(R) = ρ(R, R) = ρr,r(0)

(4.94)

From the normalization condition for the many-electron wavefunction within the basic domain, it follows that

	ρ(R, R) d3R = N n =	Rn	ρr,r(0) d3r = N	ρr,r(0) d3r	(4.95)
VN		Va	Va

Therefore, the electron density is normalized to the number of electrons per primitive cell,

	ρr,r(0) d3r = n	(4.96)
Va

By using the single-determinant approximation to the many-electron wavefunction, the DM can be expressed through the one-electron wavefunctions (crystalline orbitals):

ρ(R, R ) = (4.97)

ik

where the index i specifies the energy bands and ni(k) are the occupation numbers. In nonconducting crystals, the energy bands are either completely filled or empty; therefore, ni(k) are independent of k and ni = 0, 2.

The one-electron DM is invariant under any orthogonal transformation in the space of occupied states. In particular, in nonconducting crystals, we can go over from the orthonormal set of extended Bloch states ϕik(R) to the orthonormal set of localized Wannier functions (see Chap. 3):

134 4 Hartree–Fock LCAO Method for Periodic Systems
		1
Wi(R − Rn) =		-			k		exp(ikRn) ϕik(R)	(4.98)
				N
In this case, (4.97) for the DM takes the form
ρ(R, R ) =	ni Wi(R − Rm) Wi (R − Rm)							(4.99)

iRm

or	ni	Wi(r − Rm) Wi (r − Rm + Rn)
ρr,r (Rn) =			(4.100)
	i	Rm

It is well known that the Wannier functions Wi(R) vanish exponentially as |R| → ∞ in crystals with completely filled bands. Since the vectors r and r lie in the reference (zeroth) primitive unit cell, the products of the Wannier functions on the right-hand side of (4.100) fall o exponentially with increasing |Rn|. Therefore, we may expect the total lattice sum in (4.100) for the o diagonal elements ρr,r (Rn) of the DM to also vanish exponentially with increasing |Rn|. It should be noted that in metals, the DM decays according to a power law.

Under translation through a lattice vector, according to Bloch’s theorem, the crystal orbitals transform according to irreducible representations of the translation group

ϕik(R + Rn) = exp(−ikRn) ϕik(R)

(4.101)

This condition is satisfied for both the infinite crystal and the basic domain; only the sets of values of the wavevector k for which (4.101) is satisfied are di erent in these two cases. By applying Bloch’s theorem (4.101) to the wavefunctions in (4.97) for the one-electron DM of the basic domain, we obtain

				1
	ρr,r (Rn) =	ni(k) ϕik(r) ϕik(r + Rn) =			exp( ikRn)Pr,r (k)
	i	k		N	−
					k	(4.102)
	where Pr,r (k) is the density matrix in k-space, which is defined as
		Pr,r (k) = N i	ni(k) ϕik(r) ϕik(r )			(4.103)

From the familiar orthogonality relations for columns and rows of the matrices in the irreducible representations of the Abelian translation group, it follows that

1		1
	exp(ikRn) = δk,b;		exp(ikRn) = δRn ,A	(4.104)
N	Rn	N	k

where b is a reciprocal lattice vector and A is a primitive translation of the basic domain as a whole. Using (4.104), it is easy to derive an inverse relation of (4.102)

Pr,r (k) = exp(ikRn)ρr,r (Rn)	(4.105)

Rn

It follows from (4.105) that Pr,r (k) is a periodic function in the reciprocal space