Van der Waals interaction in a boron nitride bilayer

The van der Waals (vdW) interaction is one of the most fundamental physical quantities resulting from quantum fluctuation of charges. However, it remains a challenge to account for this interaction quantitatively in both theory and experiment. The vdW interaction is usually described in terms of the ${{R}^{-6}}$ attraction between two neutral and polarizable particles separated by a distance $R\gg r$, where r is the size of the particles. It is only recently that the direct measurement of the vdW interaction between two Rydberg atoms has been accomplished, which confirmed its ${{R}^{-6}}$ dependence and determined the strength [1]. By summing up the ${{R}^{-6}}$ contributions between the composing microscopic subsystems, the vdW interaction between two macroscopic bodies separated by a distance D can be determined. The distant attraction between two thin macroscopic layers is concluded to follow ${{D}^{-4}}$ according to the standard pair-wise summation method. However, it has been pointed out that the interaction for zero-gap layer systems should decay more slowly [2]. This was certified later for the bilayer system of a two-dimensional (2D) homogeneous electron gas by quantum Monte Carlo (QMC) methods [3] and for graphite, using the adiabatic-connection fluctuation-dissipation theorem in the random phase approximation [4], though a different conclusion of ${{D}^{-4.2}}$ behavior for graphite was also reported using QMC calculations [5]. For insulating layer systems, it is usually presumed that the ${{D}^{-4}}$ behavior applies.

The vdW interaction is one of the physical properties that the standard local and gradient approximation (LDA [6] and GGA [7]) of density functional theory (DFT) [8] methods fail to describe correctly [9]. In recent years, there have been many proposed vdW density functionals (vdWDF) to overcome this deficiency [10]. However, discrepancies in binding energy among these vdWDF results are usually apparent, and these vdWDF schemes have been used to study behavior beyond the equilibrium distance. An alternative approach to tackle this issue is to perform one of the currently most accurate electronic calculations for materials, i.e., QMC methods [11]. In this article, we present the results of studying the BN bilayer system by QMC methods. We calculated the binding energy and studied the behavior of the vdW interaction beyond the equilibrium distance (but not yet reaching the asymptotic region) in this prototypical insulating bilayer system. These results are also compared with those obtained from the standard DFT as well as those employing vdW density functionals. The outcome of the present QMC study would provide a benchmark for future generation of vdWDF and guidance for prospective experiments.

The QMC calculations were performed using the CASINO code [12]. The trial wave functions of Slater–Jastrow form are used, and the Slater orbitals are generated from the DFT calculations using the CASTEP [13] plane-wave-basis code. The tabulated relativistic Hartree–Fock pseudopotentials from the CASINO pseudopotential library are used to generate the initial wavefunctions for Slater orbitals [14, 15], and Casulaʼs approach [16] is used in DMC calculations to minimize the nonlocal pseudopotential error. In the orbital-generation calculations, the plane-wave energy cutoff is taken as 2000 eV. This is much higher than those of the ordinary DFT calculations in order to reduce the variance in the QMC calculations [17].

The QMC calculations start with the variational quantum Monte Carlo method to obtain the optimized parameters in the Jastrow factor, which include isotropic electron–electron, electron–nucleus, and electron–electron–nucleus terms. That is, J can be written as

$$\begin{eqnarray*}\begin{array}{rcl} J(\{{{r}_{i}}\},\{{{r}_{I}}\}) & = & \mathop{\sum }\limits_{i=1}^{N-1}\mathop{\sum }\limits_{j=i+1}^{N}u({{r}_{ij}})+\mathop{\sum }\limits_{I=1}^{{{N}_{{\rm ions}}}}\mathop{\sum }\limits_{i=1}^{N}{{\chi }_{I}}({{r}_{iI}}) \\ {} & {} & +\mathop{\sum }\limits_{I}^{{{N}_{{\rm ions}}}}\mathop{\sum }\limits_{i=1}^{N-1}\mathop{\sum }\limits_{j=i+1}^{N}{{f}_{I}}({{r}_{iI}},{{r}_{jI}},{{r}_{ij}}) \\ \end{array}\end{eqnarray*}$$

where N and N_ions are the numbers of electrons and ions, respectively, and r_ij and r_iI denote the distance between electrons i and j and also between electron i and ion I respectively [18]. The u, χ, and f terms are power expansions containing variational parameters and depend on the spins of the electrons. The optimization was performed with variance minimization [19] and followed by energy minimization [20]. The nodal surface of the wave function is constrained to that of the Slater determinant constructed from the DFT-generated orbitals, i.e., the fixed-node approximation [21]. The DFT-generated orbitals are represented numerically using a real-space splines grid in order to improve the system-size scaling of the QMC calculations [22]. The optimized wave functions are subsequently used to generate the initial ensemble of electronic configurations for the diffusion quantum Monte Carlo (DMC) calculation through which the ensemble is propagated by the imaginary-time Schrödinger equation toward the ground state [11]. The model periodic Coulomb interaction is used to reduce finite-size effects related to the long-ranged Coulomb interaction [23].

Throughout this study, the supercells consisting of the AA' stacking of two BN layers are used, with c = 30 Å as the size of the supercell perpendicular to the basal planes. The supercell with the equally spaced BN layers, i.e., D = c/2 = 15 Å, is taken as the reference system when the interlayer energy is evaluated. We expect excellent error cancellations in the interlayer interaction energy with respect to the k-point sampling in the DFT methods by using these corresponding reference systems⁴ .

We have carried out QMC calculations using both the lateral (3 × 3) and (4 × 4) supercells in order to estimate the finite-size bias of the QMC methods. The orbital-generation DFT calculations used the shifted Monkhorst–Pack (MP) [24] k-point meshes of (3 × 3 × 1) and (4 × 4 × 1) in the primitive supercell for the two QMC supercells, respectively. The lateral (3 × 3) supercell with c = 20 Å was also considered in order to assess the vacuum-size correction. The intralayer lattice constant is fixed at 2.5 Å, which is the experimental value [25] as well as the equilibrium lattice constant of the LDA calculations. Time steps of 0.003 a.u. were used for the DMC production calculations to ensure an acceptance ratio of more than 99.9%. The target population used in all the DMC calculations is 3840. A total of more than 5 million core-hours are required in order to achieve a standard deviation of less than 2.5 meV/2BN for the DMC total energy of all the various supercell sizes considered in this study.

For comparison purposes, we also carried out the LDA, Perdew–Burke–Ernzerhof (PBE) [7] GGA and vdWDF calculations. We used the vdW density functionals proposed by Dion et al (vdW-DF) [26] as well as those with a different exchange functional [27, 28] or correlation functional (vdW-DF2) [29], which are all implemented by the algorithm of Roman–Perez and Soler [30] in the plane-wave-based Vienna $ab$ $initio$ simulation program (VASP) [31, 32]. The vdW-DF and vdW-DF2 are first-principles vdW functionals, while the generations of vdW-optXXX involve empirical fitting to a benchmark dataset of weakly interacting dimers and water clusters. The interaction between ions and valence electrons is described by the projector augmented wave method [33], and the numbers of valence electrons included for B and N are three and five, respectively. The energy cutoff that determines the number of plane waves is more than 900 ${\rm eV}$ [34], and the (9 × 9 × 1) MP k-point mesh is used. Tests on higher energy cutoffs, denser MP meshes, and larger vacuum size show that the binding energy is converged to within 1 ${\rm meV}/2BN$, and the equilibrium distance stays the same.

The interaction energies of all the different DFT schemes are plotted in figure 1 and their binding energy and equilibrium distances (D₀) are presented in table 1. The GGA scheme leads to a tiny binding energy and too-large equilibrium distance, as demonstrated in previous works [29, 35]. The LDA approach, performed by using VASP and CASTEP (with different pseudopotentials as stated in the previous section), yields binding energies of around 55 meV/2BN and equilibrium distance of 3.3 Å, which is close to the interlayer distance of the experimental value of 3.33 Å for the hexagonal BN bulk [25]. The binding energy obtained from both the vdW-DF and vdW-DF2 schemes are about 40 meV/2BN larger than those of the LDA scheme. A similar result of 100 meV/2BN was published previously [35], using the earlier version of the vdW-DF scheme (denoted as vdW-DF-layered) [36]. The vdWDF results based on the vdW-DF version but with different exchange functionals (optPBE, optB88, and optB86b) [27] lead to an extensive range in binding energy, i.e. from 116–127 meV/2BN. However, the largest binding energy of 152 meV/2BN has also been reported previously [37] by using a different vdWDF, i.e., vdW-TS [38]. It will be demonstrated in the following section that the two first-principles vdW functions, i.e., vdW-DF and vdW-Df2, give binding energies closest to the QMC results.

Zoom In Zoom Out Reset image size

Table 1.  The binding energies E(D₀) and equilibrium distances D₀ of the BN bilayer obtained from different DFT schemes (details in the text) in this study and previous publications (the last two lines).

Schemes	E(D0) (meV/2BN)	${{{\rm D}}_{0}}(\overset{\circ}{A} )$	Codes
PBE [7]	3	4.3	VASP
LDA [6]	55	3.3	VASP
vdW-DF [26]	95	3.6	VASP
vdW-DF2 [29]	91	3.5	VASP
vdW_optPBE [27]	116	3.5	VASP
vdW_optB88 [27]	126	3.4	VASP
vdW_optB86b [27]	127	3.3	VASP
PBE [7]	3	4.4	CASTEP
LDA [6]	57	3.3	CASTEP
vdW-DF-layered [36]	100	3.6	[35]
vdW-TS [38]	152	3.4	[37]

The interaction energy calculated by the DMC method is plotted in figure 2. The binding energy close to the equilibrium interlayer distance obtained by the DMC method are listed in table 2. Due to the enormous amount of CPU time required in the QMC calculations, only selected data points at quarters of integer Ds (in unit of Å) are included in the QMC study. The lowest-energy interlayer distance was found to be either 3.25 Å or 3.50 Å. The binding energy of the larger supercell (4 × 4) was found to be around 72 meV/2BN while that of the smaller cell around 66 meV/2BN. Using the similar procedure as used previously for a Ne solid [39], the binding energy extrapolated to the system with an infinite lateral size is 79.3 meV/2BN. Extrapolation for the binding energy with an infinite vacuum between two bilayer systems of neighboring supercells can be estimated by assuming the reference energy at D = c/2 decreases as D⁻⁴. This leads to a correction of 2 meV/2BN in binding energy for those evaluated with the c = 30 Å supercell. For a shorter range behavior at D = c/2, the correction would be smaller. Inclusion of the vacuum-size correction would lead to a binding energy of 81.3 meV/2BN. It is larger than the LDA result and smaller than all the vdW_opt values by more than 20 meV/2BN. The DFT schemes closest to this value are the two first-principles vdW functionals, i.e., vdW-DF and vdW-DF2. All the results in table 2 used the PBE-generated orbitals except the one denoted with an asterisk, which used LDA-generated orbitals. However, we see that these two different initial many-electron wavefunctions lead to about the same DMC evaluated binding energies.

Zoom In Zoom Out Reset image size

In addition to the binding energy, another desirable property of the BN bilayer system is its asymptotic behavior at large interlayer separations. Nevertheless, a full convergence is unreachable in practice at this moment. As shown earlier for a graphene bilayer [4], the interlayer separation must be larger than 9 Å to see a correct asymptotic behavior. In the present system, the interaction is already smaller than the numerical error at the 7 Å interlayer separation. However, an analysis for the behavior beyond equilibrium distance (bD₀) predicted by the DMC and different DFT schemes would still shed light on the long-range behavior of this system. Note that in the conventional two-term description (one repulsive and the other attractive term) of the interaction, the attractive interaction becomes dominant already at interlayer separations around 1.3 D₀, as discussed in the following. For the 6–12 Lennard–Jones potential, the repulsive energy is about one tenth, in magnitude, of the attractive energy at D ≈ 1.3D₀. For the 4–10 potential (resulting from the summation of the 6–12 potential for 2D bilayer systems), a similar condition can be reached at D ≈ 1.25D₀. If the equilibrium distance is taken as 3.50 Å, the long-range attractive behavior is dominant for D beyond 4.55 Å and 4.375 Å, respectively, for the 6–12 and 4–10 potentials. For the three fitted Buckingham-like potentials shown in figure 2, the repulsive energies are already less than one seventeenth of the attractive energy at D = 4.5 Å, in magnitude.

The behavior beyond the equilibrium separation for the different DFT schemes (figure 1) can be sorted into four groups, i.e., PBE, LDA, vdW-DF2 and that including vdW-DF, vdW_optPBE, vdW_optB88, and vdW_optB86b. For the 4th group, we will take the vdW_optPBE as the representative scheme for later comparison with those by the QMC methods. The bD₀ behavior for these four DFT schemes are fitted with both the power law and exponential form. We found that, overall, the exponential form fitted better for these four DFT schemes. However, the power-law fitting provides us with a direct comparison to the D⁻⁴ law. All these interaction energies are plotted on a logarithmic scale in order to better display their power-law behaviors (figure 3). The bD₀ features for both the PBE and LDA results are extremely short-range with the D⁻¹⁰ dependence. The vdW-DF2 and vdW_optPBE schemes result in the D⁻⁶ and D$^{-5.3}$ dependence, respectively, a considerably longer-range behavior as compared to the LDA scheme.

Zoom In Zoom Out Reset image size

For the DMC obtained results, the behavior beyond equilibrium distance refers to the interaction between 4 Å and 7 Å interlayer separation which includes nine data points for DMC(3 × 3) and five data points for both DMC(4 × 4), and DMC(4 × 4)*. The best fitted exponents (using the least absolute deviation fitting) are −4.50, −4.48, and −4.49 for DMC(3 × 3), DMC(4 × 4), and DMC(4 × 4)*, respectively. The average residues are all smaller than the error of the DMC total energy, i.e. 2.2, 0.9, and 2.0 meV/2BN for DMC(3 × 3), DMC(4 × 4) and DMC(4 × 4)* respectively. The bD₀ DMC data can therefore be described as being best fitted with the D$^{-4.5}$ dependence. The previous study for a graphene bilayer [4] showed that the exponent can be increased (decreased in magnitude) by as much as 1.2 for interaction beyond 9 Å when compared to those of smaller separations (but beyond equilibrium distance). To what extent this conclusion can be applied to the present system is not apparent, as graphene is a zero-gap layer, while BN has a wide band gap. However, the presently obtained D$^{-4.5}$ dependence for the behavior beyond equilibrium separation is shown to be longer-range than all the available DFT schemes.

In conclusion, we have carried out QMC calculations to determine both the binding energy and distance-dependent power law of the interaction beyond the equilibrium separation in the BN bilayer. The binding energy obtained by the DMC method, i.e., 81 meV/2BN, was found to be larger than the LDA value and smaller than all the considered and previously published vdWDF ones by sizable quantities. We determined the power law of the BN interlayer interaction beyond equilibrium distance as D$^{-4.5}$ , while those of the available DFT schemes all fall off faster than this dependence.

We thank N D Drummond for helpful discussions. This work was supported by the Ministry of Science and Technology of Taiwan (under contract numbers 101-2112-M-006-012-MY3 and 102-2112-M-001-023-MY3). Part of the calculations were performed at the National Center for High-Performance Computing (NCHC). Computing resources from the Academia Sinica Computing Center of Taiwan are acknowledged. We also thank the National Center for Theoretical Sciences (NCTS) for support through the Computational Material Research (CMR) focus group.