Accurate Hartree-Fock energy of extended systems using large Gaussian basis sets

Accurate Hartree-Fock energy of extended systems using large Gaussian basis sets

Joachim Paier Department of Chemistry, Rice University, Houston, Texas 7705, USA    Cristian V. Diaconu Department of Chemistry, Rice University, Houston, Texas 7705, USA    Gustavo E. Scuseria Department of Chemistry, Rice University, Houston, Texas 7705, USA    Manuel Guidon Physical Chemistry Institute, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland    Joost VandeVondele Physical Chemistry Institute, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland    Jürg Hutter Physical Chemistry Institute, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland
July 9, 2019

Calculating highly accurate thermochemical properties of condensed matter via wave function-based approaches (such as e.g. Hartree-Fock or hybrid functionals) has recently attracted much interest. We here present two strategies providing accurate Hartree-Fock 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 Hartree-Fock exchange using a truncated Coulomb operator and an extrapolation toward the full-range Hartree-Fock limit of a Padé fit to a series of short-range screened Hartree-Fock calculations. These two techniques agreed to significant precision. We also present the Hartree-Fock cohesive energy of LiH (converged to within sub-meV) at the experimental equilibrium volume as well as the Hartree-Fock equilibrium lattice constant and bulk modulus.

61.50.Ah, 71.15.-m, 71.15.Nc, 71.15.Ap

I Introduction

The high accuracy/cost ratio of Kohn-Sham density functional theoryKohn and Sham (1965); Dreizler and Gross (1990) (KS-DFT) has been exhaustively demonstrated in the literature. In its early days, KS-DFT 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 exchange-correlation (XC) functional and the introduction of nonlocal Hartree-Fock (HF) exchange in hybrid functionalsBecke (1992, 1993) paved the way for reasonably accurate applications to molecules as well.

Within the framework of KS-DFT 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 high-temperature 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 KS-DFT lies in the accuracy of the applied XC functional.

Discussing examples for some shortcomings of KS-DFT, it is well known that standard local and semilocal approximations to the XC functional do not yield accurate results for quasiparticle band-gaps 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 KS-DFT have stimulated some DFT groups to use wave function-based methods to benchmark or correct DFT results. Note that another large community of researchers directly apply wave function-based 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 function-based and semilocal DFT methods. We refer the reader to a recent review Janesko et al. (2009) of so-called screened hybrid functionals, as e.g. the Heyd-Scuseria-ErnzerhofHeyd 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 long-range 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 quantum-chemistry community and the steady increase of computers’ efficiency induced a drive to conceive and implement even more involved wave function-based techniques, as e.g. local second-order Møller-Plesset 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 function-based 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 post-HF 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 so-called 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 all-electron 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 augmented-plane 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 sub-meV 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 full-range HFX based on the short-range (SR) HFX implementation Izmaylov et al. (2006) in the GAUSSIAN suite of programs as well as the method applied for the extrapolation of the SR-HFX 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 quasi-optimal basis sets (pc-[0-4]) that rapidly converge to the HF and DFT basis set limit. The pc-3 basis set gives atomization energies with a mean error smaller than 1 kJ/mol. For H and Li the pc-3 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 non-standard 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 densely-packed 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 Hartree-Fock 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 pc-4-like 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
Table 1: Details for the adopted basis sets for the compositions 8s3p2d1f/6s3p2d1f and 13s6p2d1f/11s5p2d1f of Hydrogen and Lithium respectively. Shown are angular momentum, Gaussian exponent and corresponding contraction coefficients.

ii.2 Extrapolation of SR-HFX to full range (GAUSSIAN)

All Gaussian calculations presented in this work are based on a very efficient implementation of the SR-HFX energy exploiting a distance based screening protocol.Izmaylov et al. (2006) Using local basis functions it is convenient to express the HFX energy for closed-shell as


where are density matrix elements and


are the four-center 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 square-root of the band gap of the system of question.Kohn (1959); Cloizeaux (1964); Kohn (1995); Goedecker (1999) This is the basic motivation behind SR-HF 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


where the long-range (LR) and short-range (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 SR-HF 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 full-range 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 self-consistent-field (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


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 least-squares 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 full-range 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 SR-HF 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 Å (seven-points-fit).

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 real-space 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 long-range tail of the Coulomb interaction (e.g. in the limit, see Eq. 3, in long-range 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 SR-HFX 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 long-range HF calculations are likely to be computationally prohibitive except for high-exponent 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 long-range 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 SR-HFX 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 Hartree-Fock 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 in-core 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


This operator was suggested by Spencer and Alavi Spencer and Alavi (2008) to obtain rapid convergence for the Hartree-Fock 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 Hartree-Fock 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


where are the wave functions at the point. In the Gaussian basis set employed, the exchange energy per cell is thus obtained from


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 full-range 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 least-squares fits obtained using rational polynomials up to order seven. In addition, correlation coefficient , root-mean-square-deviation (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 (1-r)(a) RMSD(b) RMSD%(c) E(HF) [a.u.]
6.7e-05 0.0185955 0.3654035 -32.3526129
5.7e-08 0.0005666 0.0111347 -32.2472782
5.5e-07 0.0018244 0.0358500 -32.2521383
2.0e-10 0.0000364 0.0007156 -32.2585728
7.9e-10 0.0000759 0.0014909 -32.2588628
1.5e-09 0.0001109 0.0021803 -32.2576031
6.2e-15 0.0000002 0.0000046 -32.2581712
  • r: correlation coefficient (see text for details).

  • RMSD: root mean square deviation.

  • RMSD%: normalized root mean square deviation.

Table 2: Results for the Padé fits to the 31 SR-HF energies [a.u.] of LiH at experimental lattice constant (4.084 Å) for a cell containing four LiH ion pairs. The first column shows the order of the Padé polynomials representing them by the order of the polynomial of the numerator and denominator, respectively (in squared brackets). The extrapolated total HF energy for the cell is given in Hartree atomic units.
#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
Table 3: Validation of the [7/7] Padé fit. The first column gives the number of data points (i.e. SR-HF energies) included for a [7/7] Padé fit in order to predict the SR-HF energy corresponding to the value given in the second column. Deviations between calculated and fitted SR-HF energies are given in the fifth and sixth column, respectively.
Figure 1: (Color online) Convergence of the SR-HF energy [a.u.] of a LiH unit cell containing four LiH ion pairs at experimental lattice constant (4.084 Å) with decreasing screening parameter involved in the short-range Coulomb interaction. Gaussian results for each are represented by crosses. The line shows the [7/7] Padé fit to the numerical data (see text for details). The inset gives SR-HF energies for a.u. as well as the extrapolated value for in Hartree atomic units.

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.

Figure 2: (Color online) Murnaghan equation of state (red/gray line) for LiH obtained using the HF approximation. Each of the seven points corresponds to the extrapolated least-squares fit of 31 screened HF energies to a Padé approximant of order [7/7] (see text for details).
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.

Table 4: Summary of total HF energies per formula unit, HF cohesive energies, equilibrium lattice constants and bulk moduli of LiH obtained using GAUSSIAN and CP2K. For comparison purpose, results found in the literature are included.

iii.2 HF lattice constant and bulk modulus with GAUSSIAN

Fig. 2 shows the seven-points-fit 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 SR-HF 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

Table 5: Results obtained with CP2K and the truncated Coulomb operator for unit cells that are a multiple of the cubic unit cell (4.084 Å). The columns show the size of the unit cell, the range of the truncated Coulomb operator (), the Hartree-Fock energy per four LiH ion pairs, the H atom energy, the Li atom energy, and the cohesive energy (), respectively.

Iv Conclusions

The Hartree-Fock energy of solid LiH has been calculated using large Gaussian basis sets. Two different approaches, extrapolation of a Padé fit to a series of SR-HFX 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 post-Hartree-Fock methods.


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 Hartree-Fock 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/2007-2013) under grant agreement no RI-211528. This work was supported by DOE DE-FG02-04ER15523, and the Welch Foundation C-0036.


  • 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. Zicovich-Wilson, 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,
  • (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).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description