High-Precision Calculation of Hartree-Fock Energy of Crystals M. J. GILLAN, 1,2,3 D. ALFÈ, 1,2,3,4 S. DE GIRONCOLI, 5 F. R. MANBY 6 1 Materials Simulation Laboratory, UCL, Gower Street, London WC1E 6BT, United Kingdom 2 London Centre for Nanotechnology, UCL, Gordon Street, London WC1H 0AH, United Kingdom 3 Department of Physics and Astronomy, UCL, Gower Street, London WC1E 6BT, United Kingdom 4 Department of Earth Sciences, UCL, Gower Street, London WC1E 6BT, United Kingdom 5 INFM-DEMOCRITOS and SISSA, via Beirut 2, 34014 Trieste, Italy 6 School of Chemistry, University of Bristol, Bristol BS8 1TS, United Kingdom Received 24 January 2008; Accepted 16 April 2008 DOI 10.1002/jcc.21033 Published online 5 June 2008 in Wiley InterScience (www.interscience.wiley.com). Abstract: When using quantum chemistry techniques to calculate the energetics of bulk crystals, there is a need to calculate the Hartree-Fock (HF) energy of the crystal at the basis-set limit. We describe a strategy for achieving this, which exploits the fact that the HF energy of crystals can now be calculated using pseudopotentials and plane-wave basis sets, an approach that permits basis-set convergence to arbitrary precision. The errors due to the use of pseudopotentials are then computed from the difference of all-electron and pseudopotential total energies of atomic clusters, extrapolated to the bulk-crystal limit. The strategy is tested for the case of the LiH crystal, and it is shown that the HF cohesive energy can be converged with respect to all technical parameters to a precision approaching 0.1 mE h per atom. This cohesive energy and the resulting HF value of the equilibrium lattice parameter are compared with literature values obtained using Gaussian basis sets. © 2008 Wiley Periodicals, Inc. J Comput Chem 29: 2098–2106, 2008 Key words: Hartree-Fock; crystal; solid; high accuracy; plane wave; pseudopotential Introduction Density functional theory (DFT) has dominated computational condensed-matter science for many years. 1 The reasons are well known: its mild scaling with number of atoms makes it possible to treat large complex systems; it allows basis-set convergence to be achieved to any desired tolerance; and atomic forces can be calculated at almost no extra cost, so that high-temperature dynam- ical and thermodynamic properties can be calculated by molecular dynamics simulation. However, the quantitative accuracy of DFT is often inadequate (see e.g. ref. 2), and there is currently no known way of systematically improving the description of electron correla- tion within DFT. These problems with DFT have stimulated efforts to apply wavefunction-based methods to condensed matter. These efforts date back many years, and include, for example, the incre- mental methods developed by Stoll, Fulde, and co-workers. 3–7 In the incremental approach, the total energy is divided into Hartree- Fock and correlation parts, and the correlation energy is decomposed according to a many-body expansion into single-atom, atom-pair, atom-triplet, etc... terms. The approach has been successfully used to calculate the cohesive energy, equilibrium lattice parameter, and other properties of bulk crystals. 7 Recently, there have also been efforts to implement periodic versions of correlated quantum chemistry techniques. 8–11 We have recently reported 12 an alternative way of using wavefunction-based techniques to calculate the cohesive energy. This approach was designed particularly for ionic crystals such as LiH and MgO, though we believe it is more widely applicable. The total energy is separated, as usual, into Hartree-Fock and cor- relation energies, but the latter is then further separated into the correlation energy of the appropriate number of MX molecules (M = cation, X = anion), and the remainder, which we refer to as the “correlation residual”. This residual is then extracted by analysing its systematic variation over a large hierarchy of free or embedded clusters. (Details are given in our published paper. 12 ) Taking LiH as an example, we showed that this approach is capable of giving the total cohesive energy within ∼ 1mE h (∼30 meV) per Li-H pair, and comparison with experiment confirmed this. However, it has become clear from our more recent work (unpublished) that in order to attain this level of precision, one of the major challenges is to calculate the Correspondence to: M. J. Gillan; e-mail: m.gillan@ucl.ac.uk © 2008 Wiley Periodicals, Inc.