Accurate HartreeFock energy of extended systems using large Gaussian basis sets
Abstract
Calculating highly accurate thermochemical properties of condensed matter via wave functionbased approaches (such as e.g. HartreeFock or hybrid functionals) has recently attracted much interest. We here present two strategies providing accurate HartreeFock energies for solid LiH in a large Gaussian basis set and applying periodic boundary conditions. The total energies were obtained using two different approaches, namely a supercell evaluation of HartreeFock exchange using a truncated Coulomb operator and an extrapolation toward the fullrange HartreeFock limit of a Padé fit to a series of shortrange screened HartreeFock calculations. These two techniques agreed to significant precision. We also present the HartreeFock cohesive energy of LiH (converged to within submeV) at the experimental equilibrium volume as well as the HartreeFock equilibrium lattice constant and bulk modulus.
pacs:
61.50.Ah, 71.15.m, 71.15.Nc, 71.15.ApI Introduction
The high accuracy/cost ratio of KohnSham density functional theoryKohn and Sham (1965); Dreizler and Gross (1990) (KSDFT) has been exhaustively demonstrated in the literature. In its early days, KSDFT using the local density approximationKohn and Sham (1965) was almost exclusively applied by the solid state community. However, the advent of generalized gradient approximations (GGAs, see e.g. Refs. Langreth and Perdew, 1980; Perdew et al., 1992, 1996) to the exchangecorrelation (XC) functional and the introduction of nonlocal HartreeFock (HF) exchange in hybrid functionalsBecke (1992, 1993) paved the way for reasonably accurate applications to molecules as well.
Within the framework of KSDFT it is relatively easy to achieve basis set convergence, and atomic forces can be calculated at little extra computational cost. This is of paramount importance in the calculation of hightemperature dynamical and thermodynamic properties by molecular dynamics simulations. In particular, DFT statistical mechanics for both bulk materials and for surface processes is routinely feasible (see e.g. Ref. Marx and Hutter, 2009, and references therein). The principle limitation of KSDFT lies in the accuracy of the applied XC functional.
Discussing examples for some shortcomings of KSDFT, it is well known that standard local and semilocal approximations to the XC functional do not yield accurate results for quasiparticle bandgaps of semiconductors and insulators. Generally, they do not predict the correct adsorption sites and adsorption energies of molecules on metallic surfaces (for details see e.g. Ref. Stroppa and Kresse, 2008, and references therein). Today’s DFT practitioner is confronted with these shortcomings when choosing an XC functional for a specific application. Each contemporary density functional has its relative merits but at the same time drawbacks, which might impede finding an appropriate XC functional. Different density functionals have different merits and demerits, an unsatisfactory situation. These inadequacies of semilocal KSDFT have stimulated some DFT groups to use wave functionbased methods to benchmark or correct DFT results. Note that another large community of researchers directly apply wave functionbased techniques to materials science problems (see Ref. Paulus, 2006, and references therein).
Dramatic improvements for many properties of molecules as well as solids can be achieved by mixing a fraction of nonlocal HF exchange (HFX) to the remaining part of semilocal DFT exchange. Since these functionals do not only depend on the electron density alone, but also on the KS single particle wave functions, i.e. the orbitals, they are called hybrid functionals. Therefore, these hybrid functionals can be seen as “mixed” wave functionbased and semilocal DFT methods. We refer the reader to a recent review Janesko et al. (2009) of socalled screened hybrid functionals, as e.g. the HeydScuseriaErnzerhofHeyd et al. (2003, 2006) (HSE) functional, which was proposed to extend the successes of hybrid functionals into condensed matter, by avoiding the problematic effects of longrange HFX (see Ref. Janesko et al., 2009, and references therein).
Besides the successes of screened HFX applied to condensed matter, the numerous methodological and algorithmic developments in the quantumchemistry community and the steady increase of computers’ efficiency induced a drive to conceive and implement even more involved wave functionbased techniques, as e.g. local secondorder MøllerPlesset perturbation theory (MP2) Pisani et al. (2005, 2008); Casassa et al. (2007), (resolution of the identity) atomic orbital Laplace transformed MP2 Ayala and Scuseria (2000); Ayala et al. (2001); Izmaylov and Scuseria (2008) and canonical MP2 Sun and Bartlett (1996); Marsman et al. (2009) for (infinitely) extended systems of various dimensionality and applied basis functions. Furthermore, recent reports in the literature on ab initio molecular dynamics on condensed matterGuidon et al. (2008) employing the HSE screened hybrid functional illustrate that, depending on implementation details, basis set and system, wave functionbased techniques are also applicable to statistical mechanics calculations. These successful applications of wave function methods to large systems show that they are able to tackle materials science problems with possibly much better accuracy than conventional density functionals.
Recently published HF and postHF calculations on crystalline LiH have attracted much interest in the solid state communityManby et al. (2006); Gillan et al. (2008); Nolan et al. (2009a). These calculations represent a benchmark in terms of eliminating as many inaccuracies as possible while attempting to converge toward the socalled HF limit. The approach in question employs calculations on a hierarchical series of cluster models,Manby et al. (2006); Nolan et al. (2009b) exploiting strengths and weaknesses of plane wave pseudopotentials as well as local Gaussian basis sets. Accurate evaluation of the total HF energy, as well as cohesive energy in the HF approximation employing exclusively Gaussian basis sets is desirable to bypass errors incurred by the pseudopotential approximation. Admittedly, creating an allelectron Gaussian basis set, which describes the crystal as well as the isolated atoms equally well, is challenging. Referring to the arguments of Gillan et al., Gillan et al. (2008) it is in general difficult to provide rigorous estimates how far the applied basis set is from the HF limit. However, it is reasonable to question the need for reaching the HF limit for particular materials properties, which is substantiated in the present work.
We compare total HF energies of solid LiH using two different codes employing Gaussian basis functions: (i) the Gaussian and augmentedplane wave (GAPW) Lippert et al. (1999) code CP2K/Quickstep VandeVondele et al. (2005); The CP2K developers group () and (ii) a developmental version of the GAUSSIAN suite of programs Frisch et al. (). We show that the cohesive energy of the crystal is converged to within submeV accuracy in our given large Gaussian basis set (see Tab. 1). Computational and methodological details are presented in Sec. II. Results for cohesive energies, theoretical lattice constant as well as bulk modulus are in Sec. III. Conclusions are drawn in Sec. IV.
Ii Computational details
In the following sections, we describe important computational details, such as the Gaussian basis set, the evaluation of fullrange HFX based on the shortrange (SR) HFX implementation Izmaylov et al. (2006) in the GAUSSIAN suite of programs as well as the method applied for the extrapolation of the SRHFX energy to the full range limit based on Padé approximants. Furthermore, implementation details on the direct evaluation of HFX via the CP2K/Quickstep code are presented.
ii.1 Basis set
The basis set used for this calculation has been specifically constructed for the current purpose, which is an accurate but computationally feasible HF calculation on bulk LiH. The basis constructed here is similar to the polarization consistent (pc) basis sets derived by Jensen.Jensen (2001, 2002, 2007) Jensen introduced a sequence of quasioptimal basis sets (pc[04]) that rapidly converge to the HF and DFT basis set limit. The pc3 basis set gives atomization energies with a mean error smaller than 1 kJ/mol. For H and Li the pc3 basis set has a composition 9s4p2d1f/5s4p2d1f and 14s6p2d1f/6s3p2d1f respectively, while we adopt 8s3p2d1f/6s3p2d1f and 13s6p2d1f/11s5p2d1f. However, the primitives of the basis employed here are nonstandard and optimized for the present calculations.
In a first step, we have removed primitive Gaussians with exponents smaller than , since diffuse basis functions are technically troublesome. Diffuse functions, which are needed to describe density tails in atoms or molecules, are not needed in the bulk of denselypacked solids with large band gaps as the case of LiH. Indeed, we exploit the fact that the basis functions on the lattice sites are available for the expansion of any orbitals, be it the crystal orbitals in the bulk or the atomic orbitals of the isolated atoms. This basis is thus only suited for atomic or surface calculations if ghost basis functions are left in the regular lattice positions to appropriately describe the aforementioned tails of the electron density.
In a second step, all but the core exponents have been optimized by minimizing the energy of bulk LiH subject to a restraint on the condition number of the overlap matrix. This procedure is similar to the one employed for the molecularly optimized basis sets described in Ref. VandeVondele and Hutter, 2007. In CP2K, density functionals that do not include HartreeFock exchange can be computed in a highly efficient manner, and in order to make this procedure computationally efficient, such a semilocal density functional (B88 Becke (1988)) has been employed in the optimization process. The resulting basis is well conditioned, the condition number of the overlap matrix is for bulk LiH. We have estimated the accuracy of the optimized basis by comparing to pc4like basis sets, which for this system are only feasible with local DFT, and estimate the total energy of bulk LiH (per unit of LiH) to be well within 0.001 a.u. of the basis set limit, while the basis set error on the cohesive energy is likely smaller than 0.1% (0.0001 a.u.). The details of this optimized basis set are summarized in Tab. 1.
species  l  exponent  coefficient 

H  s  0.27463675e02  1.00000000 
s  0.68559258e01  1.00000000  
s  0.17679972e01  1.00000000  
s  0.51181842e00  1.00000000  
s  0.20167548e00  1.00000000  
s  0.30797000e04  0.00023473  
0.46152000e03  0.00182450  
0.10506000e03  0.00959330  
p  0.21240865e01  1.00000000  
p  0.10736812e01  1.00000000  
p  0.56838662e00  1.00000000  
d  0.92833840e00  1.00000000  
d  0.49583000e00  1.00000000  
f  0.12073480e01  1.00000000  
Li  s  0.13360341e04  1.00000000 
s  0.44429982e03  1.00000000  
s  0.14779702e03  1.00000000  
s  0.49209451e02  1.00000000  
s  0.16428957e02  1.00000000  
s  0.55293994e01  1.00000000  
s  0.19052824e01  1.00000000  
s  0.70025874e00  1.00000000  
s  0.29958682e00  1.00000000  
s  0.16636288e00  1.00000000  
s  0.70681000e05  0.00000544  
0.13594000e05  0.00003328  
0.31004000e04  0.00019175  
p  0.15709110e01  1.00000000  
p  0.74875864e00  1.00000000  
p  0.38614089e00  1.00000000  
p  0.22620503e00  1.00000000  
p  0.28500000e02  0.00036754  
0.66400000e01  0.00322359  
d  0.77920820e00  1.00000000  
d  0.40789925e00  1.00000000  
f  0.73706300e00  1.00000000 
ii.2 Extrapolation of SRHFX to full range (GAUSSIAN)
All Gaussian calculations presented in this work are based on a very efficient implementation of the SRHFX energy exploiting a distance based screening protocol.Izmaylov et al. (2006) Using local basis functions it is convenient to express the HFX energy for closedshell as
(1) 
where are density matrix elements and
(2) 
are the fourcenter electron repulsion integrals (ERIs), represented in an atomic orbital basis. The applied interaction potential is usually equal to the Coulomb kernel .
For large gap systems, it has been shown that local single particle wave functions as well as the corresponding density matrix decay like for large , where is proportional to , the squareroot of the band gap of the system of question.Kohn (1959); Cloizeaux (1964); Kohn (1995); Goedecker (1999) This is the basic motivation behind SRHF as e.g. used in the successful HSE hybrid functional.Heyd et al. (2003, 2006) HSE is based on a screened Coulomb interaction splitting the conventional Coulomb kernel, , into
(3) 
where the longrange (LR) and shortrange (SR) parts of the interaction are described by the computationally convenient error function and its complement, respectively. The parameter in Eq. 3 determines the extent of the range separation of the Coulomb interaction.
In view of the relatively large HF band gap of LiH (10.8 eV)Baroni et al. (1995) and the fast decay of the density matrix, we will calculate the total HF energy by doing a series of SRHF calculations at different , and extrapolating to . As corroborated by numerical results shown in Sec. III, such an extrapolation of the screened HF energies of the crystal to the fullrange HF limit in the specified basis set is numerically robust and reliable.
All calculations are based on a locally modified development version of the GAUSSIAN electronic structure program.Frisch et al. () Hence, the total energies presented in Sec. III do not include any DFT contributions. Only Hartree and screened HFX energies are evaluated. The RMS convergence criterion for the density matrix in the selfconsistentfield (SCF) iteration was set to a.u., which implies an energy convergence no worse than at least a.u. (GAUSSIAN keyword: SCF=Tight). Furthermore, a mesh of points was used, which is equivalent to 6912 points and thus all calculations are sufficiently converged with respect to points. The large band gap of LiH in the HF approximation (see above) substantially helps converging the point integration.
Following ideas found in the literature,Ayala et al. (1999); Iyengar et al. (2000) we apply Padé approximants of various orders to the obtained series of screened HF energies. The actual form of the Padé approximants are the rational polynomials
(4) 
Eq. 4 represents the general expression of a Padé approximant of order . Throughout this work only diagonal rational polynomials are applied,Ayala et al. (1999); Press et al. (1992) which means that the order of the polynomial in the numerator equals the order of the polynomial in the denominator. Note that the number of parameters to be fitted is in the case of diagonal polynomials. This is the minimum number of data points, which must be included in the leastsquares fit.
For all extrapolations employed in the present work, has been chosen to lie in the interval . In order to put a higher weight to the area near to fullrange HF we decided to increment by up to and increment by up to a value that amounts to . For the remaining interval of larger values, the screening parameter was incremented by . As a consequence, each fit is based on 31 data points, representing pairs of the screening parameter and the corresponding SRHF energy.
The HF equilibrium lattice constant and bulk modulus have been obtained by fitting the volume dependence of the static lattice energy to the Murnaghan equation of state.Murnaghan (1944) The points were chosen in order to cover a range of 3% around the supposed equilibrium lattice constant of 4.108 Å (sevenpointsfit).
ii.3 HFX and periodic boundary conditions using Gaussian basis functions
In hybrid functionals, which incorporate a fraction of nonlocal HFX (Eq. 1), the decay of the Hamiltonian matrix elements (see Sec. II of Ref. Kudin and Scuseria, 2000) with distance is determined by two factors: (i) the decay behavior of the density matrix, , and (ii) the decay behavior of the ERIs. For metallic as well as insulating systems, a screened Coulomb interaction accelerates the convergence of the ERIs in realspace drastically, i.e., the number of replica cells needed for convergence is substantially decreased (see Refs. Heyd and Scuseria, 2004 and Paier et al., 2006). Fock exchange calculations involving the longrange tail of the Coulomb interaction (e.g. in the limit, see Eq. 3, in longrange corrected hybrids or in global hybrids), both the density matrix and the ERIs influence the convergence of the HFX energy.
It is a matter of fact, that due to the algebraic structure of Eq. 1, contributions to the HFX energy can be significant even far from the central cell, precluding an early truncation of the lattice sum (see Eq. 2.4 in Ref. Kudin and Scuseria, 2000). Small exponent basis functions involved in the calculation of the density matrix become important factors determining the computational workload. Calculations under periodic boundary conditions (PBC) involving SRHFX with a reasonably large value for the screening parameter (Eq. 3) are tractable for moderately diffuse Gaussian basis functions, i.e. minimal exponent equals . Conventional HF or longrange HF calculations are likely to be computationally prohibitive except for highexponent Gaussian basis sets. The relatively large and diffuse basis set used in this work (see Tab. 1) prevents calculating the HFX at or close to , i.e. the longrange limit for this particular system in the given basis. As shown by the results presented in Sec. III, a numerically stable fit to a sufficiently large series of SRHFX calculations is practicable to calculate an accurate estimate for the HF energy of extended (insulating) systems using large Gaussian basis sets. In summary, is not practical whereas a 31 point extrapolation works very well.
ii.4 Direct calculation of HFX (CP2K)
The focus of CP2K is the simulation of complex systems, with a variety of methods. Recently, the capability to perform first principles molecular dynamics simulation with density functionals including a fraction of HartreeFock exchange has been implemented and demonstrated for condensed phase systems containing a few hundred atoms Guidon et al. (2008). With this goal in mind, the implementation is massively parallel, focuses on incore calculations, uses the point only, and does not exploit molecular or crystal symmetries. The implementation was based on a minimum image (MI) convention Tymczak et al. (2005) and employed a standard 1/r Coulomb operator. The current implementation, which will be described in detail elsewhere Guidon et al. (), goes beyond the minimum image convention, and instead employs a truncated Coulomb operator which is defined as
(5) 
This operator was suggested by Spencer and Alavi Spencer and Alavi (2008) to obtain rapid convergence for the HartreeFock energy with respect to the point sampling of the exchange energy in periodic systems. Note that the use of the truncated Coulomb operator implies that the exchange energy is unconditionally convergent for all points. Furthermore, since exchange in insulators is effective on shorter range compared to the electrostatic interaction, results converge exponentially to the HartreeFock limit as is increased. In line with the results presented in Ref. Spencer and Alavi, 2008, we find that for a cubic cell with edge and , converged results of the exchange energy can be obtained using the point only. Of course, this requires that the computational cell is sufficiently large so that the point approximation is acceptable, which in turn requires that the extent of the maximally localized Wannier functions is smaller than . Consequently, the exchange energy computed in CP2K is defined as
(6) 
where are the wave functions at the point. In the Gaussian basis set employed, the exchange energy per cell is thus obtained from
(7) 
where are the indices of the basis functions in the central cell, and run over all image cells. Due to the rapid decay of the basis functions, sums over image cells and converge quickly. The sum over converges quickly and unconditionally for our choice of . Further technical details, including how to compute efficiently and accurately the required four center integrals will be presented elsewhere Guidon et al. ().
Iii Results
iii.1 HF energy of LiH at experimental volume using Padé approximants
Fig. 1 depicts the obtained series of 31 data points of screened HF energies for various values of (see Sec. II) calculated at the experimental lattice parameter. The series of calculated energies clearly converges to a certain limit with decreasing . At this point we remind the reader that for the limit the SR Coulomb kernel given in Eq. 3 approaches the fullrange operator. As shown in Fig. 1 and outlined in Sec. II, the density of data points increases significantly toward the limit. Tab. 2 presents results for several leastsquares fits obtained using rational polynomials up to order seven. In addition, correlation coefficient , rootmeansquaredeviation (RMSD) as well as relative RMSD, which is normalized to the range of observed data, i.e. calculated energies, are shown. Since is very close to , we decided to present in Tab. 2 , where a value of zero means perfect agreement between calculated data points and fit.
Apparently, the RMSD as well as relative RMSD values decrease with increasing order of the rational polynomial applied to the fit. Rational polynomials of order eight or beyond (not shown in Tab. 2) lead to unstable fits and the goodness of the fit deteriorates. According to Tab. 2 the optimal order of the Padé approximant is , which was used for all extrapolations employed in this work.
Fit  (1r)^{(a)}  RMSD^{(b)}  RMSD%^{(c)}  E(HF) [a.u.] 

6.7e05  0.0185955  0.3654035  32.3526129  
5.7e08  0.0005666  0.0111347  32.2472782  
5.5e07  0.0018244  0.0358500  32.2521383  
2.0e10  0.0000364  0.0007156  32.2585728  
7.9e10  0.0000759  0.0014909  32.2588628  
1.5e09  0.0001109  0.0021803  32.2576031  
6.2e15  0.0000002  0.0000046  32.2581712 

r: correlation coefficient (see text for details).

RMSD: root mean square deviation.

RMSD%: normalized root mean square deviation.
#dp  E [a.u.]  E [a.u.]  Error  %Error  

24  0.070  31.630779  31.630965  0.0001817  0.000574 
25  0.065  31.675185  31.675183  0.0000010  0.000003 
26  0.060  31.719533  31.719533  0.0000003  0.000001 
27  0.055  31.763998  31.763998  0.0000003  0.000001 
28  0.050  31.808571  31.808571  0.0000003  0.000001 
29  0.045  31.853244  31.853244  0.0000003  0.000001 
30  0.040  31.898007  31.898007  0.0000004  0.000001 
As a next step we had to validate the polynomial, since it is well known, that algorithms for interpolation are straightforward, whereas for extrapolation care must be taken. A plausible strategy is simply the prediction of energies for a certain value of not included in the fit. Table 3 shows a series of predictions for screened HF energies for a series of ’s starting from a.u.. The corresponding screened HF energy has been estimated based on a fit using 24 data points, where . As can be seen from Tab. 3, it is remarkable that the resulting error is only one order of magnitude larger than the applied SCF convergence criterion (see Sec. II). The error for the predicted energies is practically converged after inclusion of only one further data point and amounts to a.u.. Hence, the error incurred by the fit to the Padé approximant is much lower than the convergence threshold in the SCF procedure. By virtue of the aforementioned validations it is safe to give the total HF energy for a unit cell containing four LiH ion pairs at experimental lattice constant (4.084 Å) with a precision of five decimals in Hartree atomic units, which amounts to a.u.. The total HF energy per formula unit at experimental lattice constant is given in Tab. 4 and compared with the HF energy obtained using CP2K. Both values agree excellently to significant precision.
E(HF) [E]  [mE]  [Å]  B [GPa]  

GAUSSIAN  8.064543  4.105  32.34  
CP2K  8.064545  131.949  
CRYSTAL  129.14  4.121  28.3  
CRYSTAL  130.16  
VASP  131.7  
Gillan et al.  131.95  
Gillan et al.  131.99  4.108  32.05 

calculated at experimental lattice constant (4.084 Å).

Ref. Casassa et al., 2007.

Ref. Dovesi et al., 1984.

Ref. Marsman et al., 2009.

Ref. Gillan et al., 2008.
iii.2 HF lattice constant and bulk modulus with GAUSSIAN
Fig. 2 shows the sevenpointsfit of the obtained Padé extrapolated HF energies to the Murnaghan equation of state as outlined in Sec. II.2. The RMSD value of this fit amounts to . The resulting HF equilibrium lattice constant of LiH equals 4.105 Å and is in excellent agreement with the result obtained by Gillan et al. (see Tab. 4). The corresponding bulk modulus of LiH amounts to 32.34 GPa, which is again in very good agreement (0.9% deviation) with the results obtained by aforementioned workers. Note that bulk moduli are quite sensitive to the equilibrium volume at which they are evaluated and overall good indicators for the quality of the underlying energies at the various volumes.
iii.3 Total and cohesive energy at experimental volume with CP2K
Total energies have been computed for systems of increasing system size by explicitly repeating the cubic unit cell periodically in three dimensions. The largest cell employed is a repetition of the basic cubic cell, and contains exactly 1000 atoms. For this system, Gaussian basis functions are used for the expansion of the molecular orbitals, which makes this a computationally demanding simulation. With increasing system size, we also increase the range of the truncated Coulomb operator, in steps of 2 Å up to a maximum of 10 Å (see Tab. 5). The point approximation therefore converges quickly (exponentially) to the HF limit of this system. We thus obtain from a direct calculation, without extrapolation, an accurate estimate of the total energy per unit cell of approximately 32.258179 a.u.. The finite size error on this result is estimated to be smaller than 50 E. Furthermore, this number is in excellent agreement with the Padéextrapolated SRHF results (32.258171 a.u., Tab. 2), and thus provides numerical evidence for the quality of both approaches. Calculating the HF energy of the H atom and the Li atom with the current basis set, in periodic boundary conditions and retaining the basis functions of all other atoms in the unit cell, we can obtain a consistent estimate of the cohesive energy. In our approach, due to the fact that unrestricted calculations are needed for the atoms, these calculations are even more demanding than the bulk, and have only been performed up to a repetition of the basis unit cell. Our best estimate for the cohesive energy, obtained from just three calculations (bulk LiH, and the atoms Li, H) without extrapolation, is 131.949 mE. Also here, the finite size error is estimated to be smaller than 50 E. This number is in excellent agreement with the best estimate obtained by Gillan et al. Gillan et al. (2008) 131.95 mE.
[Å]  E(HF)[a.u.]  H[a.u.]^{(a)}  Li[a.u.]^{(b)}  [a.u.]  

4.0  32.244609  0.499957  7.428493  0.132702  
6.0  32.256844  0.499974  7.432137  0.132100  
8.0  32.258022  0.499974  7.432582  0.131949  
10.0  32.258179  N/A  N/A  N/A 

basis set limit 0.500000

basis set limit 7.432727
Iv Conclusions
The HartreeFock energy of solid LiH has been calculated using large Gaussian basis sets. Two different approaches, extrapolation of a Padé fit to a series of SRHFX calculations and direct calculation using a truncated Coulomb operator, have been found to yield total energies that agree to better than 0.1 mE. Calculations of the cohesive energy, the equilibrium lattice constant and the bulk modulus agree with the best estimates available in literature. These results show that robust and accurate calculations with nearly converged Gaussian basis sets have now become possible in the condensed phase at least for large band gap systems. These results will contribute to the growing usefulness of hybrid density functionals for condensed phase applications and opens, for these systems, the way to accurate calculations based on postHartreeFock methods.
Acknowledgments
J.V. acknowledges a fruitful collaboration with John Levesque (Cray Inc.). His extensive testing on the Cray XT5 (Jaguar) at the Oak Ridge National Laboratory (ORNL) enabled the massively parallel implementation of the HartreeFock code in CP2K. Additional computer resources have been provided by the Swiss National Supercomputing Center (CSCS), and as part of the PRACE project by the Finnish IT Center for Science (CSC). The PRACE project receives funding from the EU’s Seventh Framework Programme (FP7/20072013) under grant agreement no RI211528. This work was supported by DOE DEFG0204ER15523, and the Welch Foundation C0036.
References
 Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
 Dreizler and Gross (1990) R. M. Dreizler and E. K. U. Gross, Density Functional Theory (Springer Verlag, Berlin, 1990).
 Langreth and Perdew (1980) D. C. Langreth and J. P. Perdew, Phys. Rev. B 21, 5469 (1980).
 Perdew et al. (1992) J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992).
 Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
 Becke (1992) A. D. Becke, J. Chem. Phys. 98, 1372 (1992).
 Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
 Marx and Hutter (2009) D. Marx and J. Hutter, Ab Initio Molecular Dynamics (Cambridge University Press, Cambridge, 2009).
 Stroppa and Kresse (2008) A. Stroppa and G. Kresse, New J. Phys. 10, 063020 (2008).
 Paulus (2006) B. Paulus, Phys. Rep. 428, 1 (2006).
 Janesko et al. (2009) B. G. Janesko, T. M. Henderson, and G. E. Scuseria, Phys. Chem. Chem. Phys. 11, 443 (2009).
 Heyd et al. (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
 Heyd et al. (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124, 219906 (2006), erratum.
 Pisani et al. (2005) C. Pisani, M. Busso, G. Capecchi, S. Casassa, R. Dovesi, L. Maschio, C. ZicovichWilson, and M. Schütz, J. Chem. Phys. 122, 094113 (2005).
 Pisani et al. (2008) C. Pisani, L. Maschio, S. Casassa, M. Halo, M. Schütz, and D. Usyvat, J. Comput. Chem. 29, 2113 (2008).
 Casassa et al. (2007) S. Casassa, M. Halo, L. Maschio, C. Roetti, and C. Pisani, Theor. Chem. Acc. 117, 781 (2007).
 Ayala and Scuseria (2000) P. Y. Ayala and G. E. Scuseria, J. Comput. Chem. 21, 1524 (2000).
 Ayala et al. (2001) P. Y. Ayala, K. N. Kudin, and G. E. Scuseria, J. Chem. Phys. 115, 9698 (2001).
 Izmaylov and Scuseria (2008) A. Izmaylov and G. Scuseria, Phys. Chem. Chem. Phys. 10, 3421 (2008).
 Sun and Bartlett (1996) J. Sun and R. J. Bartlett, J. Chem. Phys. 104, 8553 (1996).
 Marsman et al. (2009) M. Marsman, A. Grüneis, J. Paier, and G. Kresse, J. Chem. Phys. 130, 184103 (2009).
 Guidon et al. (2008) M. Guidon, F. Schiffmann, J. Hutter, and J. VandeVondele, J. Chem. Phys. 128, 214104 (2008).
 Manby et al. (2006) F. Manby, D. Alfè, and M. Gillan, Phys. Chem. Chem. Phys. 8, 5178 (2006).
 Gillan et al. (2008) M. Gillan, D. Alfè, S. Gironcoli, and F. Manby, J. Comp. Chem. 29, 2098 (2008).
 Nolan et al. (2009a) S. J. Nolan, M. J. Gillan, D. Alfé, N. L. Allan, and F. R. Manby, accepted in Phys. Rev. B (2009a).
 Nolan et al. (2009b) S. J. Nolan, P. J. Bygrave, N. L. Allan, and F. R. Manby, J. Phys.: Condensed Matter (2009b), in press.
 Lippert et al. (1999) G. Lippert, J. Hutter, and M. Parrinello, Theor. Chem. Acc. 103, 124 (1999).
 VandeVondele et al. (2005) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005).
 (29) The CP2K developers group, http://cp2k.berlios.de/(2008).
 (30) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, et al., Gaussian development version Revision H.1, Gaussian Inc. Wallingford CT 2009.
 Izmaylov et al. (2006) A. Izmaylov, G. Scuseria, and M. Frisch, J. Chem. Phys. 125, 104103 (2006).
 Jensen (2001) F. Jensen, J. Chem. Phys. 115, 9113 (2001).
 Jensen (2002) F. Jensen, J. Chem. Phys. 116, 7372 (2002).
 Jensen (2007) F. Jensen, J. Phys. Chem. A 111, 11198 (2007).
 VandeVondele and Hutter (2007) J. VandeVondele and J. Hutter, J. Chem. Phys. 127, 114105 (2007).
 Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
 Kohn (1959) W. Kohn, Phys. Rev. 115, 809 (1959).
 Cloizeaux (1964) J. Cloizeaux, Phys. Rev. 135, A685 (1964).
 Kohn (1995) W. Kohn, Int. J. Quant. Chem. 56, 229 (1995).
 Goedecker (1999) S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
 Baroni et al. (1995) S. Baroni, G. P. Parravicini, and G. Pezzica, Phys. Rev. B 32, 4077 (1995).
 Ayala et al. (1999) P. Y. Ayala, G. E. Scuseria, and A. Savin, Chem. Phys. Lett. 307, 227 (1999).
 Iyengar et al. (2000) S. S. Iyengar, G. E. Scuseria, and A. Savin, Int. J. Quant. Chem. 79, 222 (2000).
 Press et al. (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and W. P. Flannery, Numerical Recipes in Fortran (Cambridge University Press, Cambridge, 1992), 2nd ed.
 Murnaghan (1944) F. D. Murnaghan, Proc. Natl. Acad. Sci. U.S.A. 30, 244 (1944).
 Kudin and Scuseria (2000) K. Kudin and G. Scuseria, Phys. Rev. B 61, 16440 (2000).
 Heyd and Scuseria (2004) J. Heyd and G. Scuseria, J. Chem. Phys. 121, 1187 (2004).
 Paier et al. (2006) J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber, and J. G. Ángyán, J. Chem. Phys. 124, 154709 (2006).
 Tymczak et al. (2005) C. J. Tymczak, V. Weber, E. Schwegler, and M. Challacombe, J. Chem. Phys. 122, 124105 (2005).
 (50) M. Guidon, J. Hutter, and J. VandeVondele, in preparation.
 Spencer and Alavi (2008) J. Spencer and A. Alavi, Phys. Rev. B 77, 193110 (2008).
 Dovesi et al. (1984) R. Dovesi, C. Ermondi, E. Ferrero, C. Pisani, and C. Roetti, Phys. Rev. B 29, 3591 (1984).