Accurate thermochemistry of covalent and ionic solids from spin-component-scaled MP2

In recent years, wavefunction-based methods, popular in the quantum chemistry community, have begun to be regularly applied to condensed-phase systems. In particular, coupled-cluster theory with single and double excitations (CCSD) and with perturbative triple excitations [CCSD(T)]—the so-called gold standard of molecular quantum chemistry—yield promising results.^1–11 However, due to their high scaling with system size N, i.e., O(N⁶) for CCSD and O(N⁷) for CCSD(T), their applicability to solids with complex unit cells is limited, and obtaining results at the thermodynamic limit and the complete basis set limit is challenging.

As a lower cost alternative to coupled-cluster theory, second-order Møller–Plesset perturbation theory (MP2) is widely used in molecular quantum chemistry and has been increasingly applied to periodic systems.^{1,2,5,6,12–23} Over a decade ago, Grüneis et al.¹⁵ demonstrated that, for simple covalent and ionic solids, ground-state structural and electronic properties predicted by MP2 are quite good and better than those predicted by DFT with the popular Perdew–Burke–Ernzerhof (PBE) exchange correlation functional.²⁴ However, it is natural to consider improvements to MP2, the limitations of which are well known in molecular quantum chemistry.

At the same cost, spin-component-scaled (SCS) MP2,²⁵ which semiempirically scales the same-spin and opposite-spin components of the correlation energy, has been shown to significantly outperform MP2 in regard to many molecular properties.^25–33 Scaled opposite-spin (SOS) MP2,³⁴ which retains and scales only the opposite-spin component of the correlation energy, shows comparably good performance^27,32–34 and can be performed with a reduced O(N⁴) scaling using density fitting and a Laplace transform of the energy denominators. Both SCS- and SOS-MP2 commonly deliver an accuracy that is better than that of CCSD and comparable to that of CCSD(T), suggesting that they could be especially promising methods for complex, insulating solids.

Here, we test periodic SCS/SOS-MP2 for the calculation of the lattice constant, bulk modulus, and cohesive energy of 12 three-dimensional covalent and ionic semiconductors and insulators, paying careful attention to the thermodynamic limit and complete basis set limit. The layout of this paper is as follows. In Sec. II, we briefly review the formalism of SCS/SOS-MP2 and provide details of our periodic MP2 calculations, basis sets, and density fitting. We also present the convergence of the HF energy and the components of the MP2 correlation energy to the thermodynamic limit and complete basis set limit using carbon diamond as an example. In Sec. III, we present the accuracy of the calculated properties for 12 solids, including a scan over the spin-scaling parameters. The conclusions are presented in Sec. IV.

The SCS-MP2 correlation energy can be split into two components,

where the spin-scaling parameters c_os = c_ss = 1 for traditional MP2. The rescaling of opposite-spin and same-spin contributions to the correlation energy can be motivated or derived in several ways,³² but it essentially originates from their different physical character: At the mean-field level of a single Slater determinant, same-spin electrons experience Pauli repulsion but opposite-spin electrons do not, suggesting that they should be treated differently at the post-mean-field level. Based on empirical fitting to calculations on molecules, the conventional SCS-MP2 parameters are (c_os, c_ss) = (1.2, 0.33)²⁵ and the conventional SOS-MP2 parameter is c_os = 1.3 (with c_ss = 0),³⁴ although reoptimization has been explored in a variety of contexts.^30,35,36

With periodic boundary conditions and N_k crystal momenta k sampled from the Brillouin zone, the spin components of the correlation energy are

where

Electron repulsion integrals are given in Mulliken (11|22) notation and, as usual, i, j refer to occupied orbitals and a, b refer to unoccupied orbitals in the periodic Hartree–Fock (HF) determinant, which we assume to be spin-restricted. The primed summation indicates conservation of crystal momentum, k_i + k_j − k_a − k_b = G, where G is a reciprocal lattice vector.

All calculations were performed with PySCF.^37,38 The Brillouin zone was sampled with a uniform Monkhorst–Pack mesh of N_k k-points that includes the Γ point. Core electrons were replaced with correlation-consistent effective-core potentials (ccECPs) and we used their corresponding correlation-consistent cc-pVXZ basis sets^39,40 (XZ). For Li and Mg, we used the large-core pseudopotentials ([He] core and [Ne] core, respectively). Gaussian density fitting with an even-tempered auxiliary basis set was used for evaluating the electron repulsion integrals.⁴¹ 

Due to the compact arrangement of atoms in crystals, linear dependencies are a common occurrence when using large atom-centered basis sets.^42–44 To avoid associated numerical problems, we removed the most diffuse primitive Gaussians from the orbital basis set of Mg (exponents $<0.05$⁠) and Li (exponents $<0.1$⁠), except in LiCl. From the even-tempered auxiliary basis set for all systems, we removed diffuse Gaussians with exponents $<0.2$⁠. Testing indicated that these modifications do not affect our results.

To assess our ability to reliably access the thermodynamic limit (TDL) and complete basis set (CBS) limit, we consider carbon diamond at its experimental crystal volume as an illustrative example. (The same behavior is seen in other materials; in  Appendix A, we present the same study of MgO, whose ionic character is quite different than the covalent character of diamond.) In Fig. 1, we show convergence to the TDL using the TZ basis set and k-point meshes of N_k = 3³–6³ for HF and N_k = 3³–6³ for MP2. To address the divergent exchange term in periodic HF, we used a Madelung constant correction,^45–47 which yields total energies and orbital energies that converge to the TDL as $Nk−1$⁠. Subsequent MP2 correlation energies converge at the same rate because the orbital pair densities appearing in the electron repulsion integrals are chargeless by orthogonality of the molecular orbitals. We extrapolate the HF energy according to the form $EHF(Nk)=EHF(∞)+aNk−1+bNk−2$ and the MP2 correlation energy according to the form $E(2)(Nk)=E(2)(∞)+aNk−1$ (numerically, we found that the inclusion of sub-leading corrections to the HF energy was most appropriate, as shown in Fig. 1, although nearly identical results are obtained with b = 0). We estimate by eye that these extrapolations give results in the TDL that are accurate to about 0.01 eV. The same-spin and opposite-spin correlation energies exhibit finite-size errors that are almost identical in magnitude, at least for this material.

In Fig. 2, we show the HF (a) and MP2 (b) equation of state (EOS) of diamond using N_k = 4³ and the DZ and TZ basis sets. For the HF EOS, we compare to results obtained using a plane-wave (PW) basis set, with a kinetic energy cutoff chosen to achieve convergence to better than 1 meV/atom. We see excellent agreement between the PW result and the TZ result (better than 0.02 eV/atom), indicating that the latter is near the CBS limit. For the MP2 energy, which is more expensive to evaluate in a large PW basis, we compare to a calculation with virtual orbitals obtained from an approximate PW resolution of the QZ basis set,⁴⁸ denoted as PW + QZ. We see more sensitivity to the basis set as expected. Comparing the TZ and PW + QZ results, the lattice constant, bulk modulus, and cohesive energy differ by 0.003 Å, 3.5 GPa, and 0.15 eV, respectively (see below for details about the calculation of these properties). A X⁻³ CBS extrapolation of the DZ and TZ correlation energies (where X = 2, 3) gives an EOS in excellent agreement with the PW + QZ result; the deviations are 0.0001 Å, 1.8 GPa, and 0.009 eV.

To summarize, in all subsequent production calculations, we use the TZ HF energy and the DZ/TZ CBS-extrapolated MP2 correlation energy and extrapolate to the TDL using results obtained with N_k = 4³–6³ (HF) and N_k = 5³–6³ (MP2). With this protocol, our largest calculations (TZ with N_k = 6³) require about 8–24 h with 12–24 cores. Because of the polynomial scaling, analogous calculations with N_k = 5³ are significantly cheaper, requiring about 1 h or less.

We studied 12 semiconductors and insulators, with diamond (C and Si), zincblende (SiC, BN, BP, AlN, and AlP), and rock salt (MgO, MgS, LiH, LiF, and LiCl) crystal structures. These solids were chosen based on their basis set availability, relatively simple crystal structures, and nonzero bandgap, as required for the application of MP2-based methods. We predicted the lattice constant a, bulk modulus B, and cohesive energy E_coh of all solids using HF, MP2, and SCS/SOS-MP2 by calculating the total energy for ten different volumes in the range of ±10% of the experimental volume. We then fitted the energies to an equation of state given by a third-order Birch–Murnaghan form. The cohesive energy is calculated at the theoretically predicted equilibrium volume and is defined with respect to the energy of isolated atoms. Atomic energies of open-shell atoms were calculated with unrestricted HF and MP2 and the basis set superposition error was accounted for by adding crystalline basis functions.

As an initial test of our methods and implementation, we have compared our calculated HF and MP2 results to those of Grüneis et al.,¹⁵ who studied 11 of the 12 solids studied here at the same levels of theory (all except MgS). At the HF level, the mean absolute deviations in the lattice parameter, bulk modulus, and cohesive energy are 0.010 Å, 2.2 GPa, and 0.09 eV, respectively. At the MP2 level, the same deviations are 0.017 Å, 3.6 GPa, and 0.05 eV. A few of the biggest discrepancies are for the lattice constant of LiF (0.039 Å with HF and 0.036 Å with MP2); the bulk moduli of diamond (5.6 GPa with HF and 16.8 GPa with MP2), SiC (5.1 GPa with MP2), and MgO (5 GPa with HF and MP2); and the cohesive energies of SiC (0.37 eV with HF) and AlP (0.25 eV with HF). Nonetheless, these mean absolute deviations provide a rough estimate of the precision of our HF and MP2 calculations due to pseudopotential, basis set, and finite-size errors, complementing the detailed study of carbon diamond in Sec. II. Agreement with experiment to higher accuracy should be viewed with caution. Our predicted HF and MP2 properties for each appear material are presented in  Appendix B, which can be directly compared to those of Grüneis et al.,¹⁵ for a more fine-grained analysis.

To evaluate the performance of HF- and MP2-based methods, we will compare to the experimental values compiled recently in Ref. 50, in which the performance of six DFT functionals was evaluated on 64 solids. We compare to experimental results that have been corrected to remove zero-point energy contributions calculated at the PBE level in the same work. For AlN and LiH, which are absent in Ref. 50, we use values from Ref. 51, which were identically corrected for zero-point energy in that work. For most of the shared materials, these two references give very similar properties. Three large outliers are the bulk moduli of BN (differing by 21.7 GPa), BP (8.5 GPa), and MgO (3.2 GPa); otherwise, lattice constants, bulk moduli, and cohesive energies agree to within 0.003 Å, 1 GPa, and 0.02 eV, respectively.

A summary of our results is given in Table I, which shows the mean absolute error (MAE) of various methods for each property compared to these experimental values; material-specific predictions for all properties are given in  Appendix B. HF predictions are unsurprisingly poor. MP2 predictions are a significant improvement, exhibiting MAEs of 0.020 Å, 5.2 GPa, and 0.23 eV for the lattice parameter, bulk modulus, and cohesive energy, respectively. SCS-MP2 predictions are a further improvement, exhibiting MAEs of 0.015 Å, 3.8 GPa, and 0.06 eV for the lattice parameter, bulk modulus, and cohesive energy, respectively; the cohesive energy is especially improved. The SOS-MP2 predictions are similarly good, exhibiting MAEs of 0.017 Å, 3.7 GPa, and 0.08 eV for the lattice parameter, bulk modulus, and cohesive energy, respectively. We therefore conclude that SCS-MP2 and SOS-MP2, using the most common spin-scaling parameters, perform remarkably well for the prediction of the structural and energetic properties of three-dimensional covalent and ionic semiconductors and insulators.

The performance of MP2-based methods can be broken down for each material, as shown in Fig. 3. Some of the apparent outliers at the MP2 level of theory are the lattice constant of LiCl, MgO, and BP; the bulk modulus of C and MgO; and the cohesive energy of C, BN, and BP. Generally, we see that spin scaling yields an increase in the lattice constant, a decrease in the bulk modulus, and a decrease in the cohesive energy. The latter, in particular, corrects the tendency of MP2 to overestimate the cohesive energy as it does for all materials except LiH.

We finally assess the extent to which further improvement is possible by a reoptimization of the spin-scaling parameters. In Fig. 4, we plot the SCS-MP2 MAE with respect to the experimental value for the lattice constant, bulk modulus, and cohesive energy as a function of the spin-scaling parameters. In addition to the MP2 and molecular SCS-MP2 parameters, we indicate the optimal SCS and SOS parameters for each property, and we see that different properties favor different combinations of spin-scaling coefficients. Moreover, many of the landscape valleys are long and narrow, indicating a range of parameters that deliver comparable performance for a given property. To identify a globally optimal set of parameters, we calculate a dimensionless average as follows:

where ⟨|ΔO|⟩ indicates the MAE of property O and ⟨·⟩_MP2 indicates the MP2 value. The weighting is such that the value of the cost function is the average error with respect to that of MP2; of course, different ways of averaging will give different optimal parameters. This average error is plotted in Fig. 4(d) and identifies the optimal SCS parameters (c_os, c_ss) = (1.24, 0.24) and optimal SOS parameter c_os = 1.36, which are quite similar to the standard values determined for molecules. The optimal parameters of the average are very similar to those of the cohesive energy, in part because the latter’s optimal parameters balance the errors of the lattice constant and the bulk modulus. Moreover, the average error of the standard SCS and SOS parameters is only marginally higher: The optimal values give 0.568 (SCS) and 0.585 (SOS) while the standard values give 0.570 (SCS) and 0.646 (SOS). In other words, spin-component scaling in any of these nearly optimal forms yields MAEs that are about 60% those of MP2.

We have studied the performance of periodic SCS- and SOS-MP2 in the TDL and CBS limit, finding excellent agreement with experimental values. As shown in Table I, the performance is significantly better than that of popular functionals for solid-state calculations, including PBE,²⁴ PBEsol,⁵² and SCAN.⁵³ The performance of SOS-MP2 is only marginally worse than that of SCS-MP2, indicating excellent promise as an especially affordable technique. Our work motivates the exploration of other SCS methods, including those for excited states, which could be valuable given the mixed success with which low-order perturbation theories can calculate bandgaps of solids.^15,54,55

Although we have carefully addressed finite-size and basis set errors, there may be residual pseudopotential errors, which are harder to remove systematically. For example, preliminary testing indicates that MP2 calculations on MgO, which here exhibited atypically large errors in the lattice constant and bulk modulus, are significantly improved through the use of a small-core pseudopotential with a [He] core. A careful study of pseudopotentials and core-correlation effects is in progress and will be presented elsewhere.

The conclusions of this work are limited to the class of materials studied, i.e., three-dimensional covalent and ionic insulators, and future work must assess the extension to surfaces, layered materials, and other weakly bound solids, such as molecular crystals. An early periodic MP2 study by Del Ben et al.¹⁹ observed good SCS-MP2 and double hybrid performance for molecular crystals, and a recent report from our group²³ found that SCS-MP2 yields a very accurate prediction of the cohesive energy of crystalline benzene, while SOS-MP2 yields a non-negligible underestimation. The present work has demonstrated good performance over a relatively wide range of spin-scaling parameters, indicating that additional systems or properties can be incorporated in the identification of optimal parameters for condensed-phase materials. Perhaps most importantly, we have found that the standard scaling parameters, which are known to deliver good performance for molecular chemistry, are almost optimal for solid-state properties. This transferability suggests that surface chemistry on insulators should be accurately described by SCS- and SOS-MP2.

This work was supported by the National Science Foundation under Grant No. CHE-1848369 (T.G.) and Grant No. OAC-1931321 (H.-Z.Y.). We acknowledge the computing resources provided by Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant No. 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract No. C090171, both awarded April 15, 2010. The Flatiron Institute is a division of the Simons Foundation.

The authors have no conflicts to disclose.

T.G. and X.W. contributed equally to this work.

Tamar Goldzak: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (lead); Writing – review & editing (supporting). Xiao Wang: Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (lead); Writing – review & editing (supporting). Hong-Zhou Ye: Investigation (supporting); Methodology (supporting); Software (supporting); Validation (lead); Writing – review & editing (supporting). Timothy C. Berkelbach: Conceptualization (equal); Data curation (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (equal); Project administration (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

In Figs. 5 and 6, we show the finite-size convergence and basis set convergence of the ionic crystal MgO, confirming the same qualitative behavior as seen in Figs. 1 and 2 for diamond.

In Tables II–IV, we provide the material-specific predictions at the HF, MP2, SCS-MP2, and SOS-MP2 levels of theory, along with their mean absolute error (MAE) and mean absolute relative error (MARE) compared to experimental values.