Faruduy Symp. Chem. SOC.,1984 19 79-84 Ab Initio Hartree-Fock Calculations for Periodic Systems BY V. R. SAUNDERS S.E.R.C. Daresbury Laboratory Daresbury Warrington WA4 4AD Received 12th December 1984 The present status of the application of Hartree-Fock theory to periodic systems using gaussian basis sets is reviewed. The most important technical aspects are discussed in a language close to that of the molecular theoretician so that the method can be seen as similar to a molecular calculation but with a number of extra features both simplifying and complicating. The level of agreement currently achieved between theory and experiment for such data as lattice parameters X-ray structure factors and electron momentum distributions (Compton profiles) is illustrated.The considerable potential for exploiting a number of the ideas stemming from periodic systems in the context of calculations on large molecules is highlighted. The application of Hartree-Fock theory to periodic systems has a number of advances in the past five years. Broadly speaking the subject’s status is similar to that of molecular Hartree-Fock theory in the early seventies. Thus it is now possible to deal with unit cells containing up to perhaps six heavy atoms using a double- zeta-quality gaussian basis set with polarization functions in favourable cases. The level of numerical accuracy is perhaps a little less than could be obtained in molecular calculations the total energy/cell being usually only correct to 10-100 phartree with concomitant errors in other properties.The largest error is not however numerical but stems from the use of finite basis sets and the Hartree-Fock approximation itself as in molecular calculations. The vast experience now gained in molecular work has had a decisive influence on the progress of the periodic-system work. For example the selection of a suitable crystal basis set is almost invariably rooted in atomic or molecular experience whilst many of the highly advanced molecular methods have migrated to periodic systems rather successfully. The present paper will argue that the time is now ripe for migration of some solid-state techniques for application in large molecules. The present work attempts to cover three aspects (1) those elements of the theory representing an extension of the molecular case (2) applications of the theory and (3) approximation schemes which have been well tested in the rather convenient context of periodic systems and which will almost certainly prove useful in calculations on very large molecular systems.THEORETICAL AND METHODOLOGICAL CONSIDERATIONS We assume a working knowledge of molecular Hartree-Fock theory and delineate the most important differences for periodic systems. Complication 1. The periodic system is of infinite extent. This leads to overlap and Fock matrices whose dimension is in principle infinite and similarly infinite summations arise when evaluating coulomb and exchange operators. This circumstance forces us to consider 79 HARTREE-FOCK CALCULATIONS FOR PERIODIC SYSTEMS Approximation 1.If the differential overlap between two gaussian basis functions is less than a pre-set tolerance then all matrix elements involving that overlap distribution are assumed zero. Typically a threshold of lop6is used in a high-quality calculation. However a problem remains; if one overlap distribution is deemed finite an infinite number may be generated by periodic repetition. This problem is solved by SimpZi$cation 1. Suppose we single out a particular unit cell call it the ‘zero cell’. Then it is sufficient to consider only overlap distributions involving at least one basis function within the zero cell since all others (where neither basis function is in the zero cell) are simply translational replicants with identical overlap and Fock matrix elements.This fact coupled with our use of approximation 1 gives rise to only a finite set of matrix elements which need to be computed to define the problem. Complication 2. Although the number of overlap distributions that need be considered explicitly has been reduced to a finite number matrix elements involving these distributions are made up of an infinite number of terms. For example the nuclear attraction operator consists of all atoms in the crystal while the coulomb and exchange matrix elements consist of double infinite summations over the entire crystal basis set. This difficulty leads to two further approximations Approximation2. Exchange operator matrix elements are evaluated using1 truncation criteria such that overlap distributions sufficiently remote from the target overlap distribution are simply ignored.Each target overlap distribution can thus be thought of as having an exchange ‘sphere’ surrounding it the size of the sphere being dependent on the nature of the overlap distribution smaller for inner-shell orbitals and larger for valence orbitals. The radius of the exchange sphere is normally surprisingly small it being rare to have to include interactions beyond fourth-nearest . neighbours. In fact the rate of convergence is considerably better by perhaps an order of magnitude than would be expected from a consideration of the largest- valued exchange integrals which have been excluded from the sum presumably because exchange interactions tend to have random phase.This quasi-localness of the exchange operator is not too surprising in view of the success of local density functional the~ry.~ Approximation 3. A coulomb matrix element represents the energy of interaction of a given target overlap distribution with the total electron distribution of the crystal while the nuclear attraction matrix elements represent the energy of interaction with respect to all nuclei in the system. The convergence of such sums is notoriously slow and recipes based on a straightforward truncation will be doomed to prohibitively long explicit summations if adequate accuracy is to be maintained. Satisfactory procedures are normally based upon a two-step analysis (a) the series is first rearranged by grouping together subsets of the total charge distribution; the purpose is to improve the convergence properties of the rearranged series whose summation index corresponds as a rule to the direct lattice vectors ordered in increasing length and (b)in the evaluation of the new series a decreasing level of accuracy is adopted with increasing index.Extreme care must be taken because connection problems may arise when passing from a zone of given accuracy to the next less precise one. In our method2 the total crystal electron density is partitioned into subsets by means of the well known Mulliken population analysis such that each shell of atomic orbitals has associated with it an electron distribution. Each shell distribution is asymptotically exponentially decaying with distance squared (assuming a gaussian basis) and so may be thought of as effectively having a finite extent.Coulomb interactions between an V. R.SAUNDERS overlap distribution and a shell charge are then evaluated according to three different schemes. (a)At short range where the penetration of the two charge distributions is appreciable the interaction is evaluated exactly as a linear combination of two-electron integrals. (b) At intermediate range each shell charge is replaced by a linear combination of point monopoles dipoles etc. (up to 1 = 6 is feasible with our current code) the poles being located at the centroid of the shell. The overlap distribution is allowed to electrostatically interact with these poles and the method can be seen as a practical outlet for the procedure of Stone5 for analysing charge distributions.Of course for the method to be economical it is necessary to be able to evaluate high-order multipole integrals and even more importantly electrostatic field integrals (involving operators of the form &,/R2z+1) efficiently. Fortunately very convenient recursive methods are now available.2 As the range increases the order of pole included may be gradually decreased until a distance is reached where only monopoles need be included so defining the border of the ‘quantum zone’. (c) Beyond this all shell charges are grouped together to form atomic charges (including the nucleus) and are said to belong to the Madelung zone. The interaction of the overlap distribution with the Madelung charges is evaluated by Ewald’s method.6 At the centroid of each overlap distribution the Madelung potential and its first gradients (with respect to distance) are evaluated whence the matrix element can be computed as a linear combination of overlap and dipole moment integrals.It is unnecessary to proceed to quadratic or higher terms as these have been shown“ to be of quite negligible importance. It may be argued that the use of the Mulliken partition of the charge is arbitrary. A ‘weighted’ Mulliken scheme has also been tried and shown2 to give negligibly different results. This invariance to the partitioning scheme is only realized if the multipole expansions are carried out to sufficiently high order. The point is that any reasonable partitioning of charge will give satisfactory results if used consistently and carefully.However the use of different partitioning algorithms in different zones will give rise to very severe zone connection problems. Simplijication 2. Consider the short-range Coulomb interaction of an overlap distribution with shell charges. Within these summations will be found sub-sets of the form J(i,j) =Z (ijlkl;m) P(k1) m where the sum is over a number of distributions (kl;m)which differ only by translation vector m,and hence have the same density matrix element P(k1). During integral evaluation it is convenient to store such quantities summed over the m index (which will normally embrace between 10 and 100 terms). The advantages are twofold (a) much less file space is consumed by the presummed integrals1 and (b)methods2.are available for evaluating such sums much more efficiently than evaluating each integral separately and then summing. Similar advantage may also be gained when evaluating sums over exchange or electrostatic field operators. Complication 3. Although the Fock and overlap matrix are now defined in terms of a finite number of distinct elements they are nonetheless infinite in dimension because of periodic replication. Fortunately this final problem may be overcome by Simplijication 3. Symmetry adaption of the basis set using the infinite number of translational symmetry operators gives rise to the Bloch basis. The Fock and overlap matrices may be evaluated at a finite set of points in k space (typically 10-150 such points are used) each k-space sub-block of the matrices being of the dimension of HARTREE-FOCK CALCULATIONS FOR PERIODIC SYSTEMS the number of basis functions per unit cell.After diagonalization the real-space density matrix may be found by numerical integration8 over the k-space density matrices. For further discussion see ref. (1) and (9). In summary the procedure outlined above is based almost entirely on real-space considerations and the severity of the approximations decay with distance in a manner which is highly regular and predictable and so can be adequately controlled. The general format of the calculation looks rather like a typical ab initio large-cluster calculation with the addition of a Madelung potential and periodic boundary conditions.The latter are crucial in allowing a much larger quantum zone than could be tolerated in a finite-cluster study since often 200-1000 atoms are involved. This means that periodic systems form the ideal test-bed for approximation schemes and in our opinion would have been invented for this purpose even if they had not occurred in nature. It is well to contrast the present method with that of Harris and Monkhorst,lO who proposed to work more directly in reciprocal space using Fourier transforms of Slater-type orbitals. Sums over real lattice points are replaced by sums over the reciprocal lattice. The principal difficulty of this method appears to be one of defining adequately accurate truncation criteria particularly when one is denied the advantage of the insight provided by being able to think in real-space terms.We may also mention here that there seems to be little advantage for the present method in the use of pseudo-potentials for the treatment of inner-shell orbitals. For example in a minimal-basis calculation on silicon (diamond structure) only ca. $ of the integrals involve inner shells. The reason is that the classification adopted in the approximation schemes delineated above makes a severe distinction between compact inner and the more diffuse valence shells. SOME ILLUSTRATIVE RESULTS We first choose to illustrate the importance of including Madelung effects particularly in ionic systems. Consider MgO in a double-zeta basis the molecule gives rise to a Mulliken oxygen charge of 0.8e.If we now perform a crystal calculation using a quantum zone of radius 25 bohr (this is about as large as can usually be afforded) but ignoring the Madelung potential the oxygen charge increases to 1.2e. After inclusion of the Madelung potential the charge separation increases to 2e and the system is fully ionic. In the calculation which did not incorporate the Madelung potential there is other evidence that all is not well for example certain Fock matrix elements which should be equal (because of periodicity) are found to differ by 0.3 hartree; this symmetry is restored (to six decimal places) upon inclusion of the Madelung field. A detailed study3 of f.c.c. lithium hydride (containing the very polarizable H anion) reveals even greater surprises because here a choice of two different unit cells produced totally different (but wrong) results if the Madelung potential was neglected (for example the total energy/cell differed by 0.2 hartree) whilst agreement was restored (to six- decimal places) if the Madelung potential was incorporated to gradient terms.If the Madelung potential included only zero-order terms only marginal improvement occurs (the energy per cell differed by 0.05 hartree). These calculations gave an optimized nearest-neighbour distance of 3.876 bohr (experimental value 3.858 bohr) and a binding energy of 0.338 hartree (experimental value 0.346 hartree). The calculated bulk modulus was within the large experimental uncertainty but it is not yet clear that the calculations are carried out with sufficient precision to allow a definitive result for such second-gradient properties.The calculated electron-density distribution which corresponded closely to a fully ionic model gave excellent V. R. SAUNDERS agreement with respect to experimental X-ray structure factor and Compton profile data. In contrast calculationsll using a simple ‘molecular simulated crystal model ’ which corresponded to a more covalent situation gave equally good agreement for the structure-factor data but rather poor agreement with the Compton profiles indicating that Compton-profile data more powerfully discriminate between theoretical models of electronic structure. A similar study’ of Li,N (space roup p6/mmm) gave rise to optimized lattice parameters a = 3.61 A and c = 3.84 x (experimental values 3.655 and 3.874 A); Li,N layers alternate with pure Li layers in the c direction with the Li atoms of the Li,N layers arranged in graphitic hexagonal structure with the N atoms at the centre of the hexagons while the Li atoms of the pure Li layers are atop the N atoms.The calculated electronic structure was found to correspond closely to the fully ionic model with the N atom being triply negatively charged and again excellent agreement with both X-ray structure-factor data and Compton-profile data was observed. Simple model calculations give good agreement for structure-factor data13 but poor agreement for the Compton profiles,14 reinforcing our view that structure-factor data are not as useful as Compton-profile data in discriminating between theoretical models.The calculated binding energy was 0.196 hartree (experimental value 0.42 hartree) an error which can be almost entirely attributed to the difference in correlation energy of the N atom and the N3- ion which may be estimated from Clementi’s tables15 to be 0.22 hartree. The calculations were also interesting because a major attempt at basis-set optimization (particularly for the N anion) was attempted in the crystalline environment. We turn now to consider briefly some results for the (SN) linear polymer,16 which is an almost planar chain and has a repeat unit of the form with the NS bonds parallel to the fibre axis being slightly shorter than those approximately perpendicular (1.593 and 1.628 A respectively).The bond angles x and y are 106.2 and 119.9’ respectively. The calculations used a minimal basis with and without d orbitals on the S atom the basis being as used by Palmer and Findlay17 in finite-cluster calculations on the same system. The polymer-optimized d-orbital exponent was 0.32 not too different from that found in the (SN) molecular system.18 The contribution of the d orbitals is found to be even more important in stabilizing the polymer than in the case of the (SN) precursor to such an extent that the chain remains unstable if d orbitals are excluded the d functions providing 0.24 hartree per cell with a sulphur atom d-orbital population of 0.29e and a total sulphur positive charge of 0.46. The d orbitals play a considerable role in reducing the ionicity of the S-N bonds (the S atom charge is +0.63 in minimal basis).Indeed in minimal basis the S atom can only achieve the necessary hypervalency by becoming ionic a tendency much reduced when the extra mechanism for bonding provided by the d functions is allowed for. The convergence of the finite-cluster calculation^^^ appears to be rather good if one considers the total energy per cell. However this does not mean that other aspects of the electronic structure are so rapidly convergent; even at the centre of the largest cluster (SN),, the electron distribution exhibits considerable differences when compared with the polymer result. One must be careful in taking even very large-scale cluster calculations as representative of the infinite system.HARTREE-FOCK CALCULATIONS FOR PERIODIC SYSTEMS CONCLUSIONS A brief and not too technical description of the techniques used in the Hartree-Fock theory of periodic systems has been given and test calculations used to illustrate the level of agreement with experiment. There are a number of hints here for large-molecule Hartree-Fock theory. (a) As a corollary of the observed slow convergence of Madelung potential effects finite-cluster calculations where one replaces a periodic system of ions with a finite array of point charges should be carried out with extreme caution taking care to use a suitable unit cell as the basis for replication. For example consider the linear chain ...A2+B2-A2+B2-A2fB2-A2+B2-... where the most obvious unit is either A2+B2-or B2-A2+.However finite-cluster calculations based on these cells will give very slow convergence because of the large dipole moment of the repeat unit.A better choice is either A+B2-A+or B-A2+B- which have zero dipole. Indeed it is possible to construct quadrupole-free repeat units (even in the general case) which would be even better. (b) The rapid convergence of exchange operator summations observed in periodic systems is almost certainly mirrored in large molecules. (c) The treatment of long-range coulomb forces in large molecules using a distributed multipole representation of the electron distribution should be tried in large molecules. (d) The procedure for efficiently summing over large numbers of interactions (simplification 2 above) is not applicable in orthodox SCF calculations because the density matrix elements vary through the SCF cycles.However in the direct SCF method,lg where one recomputes the molecular integrals at each cycle the technique should apply ufortiori. C. Pisani and R. Dovesi Znt. J. Quunturn Chem. 1980 17 501. R. Dovesi C. Pisani C. Roetti and V. R. Saunders. Phys. Rrr. B 1983 28 5781. R. Dovesi C. Ermondi E. Ferrero C. Pisani and C. Roetti Phjls. Rev. B 1984 29 3591. J. C. Slater Phys. Rev. 1951 81 385 R. Gaspar Acta Phys. Hung. 1954 3 263; W. Kohn and L. J. Sham Phys. Rep. A 1965 140 1133. A. J. Stone Chem. Phys. Lett. 1981 83 233. ti P. P. Ewald Ann. Phys. (N.Y.) 1921 64 253; M. Tosi in Solid State Physics ed. F. Seitz and D.Turnbull (Academic Press New York 1964) vol. 16 p. 1; F. E. Harris in Theoretical Chemistry Adtlances and Perspectires ed. H. Eyring and D. Henderson (Academic Press London 1975) vol. l. p. 147. V. R. Saunders in Methods in Computational Molecular Physics ed. G. H. F. Diercksen and S. Wilson (Reidel Dordrecht 1983) p. 1. * G. Gilat and L. J. Raubenheimer Phys. Rev. 1966 144 390; G. Gilat Phys. Rev. B 1973 7 891 H. J. Monkhorst and J. D. Pack Phys. Reti. B 1976 13 5188; D. J. Chadi Phys. Rev. B 1977 16 1746; J. D. Pack and H. J. Monkhorst Phys. Rez;. B 1977 16 1748. J. M. Andre L. Gouverneur and G. Leary Znt. J. Quanfurn Chem.. 1967. 1 427; 451 ;J. M. Andre J. Chem. Phys. 1969 50 1536. lo F. E. Harris and H. J. Monkhorst Phys. Rec. Lett. 1969 23 1026.l1 B. I. Ramirez W. R. McIntire and R. L. Matcha J. Chem. Phys. 1977 66,373. l2 R. Dovesi. C. Pisani F. Ricca C. Roetti and V. R. Saunders Phys. Rec. B 1984 30,972. l3 H. Schulz and K. H. Schwarz Acta Crystallogr. Sect. A 1978 34,999. P. Pattison and J. R. Schneider Acta Crystallogr. Sect. A 1980 36,390. E. Clementi. J. Chem. Phys. 1963 38,2248. l6 R. Dovesi C. Pisani C. Roetti and V. R. Saunders J. Chem. Phys. 1984 81 2839. M. H. Palmer and R. H. Findlay J. Mol. Struct. (Theochem) 1983 92 373. R. C. Haddon S. R. Wassermann F. Wudl and G. R. J. Williams J. Am. Chem. SOC.,1980 102 6687. J. Almhof K. Faegri and K. Korsell J. Comput. Chem. 1982 3 385.