Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Theory, practice, and application of the van der
Waals density functional
Ikutaro HAMADA
Department of Precision Engineering, Graduate School of Engineering, Osaka University
2-1 Yamadaoka, Suita, Osaka 565-0871, Japan

Abstract
The van der Waals density functional is a classof the density functional, which enables the description of the dispersion forces within theframework of density functional theory. In
this article, theory and practice of the van der
Waals density functional are reviewed, and applicability of a high-accuracy van der Waals
density functional to various systems is demonstrated.

1

Introduction

Density functional theory (DFT) [1, 2] within
the local density (LDA) and semilocal (GGA)
approximations has been proven to be accurate in describing covalent, ionic, and metallic bondings, but has also been known to beunable to describe the dispersion forces (or often referred to as the van der Waals forces)[3] properly. This is because the dispersion
forces originate from the instantaneous dipoleinduced dipole interaction due to the electron
fluctuation or the dynamic and nonlocal electron correlation, which are inherently missingin LDA and GGA. Wave function based theory and many body perturbation theory are
shown to describe the dispersion interaction
very accurately, but the computational cost is
prohibitively large, hindering their application
to complex systems, such as surfaces and interfaces.
There have been several proposals to correct the missing dispersion forces by adding
the term proportional to 1/r6 (and higherorder terms) to the total energy within semilocal DFT, where r is the interatomic or inter-fragment distance. This type of approach was
initiated in the field of the quantum chemistry(QC) [4], and has been established [5] and improved over the years by introducing accurate
dispersion coefficients [6] and damping functions [7], as well as the environmental effects
on the dispersion coefficients through the local
electron density [8, 9] or coordination numbers[6]. Furthermore, accuracy of such pair-wisecorrection methods have been improved by introducing the three-body term[6] or the manybody dispersion[10]. See Refs. [11, 12, 13] forrecent reviews.

The van der Waals density functional (vdWDF) [14, 15, 16, 17, 18] is a versatile methodthat allows one to describe different interactions including the dispersion forces within the
DFT framework. It is grounded on the electron
gas tradition and the many-body perturbation
theory as other constraint-based exchangecorrelation (XC) functionals, and formally has
no empirical or fitting parameters to reproduce
the results of experiments or accurate theoretical calculations. The physical foundation of
vdW-DF has been established [19], and theoretical and algorithmic improvements make ita method of choice in theoretical studies of materials and molecules.

In this article, I give a brief overview of thetheory of vdW-DF and discuss a technical issue related to the use of pseudopotentials withapplications. The Hartree atomic unit is usedthroughout unless otherwise specified.

15

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

2

Method2.1

TheoryThe XC energy, which is responsible for thequantum-mechanical electron-electron interaction, can be expressed exactly by using the adiabatic connection fluctuation dissipation theorem (ACFDT) [20, 21] as
Z 1

dλ
λ

Z ∞

i
du h λtr χ (iu)V λ − Eself ,
2π
0
0
(1)
where λ is the adiabatic parameter, whichscales the electron-electron Coulomb interaction as V λ (r, r′ ) = λ/|r − r′ |, χλ (iu) is the
reducible density-density response function at
the adiabatic parameter λ and imaginary frequency u = −iω, and Eself is the electronelectorn self-energy defined by Eself = tr [nV ]with V = V λ=1
the trace is defined
RR . Here,′as tr[AB] =drdr A(r, r′ )B(r′ , r) for arbitrary operators Â and B̂. χλ (ω) is definedby δn(ω) = χλ (ω)δvext (ω), where δn(ω) andδvext (ω) are the induced electron density andthe external potential at the frequency ω, respectively.
The vdW-DF is constructed based on
ACFDT as follows: Exc may be split as
Exc = −

0
Exc = Exc
+ Ecnl ,(2)0 and E nl are the semilocal exchangewhere Excccorrelation and the nonlocal correlation energyfunctionals, respectively, and in vdW-DF, a
tractable form of Ecnl has been developed basedon several physical constraints.
We first introduce χ̃λ , which is defined byλ (ω) with δV λ (ω) being theδn(ω) = χ̃λ (ω)δVscr
scr
screened potential. χλ (ω) and χ̃λ (ω) are related via the Dyson-like equation χλ = χ̃λ +χ̃λ V λ χλ . Then, the full-potential approximation (FPA) is introduced, which sets χ̃λ =χ̃λ=1 and is exact at the van der Waals asymptote (|r − r′ | → ∞). This enables the analytic
integration over λ and leads to Exc in FPA as
Z ∞io
du n hExc =tr ln 1 − χ̃λ=1 (iu)V− Eself .
2π
0(3)This expression has a formal similarity withthat obtained within the random phase ap-

proximation (RPA):Z ∞
du  
tr ln 1 − χ0 (iu)V − Eself ,
Exc =
2π
0
(4)
0
where χ is the irreducible density-density response function, which is connected to χλ via
a Dyson equation χλ = χ0 + χ0 V λ χλ within
RPA. Exc within FPA can also be expressed
as
Z ∞
du
Exc =
tr {ln [∇ · ε(iu)∇G]} − Eself ,2π
0(5)where ε is the dielectric function and G is the
Coulomb Green’s function, i.e., G = −4πV .
0 is approximated with
To extract Ecnl , Exc
that of homogeneous electron gas within FPA,
which is also called the internal XC energy, as
Z ∞
du
in
Exc =
tr {ln [ε(iu)]} − Eself ,(6)2π0
to obtainZ ∞dunl
tr {ln [∇ · ε(iu)∇G] − ln [ε(iu)]} .
Ec =
2π
0
(7)
nl
It is apparent that by construction, Ec vanishes in the uniform electron-gas limit.
The thus obtained Ecnl is expanded in termsof a response function S = ln ε ≈ 1 − ε−1 upto the second order as
Z ∞o
du n 2nl
Ec =tr S (iu) − [∇S(iu) · ∇G]24π
0(8)and in the plane-wave representation as
Z ∞
Z
Z
 i
du
dq
dq′ h
′ 2
Ecnl =
1
−
q̂
·
q̂
2π
(2π)3
(2π)3
0
× Sq,q′ (iu)Sq′ ,q (iu),(9)
where q̂ = q/|q| andS(ω) = S(r, r′ , ω)ZZdq′ −iq·r
dq
′ ′
=
e
Sq,q′ (ω)eiq ·r
3
3
(2π)
(2π)
(10)
This Ecnl is made tractable by applying
theplasmon-poleapproximation to Sq,q′ (ω), which is expressed as

16

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyoih
Sq,q′ (ω) = 12 S̃q,q′ (ω) + S̃−q,−q′ (ω) andZ
S̃q,q′ (ω) =×with the GGA exchange enhancement factor
for the internal XC functional given by′

dre−i(q−q )·r

ωp (r)
.
[ω + ωq (r)] −ω + ωq′ (r)Fxin (s) = 1 −(11)Here, ωp (r) = 4πn(r) is the plasma frequency,
where n(r) is the electron density, and ωq (r)is the plasmon dispersion. The specific form of
the response function S has been determinedto satisfy the following constraints: (i) thetime-reversal symmetry, (ii) the f -sum rule,(iii) cancellation of the self-energy, and (iv)charge conservation.
By writing the internal XC energy as
Z
in
(12)
Exc
= drϵinxc (r)n(r),the internal XC energy per electron ϵin
xc reads


Z
dq
1
2
ϵin
(r)
=
π
,
(13)
−
xc
(2π)3 ωq (r) q 2
where q = |q| and the second term in theright-hand side is obtained by evaluating theself-energy term explicitly. The local plasmon
dispersion is modeled byωq (r) =q2
2h[q/q0 (r)](14)

with h being the switching function which controls the behavior of the plasmon depending on
the wave vector and q0 (r), an inverse lengthscale which depends on the electron density
and its gradient |∇n(r)|. The functional formof h is determined in such a way that ωq (r) satisfies the small- and large-q cases (see below).
By imposingZ ∞3dy [h(y) − 1] = ,(15)
40
q0 is related to ϵin
xc as
ϵin
xc (r) = −3q0 (r).
4π(16)

In vdW-DF, ϵin
xc is made to be the GGA exchange plus LDA correlation functionals as
in
LDA
ϵin
(r) + ϵLDA,
xc (r) = Fx (s)ϵx
c

(17)

Zab 2
s ,
9

(18)

= −3/4πkF with kF = (3π 2 n)1/3
where ϵLDAxand s = |∇n|/kF n. By this choice, q0 canbe interpreted as kF modulated by the enLDA . In the first version ofergy ratio ϵin
xc /ϵx
vdW-DF (vdW-DF1), Zab = −0.8491 basedon the second-order gradient expansion for the
slowly varying electron gas [22, 23], while inthe second version of vdW-DF (vdW-DF2),
Zab = −1.887 based on the second-order largeN expansion for the neutral atom [17].
The switching function h for the plasmon
dispersion is determined to satisfy Eq. (15) andthe following constraints. (i) h(y) ≈ γy 2 + · · ·in the small-y limit, where γ is an arbitrary
constant to ensure h(0) = 0, which corresponds to the charge conservation of the spherical XC hole model of the internal XC energy.
(ii) h(y) → 1 in the large-y limit, which corresponds to ωq → q 2 /2 in the large q-limit. This
constraint ensures the cancellation of the divergent term of Eself in the q → 0 limit in ϵin
xc .
In Refs. [15, 17], the following form for h wasused
h(y) = 1 − exp(−γy 2 )(19)with γ = 4π/9, implying that the Gaussian
form of the XC hole is used in ϵinxc [19]. Thereis a freedom in the functional form of h and it
is not limited to Eq. (19). Indeed, Berland etal. [24] proposed a functional form, which can
improve the asymptotic behavior (and the socalled C6 coefficient) of vdW-DF . To improve
the accuracy of Ecnl as a part of the development of the third version of vdW-DF (vdWDF3), Chakraborty et al. [25] proposed thefollowing functional form
1h(y) = 1 −.
(20)
2
1 + γy + γ 2 y 4 + αy 8
See Table 1 in Ref. [25] for the parameters γ
and α. Note however, in the development of
vdW-DF3, parameters in the functionals were
optimized to the reference high-accuracy QC
results of a large molecular data set. These
parameters may vary when a plane-wave pseduopotential program is used in the optimization, depending on the pseudopotentials used.

17

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
2.0

After tedious but straightfoward calculations, one arrives at the following form of thenonlocal correlation functional
ZZ
1
nl
Ec =
drdr′ n(r)ϕ(d, d′ )n(r′ ),
(21)
2

δ=0.0
δ=0.5
δ=0.9

1.0
0.5

2

D φ (arb. unit)

1.5

where d = q0 (r)|r − r′ | and d′ = q0 (r′ )|r − r′ |.
The kernel function is given by
Z ∞
Z ∞
2
′
2
ϕ(d, d ) = 2
dbb2 W (a, b)
daa
π 0
(22)
0
′
′
× T [ν(a), ν(b), ν (a), ν (b)],where
W (a, b) =2 (3 − a2 )b sin a cos b+ (3 − b2 )a cos a sin b
2

2+ (a + b − 3) sin a sin b(23)−3ab cos a cos b]

-0.5
-1.0
0

2

4
D

6

but in practice, it results in less accurate results. By allowing crossover in the exchange
functional and letting
0
in
GGA
Exc
= Exc
+ ∆Ex = Exc+ EcLDA ,(24)

with ν(y) = y 2 /2h(y/d) and ν ′ (y) =y 2 /2h(y/d′ ) (d and d′ are defined above). Calculating the ϕ on the fly is impractical, as thetime-consuming double integral over the variables a and b is included. In practice, ϕ iscalculated and tabulated in advance as a function of D = (d + d′ )/2 and δ = (d − d′ )/(d + d′ )(0 ≤ D < ∞ and 0 ≤ |δ| < 1) and the ϕvalue is calculated on demand by interpolation. This approach works quite well, thankto the fact that ϕ is a smooth function of Dand δ (See Fig. 1). Note that ϕ should berecalculated when different switching function
h is employed, although the shapes of the hfunctions look similar, but y 2 /h(y), which isrelevant to the plasmon dispersion, differs significantly [24, 25].
Remaining task in the development of vdWDF is to determine the semilocal XC funcin as E 0 is the most consistenttional. Using Exc
xc
approach within the framework of vdW-DF,

8

Figure 1: The kernel ϕ times D2 as a functionof D at various δ’s obtained using an in-house
kernel generation code [26].

and


1
1
1
+
T (w, x, y, z) =2 w+x y+z1
×(w + y)(x + z)1+,
(w + z)(y + x)

0.0

(25)

the local correlation energy can be consistentin , and the semilocal exchangewith that of Exc
energy has been investigated.
In vdW-DF for layered materials [14] andgeneral geometry vdW-DF (vdW-DF1) [15],the revised Perdew-Burke-Ernzherof (PBE)[27] exchange (revPBE) of Zhang and Yang[28] was employed because it does not show
spurious binding of rare gas dimers from exchange only. However, it turns out that the
revPBE exchange tends to overestimate thePauli repulsion, leading to a too large separation between the fragments. In searching for more accurate exchange energy functional, Klimeš et al. [29] proposed a set of
exchange functionals, leading to optPBE-vdW
and optB88-vdW, based on the optimization
to the reference QC results. Cooper [30] proposed an exchange functional (C09) which satisfies the following constraint: (i) The gradient expansion approximation (GEA) form inthe slowly varying electron density (small-s)
limit, i.e., Fx (s) = 1 + µs2 with µ = 0.0864[31], which is mainly responsible for the covalent interaction. (ii) The revPBE behavior
in the homogeneous electron density (large-s)
limit which is relevant to the vdW bonding.
Although the µ value is not exactly the sameas that of GEA for the slowly varying electron

18

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Although the precise performance of the
vdW-DF functional depends on the combination of exchange and nonlocal correlation functionals, newly proposed functionals are shown
to be more accurate than the original one,

2.5

2.0
Fx

gas (µGEA = 10/81 [32]), this work demonstrated the importance of GEA in predicting
accurate interaction energies and atomic geometries. Murray et al. [33] showed that theexchange enhancement factor should be proportional to s2/5 in the large-s limit. Theydeveloped the revised PW86 [31] functional(PW86R), which obeys GEA in the small-s
limit and has the s2/5 behavior in the large-s
limit for the exchange enhancement factor, andpaired with a new Ecnl , to propose vdW-DF2[17]. Klimeš et al.[34] developed the optB86b
exchange functional (and optB86b-vdW) in
which the exchange enhancement factor obeys
the gradient expansion approximation (GEA)in the slowly varying density (small-s) limit
with µGEA and has the s2/5 behavior in thelarge s-limit. Berland and Hyldgaard [35] proposed an exchange functional (LV-PW86R),
which is exactly the same as that of the internal exchange-correlation functional in thesmall-s region, while it reproduces the PW86R
behavior in the large-s region. This semilocal exchange functional is the most consistent
with that of the internal XC functional (consistent exchange), and thus the resulting vdWDF functional is named vdW-DF-cx. Hamada[37] proposed an alternative one (B86R) basedon the B86b [36] exchange functional, whose
enhancement factor obeys GEA in the small-s
limit and reproduces the original B86b behavior (proportional to s2/5 ) in the large-s region.
The B86R exchange is paired with Ecnl of vdWDF2 and dubbed rev-vdW-DF2. In the latest
development of vdW-DF (vdW-DF3) [25], B88[38] and B86b [36] forms of the exchange enhancement factors are used. They are bothdesigned to obey GEA in the small-s limit,
and the remaining parameters are optimized to
the reference QC results. The proposed GGA
enhancement factors for vdW-DF are shownin Fig. 2, along with those for Perdew-BurkeErnzerhof (PBE) [27] and PBEsol [39], the latter of which also obeys GEA in the slowlyvarying density limit.

PBE
PBEsolrevPBE
optB88
optB86b
PW86R
LV-PW86R
B86R
C09

1.5

1.0
0

1

2

3

4

5
s

6

7

8

9

10

Figure 2: GGA enhancement factors (Fx ) proposed for various exchange-correlation functionals as a function of reduced density gradient (s).
and are versatile tools for accurate prediction of structures and energetics of materials,
molecules, and interfaces. See Refs. [24, 25,
40, 41, 42, 43] for assessments of the different
vdW-DF functionals on various materials andmolecules.

2.2

Pseudopotentials

In electronic structure calculations, the pseudopotential and projector augmented wave
(PAW) [44] methods are often used, in which(pseudo)potentials generated using certain XC
functionals are used. As a rule of thumb, the
XC functional used to generate pseudopotentials should be consistent with that for the
calculations of materials and molecules, but
this issue has not been discussed in vdWDF calculations. This is presumably because
the appropriate XC functionals for atomic systems have not known nor been discussed forthe vdW-DF calculations. However, it should
not be overlooked, given that the vdW-DF
employs different exchange and nonlocal correlation functionals and is not a simple dispersion correction to a semilocal XC functional. Callsen and Hamada [45] proved that
Ecnl vanishes in an isotropic system and showedthat the XC functional to be used to generate(pseudo)potentials for specific vdW-DF functional is uniquely determined. For instance,
the B86R exchange plus LDA correlation functional should be used to generate potentials

19Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo40

exchange
revPBE
PW86R
optB88
optB86b
C09
LV-PW86R
B86R

correlation
LDA
LDA
LDA
LDA
LDA
LDA
LDA

for the rev-vdW-DF2 calculations of molecular/extended systems. See Table 1 for the exchange and correlation functionals to be usedin the pseudo- and PAW-potential generationfor vdW-DF calculations.

3

Application

In the following, applications of rev-vdW-DF2
to selected systems are presented. Emphasis is put onto the comparison of the results
obtained using (psuedo)potentials generated
with different XC functionals, i.e., the PBE
and the B86R exchange plus LDA correlation(B86Rx+LDAc), the latter of which is compatible with rev-vdW-DF2. All the calculations were based on the PAW method and selfconsistent vdW-DF as implemented [23, 46] inthe Quantum-ESPRESSO [47] code. The
proper extension of vdW-DF to the spin polarized system [48, 49] was adopted when necessary. The PAW potentials were generatedusing the input files supplied in the PSLIBRARY [50].
I first calculated the interaction energies for
the S22 dataset [51], a set of 22 noncovalently bonded molecular duplexes. The calculations were performed at the S22 geometries[51] without any further optimization. Figure 3 depicts the deviations of the binding energies obtained using PBE and B86Rx+LDA
potentials, with respect to the reference results
obtained by the coupled cluster calculations
with singlet, doublet, and perturbative triplet
excitations [CCSD(T)] [52], along with the difference between the deviations obtained using

20

rev-vdW-DF2 (B86Rx+LDAc)
rev-vdW-DF2 (PBE)

0
-20
-40

∆Eb (B86Rx+LDAc) - ∆Eb (PBE) (meV)

vdW-DF functional
vdW-DF1
vdW-DF2
optB88-vdW
optB86b-vdW
vdW-DFC09x
vdW-DF-cx
rev-vdW-DF2∆Eb (meV)

Table 1: vdW-DF functional and corresponding exchange and correlation functionals to be
used in the atomic calculation and pseudo- andPAW-potential generation.

-60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

2.0
1.5
1.0
0.5
0.0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
S22 index

Figure 3: Deviations of the binding energies
obtained using rev-vdW-DF2 with different
PAW potentials with respect to the reference
CCSD(T) results (∆Eb , upper panel) and difference between the ∆Eb ’s obtained with different potentials (lower panel).
different potentials. The binding energies obtained with PBE and B86Rx+LDA potentials
are virtually the same (within the difference of
2 meV), and the mean absolute deviations are
20.5 and 20.8 meV for the former and the latter, respectively, validating the use of the PBE
potentials in the vdW-DF calculations.
I then calculated the interlayer binding energy of graphite as a function of interlayer
distance, as shown in Fig. 4. The equilibrium in-plane lattice constant is 2.46 Å, andinterlayer distance and binding energy are
3.33 Å and 59.4 meV/atom, respectively, using both potentials. The interlayer binding
energy is slightly overestimated as comparedwith quantum Monte-Carlo (56±5 meV/atom[53]), ACFDT-RPA (48 meV/atom [54]), andexperiment (52±5 meV/atom [55]) but is inreasonable agreement with them.
I also performed the calculations of a benzene crystal. Figure 5 shows the binding
energy of the benzene crystal in the Pbca
symmetry as a function of the unit cell volume. The calculated equilibrium volume is
460 Å3 with both potentials, and the binding energies are 532 and 531 meV/moleculewith B86Rx+LDAc and PBE potentials, respectively, which are in reasonable agreementwith the most recent experimental estimate of
573±23 meV/molecule [56]. Optimized cellparameters are a = 7.273 Å, b = 9.418 Å,

20

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
50

30

rev-vdW-DF2 (B86Rx+LDAc)
rev-vdW-DF2 (PBE)
ACFDT-RPA
Expt.

20
0

0
Eb (meV)Eb (meV/atom)

10

rev-vdW-DF2
PBE
ACFDT-RPA

-10
-20
-30

-50

-40
-50
-60
2.0

2.5

3.0

3.5

4.0

4.5

-100

5.0

2

Figure 4: Binding energy of graphite (Eb )as a function of interlayer distance (d) using rev-vdW-DF2 with different potentials,
along with the equilibrium energy and distancefrom ACFDT-RPA (Ref. [54]) and experiment(Ref. [55])and c = 6.723 Å, which are also in reasonable agreement with the experimental valuesof a = 7.36 Å, b = 9.37 Å, and c = 6.70 Å [57].
These results suggest that the effect of using the potentials generated using the PBE
functional is minor in vdW-DF calculations,at least with rev-vdW-DF2, although benchmark calculations on a larger set of materials
and molecules need to be done.
-470
-480

3

4

5

d (Å)

d (Å)

rev-vdW-DF2 (B86Rx+LDAc)
rev-vdW-DF2 (PBE)

Eb (meV/molecule)

-490
-500
-510

Figure 6: Binding energy of graphene onNi(111) (Eb ) as a function of distance fromthe surface (d) using rev-vdW-DF2 and PBE,
along with the ACFDT-RPA results from
Refs. [58, 59].
used but the result obtained using PBE potentials was found to be virtually identical. I also
confirmed that present result is consistent withthe one reported in Ref. [37] but with a different code. Figure. 6 shows the binding energy
as a function of graphene-Ni surface distanceby using PBE and rev-vdW-DF2 functionals(PBE potentials were used for the PBE calculations). It was found that the rev-vdWDF2 slightly overestimates the binding energies, but is able to reproduce the binding energy curve with two local minima obtained using the ACFDT-RPA [58, 59], which is suggested to be a result of competing chemicaland physical interactions. This demonstrates
the accuracy of rev-vdW-DF2 in describing theadsorption system.

-520
-530
-540
440

4
460

480

500
3
V (Å )

520

540

Figure 5: Binding energy of benzene crystal in
the P bca symmetry (Eb ) as a function of unit
cell volume (V ) obtained using rev-vdW-DF2
and different potentials.
Finally, as a critical assessment of rev-vdWDF2, I chose to calculate graphene adsorbedon Ni(111) in the atop-fcc configuration. In
this calculation, B86Rx+LDAc potentials wereSummary

I give a brief review of the theory of vdWDF, present recent theoretical advancements
and some applications based on a variant of
vdW-DF, rev-vdW-DF2. I discuss the properchoice of the XC functional to be used in
the pseudo- and PAW-potential generation for
vdW-DF calculations. It was found that theuse of PBE potentials in rev-vdW-DF2 calculations has a minor impact on the calculated
results, which may assure vdW-DF calculations in the literature. The accuracy of rev-

21

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

vdW-DF2 in describing the adsorption system
is also demonstrated, and it is anticipated thatits applicability is further broadened by introducing a fraction of the unscreened/screened
Fock exchange, i.e., nonempirical hybrid vdWDF [60, 61, 62].

Acknowledgements
The present work has been supported partly
by Grants in Aid for Scientific Researchon Innovative Areas ”Hydrogenomics” (No.
JP18H05519). Numerical calculations were inpart performed using the super computers in
Institute for Solid State Physics (ISSP), theUniversity of Tokyo.

References[1] P. Hohenberg and W. Kohn, Phys. Rev.
136, B864 (1964).
[2] W. Kohn and L. J. Sham, Phys. Rev. 140,
A1133 (1965).

[12] J. Hermann, R. A. DiStasio, and A.
Tkatchenko, Chem. Rev. 117, 4714
(2017).
[13] E. Caldeweyher and J. G. Brandenburg,
J. Phys. Condens. Matter 30, 213001(2018).
[14] H. Rydberg, et al., Phys. Rev. Lett. 91,
126402 (2003).
[15] M. Dion, et al., Phys. Rev. Lett. 92,
246401 (2004).
[16] M. Dion, et al., Phys. Rev. Lett. 95,
109902 (2005)[17] K. Lee, et al., Phys. Rev. B 82, 081101(2010).
[18] K. Berland, et al., Rep. Prog. Phys. 78,
066501 (2015).
[19] P. Hyldgaard, K. Berland and E.
Schröder, Phys. Rev. B 90, 075148 (2014).
[20] O. Gunnarsson and B. I. Lundqvist, Phys.
Rev. B 13, 4274 (1976).

[3] F. London, Z. Physik 63, 245 (1930).
[4] J. Hepburn, G. Scoles, and R. Penco,Chem. Phys. Lett. 36, 451 (1975).
[5] S. Grimme, J. Comput. Chem. 27, 1787
(2006).
[6] S. Grimme, J. Antony, S. Ehrlich, and
H. Krieg, J. Chem. Phys. 132, 154104(2010).
[7] A. D. Becke and E. R. Johnson, J. Chem.
Phys. 123, 154101 (2005).
[8] A. Tkatchenko and M. Scheffler, Phys.
Rev. Lett. 102, 073005 (2009).
[9] T. Sato and H. Nakai, J. Chem. Phys.
131, 224104 (2009).
[10] A. Tkatchenko, R. A. DiStasio, R. Car,and M. Scheffler, Phys. Rev. Lett. 108,236402 (2012).
[11] S. Grimme, A. Hansen, J. G. Brandenburg, and C. Bannwarth, Chem. Rev.
116, 5105 (2016).

[21] D. C. Langreth and J. P. Perdew, Phys.
Rev. B 15, 2884 (1977).
[22] D. C. Langreth and S. Vosko, Adv. Quantum Chem. 21, 175 (1990).
[23] T. Thonhauser, et al., Phys. Rev. B 76,
125112 (2007).
[24] K. Berland, D. Chakraborty, and T.
Thonhauser, Phys. Rev. B 99, 195418(2019).
[25] D. Chakraborty, K. Berland, and T.
Thonhauser, J. Chem. Theory Comput.
16, 5893 (2020).
[26] I. Hamada, vdwphi drsll:
A codeto calculate the kernel function of
the van der waals density functional,
https://github.com/ikuhamada/vdwphi drsll.
[27] J. P. Perdew, K. Burke, and M. Ernzerhof,
Phys. Rev. Lett. 77, 3865 (1996).
[28] Y. Zhang and W. Yang, Phys. Rev. Lett.
80, 890 (1998).

22

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[29] J. Klimeš, D. R. Bowler, and A.
Michaelides, J. Phys. Condens. Matter
22, 022201 (2009).
[30] V. R. Cooper, Phys. Rev. B 81, 161104(2010).
[31] J. P. Perdew and W. Yue, Phys. Rev. B
33, 8800 (1986).
[32] P. R. Antoniewicz and L. Kleinman, Phys.
Rev. B 31, 6779 (1985).
[33] É. D. Murray, K. Lee, and D. C. Langreth,J. Chem. Theory Comput. 5, 2754 (2009).
[34] J. Klimeš, D. R. Bowler, and A.
Michaelides, Phys. Rev. B 83, 195131(2011).
[35] K. Berland and P. Hyldgaard, Phys. Rev.
B 89, 035412 (2014).
[36] A. D. Becke, J. Chem. Phys. 85, 7184
(1986).
[37] I. Hamada, Phys. Rev. B 89, 121103(2014).
[38] A. D. Becke, Phys. Rev. A 38, 3098
(1988).

[46] R. Sabatini et al. J. Phys. Condens. Matter 24, 424209 (2012).
[47] P. Giannozzi, et al. J. Phys. Condens.
Matter 21, 395502 (2009).
[48] T. Thonhauser, et al., Phys. Rev. Lett.
115, 136402 (2015).
[49] C. M. Frostenson, et al., Electron. Struct.
4, 014001 (2022).
[50] A. Dal Corso, Comput. Mater. Sci. 95,
337 (2014).
[51] P. Jurečka, J. Sp̌oner, J. Cěrný, and P.
Hobza, Phys. Chem. Chem. Phys. 8, 1985
(2006).
[52] T. Takatani, et al., J. Chem. Phys. 132,
144104 (2010).
[53] L. Spanu, S. Sorella, and G. Galli, Phys.
Rev. Lett. 103, 196401 (2009).
[54] S. Lebégue, et al., Phys. Rev. Lett. 105,
196401 (2010).
[55] R. Zacharia, H. Ulbricht, and T. Hertel,
Phys. Rev. B 69, 155406 (2004).
[56] J. Yang, et al., Science 345, 640 (2014).

[39] J. P. Perdew, et al., Phys. Rev. Lett. 100,
136406 (2008).

[57] W. David, R. Ibberson, G. Jeffrey, and J.
Ruble, Physica B 180-181, 597 (1992).

[40] T. Björkman, J. Chem. Phys. 141, 074708(2014).

[58] F. Mittendorfer, et al., Phys. Rev. B 84,
201401 (2011).

[41] T. Rangel, et al., Phys. Rev. B 93, 115206(2016).

[59] T. Olsen and K. S. Thygesen, Phys. Rev.
B 87, 075111 (2013).

[42] F. Tran, et al., Phys. Rev. Materials 3,063602 (2019).

[60] Y. Jiao, E. Schröder, and P. Hyldgaard,
J. Chem. Phys. 148, 194115 (2018).

[43] S. Yanagisawa and I. Hamada in Theoretical Chemistry for Advanced Nanomaterials: Functional Analysis by Computation and Experiment, edited by T. Onishi(Springer Singapore, Singapore, 2020).

[61] V. Shukla, Y. Jiao, C. M. Frostenson, P.
Hyldgaard, J. Phys. Condens. Matter 34,
025902 (2022).
[62] V. Shukla et al., arXiv:2203.06682.

[44] P. E. Blöchl, Phys. Rev. B 50, 17953(1994).
[45] M. Callsen and I. Hamada, Phys. Rev. B
91, 195103 (2015).

23