A many-electron perturbation theory study of the hexagonal boron nitride bilayer system*

Abstract
In this article we explore methods to reduce the computational cost in many-electron wave function expansions including explicit correlation and compact one-electron basis sets for the virtual orbitals. These methods are applied to the calculation of the interlayer binding energy of the h-BN bilayer system. We summarize the optimized interlayer distances as well as their binding energies for various stacking faults on different levels of theory including second-order Møller-Plesset perturbation theory and the random phase approximation. Furthermore, we investigate the asymptotic behavior of the binding energy at large interlayer separation and find that it decays as D-4 in agreement with theoretical predictions, where D is the interlayer distance.



Introduction
Kohn-Sham density functional theory in the framework of approximate exchange and correlation energy density functionals is undeniably the workhorse method in ab initio calculations of real materials [1,2]. However, the lack of systematic improvability in currently available density functionals has motivated the exploration and further development of other computationally more complex manyelectron theories that inherently integrate more sophisticated electronic correlation effects such as van der Waals interactions.
Quantum chemical wave function based theories aim at a systematically improvable and very accurate approximation to the true many-electron wave function. Once obtained, these wave functions are used to calculate the physical and chemical properties of interest for atoms, molecules and solids. However, the underlying computational complexity for most of these theories is large and plagued by the high dimensionality of the Hilbert space in which the wave functions need to be optimized and approximated. Nonetheless, a large fraction of this Hilbert space is not needed and can be "down-folded", allowing for more compact wave function representations. Examples of such low-rank wave function approximations include for instance coupled cluster theory, matrix product states or stochastic full configuration interaction [3][4][5][6].
These "compactifications" have made the calculation of very accurate many-electron correlation energies for real materials possible.
This work employs a number of recently developed tools for efficient quantum chemical calculations in periodic systems based on a plane wave basis set. These tools aim at the acceleration of the convergence of the many-electron wave function expansion with respect to the virtual orbital manifold, where most of the computational cost in quantum chemical wave function based calculations originates from. We compare natural orbitals, pseudized Gaussian orbitals and explicit correlation techniques for their effectiveness and applicability to study the binding energy of the boron nitride bilayer system. It is clear that the accurate and efficient description of van der Waals interaction is a key challenge for the employed many-electron theories in this case. However, if treated correctly a number of interesting and important questions can be addressed by the employed many-electron theories, including: How large is the interlayer binding energy in h-BN and which diagrammatic techniques are needed for its accurate description? Which stacking order is the most stable at zero temperature? What is the asymptotic behavior of the interaction energy for large interlayer distances? The above questions are difficult to address by experiment and partly also by state-of-the-art ab initio theories. Currently available exchange and correlation energy functionals are often not accurate enough to allow for reliable conclusions and highly accurate diffusion Monte Carlo methods suffer from the fact that the interaction energies are small compared to the statistical errors. As such, quantum chemical wave function based methods seem to be very well-suited for this challenge: these methods are systematically improvable and do not suffer from statistical errors in general.
During the last few years hexagonal boron nitride layers have attracted a lot of interest. Some studies have focused on the interaction of single h-BN layers with molecules [7,8]. These systems have partly been investigated due to their potential applicability ranging from novel optoelectronic devices to desalination plants. From a theoretical perspective these materials are also very interesting because of their potential use as benchmark systems. Accurate and reliable benchmark data is needed for the development of more approximate and efficient theories. In contrast to molecular quantum chemistry, where accurate benchmark data for sets of weakly bound molecules are well-established, a reliable benchmark data set for weakly bound layered materials is still lacking. In this regard h-BN sheets show many advantageous features, making them well-suited for benchmarking studies. The constituent elements have few valence electrons, keeping the computational cost low. Furthermore, h-BN sheets exhibit a large band gap and are therefore computationally much cheaper and simpler to treat with high-level methods than for example graphene, which requires very dense k-point meshes for an accurate description [9].
A number of studies of the multilayer h-BN system have already been carried out using (van der Waals corrected) DFT, RPA, LMP2 theory and even quantum Monte Carlo methods [10][11][12][13][14][15][16][17][18][19][20]. These studies have mostly focused on the relative stabilities of different stacking faults, the interlayer binding energies, sliding paths and at the functional form of the binding energy for large interlayer distances. In this work we apply MP2 and coupled cluster theories for the calculation of these properties.

Theory
We now turn to the discussion of the employed theories and methods. In this work we use fully periodic implementations of (explicitly correlated) second-order Møller-Plesset perturbation and coupled cluster theory. Their implementation in the Vienna ab initio simulation package (VASP) has already been discussed elsewhere [21][22][23][24]. We note that VASP employs a plane wave basis set and the projector augmented wave method [25][26][27].
It is part of the aim of this work to compare various approaches for the treatment of the virtual orbital manifold, which is used for the construction of the excited Slater determinants in post-Hartree-Fock theories. These choices include a full canonical Hartree-Fock orbital basis, approximate natural orbitals and pseudized Gaussian type orbitals [23]. We note that the implementation of pseudized Gaussian type orbitals has recently been outlined in reference [28]. The use of canonical Hartree-Fock or DFT orbitals corresponds to the standard approach in post-HF or post-DFT calculations of extended systems using a plane wave basis set. These orbitals diagonalize the mean-field Kohn-Sham or Fock operators. However, these orbitals are not optimized in order to provide a compact and rapidly convergent virtual orbital basis for the treatment of electronic correlation. This shortcoming becomes most strongly pronounced in the calculation of molecular systems or surfaces, where supercells with large regions of vacuum have to be employed in order to minimize interactions between periodic images. Even simple atoms or small molecules require several hundred virtual orbitals in order to achieve basis set convergence.
We have already shown in reference [23] that these convergence problems can be solved efficiently by using socalled approximate natural orbitals in post-HF theories. Natural orbitals lead to a much more rapid convergence of the correlation energy with respect to the virtual orbital basis set size. Furthermore, it can be seen that natural orbitals are more localized to regions where the electronic density is large and electronic correlation effects are expected to be more important [23]. In this regard natural orbitals share some similarities with Hartree-Fock orbitals that are expanded in an atom-centered Gaussian type orbital basis set.
The approximate natural orbitals used in this work are system specific and can not be transferred to different systems, which is in contrast to the atom-centered Gaussian type orbitals basis set. The main advantage of correlation consistent Gaussian basis sets over natural orbitals is twofold, (i) they aim at a systematically improvable black box description of electronic correlation, enabling basis set extrapolation methods and (ii) their systematic truncatability also extends to chemically different systems consisting of the same atomic species, allowing for the calculation of rapidly convergent energy differences such as atomization energies. The latter is a key advantage for the efficient prediction of interaction energies between weakly interacting subsystems. For weakly interacting subsystems the dominant binding energy contribution often comes from dispersion interactions. As such it is sufficient to describe the polarizability of their constituent subsystems on the same level of accuracy. We have recently proposed a method called pseudized Gaussian type orbitals (PGTOs), which allows for employing pseudized GTOs in a planewave basis set code such as VASP [28]. In this work we will explore its accuracy for the study of the interlayer binding energy.
A large part of the total electronic correlation energy in most many-electron systems originates from the socalled electron-electron cusp condition and the surrounding correlation hole in the wave function. In the case of the uniform electron gas systems at realistic densities this contribution accounts for approximately 90% of the total correlation energy [29]. Explicitly correlated or so-called F12 theories seek to account for this contribution in the ansatz of the many-electron wave function explicitly [30,31]. This is achieved by augmenting the electron pair wave function in the employed basis of (excited) Slater determinants with terms that depend on the interelectronic distance explicitly [30,31]. In this work we will also employ a simple and efficient explicitly correlated second-order Møller-Plesset perturbation theory, which we have proposed recently in reference [24].

Computational details
For all calculations the B 2s 2 2p 1 and N 2s 2 2p 3 have been treated as valence states. The energy cutoff was set to 500 eV. As h-BN in-plane lattice constant we have employed 2.488Å, which was taken from a DFT-LDA relaxation and is in good agreement with experiment and other theoretical work [15]. All pseudized Gaussian orbital basis set calculations presented in this work employ counterpoise corrections for the basis set superposition error. To minimize the interaction between periodic images the lattice vector perpendicular to the h-BN planes was chosen to be 30Å. A summary of the different stackings that are studied in this work is provided in Figure 1. We investigate the following high-symmetry stacking faults in the h-BN bilayer system: AA, AB, AA', AB' 1 and AB' 2 . A rotation of layer A by 180 • results in A', which is equivalent to interchanging the boron and nitrogen atoms. There are always two possibilities to translate the second layer to B' 1 or B' 2 as can be seen in Figures 1d and 1e, but due to symmetry AB 1 is identical to AB 2 and will be called just AB. Whereas AB' 1 and AB' 2 differ in the atom type which are on top of each other.

Results
The results section is partitioned into three parts. The first part investigates the convergence of the calculated interlayer interactions with respect to the employed orbitals, k-point meshes and quantum chemical method. In the second part we employ MP2 theory to study the h-BN bilayer system in and around the equilibrium interlayer distance for different stackings. The third part investigates the asymptotic behavior of the binding energy for large interlayer distances.

Convergence tests
In this section we explore the convergence of the interlayer binding energies with respect to the computational parameters and methods. Therefore we restrict most of these comparisons to a rather coarse 3 × 3 × 1 k-point mesh and 64 approximate natural orbitals per k-point only, unless noted otherwise. This computational setting is good enough to benchmark the relative accuracy of different quantum chemical wave function based methods. Furthermore, this subsection will consider the AA' stacking fault only. For a better comparison, the interaction energy curves shown in Figures 3-5 have been aligned at the interlayer distance of 4.9Å. Figure 2 shows the MP2 bilayer interaction energy as a function of the interlayer distance for different k-point meshes. The interaction energy changes by less than 5 meV/(2BN) if the k-mesh is changed from 8 × 8 × 1 to 10 × 10 × 1. However, we note that even a 6 × 6 × 1 k-point mesh is not good enough if one aims at an accuracy for the interaction energy of better than 10 meV/(2BN). The calculations employ a pseudized aug-cc-pVTZ basis set and corresponding results have been converged with respect to the plane-wave cutoff energies.
As a further convergence check of our calculations we compare different choices for the virtual orbital manifold in the MP2 calculations. Figure 3 compares the interlayer binding energy as a function of distance using a complete basis set for the given cutoff energy, approximate natural orbitals and the pseudized aug-cc-pVTZ basis set. Our findings show that 224 natural orbitals are sufficient to achieve almost perfect agreement with results obtained  using a complete basis of several thousand orbitals. Furthermore, we note that approximate natural orbitals account accurately for the interlayer interaction energy even if truncated at relatively small numbers of about 96 natural orbitals. However, the procedure for calculating the natural orbitals is computationally expensive and can at the moment not be applied to very large systems. This is mostly related to technical issues in our current implementation and a solution to this problem is currently being developed. Nonetheless we conclude from the results shown in Figure 3 that a pseudized aug-cc-pVTZ basis set is accurate enough to capture the interlayer interaction energy to within 5% relative error. Another issue we seek to address in this work is the role of cusp related contributions to the interlayer binding energies. Figure 4 shows MP2 and pMP2-F12 binding energy curves using 32, 64 and 96 natural orbitals. The inset of Figure 4 shows the differences of the calculated (MP2/pMP2-F12) energies to the reference MP2 results obtained using all orbitals. We note that pMP2-F12 energies converge significantly faster with respect to the number of orbitals compared to the plain MP2 energies. Already 64 natural orbitals suffice to achieve non-parallelity errors smaller than 1 meV/atom. Furthermore we note that the pMP2-F12 energy calculated using 96 natural orbitals differs from the reference MP2 results by at most 1 meV/atom. Strictly speaking the agreement between pMP2-F12 and the reference MP2 in the complete basis set limit should be perfect. The remaining discrepancy can originate from both the reference and the pMP2-F12 energies. In passing we note that the projectors employed in pMP2-F12 theory do not include contributions from the frozen core states, which would in principle be needed for a perfect agreement [24]. However, for the current purpose we regard the achieved accuracy to be sufficient. We stress that the convergence with respect to the k-point mesh plays a much more important role. As a further test for the accuracy of the predicted MP2 interaction energies we compare MP2 to coupled cluster theory binding energies. The CCSD and CCSD(T) calculations are computationally much more demanding than MP2 calculations. Figure 5 shows that MP2 binding energy curves agree very well with CCSD(T) theory findings. MP2 predicts a slightly shorter interlayer equilibrium distance than CCSD(T). Moreover MP2 predicts a minimum which is slightly below the CCSD(T) minimum, indicating a tendency to overestimate the binding energy. We also note that the asymptotic behavior of MP2 and CCSD(T) for larger interlayer distances is almost identical. The agreement between CCSD and CCSD(T) is much worse than the agreement between MP2 and CCSD(T). CCSD strongly underestimates the binding energies and overestimates the interlayer equilibrium distance compared to CCSD(T). This indicates that MP2 is only fortuitously good for this system. However, the comparison between MP2 and CCSD(T) reveals that MP2 is likely to overestimate the interlayer binding energy by approximately 10%.
From the above convergence studies we conclude that an 8 × 8 × 1 k-point mesh and the pseudized aug-cc-pVTZ basis set is sufficient to achieve MP2 results in and around the equilibrium for the interlayer interaction energy that are converged to within the accuracy of the MP2 method. Figure 6 compares MP2 interlayer binding energy curves for different stacking faults. Our results show that AA' and AB correspond to the most stable stacking faults, deviating by less than 1 meV/atom in their binding energies. In the order of decreasing binding energies we find the other metastable stacking orders to be AB' 1 , AB' 2 Table 1. Equilibrium interlayer distance, binding energies for different stackings and energy difference between the respective and most stable stacking (AB) in meV per atom. These properties were calculated using a pseudized aug-cc-pVTZ basis set and an 8 × 8 × 1 k-point mesh. and AA. We note that our MP2 calculations predict AB to be approximately 0.25 meV/atom more stable than AA', which is in contrast to the experimentally observed more stable AA' stacking order [32,33]. However, it seems unlikely that our predictions have such a high accuracy and we will consider AA' and AB to be degenerate within the accuracy of our results. The predicted interlayer equilibrium distances for the different stackings are in the order of increasing distances AB, AA', AB' 1 , AB' 2 and AA, following the same order as for the decreasing interlayer binding energies. Table 1 summarizes the calculated equilibrium interlayer separations, binding energies, and differences in the relative binding energy for the different stackings with respect to AB. The equilibrium interlayer distances have been obtained by fitting the curves shown in Figure 6 around the equilibrium to αe −γD − β/D 4 , where D is the interlayer distance and α, β, γ are fitting parameters [16].

Equilibrium region
Overall, our findings are in good agreement with previous LMP2 findings reported in the literature [14,15]. Close inspection reveals, however, small differences as discussed in the following. We find interlayer separations for the AA' and AB stacking faults of 3.2Å and 3.19Å, respectively. These findings are in good agreement with experimental interlayer distances of 3.25±0.1Å [33]. However, our MP2 distances are significantly smaller than the corresponding LMP2 distances of 3.31Å reported in reference [15]. The comparison between MP2 and the more accurate CCSD(T) theory shown in Figure 5 supports our finding that MP2 should underestimate the interlayer distance compared to experiment, especially if zero-point vibrational effects are not considered. In passing we note that van der Waals corrected DFT calculations also predict consistently too large interlayer distances compared to experiment. For the AA' stacking the interlayer distance is predicted to be 3.37Å (PBE+MBD), 3.4Å (vdW-TS), 3.4Å (vdW optB88), 3.5Å (vdW optPBE), 3.6Å (vdW-DF-layered) as reported in references [14], [34], [35], [35] and [12], respectively. We now turn to the comparison of the relative binding energies. The LMP2 interlayer binding energies relative to AA' as reported in reference [15] correspond to 2.21, 8.25 and 9.86 meV/atom for AB' 1 , AB' 2 and AA, respectively. This is in good agreement with our MP2 results of 2, 8.25 and 8.75 meV/atom for AB' 1 , AB' 2 and AA, respectively, as summarized in the last column of Table 1. We note that LMP2 predicts the AA' stacking to be more stable than the AB by 0.12 meV/atom, which is in contrast to our findings [15]. The absolute LMP2 binding energies are not reported in reference [15]. However, we find reasonable agreement between our MP2 interlayer binding energy for the AA' stacking corresponding to 29.25 meV/atom and the PBE+MBD findings of 24.7 meV/atom [14]. Again we stress that MP2 is likely to overbind compared to experiment because MP2 also overbinds compared to CCSD(T) theory as shown in Figure 5.
To further elucidate the reliability of the MP2 results we have computed the h-BN interlayer binding energies using the (direct) random-phase phase approximation (RPA) and identical compuational settings. The RPA calculations predict AB to be 0.25 meV/atom more stable than the AA' stacking fault, confirming our MP2 findings. However, the absolute RPA interlayer binding energy and equilibrium distance is only 75 meV/(2BN) and 3.3Å, respectively. RPA suffers in general from a tendency to underestimate binding energies and overestimate bond lengths especially for weakly interacting systems [36,37]. From this and the comparison between MP2 and CCSD(T) we can conclude that the correct interlayer binding energy should be between 75 meV/(2BN) and 118 meV/(2BN), whereas the correct interlayer distance should be between 3.2Å and 3.3Å.
As for the comparison between MP2 and QMC results reported in reference [16] we note that the QMC results have been carried out using rather small supercells corresponding to a 4 × 4 × 1 k-point mesh only. Even though our MP2 results agree well with the QMC findings if we employ the 4 × 4 × 1 k-point mesh, the binding energies are not converged for this sampling density as shown in Figure 2.

Asymptotics of dispersion interaction
The final question we seek to address is the aymptotics of the binding energy for large interlayer distances. This question is of great interest because it was shown that the asymptotics is strongly influenced by the dimensionality and the metallic or insulating character of the different fragments at large distances [9,38,39].
In a recent QMC study it was shown that the binding energy between two h-BN sheets exhibits an asymptotic behavior of D −4.5 [16], whereas analytical theory suggests D −4 [38]. However, due to the large computational cost the QMC studies have been limited to rather small 3 × 3 and 4×4 supercells. In this work we employ MP2 theory to investigate the asymptotic behavior of the bilayer binding energy. Figure 7 shows the h-BN interlayer binding energy as a function of the interlayer distance on a log-log scale for k-point meshes ranging from 4×4×1 to 10×10×1. The fits are depicted by the dashed and dotted lines assuming a D −4 and D −4.5 asymptotics, respectively where D refers to the interlayer distance. We have fitted to MP2 binding energies calculated using 4 × 4 × 1 and 10 × 10 × 1 k-point meshes. The comparison between the 4 × 4 × 1 k-point 4 mesh results and the fitted D −4 and D −4.5 asymptots reveals that for larger distances the binding energy decays even faster than D −4.5 as D increases. However, Figure 7 illustrates convincingly that the slope for denser k-point meshes becomes less steep and approaches a D −4 behavior. We note that even the MP2 binding energies obtained using the 10 × 10 × 1 k-point mesh start to deviate from the D −4 fit for D ≥ 4.9Å, although it seems probable that this discrepancy would become smaller if denser k-point meshes would be employed.

Conclusions
In conclusion we have shown that MP2 theory allows for a relatively accurate description of the binding energies in the h-BN bilayer system for a wide range of interlayer distances and different stacking faults. Our MP2 results predict the following (metastable) stackings in the order of decreasing binding energies AB, AA', AB' 1 , AB' 2 and AA. The binding energies of the metastable stackings relative to the most stable binding energy is in good agreement with previously reported LMP2 results. However, in contrast to previously reported findings, our MP2 and RPA results indicate that the AB stacking fault is more stable than the AA' for the bilayer system [14,40]. The energy difference between AB and AA' is only 0.25 meV/atom and it seems probable that this is smaller than the accuracy of our calculations. We note that all van der Waals corrected density functionals reported in the literature seem to predict too large interlayer distances compared to our MP2 (3.2Å) and RPA (3.3Å) findings that are in good agreement with experiment (3.25 ± 0.1Å).
We have shown that approximate natural orbitals, pseudized Gaussian orbitals and explicitly correlated methods help to significantly improve the accuracy and rate of convergence of the binding energies with respect to the employed size of the basis set. This indicates that the interlayer binding energies and the polarizabilities of the individual sheets are also sensitive to cusp related electronic correlation effects.
Finally, we have addressed the aysmptotics of the interlayer binding energy using MP2 theory and very dense k-point meshes. Our findings conclusively show that the binding energy decays as D −4 in the investigated range of distances, which is in good agreement with other theoretical predictions.