QMMS Lecture #8 (Venables/Heggie)

Notes for QMMS Lecture 8 (Venables/Heggie)

Lecture notes by John A. Venables and Edward Hernández. Lecture scheduled for 6 Dec 05. Latest version 4th November 2005.

The references for this lecture are here. Note that this lecture needs the Symbol font enabled on your browser.


8. Matrix (N-cubed) and order-N methods

8.1 The LCAO bands of silicon: an 8x8 matrix

In the computation sessions so far, we have been building up various one-dimensional (1D) problems, starting from particles in isolated (infinitely deep) wells, to finite wells, and on to coupled wells as models of solids. The latter models have illustrated Bloch's theorem, and allow us to calculate the band structure, E(k) of these idealised solids. Now we are interested in making the transition to a real, 3D, solid with a specific (diamond cubic) crystal structure, and we wish to calculate E(k), where the energy E is a function of both the magnitude and direction of the wave vector (k). The energy E is thus a graph in reciprocal space, where typically one plots E(k) along particular (relatively high symmetry) directions. Conventions include capital letters (e.g. K, L, X) for points on the Brillouin zone boundary (BZB) and greek capital letters (e.g. D, S, L) for particular directions in the first BZ, whose centre is at G. A handout was issued to make this notation clear. If anyone is still unclear about the crystal (real space) structures of semiconductors, nice Java applets which can visualise them in 3D are on a semiconductor web site run by Chu Wie at the University of Buffalo. Doubtless there is much more good stuff on this site too...

So, one can ask what is the smallest size of matrix which can attempt to model the band structure of diamond-structure semiconductors. The answer is 8x8, since there are two independent atoms per unit cell (at [0,0,0] and [1/4,1/4,1/4] in the f.c.c arrangement), and we need a minimal basis set of four states per atom to model the valence band. If we consider these as two s- and two p-electrons, we will not stabilise the diamond structure; but, as explained in the last lecture, if we consider one s- and three p- electrons, then we do. These electrons can be thought of (notionally, maybe) as forming sp3 hybrids; the key point is that we need to get back more energy, from overlap (of the hybrids), than we have to give up by promoting one s-electron to a p-state. In any case we need at least an 8x8 matrix to solve this problem by matrix methods. Note that this does not include any d-states, which have been shown to be important in getting the calculated band structure to agree with experiment. Silicon has an indirect gap (see problem 7.4.1), and the position of the minimum in the conduction band turns out to be a sensitive test of the involvement of d-states. These points are all too detailed for our present purposes.

There are two elements to the calculation of matrix elements involving only nearest neighbour overlap. The first is the size of the on-site and overlap integrals, which in principle are exactly the same as in the diatomic molecules discussed in lecture 1. However, now we have to take the 3D nature of the wavefunctions and the direction of bonding into account. This is summarised by Sutton on page 122, and by Pettifor on page 200-201. The largest contribution to the overlap is between hybrids pointing towards each other along the bond, and this is given by


b, or h = 1/4(sss - 2Ö 3sps - 3pps)         (8.1)

Here we need to recognise s-bonds, as in sss or sometimes Vsss, and also p-bonds. There are actually 4 of these, the other three being Vsps, Vpps and Vppp, in the notation used by Yu & Cardona (1996). The notation can be confusing here: Pettifor uses f for the hybrid and h for the bond integral, whereas Sutton uses h for the hybrid and b for the bond integral as in equation 8.1 above. In addition the promotion energy is called D, with


D = 1/4(Ep - Es)         (8.2)

If hybrids are used as the basis set, then D is the on-site matrix element between different hybrids on the same atom, and b is the matrix element between the main overlapping hybrids on the two different atoms. A gross approximation, which nevetheless retains some of the essential features, is to neglect all other overlaps; this is known as the Weaire-Thorpe model, which depends on just these two parameters D and b.

8.2 A detailed tight-binding calculation for silicon

If the first part of the calculation was the size of the bond integrals, the second is the disposition of the nearest neighbours at vector positions d, which gives rise to sums of phase factors of the form exp(ik.d). Thus the matrix elements are different, via this factor, for each value of k. To calculate the band structure (of silicon) we have to have all these features in place. The details were in the following computation session, and the method is described in outline below.

The classic tight binding calculation, using the minimal basis set needed to describe silicon, was developed in a series of papers by Chadi and Cohen, and is described in detail by Yu and Cardona (1996, section 2.7, pages 78-91).

The first decision to make is whether to use the sp3 basis set or the s, px, py, pz set. It turns out that unless one is going to make further approximations, such as neglecting small overlaps, there is no advantage in the hybrid set, and the latter set is just fine. Then the matrix elements are proportional to Vss, Vsp, Vxx and Vyy, in the notation used by Yu & Cardona (1996, equation 2.80, page 84). Each of these are linear combinations of the 4 V's given previously, as shown below in equation 8.3. These details can be explored further via problem 7.4.2, or equivalently by reverse engineering Edward's computer program for the silicon band structure. In addition, Jony Hudson did a project on this topic in 2000, and produced a mathematica program which can be seen here, complete with animations. In the next project of this type Fridrik Magnus has produced a MatLab program to explore the differences between C(diamond), Si and Ge: why are the band structures of these supposedly similar materials different?


Vss = 4Vsss; Vsp = -(4/Ö 3)Vsps;
3Vxx = 4Vpps + 8Vppp and
3Vyy = 4Vpps - 4Vppp         (8.3)

The sums over phase factors, called g1 to g4, are given by Yu & Cardona (1996, equations 2.81-82, page 85). They have a structure very like the hybrids themselves, reflecting the four neighbours and whether the phases are positive or negative, as in equation 7.2. But, as a warning, take care to get the signs right if you actually use these equations, don't take my word for it.

Now, we can begin to construct the 8x8 matrix itself, in terms of 8 columns labelled S1, S2, X1, Y1, Z1, X2, Y2, Z2 and equivalently 8 rows, as in Table 8.1 below. Here 1,2 represent the two non-equivalent atoms, S is an s-state, and X,Y,Z are p-states oriented along the corresponding axes. Since we are using atomic s- and p-orbitals for the basis, you should be able to fill in which matrix elements are zero by orthogonality. The diagonal elements should also be clear. The matrix, as all Hamiltonian matrices, is Hermitian, so that gives a further clue. The rest, however, is in the details given below the table, which you can find in Yu and Cardona (1996, Table 2.25, page 85). In Edward's program, which was implemented from this reference, specific directions for the k-values have been chosen so as to produce the equivalent of Yu and Cardona's figure 2.24. This is not a trivial exercise- Jony Hudson spent a good deal of time on this during his project in 2000.

Table 8.1: 8x8 Matrix for the s- and p-bands of the diamond structure within tight binding
S1
S2
X1
Y1
Z1
X2
Y2
Z2
 
A
C
0
0
0
D
E
F
S1
 
A
-D*
-E*
-F*
0
0
0
S2
   
B
0
0
G
K
H
X1
     
B
0
K
G
H
Y1
       
B
H
J
G
Z1
         
B
0
0
X2
           
B
0
Y2
             
B
Z2

Here A = Es - Ek; B = Ep - Ek; C = Vssg1; D = Vspg2; E = Vspg3; F = Vspg4; G = Vxxg1; H = Vxyg2; J = Vxyg3 and K = Vxyg4. From Yu and Cardona (1996, Table 2.25).

The Vss etc., values to be used in this program were given by Qian and Chadi (1987), as discussed by Sutton, pages 124-131. These values (in eV) are D = 1.61, (sss) = -1.9375, (sps) = 1.745, (pps) = 3.050, and (ppp) = -1.075. (Don't be overawed by the apparent accuracy here, the .025's must result from dividing by 4). Thus such a calculation represents in practice a five-parameter fit to the experimental band structure of silicon. Different values can be used to model the band structures of diamond or germanium, with appropriate values given by Yu and Cardona (1996, Table 2.26, page 86), based on Chadi and Cohen's work. The figures and discussion which follow indicate, both qualitatively and quantitatively, some of the differences between these solids. Here, we have set out on a walk through the foothills of this topic; there are a lot more details to explore if we had the time or inclination to put on our climbing boots, and the stamina to hold on tight and keep going.

8.3 Large scale cluster and "order-N" calculations

The calculation described in the last section is a 'solid state' type calculation, in that the periodicity of the lattice is implied, and built into the lattice sums which yield the various g's. Yet the structure of using a matrix really doesn't build this in; we got down to an 8x8 matrix by using translational and other crystal symmetry to say that there are only 2 non-equivalent atoms. If we had been prepared to allow a large number, N, of non-equivalent atoms, then the symmetry would not have been needed. But the entity we would then be studying is a (large) cluster. With M basis states associated with each atom, the size of the matrix now becomes the total number of electronic states; i.e. it rapidly becomes huge, of order (NM)x(NM).

When computers had small memories, it was essential to use symmetry arguments to the maximum extent to reduce the size of the matrix, and hence the time of the calculation. It still is a good idea, and at the limit of what is possible, this will always be important. However, as speeds and memory have increased dramatically, symmetry (or equivalently group theory) arguments have become less important, and small clusters can now be routinely solved very quickly. The development of these cluster codes was the reason why John Pople received his share of the 1998 Nobel prize for Chemistry - the other half going to Walter Kohn for density functional methods. Several cluster and DFT codes are available locally on our bfg computer, and the Theoretical Chemistry group is linked into various consortia which use and develop such codes, e.g. AIMPRO at Exeter and Newcastle.

As one can appreciate from this outline description, actually doing a real calculation is computationally very intensive, and it also has to be repeated for different ionic positions R, and relaxed to the minimum energy configuration in the most efficient way. Individual orbitals are often expanded in 'Gaussian orbitals', or in plane waves with a large number M of independent coefficients ck of the various k-vectors which have to be computed as the calculation proceeds. For N electrons, the number of computing operations to diagonalise a matrix scales as (MN)3. However, many of the operations required for each k are identical, so the code can be written for implemention on parallel (super)computers. Now, systems with N ~several 100’s and typical M~1000 can be tackled. The virtue of pseudo-potential calculations is that they reduce the number of electrons per atomic site, at the cost of increased complexity of the ionic (external) potential. There is a strong impetus to reduce the cubic power law to something lower, which one can summarise by exclaiming: O(N) methods are in! This is a major area for development, with even small advances being highly prized.

However, no actual methods are as good as O(N), the best perhaps scaling as Nln(N), and one also needs to look critically at the multiplying constants. Many careers have been spent trying to crack these highly technical conceptual and computational problems. A feel for how this is going can be gauged by consulting Goringe et al. (1997), Turchi et al. (1998), or Gödecker (1999), and doubtless yet more in the period 2000-05. If you do that, you will know more than we do: so we're done for this year. See you sometime in the next.


Return to timetable or to course home page.