Rehabilitation of PBE-GGA for Layered Materials
aa r X i v : . [ c ond - m a t . m t r l - s c i ] D ec Rehabilitation of PBE-GGA for Layered Materials
Haowei Peng ∗ and John P. Perdew
1, 2 Department of Physics, Temple University, Philadelphia, PA 19122, USA Department of Chemistry, Temple University, Philadelphia, PA 19122, USA
The structural and energetic properties of layered materials propose a challenge to density func-tional theory with common semilocal approximations to the exchange-correlation. By combining themost-widely used semilocal generalized gradient approximation (GGA), Perdew–Burke–Ernzerhof(PBE), with the revised Vydrov–Van Voorhis non-local correlation functional (rVV10), both ex-cellent structural and energetic properties of 28 layered materials were recovered with a judiciousparameter selection. We term the resulting functional as PBE+rVV10L with “L” denoting for lay-ered materials. Such combination is not new, and involves only refitting a single global parameter,however the resulting excellency suggests such corrected PBE for many aspects of theoretical stud-ies on layered materials. For comparison, we also present the results for PBE+rVV10 where theparameter is determined by the 22 interaction energies between molecules.
PACS numbers: 31.15.E-, 71.15.Mb, 71.15.Nc, 68.43.Bc
I. INTRODUCTION
Plenty of interests have been attracted by the two-dimensional (2D) materials and their parent layeredmaterials [1–4] since the experimentally realization ofgraphene [5]. Layered materials presented a huge chal-lenge for density functional theory (DFT) [6], the cur-rent work-horse first-principles method. The difficultycomes from the coexistence of inter–layer van der Waals(vdW) interaction and strong chemical bonding withinthe layer. The vdW-correction needed for layered mate-rials is usually weaker than that for molecular systems.Hence, most of vdW density functionals which are goodfor molecular systems [7–15] overbind layered materi-als significantly [16–19]. It is even not easy to find adispersion-corrected generalized gradient approximation(GGA), which is able to predict well for both the geomet-ric and energetic properties, i.e., the intra-layer latticeconstant a , inter-layer lattice constant c , and inter-layerbinding energy E b [18, 19].We have found a solution on the meta-GGA level,where we combined the strongly constrained appro-priately normed (SCAN) [20] meta-GGA and the re-vised Vydrov–Van Voorhis non-local correlation func-tional (rVV10) [14, 15] with one parameter adjusted tothe Ar binding curve. The so-termed SCAN+rVV10[21] functional not only gives the best description forlayered materials, but also excellently describes solids,molecular systems, adsorption of benzene on metal sur-faces, and hence we expected it to be a versatile vdWdensity functional. One important conceptual feature ofSCAN+rVV10, compared to other popular vdW densityfunctionals [22], is that we deliberately combine the non-local vdW correlation functional with a semilocal func-tional which already includes certain amount of inter-mediate range vdW binding from the exchange. Instead, ∗ [email protected]. The implementation of r VV10 in VASP will be available within the next version. previous work was usually pursuing a vdW-free exchangefunctional.Inspired by this new concept, we revisit the request ofa GGA-based vdW density functional for layered ma-terials, and end up with a solution by combining themost commonly used Perdew-Burke-Ernzerhof (PBE)[23] GGA and rVV10. The resulting functional, termedas PBE+rVV10L with the “L” denoting for layeredmaterials, achieves similar accuracy as SCAN+rVV10.PBE+rVV10L is not as versatile as its meta-GGA coun-terpart SCAN+rVV10, however it is noticeably cheaperin computation, and numerically more stable thanks tothe much simpler mathematical form of PBE. Besides,the rVV10 and PBE have been implemented in manyab-inito codes, and hence PBE+rVV10L provides a veryhandy solution for many problems related to layered ma-terials. PBE+rVV10L is even better than the AM05-VV10sol functional [14, 18, 24], which is constructed ina similar way as here but with an additional parameteradjusted (unfortunately a physically-sound justificationhas not been provided yet). Combining PBE and rVV10(or VV10) is new in the condensed matter physics com-munity, but is not in the quantum chemistry commu-nity where the PBE+VV10 has already been proposedfor molecular systems [25]. In this work, we will re-port the benchmarking results of this newly proposedPBE+rVV10L, and compared with the PBE+rVV10where the parameter is adjusted to the interaction en-ergies of 22 molecular complexes (S22) [26, 27] as theoriginal VV10 [14] and rVV10 [15]. All calculations inthis work were performed with the projector augmentedwave (PAW) method [28] as implemented in the VASPcode (version 5.4.1) [29–31]. For more details, we referto the Appendix in Ref. [21].The r VV10 [14, 15] nonlocal correlation functionaltakes a similar form as the popular family of Rutgers-Chalmers vdW-DFs [7–12], E nlc = Z d r n ( r )[ ¯ h Z d r ′ Φ( r , r ′ ) n ( r ′ ) + β ] , (1) TABLE I. Layer–layer binding energy E b in meV/˚A , inter–layer lattice constant c in ˚A, and intra–layer lattice constant a in ˚A for 28 layered materials from SCAN+rVV10 [21] andPBE+rVV10L. The mean error (ME) and mean absolute er-ror (MAE) are also given in the same units, and the mean rel-ative error (MRE) and mean absolute relative errors (MARE)are given in percentage. The reference values for E b are fromRPA calculations, and from experiments for c and a [17, 19]. Reference SCAN+ r VV10 PBE+ r VV10L E b c a E b c a E b c a TiS h -BN 14.49 6.69 2.51 18.45 6.48 2.50 14.43 6.85 2.51PbO 20.25 5.00 3.96 22.93 4.81 3.98 17.95 5.08 4.04ME 0.20 0.02 -0.01 -0.02 0.15 0.00MAE 1.48 0.08 0.02 1.74 0.15 0.02MRE 1.7 0.2 -0.4 0.2 1.8 0.0MARE 7.7 1.2 0.5 8.9 1.9 0.7 β vanishes for the Rutgers-Chalmers vdW-DFs, and thetotal exchange correlation functional reads E xc = E xc + E nlc . (2)Here n ( r ) is the electron density, and Φ( r , r ′ ) is the ker-nel describing the density-density interactions, E xc is thecompanying semilocal exchange correlation. To ensurezero E nlc for the uniform electron gas, β = ( b ) inHartree is required. Two empirical parameters C and b appear in the kernel Φ( r , r ′ ): C chosen for accurate − C /R vdW interactions between molecules at largeseparation R , and b controlling the damping of E nlc atshort range.In the original form for both VV10 and rVV10 [14, 15],the nonlocal correlation was proposed to combine withthe semilocal exchange-correlation functional [23, 32] E xc = E rP W x + E P BEc , partly because of the rPW86exchange is nearly vdW-free [32]. For a semilocal E xc , C = 0 . b parameter was determined as 5 . . r VV10. Increasing C or b generally results insmaller vdW correction. Keeping the original semilocalpart, b = 9 .
15 is required to fit the binding energies of26 layered materials for both VV10 [18] and rVV10 [21],implying that layered materials require weaker vdW cor-rection than molecular complexes. However such fitted b value leads to much worse performance for both theintra- and inter-layer lattice constants in layered materi-als, and also for solids [18]. Besides, b = 9 . E xc to the SCAN meta-GGA[20] results in the versatile SCAN+rVV10 with b = 15 . b = 14 . E xc to the AM05 form [24] results in the AM05-VV10solwith b = 10 .
25 and C = 10 − [18], which works wellfor layered materials but not for S22. Note the practi-cally zero value for the C parameter is required by thefitting but not by chemistry. For PBE+VV10, b = 6 . b = 10 . b = 6 . b value for PBE+rVV10 is slightly larger than thatfor the original rVV10, in accordance with the slightlymore vdW binding from the PBE exchange than thatfrom the rPW86 exchange [32]. This difference is re-lated to that the exchange enhancement factor of PBEis bounded by the Lieb-Oxford constraintas as the re-duced density gradient s increases [23], while the rPW86enhancement factor diverges as s . [32]. The referencebinding energies of the layered materials are not fromexperiments, but from adiabatic-connection fluctuation-dissipation theorem within the random-phase approxi-mation (RPA) [35–37], which are yet the best availablechoice.The lattice constants, both intra-layer a and inter-layer c , and the layer–layer binding energy E b are the most fun-damental quantities when one embarks on first-principlescomputation of layered materials. For the benchmarking,we use the binding energies from RPA, and lattice con-stants from experiments as references [17, 19]. Until now,SCAN+rVV10 is the only one which can simultaneouslypredict with the mean absolute relative error <
10% for E b , <
2% for c , and <
1% for a [21]. In Table I, wecompare the results from PBE+rVV10L to those fromSCAN+rVV10, and the reference values. PBE+rVV10Lachieves excellent accuracy for both the geometrical andenergetic properties by adjusting only one parameter,which is not trivial as discussed above. PBE+rVV10L ac-tually is comparable with SCAN+rVV10 for this class ofmaterials, with a slightly overestimated layer-layer spac-ing. Nevertheless, PBE+rVV10L is another member ofthe “10-2-1” club for layered materials. Considering theextremely simply mathematical form of both PBE andrVV10, and their widely availability in many scientificcodes, PBE+rVV10L can be a very handy theoreticaltool for layered materials. It can be also used to preparereasonably good initial (relaxed) structure and orbitalsto accelerate the convergence of following SCAN+rVV10calculations. FIG. 1. (Color online) Box-plots for the absolute relativeerrors and relative errors of the inter-layer binding ener-gies, inter- and intra-layer lattice constants ( c and a ) fromSCAN+rVV10 [21], PBE+rVV10, and PBE+rVV10L, for 28layered materials. The reference values are from RPA for thebinding energy, and from experiment for the lattice constants[17, 19]. The Tukey box-plot used here summarizes the overalldistribution of a set of data points: The bottom and the topof the box are the first (Q1) and third (Q3) quartiles (25%of data points lies below Q1, and another 25% above Q3);The band inside the box denotes the median; The circles ifany denote outliers which lie further than 1.5 ∗| Q − Q | awayfrom the box; The vertical line extends from the minimum tothe maximum, except for the outliers. Besides, we also denotethe mean value with a shape inside the box. Abs. Rel. Err. for E b Rel. Err. for E b Abs. Rel. Err. for c Rel. Err. for c Abs. Rel. Err. for a Rel. Err. for a P e r c en t age SCAN+rVV10 PBE+rVV10 PBE+rVV10L
In Fig. 1, we summarise the absolute relative errors and relative errors for a , c and E b from SCAN + rV V P BE + rV V
10, and
P BE + rV V L . With a smaller b parameter, PBE+rVV10 overbinds the layered materi-als by about 50%, similar to the original rVV10, but thelattice constants are still reasonably accurate. Therefore,a b parameter between 6 . . FIG. 2. (Color online) Box-plots for the absolute relativeerrors and relative errors of the atomization energies E a ,and the lattice equilibrium volumes V , from RPA, PBE,PBE+rVV10, and PBE+ r VV10L for 50 solids, with respectto the experimental values. The RPA, PBE, and experimen-tal values (after the zero-point correction) are from Refs. 38and 39.
Abs. Rel. Err. for E a Rel. Err. for E a Abs. Rel. Err. for V Rel. Err. for V P e r c en t age RPA PBE PBE+rVV10 PBE+rVV10L
For solid systems, we further benchmark the perfor-mance of PBE+rVV10L and PBE+rVV10 for 50 solidscompiled as in Ref. 21, which includes ( i ) 13 group–IVand III–V semiconductors, ( ii ) 5 insulators, ( iii ) 8 main-group metals, ( iv ) 3 ferromagnetic transition metals Fe,Co, and Ni, and ( v ) 21 other transition metals for whichnon-spin-polarized calculations were performed [38, 39].We compared the mean relative errors and mean abso-lute errors for atomization energies and lattice volumes.The rVV10 correction decreases the mean absolute rela-tive error for the atomization energies slightly for PBE,and the otherwise underestimated atomization energiesare slightly overestimated now with both PBE+rVV10and PBE+rVV10L. However, atomization energy maynot be a good choice to assess a semilocal functional[40]. It is well-known that PBE overestimate the lat-tice volume, with a mean absolute relative error over3% as shown in Fig. 2. The attractive vdW correc-tion slightly remedies this systematic overestimation byabout 1%. This has been already known in recent worksof Tao et al. [41, 42]. Overall, the structure and ener-getic properties for solids are not skewed by the rVV10correction with both PBE+rVV10 and PBE+rVV10L. FIG. 3. (Color online) Box-plots for the absolute relativeerrors and relative errors of the interaction energies from r VV10, SCAN+rVV10, PBE+rVV10, PBE+rVV10L, withrespect to the CCSD(T) results [26, 27], for the molecularcomplexes in the S22 dataset.
Abs. Rel. Err. Rel. Err. P e r c en t age rVV10 SCAN+rVV10 PBE+rVV10 PBE+rVV10L Fig. 3 shows the results of PBE+rVV10 andPBE+rVV10L for S22 which includes seven hydrogen-bonded, eight dispersion-bound, and seven mixed com-plexes, compared to the rVV10, SCAN+rVV10, and theCCSD(T) reference [26, 27]. Similar to PBE+VV10[25], the fitting of PBE+rVV10 is less accurate thanthe original rVV10 with the rPW86 exchange, and themean absolute relative error of 6% (4.5% for rVV10).PBE+rVV10L with the b parameter fitted to layered ma-terials significantly underbinds with a mean absolute rel-ative error of 24% and a mean relative error of 22%. Nev-ertheless, PBE+rVV10L is noticeably better than AM05- VV10sol, whose mean absolute relative error is 36% [18].Besides, PBE+rVV10L performs very well for the sevenhydrogen-bonding complexes with a mean absolute rela-tive error of only 2%. This indicates that PBE+rVV10Lshould be better than PBE+rVV10 for structural prop-erties of water [33].At last, we benchmarked the performance ofPBE+rVV10 and PBE+rVV10L with the adsorption ofBenzene ring on Cu, Ag and Au (111) surfaces, whichhave been widely studied [43–48]. The SCAN+rVV10results [21] are chosen for reference, which agrees verywell with available experiments [49–54]. In these sys-tems, PBE+rVV10L is slightly worse than PBE+rVV10,underestimating the binding energy ∆ E by about 0.2 eVand overestimating the distance ∆ z between benzene and TABLE II. Adsorption energy E ad and distance ∆ z be-tween benzene and the (111) surface of Cu, Ag, and Aufrom PBE+rVV10 and PBE+ r VV10L, compared with theSCAN+rVV10 results [21]. The data for the lowest-energyhcp30 ◦ configuration [43] is shown.PBE+rVV10 PBE+ r VV10L SCAN+rVV10 E ad (eV) ∆ z (˚A) E ad (eV) ∆ z (˚A) E ad (eV) ∆ z (˚A)Cu 0 .
84 2 .
88 0 .
52 3 .
05 0 .
74 2 . .
74 3 .
02 0 .
46 3 .
18 0 .
68 3 . .
82 3 .
04 0 .
51 3 .
20 0 .
73 3 . metal surface by about 0.14 ˚A. PBE+rVV10 is betterthan PBE+rVV10L, and only overbinds comparing toSCAN+rVV10 very slightly.In conclusion, we provide here two set of parameters forthe combination between PBE and rVV10. For systemsinvolving molecules, the PBE+rVV10, where b = 6 . b = 10 .
0, achieves the accuracy of the bestdispersion-corrected, semilocal density functional, andhence is highly recommended. Values between thesetwo may also be employed for specific problems. ThePBE+VV10L and PBE+VV10 is not as versatile as themeta-GGA-level SCAN+rVV10, but they are very handyand computationally high-efficiency alternatives. [1] Q. Wang, K. Kalantar-Zadeh, A. Kis, J. Coleman, andM. S. Strano, Nat. Nanotechnol. , 699 (2012).[2] M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh,and H. Zhang, Nat. Chem. , 263 (2013).[3] M. Xu, T. Liang, M. Shi, and H. Chen,Chem. Rev. , 3766 (2013).[4] S. Butler, S. Hollen, L. Cao, and Y. Cui,ACS Nano , 2898 (2013).[5] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 666 (2004).[6] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). [7] K. Lee, E. Murray, L. Kong, B. I. Lundqvist,D. C. Langreth, and ´E. D. Murray,Phys. Rev. B , 081101 (2010).[8] V. R. Cooper, Phys. Rev. B , 161104 (2010).[9] J. Klimeˇs, D. R. Bowler, and A. Michaelides,J. Phys.: Condens. Matter. , 022201 (2010).[10] J. Klimeˇs, D. R. Bowler, and A. Michaelides,Phys. Rev. B , 195131 (2011).[11] I. Hamada, Phys. Rev. B , 121103 (2014).[12] K. Berland and P. Hyldgaard,Phys. Rev. B , 035412 (2014).[13] O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. , 063004 (2009).[14] O. A. Vydrov and T. Van Voorhis,J. Chem. Phys. , 244103 (2010).[15] R. Sabatini, T. Gorni, and S. De Gironcoli,Phys. Rev. B , 041108 (2013).[16] T. Bj¨orkman, A. Gulans, A. V.Krasheninnikov, and R. M. Nieminen,J. Phys.: Condens. Matter. , 424218 (2012).[17] T. Bj¨orkman, A. Gulans, A. V. Krasheninnikov, andR. M. Nieminen, Phys. Rev. Lett. , 235502 (2012).[18] T. Bj¨orkman, Phys. Rev. B , 165109 (2012).[19] T. Bj¨orkman, J. Chem. Phys. , 074708 (2014).[20] J. Sun, A. Ruzsinszky, and J. P. Perdew,Phys. Rev. Lett. , 036402 (2015).[21] H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun,Phys. Rev. X , 041005 (2016).[22] K. Berland, V. R. Cooper, K. Lee, E. Schr¨oder,T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist,Reports Prog. Phys. , 066501 (2015).[23] J. P. Perdew, K. Burke, and M. Ernzerhof,Phys. Rev. Lett. , 3865 (1996).[24] R. Armiento and A. E. Mattsson,Phys. Rev. B , 1 (2005).[25] J. Arag´o, E. Ort´ı, and J. C. Sancho-Garc´ıa,J. Chem. Theory Comput. , 3437 (2013).[26] R. Podeszwa, K. Patkowski, and K. Szalewicz,Phys. Chem. Chem. Phys. , 5974 (2010).[27] T. Takatani, E. G. Hohenstein, M. Malagoli,M. S. Marshall, and C. D. Sherrill,J. Chem. Phys. , 144104 (2010).[28] P. Bl¨ochl, Phys. Rev. B , 17953 (1994).[29] G. Kresse and J. Hafner, Phys. Rev. B , 14251 (1994).[30] G. Kresse and J. Furthm¨uller,Phys. Rev. B , 11169 (1996).[31] G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999).[32] ´E. D. Murray, K. Lee, and D. C. Langreth,J. Chem. Theory Comput. , 2754 (2009).[33] G. Miceli, S. de Gironcoli, and A. Pasquarello,J. Chem. Phys. , 034501 (2015).[34] J. G. Brandenburg, J. E. Bates, A. Ruzsinszky, J. Sun,and J. P. Perdew, Phys. Rev. B , 115144 (2016).[35] J. Harl and G. Kresse, Phys. Rev. Lett. , 056401 (2009).[36] J. Harl, L. Schimka, and G. Kresse,Phys. Rev. B , 115126 (2010).[37] H. Eshuis, J. E. Bates, and F. Furche,Theor. Chem. Acc. , 1084 (2012).[38] J. Harl, L. Schimka, and G. Kresse,Phys. Rev. B , 115126 (2010).[39] L. Schimka, R. Gaudoin, J. Klimeˇs, M. Marsman, andG. Kresse, Phys. Rev. B , 214102 (2013).[40] J. P. Perdew, J. Sun, A. Ruzsin-szky, P. D. Mezei, and G. I. Csonka,Period. Polytech. Chem. Eng. , 2 (2015).[41] J. Tao, J. P. Perdew, and A. Ruzsinszky,Phys. Rev. B , 233102 (2010).[42] J. Tao, F. Zheng, F. Wang, S. Liu, J. P. Perdew, andA. M. Rappe, (in preparation).[43] W. Liu, V. G. Ruiz, G. X. Zhang, B. Santra,X. Ren, M. Scheffler, and A. Tkatchenko,New J. Phys. , 053046 (2013).[44] A. Bili´c, J. R. Reimers, N. S. Hush, R. C. Hoft, andM. J. Ford, J. Chem. Theory Comput. , 1093 (2006).[45] K. Toyoda, Y. Nakano, I. Hamada, K. Lee, S. Yanagi-sawa, and Y. Morikawa, Surf. Sci. , 2912 (2009).[46] H. Yildirim, T. Greber, and A. Kara,J. Phys. Chem. C , 20572 (2013).[47] W. Reckien, M. Eggers, and T. Bredow,Beilstein J. Org. Chem. , 1775 (2014).[48] D. J. Carter and A. L. Rohl,J. Comput. Chem. , 2263 (2014).[49] W. Liu, F. Maaß, M. Willenbockel, C. Bronner,M. Schulze, S. Soubatch, F. S. Tautz, P. Tegeder, andA. Tkatchenko, Phys. Rev. Lett. , 036104 (2015).[50] C. T. Campbell and J. R. V. Sellers,J. Am. Chem. Soc. , 18109 (2012).[51] M. Xi, M. X. Yang, S. K. Jo, B. E. Bent, and P. Stevens,J. Chem. Phys. , 9122 (1994).[52] X.-L. Zhou, M. Castro, and J. White,Surf. Sci. , 215 (1990).[53] D. Syomin, J. Kim, B. E. Koel, and G. B. Ellison,J. Phys. Chem. B , 8387 (2001).[54] L. Ferrighi, G. K. H. Madsen, and B. Hammer,J. Chem. Phys.135