Introduction
Bulk hBN (30-50nm thick, for example) is often used in modern research to make heterostructures and its lattice constants are very useful information to researchers in considering lattice mismatch and strains. Hereby density functional theory (DFT) calculations are performed to numerically solve for the optimal lattice constants (a and c) of hBN, which are then compared with the experimental values. The experimental optimal values for a/b axis are 2.502Å, and for c it is 6.617Å [1]. The specific code package used to perform the calculations is the Cambridge Serial Total Energy Package (CASTEP), which uses a plane-wave basis set. The functional used to perform the calculations is GGA-PBE. hBN is not a metal and does not have spin-orbit coupling or spin-polarization. On the fly generated (OTFG) ultrasoft pseudopotential is used for the calculations.
An illustration of the hBN unit cell
Calculation Results and Discussions
The convergence test is initially performed using a=2.50Å and c=6.62Å and then calculations with large step size are performed. Surprisingly, the energy as a function of c is minimized at around 8.55Å, which is much larger than the experimental value (6.62Å).
It turns out that it is not considering the dispersion correction that caused this problem. Noncovalent forces such as van der Waals interaction play a crucial role in the formation of layered materials such as hBN. However, the noncovalent forces are not easy to be accounted for in calculations and are ignored in the local density approximation (LDA) and generalized gradient approximation (GGA) functionals. Therefore when performing calculations without dispersion corrections, the interlayer interaction is weaker due to the absence of van der Waals interaction and the preferred value of c is larger. Luckily, many semiempirical solutions have been developed regarding this problem. To better demonstrate the importance of the dispersion corrections, here the calculation results are reported in two parts, with the first part not using the TS method (one semiempirical dispersion correction that works well for GGA-PBE functional) and having preferred c around 8.55Å as result and the second part using the TS method and having preferred c close to the experimental value as result.
I. Finding the optimal value of c without dispersion corrections
Since the results from calculations without dispersion corrections are not precise in the first place, the report of the convergence test is skipped here. A k-point mesh of 9 by 9 by 4 and an energy cutoff of 600 eV should be fine enough for the purpose of demonstrations. With a=2.50Å and a step size of 0.3Å, the optimal value of c is found to be around 8.55Å with an error range of 0.15Å (since the smallest two energies are with c=8.4Å and c=8.7Å, the optimal value of c must be within the interval from 8.4Å to 8.7Å), which is a very coarse calculation in terms of error range, but it is good enough to show that the optimal value of c without the consideration of dispersion corrections is much larger than the experimental value.
The total energy versus the value of c
II. Finding the optimal values of a and c with dispersion corrections
II. i. Convergence test
Using a=2.50Å and c=6.62Å, the convergence test is performed on k-points and the energy cutoff. The TS method is used for dispersion correction. SCF tolerance is chosen to be medium. The valance electron configuration of B is 2s2 2p1, with partial core correction Rc = 0.838 Bohr radius. The valance electron configuration of N is 2s2 2p3, with partial core correction Rc = 0.769 Bohr radius.
The convergence test on k-points is first performed. An energy cutoff of 600eV should be enough for the purpose of this test. It turns out that within the precision of 0.0001eV, the energy converges very quickly and the choice of the k-points mesh can be rather arbitrary, with the choices of 8/9/10 by 8/9/10 by 3/4/5. A mesh of 9 by 9 by 4 is chosen for later calculations.
Then a proper energy cutoff needs to be decided. A series of calculations is performed, with the cutoff energy begins at 400eV and increases with steps of 50eV. Since the energy monotonically decreases with the increase of the cutoff energy, the total energy will never converge. But a cutoff that is good enough can be determined by calculating the difference between each total energy. And the magnitude of the difference in total energy monotonically decreases with the increase of the energy cutoff. In this way, the change in energy versus the cutoff energy can be plotted and 800eV is good enough for later calculations since the change in total energy is only -0.0004eV when cutoff energy increases from 750eV to 800eV.
The change in total energy due to the increase of cutoff energy versus the cutoff energy
II. ii. Finding the optimal value of c coarsely
With the value of a/b axis fixed at 2.50Å and a step size of 0.20Å for c, the optimal value of c is coarsely found to be around 6.70Å with an error range of 0.10Å, which is much closer to the experimental value (6.617Å) compared with the optimal value from the calculations without the consideration of the dispersion corrections (8.55Å). The plot of the total energy versus the value of c is as follows:
The total energy E (eV) versus the value of c (Å). The smallest two energies are found at c=6.60Å and c=6.80Å. Therefore the optimal c must be within that interval.
II. iii. Finding the optimal values of a/b and c finely
The step size in a/b direction and in c direction is 0.02Å and 0.05Å respectively.
| 6.50 | 6.55 | 6.60 | 6.65 | 6.70 | |
|---|---|---|---|---|---|
| 2.46 | -740.1006 | -740.1011 | -740.1012 | -740.1006 | -740.0996 |
| 2.48 | -740.1573 | -740.1580 | -740.1581 | -740.1578 | -740.1568 |
| 2.50 | -740.1807 | -740.1815 | -740.1818 | -740.1816 | -740.1808 |
| 2.52 | -740.1729 | -740.1738 | -740.1743 | -740.1741 | -740.1735 |
| 2.54 | -740.1358 | -740.1368 | -740.1373 | -740.1373 | -740.1368 |
The values in the first column are the values of a/b used in the corresponding calculations, and the values in the first row are the values of c used in the corresponding calculations. The lower right 5 by 5 corner of the table are the corresponding calculated energy results.
It can be seen from the table that for each different c, the optimal value of a/b is always 2.51Å plus or minus 0.01Å. This agrees really well with the experimental value (2.502Å), with only 0.32% off.
However, for different values of a/b, the optimal c is different. For a/b=2.46Å and 2.48Å, the optimal c is 6.575Å plus or minus 0.025Å, which is 0.63% off the experimental value (6.617Å). And for a/b=2.50Å, 2.52Å, and 2.54Å, the optimal c is 6.625Å plus or minus 0.025Å, which is 0.12% off the experimental value (6.617Å).
0.63% and 0.12% are both reasonably good. But it is not hard to notice that the total energy is more sensitive to the change in a/b axis than the change in c axis, with ΔE due to 1% change in a/b is in the order of 0.01eV whereas ΔE due to 1% change in c is in the order of 0.001eV. A possible explanation is that the semiempirical dispersion correction potential is proportional to 
The energy versus a/b at different values of c. It is easy to read from the plot that the change in energy due to a/b is much larger than the change in energy due to c which is barely readable from the plot.
Conclusion
The optimal lattice constants given by the CASTEP calculations are a/b=2.51Å plus or minus 0.01Å which is 0.32% off the experimental value (2.502Å), and c=6.625Å plus or minus 0.025Å which is 0.12% off the experimental value (6.617Å). The numerical results agree well with the experimental values and provide a good demonstration of how powerful DFT calculations can be.
References
- http://www.hqgraphene.com/h-BN-CAN1.php
