Van der Waals Density Functional Theory vdW-DFq for Semihard Materials

Abstract

There are a large number of materials with mild stiffness, which are not as soft as tissues and not as strong as metals. These semihard materials include energetic materials, molecular crystals, layered materials, and van der Waals crystals. The integrity and mechanical stability are mainly determined by the interactions between instantaneously induced dipoles, the so called London dispersion force or van der Waals force. It is challenging to accurately model the structural and mechanical properties of these semihard materials in the frame of density functional theory where the non-local correlation functionals are not well known. Here, we propose a van der Waals density functional named vdW-DFq to accurately model the density and geometry of semihard materials. Using $\beta$-cyclotetramethylene tetranitramine as a prototype, we adjust the enhancement factor of the exchange energy functional with generalized gradient approximations. We find this method to be simple and robust over a wide tuning range when calibrating the functional on-demand with experimental data. With a calibrated value $q=1.05$, the proposed vdW-DFq method shows good performance in predicting the geometries of 11 common energetic material molecular crystals and three typical layered van der Waals crystals. This success could be attributed to the similar electronic charge density gradients, suggesting a wide use in modeling semihard materials. This method could be useful in developing non-empirical density functional theories for semihard and soft materials.

1. Introduction

Semihard materials generally consist of molecules. Energetic materials are a common class of semihard materials with light composites of elements C, H, O, and N, of which the molecular sizes are moderate (ranging from 10 to 1000 atoms) leading to fast chemical kinetics so that the chemical reactions during decomposition or detonation occurs in a few picoseconds. To gain high energy densities, these molecules are mostly bonded via intermolecular London dispersion forces to form molecular crystals (Figure 1a–k), which belong to semihard materials since they have mild hardness. The mass density is closely related to the energy density and detonation performance. Therefore, the accurate prediction of the mass density within 3% of experimental values is critical in predicting the reaction rate and energy density. In general, Density Functional Theory (DFT) [1] is a reliable and popular tool to predict material properties and material designs. However, the London dispersion force is poorly described with conventional semi-local approximations to DFT [2] since it is essentially a nonlocal electron correlation effect. As a result, an accurate description of the dispersion force is essential in DFT modeling of energetic materials.

There are extensive studies to improve the modeling of dispersion, or van der Waals (vdW) interactions [2,3], such as DFT-D [4,5,6,7,8,9], vdW-TS [10,11], DFT-vdW${}^{surf}$ [12], and vdW-DF [13,14,15,16] methods. The DFT-D1 and DFT-D2 methods are classified as the 1st rung of Jacob’s Ladder for DFT-based dispersion correction schemes according to Klimes and Michaelides [3]. The DFT-D3 and vdW-TS methods belong to the next rung that has environment dependent $C_6$ corrections. Non-local density functionals like DFT-DF and DFT-DF2 methods sit in the 3rd rung [3], of which the overall accuracy increases, as well as the computational demand. Recently, a few improvements in vdW-DF methods have been proposed, including optPBE-vdW, optB88-vdW, optB86b-vdW [17,18], vdW-DF2-C09 [19], vdW-DF-cx [20], rev-vdW-DF2 [21], PBE+rVV10 [22], PBEsol+rVV10 [23], SCAN+rVV10 [24], and SG4+rVV10 [25]. These methods are designed to minimize the error of the binding energies and interaction energy curves in a certain collection of materials, such as the S22 data set [26], aiming for a general use for all materials. The accuracy of these models on energetic materials like HMX has not been assessed yet. We found that the prediction of the mass density in HMX is still not satisfied, as opposed to the claimed success in many other materials. Here we propose a vdW-DF-based method (we denote it as vdW-DFq) using a single parameter to tune the exchange functional so that the target mass density can be achieved. We are looking for a reliable model for a particular material, instead of a general formula for all materials.

The cyclotetramethylene tetranitramine (HMX) shown in Figure 1a–c is an important secondary explosive. The P2${}_1$/c monoclinic crystal structure, namely $\beta$-phase, (Figure 1a) is thermodynamically the most stable phase of HMX under ambient conditions. $\beta$-HMX has a high detonation velocity and high energy density. HMX is most commonly used in military and industrial applications including automobile air bags, rocket propellants, and polymer-bonded explosives. Therefore, $\beta$-HMX has been extensively studied [27,28,29]. We focus on $\beta$-HMX as a representative of semihard materials because of its moderate bulk modulus (around 15 GPa).

2. Methodology

Within the DFT frame work, the van der Waals density functional (vdW-DF) [13] is good for both covalent and van der Waals (vdW) interactions in a seamless fashion. The exchange correlation energy $E_{xc}$ as sum of exchange energy $E_x$ and correlation energy $E_c$ is expressed as

$$E_{xc}=E_x+E_c=E_x^{GGA}+E_c^{LDA}+E_c^{NL}$$

(1)

where $E_x^{GGA}$ is the exchange energy within the generalized gradient approximation (GGA) [30], $E_c^{LDA}$ is the short-range correlation energy within the local density approximation (LDA) [31], and $E_c^{NL}$ is the non-local correlation energy is long-range beyond LDA. $E_c^{NL}$ could be expressed as

$$E_c^{NL}=\frac{1}{2}\int \int d\mathbf{r}d{\mathbf{r}}^{\prime }n\left(\mathbf{r}\right)\varphi (d,d^{\prime })n\left({\mathbf{r}}^{\prime }\right),$$

(2)

where $n\left(\mathbf{r}\right)$ is the electron density at the point $\mathbf{r}$, and $\varphi (d,d^{\prime })$ is the vdW kernel [13], $d=q_0\left(\mathbf{r}\right)|r-r^{\prime }\left|\right)$, $d=q_0\left({\mathbf{r}}^{\prime }\right)|r-r^{\prime }\left|\right)$, $q_0$ is a function of $n\left(\mathbf{r}\right)$ and its gradient $|\nabla n(\mathbf{r}\left)\right|$, which determines the long-range asymptote as well as the short-range damping of the vdW kernel. The choice of LDA for the short-range correlation is to avoid possible double counting of the gradient on $n\left(\mathbf{r}\right)$ and $|\nabla n(\mathbf{r}\left)\right|$. Ultimately, this choice implicitly implies that the only non-local correlation term is due to vdW interactions, where the vdW kernel is actually derived in a way to model dispersion. The computational cost for vdW-DF calculations is marginal as compared with orbital dependent functionals and recent developments in efficient algorithms [32,33,34]. As a result, vdW-DF calculations are feasible in large calculations at a computational cost comparable with that of conventional semi-local DFT, at about only 30% additional computation time from the authors’ experience.

The general form of the GGA exchange energy functional is given as

$$E_x^{GGA}=\int d\mathbf{r}n\left(\mathbf{r}\right){ϵ}_x^{unif}\left[n\left(\mathbf{r}\right)\right]F_X^{GGA}\left[s\left(\mathbf{r}\right)\right],$$

(3)

where ${ϵ}_x^{unif}\left(n\right)=-3k_F/4\pi$ with $k_F={\left(3{\pi }^{2}n\right)}^{\frac{1}{3}}$ is the exchange energy for the uniform electron gas, $F_X^{GGA}\left(s\right)$ is the GGA enhancement factor, which distinguishes one GGA from another. The variable s in $F_X^{GGA}\left(s\right)$ is the dimensionless reduced gradient, defined as

$$s=\frac{|\nabla n|}{2{\left(3{\pi }^{2}n\right)}^{1/3}n}.$$

(4)

The reduced gradient expresses how fast the density varies on the scale of the local Fermi wavelength $k_F$. When s is set to zero, the GGA functional reduces to LDA. There are quite a few proposals of the $E_x^{GGA}$ for vdW-DF [17,18,19,32].

The nonlocal vdW-DF consistently overestimates the inter-fragment distance, although it improves the description of van der Waals interactions over conventional semi-local DFT. The original vdW-DF uses the revised Perdew-Burke-Ernzerhof (revPBE) [35] exchange functional for weakly interacting molecular pairs. However, the revPBE exchange correction is too repulsive at small separations [36]. Therefore, a few efforts are devoted to develop a better exchange functional for the nonlocal correlation functional. Klimeš et al. [17,18] have introduced several exchange functionals with accurate binding energies of molecular duplexes in the S22 data set [37]. Cooper constructed an exchange functional (named C09x) based on revPBE [19]. The second version of vdW-DF (vdW-DF2) uses the revised PW86 (rPW86) exchange functional [14]. Additionally, there are a few other nonlocal correlation functionals [38,39,40,41].

Recently a vdW-DF optimized semi-empirically based on the Bayesian error estimation named BEEF was proposed [42]. Hamada also proposed an exchange functional for the second version of vdW-DF (vdW-DF2). The so called rev-vdW-DF2 [21] is based on the GGA exchange functional of optB86b [18] with a revised q value of 0.7114 (as named B86R). A very recent study suggests that rev-vdW-DF2 is very accurate for weakly bound solids overall compared to other methods. However, it is still not satisfactory on our demanded accurancy in geometry. Nevertheless, the performance of the improved functionals depend on the system, and for the applications to a variety of systems including surfaces and interfaces, which consist of materials with different bonding natures, accuracy of vdW-DF has yet to be improved.

It is worth noting that all these methods are targeting general materials and all the properties, which is definitely a non-trivial work because the exchange enhancement factor is hard to satisfy in all systems where the reduced density gradient is different from case to case. In some cases, one needs a well tailored density functional only for a particular property of interest. For example, the mass density is the primary concern in prediction of the detonation performance of energetic materials. Thus, a relatively simple, accurate, and tunable density functional theory is indispensable.

To obtain an accurate exchange functional for vdW-DF, the following properties are considered here. First, $E_x^{GGA}$ should recover the second-order gradient expansion approximation (GEA) [43] at the slowly varying density limit as

$$F_X^{GEA}\left(s\right)=1+\mu s^2,$$

(5)

with $\mu =10/81$. This property plays an important role in predicting equilibrium geometries for solids and surfaces [44]. This constraint has been applied in other exchange functionals [18,20,45]. It is worth pointing out that C09 employs $\mu =0.0864$. Second, at large gradient limit, it proposed that $F_x$ should have $s^{2/5}$ dependence so that it can avoid the spurious binding from exchange only. So far only Perdew and Wang (PW86) [46] and Becke (B86b) [47] exchange functionals fulfill this reqirement [45].

Murray et al. [45] re-parameterized PW86 to satisfy these known limits (PW86R), which leads to vdW-DF2. However, it turns out that the exchange functional is repulsive in solids and adsorption systems, leading to overestimated equilibrium distances [48].

Another efforts is that Hamada proposed a revised B86b exchange functional (B86R), which improves the description of the attractive van der Waals interactions near the equilibrium over the original vdW-DF. The enhancement factor for the B86b exchange functional is given by

$$\begin{array}{cc}\hfill F_X\left(s\right)& =1+f_q{\mu }_0{s}^2,\end{array}$$

(6)

$$\begin{array}{cc}\hfill f_q& =\frac{1}{{(1+{\mu }_0{s}^2/q)}^{4/5}}\end{array}$$

(7)

When $s\to 0$, the factor $f_q=1$ and the first condition is fulfilled. The value of q is arbitrarily determined by a least square fit to $F_X$ of the original B86b at $8<s<10$.

The first order derivative of the exchange enhancement factor with respect to the reduced density gradient $\frac{d}{ds}{F}_X$ is important for binding separations [20], which can be obtained directly, as

$$\frac{d}{ds}{F}_X\left(s\right)=2{\mu }_0{f}_{q}s-1.8{\mu }_0^2{f}_q^{9/4}/q{s}^3.$$

(8)

We can determine q by reproducing the B86b asymptote at $s\to \infty$. However, the results are still questionable in improving accuracy. Thus Hamada have obtained $q=0.7114$ by the least square fit. Using the same principle, Berland and Hyldgaard [20] introduce an exchange functional, namely LV-PW86r, for the original vdW-DF which leads to vdW-DF-cx.

Here we propose the vdW-DFq method to tune the q for the particular system as a calibration to the experiments. For instance, we find a very good agreement with experiment in predicting the mass density of $\beta$-HMX using $q=1.05$. Figure 2a,b illustrates the enhancement factor of B86q with the $q=1.05$ exchange energy functional along with ten other corresponding ones in the low-gradient (panel a) and high-gradient domains (panel b). Their derivatives are displayed in the lower part of the figure.

It is desirable that the exchange energy functional fullfills the gradient at the slowly varying density limit, in addition to the revPBE asymptote at the large reduced gradient. Therefore, we roughly divide the range of the reduced density gradient into two domains: 0–3 for the slowly varying density domain (low-gradient domain) and 3–10 for the large reduced gradient domain (high-gradient domain). Plots of the enhancement factors provide a way to visualize the s dependence of the GGA.

A universal accurate function of the enhancement factor might be hard find. However, it is practical to tune the enhancement factor for a particular kind of materials that have similar electronic charge density gradients s. To achieve that, a function of enhancement factor that can be tuned in a large range is expected to explore the more $F_X-s$ and $\frac{d}{ds}{F}_X-s$ space. We find that the enhancement factor in form of Equation (5) are surprisingly good for tuning enhancement factors by varying q values. One can see in Figure 3 that both the enhancement factor and its derivative vary over a wide range when q is changed from 0.1 to 100.0. This illustrates that q is an effective variable to tune the exchange energy functional and binding separations.

3. Results and Analysis

The proposed vdW-DFq method is firstly tested for $\beta$-HMX, whose primitive cell is shown in Figure 1a. The geometry optimization of the primitive unit cell starts with the experimental configuration [49]. The atoms are then relaxed to the configuration with the minimum total energy. Our optimized crystal structure has a space group of $P{2}_1/c$. The optimized lattice constants using various methods are summarized in

Table 1. Geometry. Lattice constants $a,b,c$, lattice angle $\beta$, and volume of the unit cell V from various methods along with the relative standard deviation (RSD%) compared with experiments.

	a(Å)	b(Å)	c(Å)	$\mathit{\beta }$	V(Å${}^3$)	RSD
Expt. ${}^a$	6.537	11.054	8.7018	124.443	518.558
PBE ${}^b$	6.673	11.312	8.894	124.395	553.99	6.8
PBE-D2 ${}^b$	6.542	10.842	8.745	124.41	511.73	−1.4
PBE-D3 ${}^c$	6.541	10.894	8.748	124.38	514.488	−0.8
PBE-TS ${}^c$	6.530	10.985	8.770	124.572	517.992	−0.1
RPBE-D3 ${}^c$	6.580	11.061	8.762	124.357	526.501	1.5
PBEsol ${}^c$	6.603	10.928	8.709	124.002	520.949	0.5
vdW-DF	6.653	10.096	8.901	124.756	539.880	4.1
vdW-DF2	6.575	10.946	8.765	124.579	519.387	0.2
optB88-vdW	6.476	10.662	8.680	124.735	492.535	−5.0
optB86b-vdW	6.481	10.635	8.700	124.749	492.710	−5.0
optPBE-vdW	6.550	10.877	8.760	124.703	513.088	−1.1
vdW-DF-C09	6.435	10.494	8.650	124.807	479.676	−7.5
vdW-DF-cx	6.556	10.787	8.779	124.557	511.324	−1.4
rev-vdW-DF2 ${}^d$	6.504	10.696	8.685	124.670	496.905	−4.2
vdW-DFq e	6.569	10.921	8.737	124.316	517.742	−0.16
optB86q ${}^f$	6.922	11.773	8.985	121.986	620.991	19.8
vdW-DFk’ ${}^g$	6.486	10.636	8.694	124.713	493.013	−4.9

${}^a$ Ref. [50]; ${}^b$ Ref. [51]; ${}^c$ Ref. [64]; ${}^d$$q=0.7114$; ^e$q=1.05$; ${}^f$ optB86q-GGA calculations only (without vdW-DF2 corrections); ${}^g$ using $Z_{ab}=-0.8491$.

. The relative standard deviation (RSD) or the error of the predicted volumes to the experiment are also listed as a measurement of the accuracy of these methods.

The experimental values are taken from the report of Herrmann et al. [50] in which they measured the crystal structures under various temperatures using X-ray diffraction. The standard DFT calculations (denoted as PBE hereafter) without van der Waals corrections agree well with the previous theoretical prediction [27]. It shows that standard DFT calculations give poor predictions. For example, the volume of the unit cell is 6.8% larger than the experimental value. The volume of $\beta$-HMX from rung 1 vdW methods (PBE-D2), rung 2 vdW methods (PBE-D3, PBE-TS, RPBE-D3), and PBEsol calculations [51] has a significant improvement over standard DFT calculations. They have a RSD within 2% compared to the experimental value.

The van der Waals density functional methods are in rung 3. The explicitly tested methods are original vdW-DF [13], vdW-DF2 [14], optPBE-vdW [18], optB88-vdW [18], optB86b-vdW [18], vdW-DF2-C09 [19], LV-PW86r, vdW-DF-cx [20], rev-vdW-DF2 [18], screened-vdW [52], and vdW-DFq. It is worth pointing out that the gradient coefficient $Z_{ab}$ is −0.8491 and −1.887 in vdW-DF and vdW-DF2 methods, respectively [18]. One might expect that the more computing intensive vdW-DF models including the latest rev-vdW-DF2 method have higher accuracy as rung 3 vdW methods compared to rung 0, 1 and 2 methods. However, it turns out that the accuracy some of the vdW-DF functionals are even worse than standard PBE calculations as shown in Table 1. This implies that the exchange energy functional plays an important role in determining the geometry of semihard materials in vdW-DF methods. Only a carefully selected exchange energy functional together with the vdW-DF methods can predict the volume accurately.

It is also seen from the results listed in Table 1 that the predicted volume scatters over a wide range. One might want a way to monotonically tune the volume to that of the experimental volume. Using the one-parameter vdW-DFq method proposed here, one can smoothly tune the system value to the target one. It is worth noting that the only difference between the rev-vdW-DF2 and the vdW-DFq methods is the q value, which is fixed as 0.7114 in rev-vdW-DF2, but tunable in vdW-DFq from system to system, where a value of $q=1.05$ is good for $\beta$-HMX in this work, as illustrated in Figure 4. In addition, we have carried out the calculation using optB86q-GGA without vdW-DF2 corrections to examine the effect of vdW-DF2 corrections to the system’s volume. We find that the volume will be 19.75% larger than the experiment, clearly showing that vdW-DF2 corrections are important.

It is seen from Figure 4 that the lattice constants, lattice angle, and the volume can be tuned continuously at a very large range. For example, the volume varies from 82% to 200% when q changed from 0.1 to 100.0. Within the vicinity of the $q=1.0$, the volume varies linearly with respect to q. $q=1.05$ gives the best match to the experimental volume with RSD of −0.16%. As a comparison, the rev-vdW-DF2 method [18], the latest vdW-DF series method, predicted the volume as 496.905 Å${}^3$, 4.2% smaller than the experiment, which falls out of the satisfactory scope of 3%. At large deviation from 1.0, the relationship between the volume and the variable q is nonlinear as expected.

The large tunable range and linear behavior in the vicinity of the q corresponding to the experimental volume make the calibration of the vdW-DF much easier. This also indicates that our method is simple and robust in modeling various semihard materials where the vdW interactions are critical but the analytical formula is unknown a priori. We focus on the geometry study because the density is our main concern. One could tune q for other properties, such as for binding and cohesive energies. One needs to keep in mind that the tuning parameter q is not unique: there are quite a few parameters that affect the exchange energies. For example, one can also vary the gradient coefficient $Z_{ab}$ to modify the exchange energy functionals and it could tune the system’s volume as well. Our test shows in Table 1 that when $Z_{ab}$ = −0.8491 (vdW-DFk’) is used instead of −1.887 in vdW-DFq, the system volume will have a −4.9% change.

In the proposed vdW-DFq method, the $q=1.05$ is calibrated to $\beta$-HMX with the primitive unit cell volume. A general concern rises about the predictive power of this method. To that end, we further study the other semihard materials utilizing vdW-DFq method. Additional ten common structures of high explosives are investigated as $\alpha$ and $\delta$-HMX [53], $\alpha$ [54], $\beta$ [55], and $\gamma$-RDX [56], $\alpha$ (Monoclinic) and $\beta$ (Orthorhombic)-TNT [57], PETN [58], Tetryl [59], and TATB [60]. We focus on the mass density as it is the predominant factor for detonation density and pressure. The results of the lattice parameters, volume, mass density, and density RSD of these 10 structures are summarized in Table S1 in the Supplementary Information, compared with the predictions from vdW-DF2 and rev-vdW-DF2 methods as well as experiments. The relative standard deviations of mass densities are referred to the corresponding experiments. Particularly, we plotted the mass density RSD in Figure 5.

It shows that the recent developed rev-vdW-DF2 method is worse than conventional vdW-DF2 method in predicting the mass density of the energetic materials, and the proposed vdW-DFq method has the best performance of the examined three vdW-DF methods. The root mean square error of densities predicted from the vdW-DF2, rev-vdW-DF2, and vdW-DFq methods for the examined 11 molecular crystals is 1.7%, 3.7%, and 1.4%, respectively. These results indicate that vdW-DFq method with $q=1.05$ might be a good first-principles approach to study energetic materials. This success could be attributed to the similar electronic charge density gradients, because all these materials are composed of similar elements of carbon, nitrogen, oxygen, and hydrogen, and all have similar mass densities around 1.9 g/cm${}^3$.

In addition to energetic materials, we further study the mass density of the layered van der Waals crystals [61,62]. Graphite, hexagonal boron nitride (h-BN), and molybdenum disulfide (MoS${}_2$) are typical layered crystals. Their lattice parameters, volume, mass density, and density RSD are summarized in Table S1 in the Supplementary Information.

As shown in Figure 5, the vdW-DFq method prediction also has good agreements with experimental measurements. The performance of the vdW-DFq method is better than the vdW-DF2 and rev-vdW-DF2 methods in the layered van der Waals crystals. For the total 14 structures, the root mean square error of densities is 5.52%, 3.46%, and 1.37% for the vdW-DF2, rev-vdW-DF2, and vdW-DFq methods, respectively. The better performance of the vdW-DFq method over the vdW-DF2 and rev-vdW-DF2 methods in predicting the geometries of the 14 representative molecular crystals implies that the proposed vdW-DFq method could also have good performance in modeling other semihard materials. It is reasonable because the q value is calibrated and thus “localized” to semihard materials, as opposed to the other two methods which aim at a universal method for all materials.

It is worth mentioning that graphite, h-BN, and MoS${}_2$ are actually highly anisotropic materials. The success of our vdW-DFq method also suggests that, in some cases, finding correct equilibrium configurations does not necessarily imply finding the correct dispersion interactions [63], but a compromise of various approximations, which, however, vary from case to case. From this aspect, it evidences the essential role of q that is tunable to optimize such an accommodation. It might be expected that vdW-DFq method is a handy tool for various vdW systems on-demand.

4. Conclusions

In summary, we assess the performance of van der Waals (vdW) density functionals in predicting the geometry of semihard materials. We propose a one-parameter empirical van der Waals density functional named vdW-DFq to continuously tune the lattice constants by adjusting the enhancement factor of the exchange energy functionals. We illustrate that the lattice constants and volumes can be tuned over a wide range, and with a linear behavior in the vicinity of the q corresponding to the experimental volume. The application on $\beta$-HMX shows that our method is simple and robust in modeling various semihard materials. By one-parameter tuning, the vdW density functional can be well calibrated to the experimental value. This method provides a controllable approach to accurately model the mechanical, electrical, chemical, and other properties within the framework of van der Waals density functionals. Further study on the 10 energetic material molecular crystals and 3 layered van der Waals crystals shows that the proposed vdW-DFq method has better performance than the vdW-DF2 and rev-vdW-DF2 methods in predicting the geometries, which indicates that the vdW-DFq method could be good at modeling semihard materials. Our investigation might also inspire computational physicists to develop a better non-empirical vdW density functional. Nevertheless, our vdW-DFq method is a on-demand handy tool for various vdW systems by tuning q value.

5. Computational Method

We used a conventional unit cell with two HMX molecules. There are 56 atoms in total in one unit cell. The periodic boundary conditions are applied in three normal directions. The first-principles calculations based on density-functional theory (DFT) were carried out with the Vienna Ab-initio Simulation Package (VASP) [65]. The Kohn-Sham Density Functional Theory (KS-DFT) [1] is employed with the generalized gradient approximation as parameterized by Perdew, Burke, and Ernzerhof (PBE) for exchange-correlation functionals [30]. We have examined two additional exchange-correlation functionals: revised Perdew-Burke-Ernzerhof (RPBE) [66], and Perdew-Burke-Ernzerhof revised for solids (PBEsol) [44].

The electrons explicitly included in the calculations are the $1s^1$ for hydrogen atoms, $2s^{2}2p^2$ for carbon atoms, $2s^{2}2p^3$ for nitrogen atoms, and $2s^{2}2p^4$ for oxygen atoms. The core electrons are discribed by the projector augmented wave (PAW) and pseudo-potential approach [67,68]. A energy cutoff of 520 eV is applied in this study. All the vdW-DF calculations were done self-consistently using the efficient algorithm of Román-Pérez and Soler [33], implemented in VASP by Kliměs et al. [18]. It is worth noting that it is not evident a priori that self consistency should necessarily lead to higher accuracy: although self-consistent, vdW-DF and related methods still imply a number of approximations which may be way larger than self-consistency vdW effects. One of these is the lack of long-range many body effects which may amount to more than 10% of the total dispersion energy (see [69,70,71], and can lead to unexpected features in highly anisotropic systems [72,73].

The criterion electronic relaxation is 10${}^{-6}$ eV. The optimized atomic geometry was achieved when the forces on each atom is smaller than 0.001 eV/Å. The irreducible Brillouin Zone was sampled with a $5\times 3\times 4$ Gamma-centered k-mesh. The other 13 structures are optimized using similar parameters in DFT calculations.