WWW.C-CHEM.ORG
FULL PAPER
Local Response Dispersion Method in Periodic Systems: Implementation and Assessment Yasuhiro Ikabata,[a] Yusuke Tsukamoto,[a] Yutaka Imamura,[a]† and Hiromi Nakai*[a,b,c,d] We report the implementation of the local response dispersion (LRD) method in an electronic structure program package aimed at periodic systems and an assessment combined with the Perdew–Burke–Ernzerhof (PBE) functional and its revised version (revPBE). The real-space numerical integration was implemented and performed exploiting the electron distribution given by the plane-wave basis set. The dispersion-
corrected density functionals revPBE1LRD was found to be suitable for reproducing energetics, structures, and electron distributions in simple substances, molecular crystals, and C 2014 Wiley Periodicals, Inc. physical adsorptions. V
Introduction
sion (MBD) method[34] and investigated the relation between pairwise methods and the MBD method.[35] Sato and Nakai proposed the local response dispersion (LRD) method in 2009 as a density-dependent dispersion correction scheme.[36] This method calculates dispersion coefficients and energies on the basis of effective atomic polarizabilities computed using the electron density and its gradient. Although the LRD method is one of the pairwise correction methods, the dispersion energy is given as a density functional. Similar to conventional XC functionals and nonlocal correlation functionals, the LRD method does not depend on any experimental values explicitly. Furthermore, the extension of the LRD method to multicenter interactions was achieved,[37] followed by the formulation of the self-consistent field (SCF) treatment and the analytic energy gradient.[38]
Density functional theory (DFT) is a standard electronic structure methodology used in periodic systems. From the establishment of Kohn–Sham DFT[1,2] to the present, the local density approximation (LDA)[3,4] to the exchange-correlation (XC) functional has been employed. The generalized gradient approximation (GGA)[5–8] and meta-GGA[9,10] were proposed partly to include the nonlocal XC effect. Hybrid functionals,[11–13] which are theoretically justified by the generalized Kohn–Sham scheme,[14] handle further nonlocal exchange by mixing the Hartree–Fock exchange. Although these improvements have increased the accuracy and availability of DFT, it is well known that this is not the case with noncovalent interactions. LDA evaluates the van der Waals force qualitatively, although the attractive force originates in the inaccuracy of the exchange functional.[15] To make the matter worse, GGA and most later XC functionals cannot describe van der Waals interaction even qualitatively.[16] This failure originates in the inaccurate description of dispersion interactions. The development of a dispersion correction technique is an increasingly important topic in DFT.[17,18] Because the dispersion interaction is derived from the nonlocal correlation effect, several nonlocal correlation functionals have been proposed. One example, proposed by Dobson and Dinte,[19] is equivalent to the Andersson–Langreth–Lundqvist (ALL) functional.[20] More sophisticated nonlocal correlation functionals developed later were included in vdW-DF series[21–23] and XC functionals by Vydrov and Van Voorhis.[24,25] Conversely, a number of pairwise correction methods that compute dispersion energy as a sum of atom–atom contributions, have been proposed. The DFT-D method[26–30] has a conceptual simplicity, and various combinations with XC functionals have been examined.[31] Instead of predetermined dispersion coefficients, some methods introduced a dependency on the electronic state, for example, the exchange-dipole moment (XDM) scheme[32] and the vdW-TS method.[33] Recently, Tkatchenko et al. proposed the many-body disper-
DOI: 10.1002/jcc.23807
[a] Y. Ikabata, Y. Tsukamoto, Y. Imamura, H. Nakai Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo, 1698555, Japan E-mail: [email protected] [b] H. Nakai Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo, 169-8555, Japan [c] H. Nakai CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan [d] H. Nakai Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Katsura, Kyoto, 615–8520, Japan † Present address: Advanced Institute for Computational Science, RIKEN, 7-1-26, Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan. Contract grant sponsor: Ministry of Education, Culture, Sports, Science and Technology [MEXT; Strategic Programs for Innovative Research (SPIRE) and the Computational Materials Science Initiative (CMSI)]; Contract grant sponsor: Elements Strategy Initiative for Catalyst & Batteries (ESICB) [MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012)]; Contract grant sponsor: Japan Science and Technology Agency [JST; Core Research for Evolutional Science and Technology (CREST) program “Theoretical Design of Materials with Innovative Functions Based on Relativistic Electronic Theory”] C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
1
FULL PAPER
WWW.C-CHEM.ORG
Furthermore, the combination of the LRD method with the long-range corrected (LC) functional[39] was found to be highly accurate not only for conventional molecular dimers[40] but also for open-shell[41] and excited-state molecular complexes.[42] These developments and assessments treated isolated molecular complexes with Gaussian basis functions. To treat molecular solids and physical adsorption, the extension of the LRD method to periodic systems is of vital importance. In this study, the LRD method is implemented in a package based on the periodic boundary condition and the plane-wave basis functions. We here combine the LRD method with the Perdew–Burke–Ernzerhof functional (PBE)[7] and its revised version (revPBE),[8] which are representatives of pure GGA functionals. Simple substances, molecular crystals, and physical adsorptions of n-butane on transition metals are considered in the assessment of these dispersion-corrected functionals.
jrqj 6q
pffiffiffiffiffiffiffiffi 4pq : kF
(3)
Equation (3) is not implemented in the LRD method with Gaussian basis functions while it is adopted in practical calculations of the ALL functional.[20] The present study requires a cutoff because the plane-wave basis generates small electron densities in the region far from the nucleus, which sometimes causes serious errors in the polarizability. Becke’s partition function[47] was implemented as wa and was computed considering all atoms in 3 3 3 3 3 unit cells. On the basis of the polarizability and geometry, the nthorder (n 5 6, 8, and 10) interatomic dispersion coefficient is calculated as Cnab 5
ð1 1 X X X ab duaal m ;l m ’ ðiuÞabl m ;l m ’ ðiuÞ; Sl1 m1 ;l2 m2 Sab l1 m1 ’;l2 m2 ’ 1 1 1 1 2 2 2 2 2p l l m1 m2 0 1 2 m1 ’m2 ’
(4)
Theory and Implementation We implemented the LRD method in the Quantum ESPRESSO package,[43] following the first paper[36] in which pairwise dispersion energy was computed in the post-SCF manner. Our previous study indicates that the SCF treatment brings minor effect.[38] The procedure to calculate dispersion correction energy is explained in this section. The energy gradient required for the lattice structure optimization is calculated by simply differentiating the LRD energy with respect to nuclear coordinates. The differentiation of total energy with respect to density is neglected. First, atomic polarizabilities are computed using the electron density and its gradient, obtained by the conventional DFT calculation. The polarizability is obtained from the intermolecular dispersion energy by the second-order perturbation theory, the multicenter-multipole expansion of the Coulomb operator,[44] and the local response approximation to the density response function.[19] As a result, the polarizability of atom a at an imaginary frequency iu is given by ð qðrÞ aalm;l’m’ ðiuÞ5 drwa2 ðrÞ 2 rRlm ðr2Ra Þ rRl’m’ ðr2Ra Þ; (1) x0 ðrÞ1u2 where wa is the partition function of atom a, q is the electron density, and R is the solid harmonics. l and m are integers with the restrictions l 1 and – l m l. l 6¼ l0 terms are neglected in the following because of their minute numerical impact. x0 in eq. (1) is a dispersion relation of a dielectric medium in the long-wavelength limit. The LRD method adopts following gradient-corrected form[45]: k2 ð11ks2 Þ2 x0 ðrÞ5 F ; (2) 3 2 1=3 with the local Fermi wave vector kF 5ð3p qÞ , an empirical parameter k, and the reduced density gradient s5jrqj=ð2kF qÞ. The polarizabilities of atoms in one unit cell are numerically integrated using the Cartesian grid. We consider 3 3 3 3 3 unit cells for real-space integration with a cutoff criteria proposed by Rapcewicz and Ashcroft,[46] 2
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
where Sab l1 m1 ;l2 m2 is the angular factor of the multipole expansion defined by l1 Sab l1 m1 ;l2 m2 5ð21Þ
ð2l1 12l2 11Þ! 1=2 dL;l1 1l2 3 ð2l1 Þ!ð2l2 Þ!
l1 m1
l2
L
m2
M
! ^ ab Þ: CL;M ðR
(5) The 2 3 3 matrix is Wigner’s 3-j symbol, and CL;M is the spherical harmonics with Racah’s normalization factor evaluated at ^ ab 5Rab =Rab. Note that the the polar angles of the unit vector R order is related to l1 and l2 by n 5 2(l1 1 l2 1 1). The integration with respect to the imaginary frequency in eq. (4) is possible pffiffiffiffiffiffiffiffiffiffiffi by transforming the variable u to v5u= 11u2 and applying Gauss–Chevyshev’s rule. The modified atomic polarizability is given by a alm;lm’ ðvk Þ ð12vk 2 Þ21=2 aalm;lm’ ðvk Þ ð ð12v2 Þ1=2 qðrÞ 5 drwa2 ðrÞ 3rRlm ðr2Ra ÞrRlm’ ðr2Ra Þ: ð12v 2 Þx20 ðrÞ1v 2
(6)
The dispersion coefficient is computed using the modified atomic polarizabilities: Cnab 5
N X 1 X X X ab a al m ;l m ’ ðvk Þ a bl m ;l m ’ ðvk Þ; Sl1 m1 ;l2 m2 Sab l1 m1 ’;l2 m2 ’ 1 1 1 1 2 2 2 2 4N l l m1 m2 k51 1 2 m1 ’m2 ’
(7) with N quadrature points at vk 5cos ½ð2k21Þp=4N. Twelve points were confirmed to be sufficient to maintain numerical accuracy. Finally, the dispersion correction energy for the unit cell is calculated. ELRD 52
1 X X ab 2n ðnÞ C R f ðRab Þ 2 n56;8;10 a;b n ab damp
(8)
We consider all atomic pairs in one unit cell and atoms in other unit cells containing at least one atom that is within 20 WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Bohr from the atom in the focused unit cell. The damping function is " ðnÞ fdamp 5exp
26 # Rab 2m ; n52m14; R
(9)
where R is interpreted as the sum of van der Waals radii with respect to atoms a and b. Assuming a linear relationship between atomic radius and polarizability, R is written by R5j
n 1=3 b 1=3 o 1 aeff ðidÞ 1R0 ; aaeff ðidÞ
(10)
with two empirical parameters j and R0. a a is the effective static polarizability,[36] defined by ð aaeff ðidÞ5 drwa ðrÞ
qðrÞ : x20 ðrÞ1d2
(11)
Note that d is set to 0.01 a.u. to maintain numerical stability.
Illustrative Applications Simple substances First, we calculated the atomization energies of simple substances. The number of elements treated in this subsection is 47. We excluded elements that form molecules in the most stable solid phase. Our aim in this subsection is to investigate the effect of the dispersion correction in metal- or covalentbonded periodic systems after fitting empirical parameters to rare-gas elements. The lattice parameters were set to experimental values.[48] The atomization energy is given as the difference in atomic energies between the solid and dissociated states. As for the dissociated state, we set one atom in a large unit cell so that the nearest interatomic distance is 20 Bohr. Ultrasoft pseudopotentials[49,50] with a plane-wave cutoff of 100 Ry and charge density cutoff of 800 Ry were adopted for some elements. The other elements were calculated using normconserving pseudopotentials[51] with a plane-wave cutoff of 200 Ry and charge density cutoff of 800 Ry. The (12 3 12 3 12) k-point grid was introduced based on the Monkhorst– Pack scheme. Before the LRD calculation of simple substances, we optimized three empirical parameters in eqs. (2) and (10). Similar to our previous work,[35] k was determined to minimize the mean absolute percent deviation (MAPD) of experimental C6 coefficients of rare-gas (Ne, Ar, Kr, and Xe) dimers. The optimized values were k 5 0.215 and 0.211 a.u. for the PBE and revPBE functionals, giving 2.92 and 3.11% of the MAPD, respectively. In principle, two empirical parameters j and R0 in eq. (10) should be fitted to reproduce the dispersion contribution in cohesive energies. However, highly accurate dispersion energy is difficult to evaluate for the periodic system. In this study, j and R0 were fitted using R values of rare-gas solids that reproduce experimental cohesive energies.[48] This treatment cor-
rects the error from the lack of dispersion as well as the inaccuracy of the exchange-repulsion interaction, which is caused by the approximated exchange functional.[15] Rare gases are appropriate as a training set for two reasons. First, the dominant attractive force of rare-gas elements is the dispersion force, and other components such as electrostatic and inductive interactions are negligible. Second, the effective radii given by eq. (10) are same in all atomic pairs in rare-gas solids, which enables the fitting of j and R0 for periodic systems. In our previous reports for isolated systems, the damping parameters for the LC functional were fit to reproduce the interaction of six rare-gas diatomics from He2 to Ar2. For pure functionals, it is difficult to reproduce weakly bound systems including He because of inaccuracy of exchange functional.[15] Therefore, we adjusted the damping parameters based on Ne (only for revPBE), Ar, Kr, and Xe. As for the PBE functional, the least square fitting determined j 5 0.4981 a.u. and R0 5 4.2060 bohr with R2 5 0.9968. The optimized damping parameters for the revPBE functional were j 5 0.4740 a.u. and R0 5 3.5348 bohr determined with R2 5 0.9976. These parameters were fixed throughout this study. The calculated atomization energies of simple substances were compared with the experimental values.[48] All numerical results are summarized in Supporting Information. The experimental atomization energies were corrected by zero-point energies based on the following formula: 9 EZPE 5 kB uD ; 8
(12)
where kB is the Boltzmann constant and hD is the Debye temperature of each element. The results are shown in Figure 1 and Table 1 as relative and mean deviations (MDs), respectively. Figure 1 shows that the dispersion-uncorrected PBE and revPBE functionals underestimated the atomization energies not only of rare-gas solids but also of other elements. Although the revPBE functional is known for reproducing the atomization energies of small molecules,[8] the lack of long-range correlation causes the significant underestimation in solids. According to Table 1, the LRD method decreased the deviations of atomization energies with respect to transition metals as well as rare-gas elements. In addition, smaller improvements were confirmed in the case of main group elements. The present results show the suitability of dispersion correction methods in metal- or covalent-bonded systems. We also investigated the effect of the dispersion correction in properties. In the present implementation, the LRD energy is computed in the post-SCF manner, meaning that the dispersion potential derived from the LRD method is not included in the Kohn–Sham equation. Therefore, the dispersion effect in electronic properties is only caused by geometric difference. We optimized lattice structures of several simple substances to investigate the dispersion effect in band gaps, spin magnetic moments, and bulk moduli. The bulk moduli were calculated by numerical fitting to the third-order Birch–Murnaghan isothermal equation of state.[52] Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Relative errors in cohesive energies of simple substances using (a) PBE and (b) revPBE functionals with and without the LRD method. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
8" 9 " 2 #2 " 2 #3 23 #= 9V0 B0 < V0 3 V0 3 V0 0 E ðV Þ5E0 1 21 B 0 1 21 624 ; 16 : V V V
(13) where E, V, V0, and B0 represent the total energy, lattice volume, equilibrium lattice volume, and bulk modulus, respectively. As shown in Table 2, the LRD method did not affect the values of properties largely, except for the bulk moduli of raregas elements. No obvious improvement or deterioration of the agreement with experimental values was found. In contrast, the bulk moduli of rare-gas elements were significantly
Table 1. Mean deviation (MD), mean absolute deviation (MAD), mean percentage deviation (MPD), and mean absolute percentage deviation (MAPD) of atomization energies of simple substances with experimental lattice constants. PBE Main group elements
Transition elements
Rare gases
Total
MD MAD MPD MAPD MD MAD MPD MAPD MD MAD MPD MAPD MD MAD MPD MAPD
20.074 0.189 24.7% 9.7% 20.343 0.456 27.7% 10.1% 20.103 0.103 286.0% 86.0% 20.219 0.324 213.2% 16.4%
PBE1LRD
revPBE
0.122 20.287 0.177 0.334 5.6% 214.4% 8.5% 16.5% 20.147 20.800 0.405 0.808 23.1% 217.9% 8.1% 18.1% 20.003 20.198 0.008 0.198 6.2% 2185.9% 14.5% 185.9% 20.032 20.552 0.284 0.575 1.0% 230.8% 8.8% 31.7%
Molecular crystals We evaluated the cohesive energies of the X23 benchmark set,[60] which contains reference electronic cohesive energies of 23 molecular crystals. Lattice parameters and atomic Table 2. Properties of simple substances calculated using optimized lattice structures.
revPBE1LRD 0.046 0.202 3.4% 9.1% 20.492 0.560 210.6% 12.2% 0.001 0.005 0.7% 4.6% 20.244 0.375 24.3% 10.4%
MD and MAD are given in eV.
4
changed. Both the PBE and revPBE functionals cannot reproduce the bulk moduli even qualitatively. By applying the LRD method, the bulk moduli increased in the order of Ne, Ar, and Kr, which is in accordance with experimental results. The result indicates that the properties related to the shape of the potential energy surface are expected to be described reasonably if the dispersion correction is applied to the DFT calculation of dispersion-dominated periodic systems.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
Band gap C Si Ge Magnetic moment Fe Ni Co Bulk modulus C Si Ge Ne Ar Kr
PBE
PBE1LRD
revPBE
revPBE1LRD
Expt.
4.04 1.21 1.09
4.05 1.20 1.10
4.05 1.23 1.08
4.05 1.19 1.09
5.5[48] 1.11[53] 0.67[53]
2.37 0.65 1.80
2.32 0.64 1.79
2.46 0.66 1.82
2.40 0.65 1.81
2.2[54] 0.6[54] 1.7[54]
430 86.2 54.5 1.28 0.75 0.79
423 83.2 58.5 2.01 2.14 2.41
420 87.5 61.1 2.02 2.91 3.39
443[55] 99.2[56] 75.8[56] 1.17[57] 2.7[58] 3.6[59]
422 84.3 56.1 n/a n/a n/a
Band gap, magnetic moment, bulk modulus are given in eV, bohr magneton, and GPa, respectively.
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 3. Cohesive energies of the X23 benchmark set (kcal/mol). PBE Hydrogen bonded Acetic acid Ammonia Cyanamide Cytosine Ethylcarbamate Formamide Imidazole Oxalic acid (a) Oxalic acid (b) Succinic acid Uracil Urea VdW bonding 1,4-cyclohexanedione Adamantane Anthracene Benzene CO2 Hexamine Naphthalene Pyrazine Pyrazole Triazine Trioxane Overall MD MAD MPD MAPD
PBE1LRD
revPBE
revPBE1LRD
Ref.
48.1 28.2 70.2 99.8 47.6 58.9 53.4 83.6 86.8 82.5 84.0 78.9
(224.7) (29.0) (29.5) (270.0) (238.7) (220.3) (233.4) (212.7) (29.3) (247.8) (251.7) (223.6)
83.2 45.3 97.9 166.3 95.7 86.3 98.0 126.4 127.6 144.2 144.4 112.8
(10.4) (8.1) (18.2) (23.5) (9.4) (7.1) (11.2) (30.1) (31.5) (13.9) (8.7) (10.3)
22.1 11.4 43.7 49.3 18.1 33.8 24.4 44.1 48.3 37.9 37.7 48.5
(250.7) (225.8) (236.0) (2120.5) (268.2) (245.4) (262.4) (252.2) (247.8) (292.4) (298.0) (254.0)
70.2 37.0 88.4 148.5 82.0 75.2 90.5 107.1 108.8 123.4 125.7 96.4
(22.6) (20.2) (8.7) (221.3) (24.3) (24.0) (3.7) (10.8) (12.7) (26.9) (210.0) (26.1)
72.8 37.2 79.7 169.8 86.3 79.2 86.8 96.3 96.1 130.3 135.7 102.5
35.1 24.3 3.2 5.1 9.9 15.9 3.1 17.1 43.0 14.8 19.4
(253.5) (273.7) (2109.5) (246.6) (218.5) (270.3) (278.6) (244.2) (234.7) (246.9) (247.0)
102.2 90.3 121.7 61.8 26.0 98.9 91.4 70.9 84.3 66.2 68.5
(13.6) (20.9) (9.0) (10.1) (22.4) (12.7) (9.7) (9.6) (6.6) (4.5) (2.1)
1.8 236.6 247.6 219.8 20.5 215.2 235.1 214.4 18.4 213.1 24.0
(286.8) (2106.0) (2160.3) (271.5) (228.9) (2101.4) (2116.8) (275.7) (259.3) (274.8) (270.4)
89.7 76.4 117.2 59.2 23.9 90.2 88.1 64.2 77.3 60.1 60.1
(1.1) (7.0) (4.5) (7.5) (24.5) (4.0) (6.4) (2.9) (20.4) (21.6) (26.3)
88.6 69.4 112.7 51.7 28.4 86.2 81.7 61.3 77.7 61.7 66.4
242.4 42.4 252.2% 52.2%
10.9 11.5 13.3% 14.2%
0.1 6.0 0.6% 7.0%
274.1 74.1 290.6% 90.6%
Differences from reference values are shown in parentheses.
locations in the solid state were optimized using the LRD method combined with the PBE or revPBE functional. The dissociated state was modeled by one molecule in the center of a cubic unit cell with a cell length of 30 Bohr. The ultrasoft pseudopotentials determined by the Rabe–Rappe–Kaxiras– Joannopoulos (RRKJ) method[50] were adopted for all elements. The plane-wave cutoff and charge density cutoff were 49 and 400 Ry, respectively. The k-point sampling was performed in the same manner as the benchmark study of 21 molecular crystals (C21 set) by Otero-de-la-Roza and Johnson.[61] The calculated and reference cohesive energies are summarized in Table 3. The cohesive energies without a dispersion correction were evaluated using the lattice structure optimized using the LRD method. Although the sign of the cohesive energies of hydrogen-bonded systems is reproduced without any dispersion correction, the application of the LRD method significantly reduced deviations. Considering that the dispersion-uncorrected GGA functional can evaluate a large portion of the interaction energy of hydrogen-bonded molecular dimers,[17] the result indicates the importance of dispersion correction for molecular crystals. As expected, the LRD method reduced errors in the cohesive energies of vdW bonding systems. The PBE1LRD approach somewhat overestimates the cohesive energies on average. The overbinding may be related to that in the cohesive energy of neon shown in Figure 1, where PBE1LRD overestimates the cohesive energy due to the
attractive nature without any dispersion correction. The ME of revPBE1LRD is nearly zero, indicating that three empirical parameters adjusted for the cohesive energies of rare gases are transferrable to molecular crystals. The overall mean absolute deviation (MAD) of the revPBE functional is reduced from 74.14 to 5.97 kJ/mol, and the overall MAPD is reduced from 90.60 to 6.97%. This accuracy is comparable to the MBD method combined with the PBE functional: 5.92 kJ/mol (8.05%).[60] We compare the accuracies of PBE1LRD and revPBE1LRD with those of other dispersion-corrected functionals in terms of cohesive energies, optimized lattice parameters, and lattice volumes. Here, we consider the C21 benchmark set,[61] which
Table 4. Mean and mean absolute deviations from experimental lattice parameters of cohesive energies (kJ/mol), optimized cell lengths (A˚), angles ( ), and volumes (A˚3) of the C21 benchmark set. Cohesive energy
Length
Angles
Volume
MD
MAD
MD
MAD
MD
MAD
MD
MAD
PBE1LRD 10.7 PBE-D[61] 8.2 PBE-TS[61] 16.8 PBE-XDM[61] 0.2 revPBE1LRD 0.2 vdW-DF[61] 8.8
11.3 9.0 17.0 5.4 6.0 10.2
0.00 20.07 0.04 0.16 0.09 0.29
0.11 0.11 0.10 0.20 0.14 0.31
20.22 0.05 20.06 20.15 20.05 20.06
0.24 0.19 0.14 0.29 0.12 0.11
4.2 211.7 7.3 31.2 11.7 48.0
8.6 14.4 10.1 31.2 14.0 48.0
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 2. Adsorption structures and unit cells of n-butane on (a) M (111) (M 5 Ag, Au, Pd, and Pt) and (b) Cu (111) surfaces. (c) Side view of the adsorbed structure with the Hmetal interaction shown by the dashed line.
consist 21 molecular crystals also included in the X23 set. All results of our calculations are summarized in Supporting Information. The deviations of the calculated results are summarized in Table 4 as well as previously reported values. In addition to PBE-based methods, vdW-DF is listed because the exchange part of the revPBE functional is included. As for cohesive energies, PBE1LRD and revPBE1LRD showed the same tendency as the X23 set. PBE-XDM and revPBE1LRD gave nearly zero MDs and are in a good agreement with the reference. The other methods resulted in overestimation and larger MADs. In terms of structural parameters, PBE1LRD reproduced them well while revPBE1LRD slightly overestimated the lattice lengths and volumes. This originates in the strong repulsive character of the revPBE exchange at short internuclear distances and remarkably rapid decay at long distances.[15] The overestimation is not so large, smaller than that of PBE-XDM and vdW-DF. We speculate that the GGA correlation, which is not included in vdW-DF, results in decreased deviations to some extent. The lattice angles were well reproduced by revPBE1LRD and vdW-DF. The present numerical assessment of cohesive energies and structural parameters support the usefulness of the LRD method combined with the revPBE functional for molecular crystals. Physical adsorption In the present subsection, n-butane on a transition metal surface is chosen as an example of physical adsorption. This system was previously studied by Morikawa and coworkers,[62,63] in which adsorptions on Cu(100), Cu(111), Au(111), and Pt(111) 6
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
surfaces were calculated using the PBE functional[7] and the vdW-DF series.[21,22] Here, we apply the revPBE1LRD method to adsorptions on the (111) surfaces of Pd and Ag in addition to those of Cu, Au, and Pt. The unit cell of the adsorption structure is shown in Figure 2a. We introduced a slab model, consisting of four layers of metal atoms, to represent the (111) surface. Interatomic distances of metal atoms werepset ffiffiffi to experimentally determined pffiffiffi bulk structures. The 73 3 surface unit cells were adopted to represent the fully adsorbed state. The intersurface distance ˚ . The geometry of n-butane was optimized was set to 20 A with fixed layers and lattice constants. However, an appropriate adsorption structure was not found in the case of Cu. We speculate that this is a result of intermolecular repulsion due pffiffiffi to the small surface unit cell. Instead, we used the 732 surface unit cell shown in Figure 2b for Cu. The ultrasoft
Table 5. Adsorption energies (eV) and distances (A˚) of n-butane on the (111) surface of transition metals. Adsorption energy [63]
revPBE1LRD vdW-DF2 Pd Pt Cu Ag Au
0.60 0.58 0.55 0.49 0.48
– 0.61 0.37 – 0.52
Adsorption distance Expt. – 0.53[64] 0.51[65] – 0.42[66]
revPBE1LRD vdW-DF2[63] Est.[a] 3.41 3.59 3.53 3.69 3.79
– 3.88 3.95 – 3.97
– 3.56 – – 3.52
[a] Reported in Ref. [62], where the distances are fit to reproduce the red-shifts of the CAH stretching mode and the work function changes due to the adsorption.
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 3. Difference in electron densities along the CAHmetal interactions given by the revPBE1LRD method.
pseudopotentials determined by the RRKJ method[50] were adopted for all elements. The (4 3 6 3 1) k-point grid was introduced based on the Monkhorst–Pack scheme. In Table 5, the calculated adsorption energies and distances are summarized together with the results of the previous studies.[62,63] Here, the adsorption distance is defined as the average of distances between each carbon atom and the plane including the nuclei of surface atoms. The adsorption energies determined using the revPBE1LRD method are in good agreement with the experimental values. Compared to vdW-DF2,[22] which consists of the refitted Perdew-Wang 86 exchange,[67] LDA correlation,[4] and nonlocal correlation functionals,[22] the present scheme provided shorter adsorption distances and smaller differences from the distances fitted to experimental
results. The adsorption energies in the case of Ag and Au are slightly smaller than those of the other metals. This is consistent with the longer adsorption distances of Ag and Au than those of the other metals. In addition to adsorption energies and structures, electron transfer between the CAH bonds and the metal layers is worth investigating. According to the previous report,[68] this is interpreted as a result of mixing the CAH bonding and antibonding orbitals and the d band of transition metals. In principle, the LRD method is suitable for charge-transferred systems because this method is dependent on the electron distribution obtained as a solution of the DFT calculation. Previously, we confirmed that the electronic excitation changes atomic polarizabilities given by the LRD method reasonably.[42] Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
7
FULL PAPER
WWW.C-CHEM.ORG
Figure 3 visualizes the electron transfer obtained by the revPBE1LRD method. All systems clearly show the electron transfer from the atomic region of hydrogen and metal to the intermediate region. A larger transfer is confirmed in the case of Pd and Pt. This is because Pd and Pt have larger vacancies in the d band than Cu, Ag, and Au. As for the noble metals, the charge transfer is larger in the order of Cu > Ag > Au. This is consistent with the reactivity of the noble metals. These tendencies for difference densities may cause the difference in the adsorption energies. The results indicate that the revPBE1LRD method can provide not only accurate adsorption energies but also reasonable adsorption structures and electron transfers.
Conclusion We implemented the LRD method in the Quantum ESPRESSO and assessed the combination with the PBE and revPBE functionals. This is based on the feasibility of the LRD method, which does not depend on the electronic basis functions such as Gaussian and plane-wave basis functions but the electron density. Numerical applications confirmed the accuracy of the LRD method even in the case of periodic systems. The LRD method improved the underestimation of cohesive energies not only in rare gases but also in strongly bonded elements. The LRD method also improved the accuracy of the cohesive energies and lattice parameters of molecular crystals. Furthermore, we confirmed that not only the physical adsorption energy but also the electron transfer in n-butane on metal surfaces were calculated reasonably.
Acknowledgment Some of the present calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). Keywords: density functional theory dispersion force periodic boundary condition plane-wave basis set physical adsorption
How to cite this article: Y. Ikabata Y. Tsukamoto Y. Imamura, H. Nakai. J. Comput. Chem. 2014, DOI: 10.1002/jcc.23807
]
Additional Supporting Information may be found in the online version of this article.
[1] [2] [3] [4] [5] [6] [7] [8] [9]
P. Hohenberg, W. Kohn, Phys. Rev. 1964, 136, B864. W. Kohn, L. J. Sham, Phys. Rev. 1965, 140, A1133. J. C. Slater, Phys. Rev. 1951, 81, 385. S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. A. D. Becke, Phys. Rev. A 1988, 38, 3098. C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. Y. Zhang, W. Yang, Phys. Rev. Lett. 1998, 80, 890. J. P. Perdew, S. Kurth, A. Zupan, P. Blaha, Phys. Rev. Lett. 1999, 82, 2544. [10] J. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Phys. Rev. Lett. 2003, 91, 146401.
8
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
[11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
[44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55]
[56] [57] [58]
A. D. Becke, J. Chem. Phys. 1993, 98, 5648. C. Adamo, V. Barone, J. Chem. Phys. 1999, 110, 6158. J. Heyd, G. E. Scuseria, M. Ernzerhof, J. Chem. Phys. 2003, 118, 8207. A. Seidl, A. G€ orling, P. Vogl, J. A. Majewski, M. Levy, Phys. Rev. B 1986, 53, 3764. F. O. Kannemann, A. D. Becke, J. Chem. Theory Comput. 2009, 5, 719. N. Marom, A. Tkatchenko, M. Rossi, V. V. Gobre, O. Hod, M. Scheffler, L. Kronik, J. Chem. Theory Comput. 2011, 7, 3944. J. Klimes, A. Michaelides, J. Chem. Phys. 2012, 137, 120901. A. J. Cohen, P. Mori-Sanchez, W. Yang, Chem. Rev. 2012, 112, 289. J. F. Dobson, B. P. Dinte, Phys. Rev. Lett. 1996, 76, 1780. Y. Andersson, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Lett. 1996, 76, 102. M. Dion, H. Rydberg, E. Schr€ oder, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Lett. 2004, 92, 246401. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, D. C. Langreth, Phys. Rev. B 2010, 82, 081101. I. Hamada, Phys. Rev. B 2014, 89, 121103. O. A. Vydrov, T. Van Voorhis, Phys. Rev. Lett. 2009, 103, 063004. O. A. Vydrov, T. Van Voorhis, J. Chem. Phys. 2010, 133, 244103. M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 2001, 114, 5149. Q. Wu, W. Yang, J. Chem. Phys. 2002, 116, 515. S. Grimme, J. Comput. Chem. 2004, 25, 1463. S. Grimme, J. Comput. Chem. 2006, 27, 1787. S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104. L. Goerigk, S. Grimme, J. Chem. Theory Comput. 2011, 7, 291. A. D. Becke, E. R. Johnson, J. Chem. Phys. 2005, 122, 154104. A. Tkatchenko, M. Scheffler, Phys. Rev. Lett. 2009, 102, 073005. A. Tkatchenko, R. A. DiStasio, Jr., R. Car, M. Scheffler, Phys. Rev. Lett. 2012, 108, 236402. A. Tkatchenko, A. Ambrosetti, R. A. DiStasio, Jr., J. Chem. Phys. 2013, 138, 074106. T. Sato, H. Nakai, J. Chem. Phys. 2009, 131, 224104. T. Sato, H. Nakai, J. Chem. Phys. 2010, 133, 194101. Y. Ikabata, T. Sato, H. Nakai, Int. J. Quantum Chem. 2013, 113, 257. H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540. R. Kar, J.-W. Song, T. Sato, K. Hirao, J. Comput. Chem. 2013, 34, 2353. Y. Ikabata, H. Nakai, Chem. Phys. Lett. 2013, 556, 386. Y. Ikabata, H. Nakai, J. Chem. Phys. 2012, 137, 124106. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, R. M. Wentzcovitch, J. Phys.: Condens. Matter 2009, 21, 395502. A. J. Stone, In Theory of Intermolecular Forces; Clarendon: Oxford, 1996. O. A. Vydrov, T. van Voorhis, J. Chem. Phys. 2009, 130, 104105. K. Rapcewicz, N. W. Ashcroft, Phys. Rev. B 1991, 44, 4032. A. D. Becke, J. Chem. Phys. 1988, 88, 2547. C. Kittel, In Introduction to Solid State Physics, 8th ed.; Wiley: New York, 2005. D. Vanderbilt, Phys. Rev. B 1990, 41, 7892. A. M. Rappe, K. M. Rabe, E. Kaxiras, J. D. Joannopoulos, Phys. Rev. B 1990, 41, 1227. N. Troullier, J. L. Martins, Phys. Rev. B 1991, 43, 1993. F. Birch, Phys. Rev. 1947, 71, 809. B. G. Streetman, S. Banerjee, In Solid State Electronic Devices, 5th ed.; Prentice Hall: New Jersey, 2000. R. M. Bozorth, In Ferromagnetism; D. Van Nostrand Company, Inc.: Princeton, NJ, 1956. M. E. Levinshtein, S. L. Rumyantsev, In Handbook Series on Semiconductor Parameters; M. Levinshtein, S. Rumyantsev, M. Shur, Eds.; Vol. 1; World Scientific: Singapore, 1996. K. H. Hellwege, In Landolt-Bornstein New Series, Group III: Condensed Matter, Vol. 17a; Springer: Berlin, 1966. D. N. Batchelder, D. L. Losee, R. O. Simmons, Phys. Rev. 1967, 162, 767. O. G. Peterson, D. N. Batchelder, R. O. Simmons, Phys. Rev. 1966, 150, 703.
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
[59] [60] [61] [62] [63] [64]
D. L. Losee, R. O. Simmons, Phys. Rev. 1968, 172, 944. A. M. Reilly, A. Tkatchenko, J. Chem. Phys. 2013, 139, 024705. A. Otero-de-la-Roza, E. R. Johnson, J. Chem. Phys. 2012, 137, 054103. Y. Morikawa, H. Ishii, K. Seki, Phys. Rev. B 2004, 69, 041403. K. Lee, Y. Morikawa, D. C. Langreth, Phys. Rev. B 2010, 82, 155461. S. L. Tait, Z. Dohnalek, C. T. Campbell, B. D. Kay, J. Chem. Phys. 2006, 125, 234308. [65] R. Z. Lei, A. J. Gellman, B. E. Koel, Surf. Sci. 2004, 554, 125. [66] S. M. Wetterer, D. J. Lavrich, T. Cummings, S. L. Bernasek, G. Scoles, J. Phys. Chem. B 1998, 102, 9266.
FULL PAPER
[67] E. D. Murray, K. Lee, D. C. Langreth, J. Chem. Theory Comput. 2009, 5, 2754. € om, L. Triguero, M. Nyberg, H. Ogasawara, L. G. M. Pettersson, [68] H. Ostr€ A. Nilsson, Phys. Rev. Lett. 2004, 91, 046102.
Received: 26 July 2014 Revised: 6 November 2014 Accepted: 8 November 2014 Published online on 00 Month 2014
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
9
FULL PAPER
Local Response Dispersion Method in Periodic Systems: Implementation and Assessment Yasuhiro Ikabata,[a] Yusuke Tsukamoto,[a] Yutaka Imamura,[a]† and Hiromi Nakai*[a,b,c,d] We report the implementation of the local response dispersion (LRD) method in an electronic structure program package aimed at periodic systems and an assessment combined with the Perdew–Burke–Ernzerhof (PBE) functional and its revised version (revPBE). The real-space numerical integration was implemented and performed exploiting the electron distribution given by the plane-wave basis set. The dispersion-
corrected density functionals revPBE1LRD was found to be suitable for reproducing energetics, structures, and electron distributions in simple substances, molecular crystals, and C 2014 Wiley Periodicals, Inc. physical adsorptions. V
Introduction
sion (MBD) method[34] and investigated the relation between pairwise methods and the MBD method.[35] Sato and Nakai proposed the local response dispersion (LRD) method in 2009 as a density-dependent dispersion correction scheme.[36] This method calculates dispersion coefficients and energies on the basis of effective atomic polarizabilities computed using the electron density and its gradient. Although the LRD method is one of the pairwise correction methods, the dispersion energy is given as a density functional. Similar to conventional XC functionals and nonlocal correlation functionals, the LRD method does not depend on any experimental values explicitly. Furthermore, the extension of the LRD method to multicenter interactions was achieved,[37] followed by the formulation of the self-consistent field (SCF) treatment and the analytic energy gradient.[38]
Density functional theory (DFT) is a standard electronic structure methodology used in periodic systems. From the establishment of Kohn–Sham DFT[1,2] to the present, the local density approximation (LDA)[3,4] to the exchange-correlation (XC) functional has been employed. The generalized gradient approximation (GGA)[5–8] and meta-GGA[9,10] were proposed partly to include the nonlocal XC effect. Hybrid functionals,[11–13] which are theoretically justified by the generalized Kohn–Sham scheme,[14] handle further nonlocal exchange by mixing the Hartree–Fock exchange. Although these improvements have increased the accuracy and availability of DFT, it is well known that this is not the case with noncovalent interactions. LDA evaluates the van der Waals force qualitatively, although the attractive force originates in the inaccuracy of the exchange functional.[15] To make the matter worse, GGA and most later XC functionals cannot describe van der Waals interaction even qualitatively.[16] This failure originates in the inaccurate description of dispersion interactions. The development of a dispersion correction technique is an increasingly important topic in DFT.[17,18] Because the dispersion interaction is derived from the nonlocal correlation effect, several nonlocal correlation functionals have been proposed. One example, proposed by Dobson and Dinte,[19] is equivalent to the Andersson–Langreth–Lundqvist (ALL) functional.[20] More sophisticated nonlocal correlation functionals developed later were included in vdW-DF series[21–23] and XC functionals by Vydrov and Van Voorhis.[24,25] Conversely, a number of pairwise correction methods that compute dispersion energy as a sum of atom–atom contributions, have been proposed. The DFT-D method[26–30] has a conceptual simplicity, and various combinations with XC functionals have been examined.[31] Instead of predetermined dispersion coefficients, some methods introduced a dependency on the electronic state, for example, the exchange-dipole moment (XDM) scheme[32] and the vdW-TS method.[33] Recently, Tkatchenko et al. proposed the many-body disper-
DOI: 10.1002/jcc.23807
[a] Y. Ikabata, Y. Tsukamoto, Y. Imamura, H. Nakai Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo, 1698555, Japan E-mail: [email protected] [b] H. Nakai Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo, 169-8555, Japan [c] H. Nakai CREST, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan [d] H. Nakai Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Katsura, Kyoto, 615–8520, Japan † Present address: Advanced Institute for Computational Science, RIKEN, 7-1-26, Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan. Contract grant sponsor: Ministry of Education, Culture, Sports, Science and Technology [MEXT; Strategic Programs for Innovative Research (SPIRE) and the Computational Materials Science Initiative (CMSI)]; Contract grant sponsor: Elements Strategy Initiative for Catalyst & Batteries (ESICB) [MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012)]; Contract grant sponsor: Japan Science and Technology Agency [JST; Core Research for Evolutional Science and Technology (CREST) program “Theoretical Design of Materials with Innovative Functions Based on Relativistic Electronic Theory”] C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
1
FULL PAPER
WWW.C-CHEM.ORG
Furthermore, the combination of the LRD method with the long-range corrected (LC) functional[39] was found to be highly accurate not only for conventional molecular dimers[40] but also for open-shell[41] and excited-state molecular complexes.[42] These developments and assessments treated isolated molecular complexes with Gaussian basis functions. To treat molecular solids and physical adsorption, the extension of the LRD method to periodic systems is of vital importance. In this study, the LRD method is implemented in a package based on the periodic boundary condition and the plane-wave basis functions. We here combine the LRD method with the Perdew–Burke–Ernzerhof functional (PBE)[7] and its revised version (revPBE),[8] which are representatives of pure GGA functionals. Simple substances, molecular crystals, and physical adsorptions of n-butane on transition metals are considered in the assessment of these dispersion-corrected functionals.
jrqj 6q
pffiffiffiffiffiffiffiffi 4pq : kF
(3)
Equation (3) is not implemented in the LRD method with Gaussian basis functions while it is adopted in practical calculations of the ALL functional.[20] The present study requires a cutoff because the plane-wave basis generates small electron densities in the region far from the nucleus, which sometimes causes serious errors in the polarizability. Becke’s partition function[47] was implemented as wa and was computed considering all atoms in 3 3 3 3 3 unit cells. On the basis of the polarizability and geometry, the nthorder (n 5 6, 8, and 10) interatomic dispersion coefficient is calculated as Cnab 5
ð1 1 X X X ab duaal m ;l m ’ ðiuÞabl m ;l m ’ ðiuÞ; Sl1 m1 ;l2 m2 Sab l1 m1 ’;l2 m2 ’ 1 1 1 1 2 2 2 2 2p l l m1 m2 0 1 2 m1 ’m2 ’
(4)
Theory and Implementation We implemented the LRD method in the Quantum ESPRESSO package,[43] following the first paper[36] in which pairwise dispersion energy was computed in the post-SCF manner. Our previous study indicates that the SCF treatment brings minor effect.[38] The procedure to calculate dispersion correction energy is explained in this section. The energy gradient required for the lattice structure optimization is calculated by simply differentiating the LRD energy with respect to nuclear coordinates. The differentiation of total energy with respect to density is neglected. First, atomic polarizabilities are computed using the electron density and its gradient, obtained by the conventional DFT calculation. The polarizability is obtained from the intermolecular dispersion energy by the second-order perturbation theory, the multicenter-multipole expansion of the Coulomb operator,[44] and the local response approximation to the density response function.[19] As a result, the polarizability of atom a at an imaginary frequency iu is given by ð qðrÞ aalm;l’m’ ðiuÞ5 drwa2 ðrÞ 2 rRlm ðr2Ra Þ rRl’m’ ðr2Ra Þ; (1) x0 ðrÞ1u2 where wa is the partition function of atom a, q is the electron density, and R is the solid harmonics. l and m are integers with the restrictions l 1 and – l m l. l 6¼ l0 terms are neglected in the following because of their minute numerical impact. x0 in eq. (1) is a dispersion relation of a dielectric medium in the long-wavelength limit. The LRD method adopts following gradient-corrected form[45]: k2 ð11ks2 Þ2 x0 ðrÞ5 F ; (2) 3 2 1=3 with the local Fermi wave vector kF 5ð3p qÞ , an empirical parameter k, and the reduced density gradient s5jrqj=ð2kF qÞ. The polarizabilities of atoms in one unit cell are numerically integrated using the Cartesian grid. We consider 3 3 3 3 3 unit cells for real-space integration with a cutoff criteria proposed by Rapcewicz and Ashcroft,[46] 2
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
where Sab l1 m1 ;l2 m2 is the angular factor of the multipole expansion defined by l1 Sab l1 m1 ;l2 m2 5ð21Þ
ð2l1 12l2 11Þ! 1=2 dL;l1 1l2 3 ð2l1 Þ!ð2l2 Þ!
l1 m1
l2
L
m2
M
! ^ ab Þ: CL;M ðR
(5) The 2 3 3 matrix is Wigner’s 3-j symbol, and CL;M is the spherical harmonics with Racah’s normalization factor evaluated at ^ ab 5Rab =Rab. Note that the the polar angles of the unit vector R order is related to l1 and l2 by n 5 2(l1 1 l2 1 1). The integration with respect to the imaginary frequency in eq. (4) is possible pffiffiffiffiffiffiffiffiffiffiffi by transforming the variable u to v5u= 11u2 and applying Gauss–Chevyshev’s rule. The modified atomic polarizability is given by a alm;lm’ ðvk Þ ð12vk 2 Þ21=2 aalm;lm’ ðvk Þ ð ð12v2 Þ1=2 qðrÞ 5 drwa2 ðrÞ 3rRlm ðr2Ra ÞrRlm’ ðr2Ra Þ: ð12v 2 Þx20 ðrÞ1v 2
(6)
The dispersion coefficient is computed using the modified atomic polarizabilities: Cnab 5
N X 1 X X X ab a al m ;l m ’ ðvk Þ a bl m ;l m ’ ðvk Þ; Sl1 m1 ;l2 m2 Sab l1 m1 ’;l2 m2 ’ 1 1 1 1 2 2 2 2 4N l l m1 m2 k51 1 2 m1 ’m2 ’
(7) with N quadrature points at vk 5cos ½ð2k21Þp=4N. Twelve points were confirmed to be sufficient to maintain numerical accuracy. Finally, the dispersion correction energy for the unit cell is calculated. ELRD 52
1 X X ab 2n ðnÞ C R f ðRab Þ 2 n56;8;10 a;b n ab damp
(8)
We consider all atomic pairs in one unit cell and atoms in other unit cells containing at least one atom that is within 20 WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Bohr from the atom in the focused unit cell. The damping function is " ðnÞ fdamp 5exp
26 # Rab 2m ; n52m14; R
(9)
where R is interpreted as the sum of van der Waals radii with respect to atoms a and b. Assuming a linear relationship between atomic radius and polarizability, R is written by R5j
n 1=3 b 1=3 o 1 aeff ðidÞ 1R0 ; aaeff ðidÞ
(10)
with two empirical parameters j and R0. a a is the effective static polarizability,[36] defined by ð aaeff ðidÞ5 drwa ðrÞ
qðrÞ : x20 ðrÞ1d2
(11)
Note that d is set to 0.01 a.u. to maintain numerical stability.
Illustrative Applications Simple substances First, we calculated the atomization energies of simple substances. The number of elements treated in this subsection is 47. We excluded elements that form molecules in the most stable solid phase. Our aim in this subsection is to investigate the effect of the dispersion correction in metal- or covalentbonded periodic systems after fitting empirical parameters to rare-gas elements. The lattice parameters were set to experimental values.[48] The atomization energy is given as the difference in atomic energies between the solid and dissociated states. As for the dissociated state, we set one atom in a large unit cell so that the nearest interatomic distance is 20 Bohr. Ultrasoft pseudopotentials[49,50] with a plane-wave cutoff of 100 Ry and charge density cutoff of 800 Ry were adopted for some elements. The other elements were calculated using normconserving pseudopotentials[51] with a plane-wave cutoff of 200 Ry and charge density cutoff of 800 Ry. The (12 3 12 3 12) k-point grid was introduced based on the Monkhorst– Pack scheme. Before the LRD calculation of simple substances, we optimized three empirical parameters in eqs. (2) and (10). Similar to our previous work,[35] k was determined to minimize the mean absolute percent deviation (MAPD) of experimental C6 coefficients of rare-gas (Ne, Ar, Kr, and Xe) dimers. The optimized values were k 5 0.215 and 0.211 a.u. for the PBE and revPBE functionals, giving 2.92 and 3.11% of the MAPD, respectively. In principle, two empirical parameters j and R0 in eq. (10) should be fitted to reproduce the dispersion contribution in cohesive energies. However, highly accurate dispersion energy is difficult to evaluate for the periodic system. In this study, j and R0 were fitted using R values of rare-gas solids that reproduce experimental cohesive energies.[48] This treatment cor-
rects the error from the lack of dispersion as well as the inaccuracy of the exchange-repulsion interaction, which is caused by the approximated exchange functional.[15] Rare gases are appropriate as a training set for two reasons. First, the dominant attractive force of rare-gas elements is the dispersion force, and other components such as electrostatic and inductive interactions are negligible. Second, the effective radii given by eq. (10) are same in all atomic pairs in rare-gas solids, which enables the fitting of j and R0 for periodic systems. In our previous reports for isolated systems, the damping parameters for the LC functional were fit to reproduce the interaction of six rare-gas diatomics from He2 to Ar2. For pure functionals, it is difficult to reproduce weakly bound systems including He because of inaccuracy of exchange functional.[15] Therefore, we adjusted the damping parameters based on Ne (only for revPBE), Ar, Kr, and Xe. As for the PBE functional, the least square fitting determined j 5 0.4981 a.u. and R0 5 4.2060 bohr with R2 5 0.9968. The optimized damping parameters for the revPBE functional were j 5 0.4740 a.u. and R0 5 3.5348 bohr determined with R2 5 0.9976. These parameters were fixed throughout this study. The calculated atomization energies of simple substances were compared with the experimental values.[48] All numerical results are summarized in Supporting Information. The experimental atomization energies were corrected by zero-point energies based on the following formula: 9 EZPE 5 kB uD ; 8
(12)
where kB is the Boltzmann constant and hD is the Debye temperature of each element. The results are shown in Figure 1 and Table 1 as relative and mean deviations (MDs), respectively. Figure 1 shows that the dispersion-uncorrected PBE and revPBE functionals underestimated the atomization energies not only of rare-gas solids but also of other elements. Although the revPBE functional is known for reproducing the atomization energies of small molecules,[8] the lack of long-range correlation causes the significant underestimation in solids. According to Table 1, the LRD method decreased the deviations of atomization energies with respect to transition metals as well as rare-gas elements. In addition, smaller improvements were confirmed in the case of main group elements. The present results show the suitability of dispersion correction methods in metal- or covalent-bonded systems. We also investigated the effect of the dispersion correction in properties. In the present implementation, the LRD energy is computed in the post-SCF manner, meaning that the dispersion potential derived from the LRD method is not included in the Kohn–Sham equation. Therefore, the dispersion effect in electronic properties is only caused by geometric difference. We optimized lattice structures of several simple substances to investigate the dispersion effect in band gaps, spin magnetic moments, and bulk moduli. The bulk moduli were calculated by numerical fitting to the third-order Birch–Murnaghan isothermal equation of state.[52] Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Relative errors in cohesive energies of simple substances using (a) PBE and (b) revPBE functionals with and without the LRD method. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
8" 9 " 2 #2 " 2 #3 23 #= 9V0 B0 < V0 3 V0 3 V0 0 E ðV Þ5E0 1 21 B 0 1 21 624 ; 16 : V V V
(13) where E, V, V0, and B0 represent the total energy, lattice volume, equilibrium lattice volume, and bulk modulus, respectively. As shown in Table 2, the LRD method did not affect the values of properties largely, except for the bulk moduli of raregas elements. No obvious improvement or deterioration of the agreement with experimental values was found. In contrast, the bulk moduli of rare-gas elements were significantly
Table 1. Mean deviation (MD), mean absolute deviation (MAD), mean percentage deviation (MPD), and mean absolute percentage deviation (MAPD) of atomization energies of simple substances with experimental lattice constants. PBE Main group elements
Transition elements
Rare gases
Total
MD MAD MPD MAPD MD MAD MPD MAPD MD MAD MPD MAPD MD MAD MPD MAPD
20.074 0.189 24.7% 9.7% 20.343 0.456 27.7% 10.1% 20.103 0.103 286.0% 86.0% 20.219 0.324 213.2% 16.4%
PBE1LRD
revPBE
0.122 20.287 0.177 0.334 5.6% 214.4% 8.5% 16.5% 20.147 20.800 0.405 0.808 23.1% 217.9% 8.1% 18.1% 20.003 20.198 0.008 0.198 6.2% 2185.9% 14.5% 185.9% 20.032 20.552 0.284 0.575 1.0% 230.8% 8.8% 31.7%
Molecular crystals We evaluated the cohesive energies of the X23 benchmark set,[60] which contains reference electronic cohesive energies of 23 molecular crystals. Lattice parameters and atomic Table 2. Properties of simple substances calculated using optimized lattice structures.
revPBE1LRD 0.046 0.202 3.4% 9.1% 20.492 0.560 210.6% 12.2% 0.001 0.005 0.7% 4.6% 20.244 0.375 24.3% 10.4%
MD and MAD are given in eV.
4
changed. Both the PBE and revPBE functionals cannot reproduce the bulk moduli even qualitatively. By applying the LRD method, the bulk moduli increased in the order of Ne, Ar, and Kr, which is in accordance with experimental results. The result indicates that the properties related to the shape of the potential energy surface are expected to be described reasonably if the dispersion correction is applied to the DFT calculation of dispersion-dominated periodic systems.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
Band gap C Si Ge Magnetic moment Fe Ni Co Bulk modulus C Si Ge Ne Ar Kr
PBE
PBE1LRD
revPBE
revPBE1LRD
Expt.
4.04 1.21 1.09
4.05 1.20 1.10
4.05 1.23 1.08
4.05 1.19 1.09
5.5[48] 1.11[53] 0.67[53]
2.37 0.65 1.80
2.32 0.64 1.79
2.46 0.66 1.82
2.40 0.65 1.81
2.2[54] 0.6[54] 1.7[54]
430 86.2 54.5 1.28 0.75 0.79
423 83.2 58.5 2.01 2.14 2.41
420 87.5 61.1 2.02 2.91 3.39
443[55] 99.2[56] 75.8[56] 1.17[57] 2.7[58] 3.6[59]
422 84.3 56.1 n/a n/a n/a
Band gap, magnetic moment, bulk modulus are given in eV, bohr magneton, and GPa, respectively.
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 3. Cohesive energies of the X23 benchmark set (kcal/mol). PBE Hydrogen bonded Acetic acid Ammonia Cyanamide Cytosine Ethylcarbamate Formamide Imidazole Oxalic acid (a) Oxalic acid (b) Succinic acid Uracil Urea VdW bonding 1,4-cyclohexanedione Adamantane Anthracene Benzene CO2 Hexamine Naphthalene Pyrazine Pyrazole Triazine Trioxane Overall MD MAD MPD MAPD
PBE1LRD
revPBE
revPBE1LRD
Ref.
48.1 28.2 70.2 99.8 47.6 58.9 53.4 83.6 86.8 82.5 84.0 78.9
(224.7) (29.0) (29.5) (270.0) (238.7) (220.3) (233.4) (212.7) (29.3) (247.8) (251.7) (223.6)
83.2 45.3 97.9 166.3 95.7 86.3 98.0 126.4 127.6 144.2 144.4 112.8
(10.4) (8.1) (18.2) (23.5) (9.4) (7.1) (11.2) (30.1) (31.5) (13.9) (8.7) (10.3)
22.1 11.4 43.7 49.3 18.1 33.8 24.4 44.1 48.3 37.9 37.7 48.5
(250.7) (225.8) (236.0) (2120.5) (268.2) (245.4) (262.4) (252.2) (247.8) (292.4) (298.0) (254.0)
70.2 37.0 88.4 148.5 82.0 75.2 90.5 107.1 108.8 123.4 125.7 96.4
(22.6) (20.2) (8.7) (221.3) (24.3) (24.0) (3.7) (10.8) (12.7) (26.9) (210.0) (26.1)
72.8 37.2 79.7 169.8 86.3 79.2 86.8 96.3 96.1 130.3 135.7 102.5
35.1 24.3 3.2 5.1 9.9 15.9 3.1 17.1 43.0 14.8 19.4
(253.5) (273.7) (2109.5) (246.6) (218.5) (270.3) (278.6) (244.2) (234.7) (246.9) (247.0)
102.2 90.3 121.7 61.8 26.0 98.9 91.4 70.9 84.3 66.2 68.5
(13.6) (20.9) (9.0) (10.1) (22.4) (12.7) (9.7) (9.6) (6.6) (4.5) (2.1)
1.8 236.6 247.6 219.8 20.5 215.2 235.1 214.4 18.4 213.1 24.0
(286.8) (2106.0) (2160.3) (271.5) (228.9) (2101.4) (2116.8) (275.7) (259.3) (274.8) (270.4)
89.7 76.4 117.2 59.2 23.9 90.2 88.1 64.2 77.3 60.1 60.1
(1.1) (7.0) (4.5) (7.5) (24.5) (4.0) (6.4) (2.9) (20.4) (21.6) (26.3)
88.6 69.4 112.7 51.7 28.4 86.2 81.7 61.3 77.7 61.7 66.4
242.4 42.4 252.2% 52.2%
10.9 11.5 13.3% 14.2%
0.1 6.0 0.6% 7.0%
274.1 74.1 290.6% 90.6%
Differences from reference values are shown in parentheses.
locations in the solid state were optimized using the LRD method combined with the PBE or revPBE functional. The dissociated state was modeled by one molecule in the center of a cubic unit cell with a cell length of 30 Bohr. The ultrasoft pseudopotentials determined by the Rabe–Rappe–Kaxiras– Joannopoulos (RRKJ) method[50] were adopted for all elements. The plane-wave cutoff and charge density cutoff were 49 and 400 Ry, respectively. The k-point sampling was performed in the same manner as the benchmark study of 21 molecular crystals (C21 set) by Otero-de-la-Roza and Johnson.[61] The calculated and reference cohesive energies are summarized in Table 3. The cohesive energies without a dispersion correction were evaluated using the lattice structure optimized using the LRD method. Although the sign of the cohesive energies of hydrogen-bonded systems is reproduced without any dispersion correction, the application of the LRD method significantly reduced deviations. Considering that the dispersion-uncorrected GGA functional can evaluate a large portion of the interaction energy of hydrogen-bonded molecular dimers,[17] the result indicates the importance of dispersion correction for molecular crystals. As expected, the LRD method reduced errors in the cohesive energies of vdW bonding systems. The PBE1LRD approach somewhat overestimates the cohesive energies on average. The overbinding may be related to that in the cohesive energy of neon shown in Figure 1, where PBE1LRD overestimates the cohesive energy due to the
attractive nature without any dispersion correction. The ME of revPBE1LRD is nearly zero, indicating that three empirical parameters adjusted for the cohesive energies of rare gases are transferrable to molecular crystals. The overall mean absolute deviation (MAD) of the revPBE functional is reduced from 74.14 to 5.97 kJ/mol, and the overall MAPD is reduced from 90.60 to 6.97%. This accuracy is comparable to the MBD method combined with the PBE functional: 5.92 kJ/mol (8.05%).[60] We compare the accuracies of PBE1LRD and revPBE1LRD with those of other dispersion-corrected functionals in terms of cohesive energies, optimized lattice parameters, and lattice volumes. Here, we consider the C21 benchmark set,[61] which
Table 4. Mean and mean absolute deviations from experimental lattice parameters of cohesive energies (kJ/mol), optimized cell lengths (A˚), angles ( ), and volumes (A˚3) of the C21 benchmark set. Cohesive energy
Length
Angles
Volume
MD
MAD
MD
MAD
MD
MAD
MD
MAD
PBE1LRD 10.7 PBE-D[61] 8.2 PBE-TS[61] 16.8 PBE-XDM[61] 0.2 revPBE1LRD 0.2 vdW-DF[61] 8.8
11.3 9.0 17.0 5.4 6.0 10.2
0.00 20.07 0.04 0.16 0.09 0.29
0.11 0.11 0.10 0.20 0.14 0.31
20.22 0.05 20.06 20.15 20.05 20.06
0.24 0.19 0.14 0.29 0.12 0.11
4.2 211.7 7.3 31.2 11.7 48.0
8.6 14.4 10.1 31.2 14.0 48.0
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 2. Adsorption structures and unit cells of n-butane on (a) M (111) (M 5 Ag, Au, Pd, and Pt) and (b) Cu (111) surfaces. (c) Side view of the adsorbed structure with the Hmetal interaction shown by the dashed line.
consist 21 molecular crystals also included in the X23 set. All results of our calculations are summarized in Supporting Information. The deviations of the calculated results are summarized in Table 4 as well as previously reported values. In addition to PBE-based methods, vdW-DF is listed because the exchange part of the revPBE functional is included. As for cohesive energies, PBE1LRD and revPBE1LRD showed the same tendency as the X23 set. PBE-XDM and revPBE1LRD gave nearly zero MDs and are in a good agreement with the reference. The other methods resulted in overestimation and larger MADs. In terms of structural parameters, PBE1LRD reproduced them well while revPBE1LRD slightly overestimated the lattice lengths and volumes. This originates in the strong repulsive character of the revPBE exchange at short internuclear distances and remarkably rapid decay at long distances.[15] The overestimation is not so large, smaller than that of PBE-XDM and vdW-DF. We speculate that the GGA correlation, which is not included in vdW-DF, results in decreased deviations to some extent. The lattice angles were well reproduced by revPBE1LRD and vdW-DF. The present numerical assessment of cohesive energies and structural parameters support the usefulness of the LRD method combined with the revPBE functional for molecular crystals. Physical adsorption In the present subsection, n-butane on a transition metal surface is chosen as an example of physical adsorption. This system was previously studied by Morikawa and coworkers,[62,63] in which adsorptions on Cu(100), Cu(111), Au(111), and Pt(111) 6
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
surfaces were calculated using the PBE functional[7] and the vdW-DF series.[21,22] Here, we apply the revPBE1LRD method to adsorptions on the (111) surfaces of Pd and Ag in addition to those of Cu, Au, and Pt. The unit cell of the adsorption structure is shown in Figure 2a. We introduced a slab model, consisting of four layers of metal atoms, to represent the (111) surface. Interatomic distances of metal atoms werepset ffiffiffi to experimentally determined pffiffiffi bulk structures. The 73 3 surface unit cells were adopted to represent the fully adsorbed state. The intersurface distance ˚ . The geometry of n-butane was optimized was set to 20 A with fixed layers and lattice constants. However, an appropriate adsorption structure was not found in the case of Cu. We speculate that this is a result of intermolecular repulsion due pffiffiffi to the small surface unit cell. Instead, we used the 732 surface unit cell shown in Figure 2b for Cu. The ultrasoft
Table 5. Adsorption energies (eV) and distances (A˚) of n-butane on the (111) surface of transition metals. Adsorption energy [63]
revPBE1LRD vdW-DF2 Pd Pt Cu Ag Au
0.60 0.58 0.55 0.49 0.48
– 0.61 0.37 – 0.52
Adsorption distance Expt. – 0.53[64] 0.51[65] – 0.42[66]
revPBE1LRD vdW-DF2[63] Est.[a] 3.41 3.59 3.53 3.69 3.79
– 3.88 3.95 – 3.97
– 3.56 – – 3.52
[a] Reported in Ref. [62], where the distances are fit to reproduce the red-shifts of the CAH stretching mode and the work function changes due to the adsorption.
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 3. Difference in electron densities along the CAHmetal interactions given by the revPBE1LRD method.
pseudopotentials determined by the RRKJ method[50] were adopted for all elements. The (4 3 6 3 1) k-point grid was introduced based on the Monkhorst–Pack scheme. In Table 5, the calculated adsorption energies and distances are summarized together with the results of the previous studies.[62,63] Here, the adsorption distance is defined as the average of distances between each carbon atom and the plane including the nuclei of surface atoms. The adsorption energies determined using the revPBE1LRD method are in good agreement with the experimental values. Compared to vdW-DF2,[22] which consists of the refitted Perdew-Wang 86 exchange,[67] LDA correlation,[4] and nonlocal correlation functionals,[22] the present scheme provided shorter adsorption distances and smaller differences from the distances fitted to experimental
results. The adsorption energies in the case of Ag and Au are slightly smaller than those of the other metals. This is consistent with the longer adsorption distances of Ag and Au than those of the other metals. In addition to adsorption energies and structures, electron transfer between the CAH bonds and the metal layers is worth investigating. According to the previous report,[68] this is interpreted as a result of mixing the CAH bonding and antibonding orbitals and the d band of transition metals. In principle, the LRD method is suitable for charge-transferred systems because this method is dependent on the electron distribution obtained as a solution of the DFT calculation. Previously, we confirmed that the electronic excitation changes atomic polarizabilities given by the LRD method reasonably.[42] Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
7
FULL PAPER
WWW.C-CHEM.ORG
Figure 3 visualizes the electron transfer obtained by the revPBE1LRD method. All systems clearly show the electron transfer from the atomic region of hydrogen and metal to the intermediate region. A larger transfer is confirmed in the case of Pd and Pt. This is because Pd and Pt have larger vacancies in the d band than Cu, Ag, and Au. As for the noble metals, the charge transfer is larger in the order of Cu > Ag > Au. This is consistent with the reactivity of the noble metals. These tendencies for difference densities may cause the difference in the adsorption energies. The results indicate that the revPBE1LRD method can provide not only accurate adsorption energies but also reasonable adsorption structures and electron transfers.
Conclusion We implemented the LRD method in the Quantum ESPRESSO and assessed the combination with the PBE and revPBE functionals. This is based on the feasibility of the LRD method, which does not depend on the electronic basis functions such as Gaussian and plane-wave basis functions but the electron density. Numerical applications confirmed the accuracy of the LRD method even in the case of periodic systems. The LRD method improved the underestimation of cohesive energies not only in rare gases but also in strongly bonded elements. The LRD method also improved the accuracy of the cohesive energies and lattice parameters of molecular crystals. Furthermore, we confirmed that not only the physical adsorption energy but also the electron transfer in n-butane on metal surfaces were calculated reasonably.
Acknowledgment Some of the present calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). Keywords: density functional theory dispersion force periodic boundary condition plane-wave basis set physical adsorption
How to cite this article: Y. Ikabata Y. Tsukamoto Y. Imamura, H. Nakai. J. Comput. Chem. 2014, DOI: 10.1002/jcc.23807
]
Additional Supporting Information may be found in the online version of this article.
[1] [2] [3] [4] [5] [6] [7] [8] [9]
P. Hohenberg, W. Kohn, Phys. Rev. 1964, 136, B864. W. Kohn, L. J. Sham, Phys. Rev. 1965, 140, A1133. J. C. Slater, Phys. Rev. 1951, 81, 385. S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. A. D. Becke, Phys. Rev. A 1988, 38, 3098. C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. Y. Zhang, W. Yang, Phys. Rev. Lett. 1998, 80, 890. J. P. Perdew, S. Kurth, A. Zupan, P. Blaha, Phys. Rev. Lett. 1999, 82, 2544. [10] J. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Phys. Rev. Lett. 2003, 91, 146401.
8
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
[11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
[44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55]
[56] [57] [58]
A. D. Becke, J. Chem. Phys. 1993, 98, 5648. C. Adamo, V. Barone, J. Chem. Phys. 1999, 110, 6158. J. Heyd, G. E. Scuseria, M. Ernzerhof, J. Chem. Phys. 2003, 118, 8207. A. Seidl, A. G€ orling, P. Vogl, J. A. Majewski, M. Levy, Phys. Rev. B 1986, 53, 3764. F. O. Kannemann, A. D. Becke, J. Chem. Theory Comput. 2009, 5, 719. N. Marom, A. Tkatchenko, M. Rossi, V. V. Gobre, O. Hod, M. Scheffler, L. Kronik, J. Chem. Theory Comput. 2011, 7, 3944. J. Klimes, A. Michaelides, J. Chem. Phys. 2012, 137, 120901. A. J. Cohen, P. Mori-Sanchez, W. Yang, Chem. Rev. 2012, 112, 289. J. F. Dobson, B. P. Dinte, Phys. Rev. Lett. 1996, 76, 1780. Y. Andersson, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Lett. 1996, 76, 102. M. Dion, H. Rydberg, E. Schr€ oder, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Lett. 2004, 92, 246401. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, D. C. Langreth, Phys. Rev. B 2010, 82, 081101. I. Hamada, Phys. Rev. B 2014, 89, 121103. O. A. Vydrov, T. Van Voorhis, Phys. Rev. Lett. 2009, 103, 063004. O. A. Vydrov, T. Van Voorhis, J. Chem. Phys. 2010, 133, 244103. M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 2001, 114, 5149. Q. Wu, W. Yang, J. Chem. Phys. 2002, 116, 515. S. Grimme, J. Comput. Chem. 2004, 25, 1463. S. Grimme, J. Comput. Chem. 2006, 27, 1787. S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104. L. Goerigk, S. Grimme, J. Chem. Theory Comput. 2011, 7, 291. A. D. Becke, E. R. Johnson, J. Chem. Phys. 2005, 122, 154104. A. Tkatchenko, M. Scheffler, Phys. Rev. Lett. 2009, 102, 073005. A. Tkatchenko, R. A. DiStasio, Jr., R. Car, M. Scheffler, Phys. Rev. Lett. 2012, 108, 236402. A. Tkatchenko, A. Ambrosetti, R. A. DiStasio, Jr., J. Chem. Phys. 2013, 138, 074106. T. Sato, H. Nakai, J. Chem. Phys. 2009, 131, 224104. T. Sato, H. Nakai, J. Chem. Phys. 2010, 133, 194101. Y. Ikabata, T. Sato, H. Nakai, Int. J. Quantum Chem. 2013, 113, 257. H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540. R. Kar, J.-W. Song, T. Sato, K. Hirao, J. Comput. Chem. 2013, 34, 2353. Y. Ikabata, H. Nakai, Chem. Phys. Lett. 2013, 556, 386. Y. Ikabata, H. Nakai, J. Chem. Phys. 2012, 137, 124106. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, R. M. Wentzcovitch, J. Phys.: Condens. Matter 2009, 21, 395502. A. J. Stone, In Theory of Intermolecular Forces; Clarendon: Oxford, 1996. O. A. Vydrov, T. van Voorhis, J. Chem. Phys. 2009, 130, 104105. K. Rapcewicz, N. W. Ashcroft, Phys. Rev. B 1991, 44, 4032. A. D. Becke, J. Chem. Phys. 1988, 88, 2547. C. Kittel, In Introduction to Solid State Physics, 8th ed.; Wiley: New York, 2005. D. Vanderbilt, Phys. Rev. B 1990, 41, 7892. A. M. Rappe, K. M. Rabe, E. Kaxiras, J. D. Joannopoulos, Phys. Rev. B 1990, 41, 1227. N. Troullier, J. L. Martins, Phys. Rev. B 1991, 43, 1993. F. Birch, Phys. Rev. 1947, 71, 809. B. G. Streetman, S. Banerjee, In Solid State Electronic Devices, 5th ed.; Prentice Hall: New Jersey, 2000. R. M. Bozorth, In Ferromagnetism; D. Van Nostrand Company, Inc.: Princeton, NJ, 1956. M. E. Levinshtein, S. L. Rumyantsev, In Handbook Series on Semiconductor Parameters; M. Levinshtein, S. Rumyantsev, M. Shur, Eds.; Vol. 1; World Scientific: Singapore, 1996. K. H. Hellwege, In Landolt-Bornstein New Series, Group III: Condensed Matter, Vol. 17a; Springer: Berlin, 1966. D. N. Batchelder, D. L. Losee, R. O. Simmons, Phys. Rev. 1967, 162, 767. O. G. Peterson, D. N. Batchelder, R. O. Simmons, Phys. Rev. 1966, 150, 703.
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
[59] [60] [61] [62] [63] [64]
D. L. Losee, R. O. Simmons, Phys. Rev. 1968, 172, 944. A. M. Reilly, A. Tkatchenko, J. Chem. Phys. 2013, 139, 024705. A. Otero-de-la-Roza, E. R. Johnson, J. Chem. Phys. 2012, 137, 054103. Y. Morikawa, H. Ishii, K. Seki, Phys. Rev. B 2004, 69, 041403. K. Lee, Y. Morikawa, D. C. Langreth, Phys. Rev. B 2010, 82, 155461. S. L. Tait, Z. Dohnalek, C. T. Campbell, B. D. Kay, J. Chem. Phys. 2006, 125, 234308. [65] R. Z. Lei, A. J. Gellman, B. E. Koel, Surf. Sci. 2004, 554, 125. [66] S. M. Wetterer, D. J. Lavrich, T. Cummings, S. L. Bernasek, G. Scoles, J. Phys. Chem. B 1998, 102, 9266.
FULL PAPER
[67] E. D. Murray, K. Lee, D. C. Langreth, J. Chem. Theory Comput. 2009, 5, 2754. € om, L. Triguero, M. Nyberg, H. Ogasawara, L. G. M. Pettersson, [68] H. Ostr€ A. Nilsson, Phys. Rev. Lett. 2004, 91, 046102.
Received: 26 July 2014 Revised: 6 November 2014 Accepted: 8 November 2014 Published online on 00 Month 2014
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23807
9
