van der Waals 力が支配する系：vdW-DF 汎関数による第一原理解析#

PHASE logo

希ガス固体、分子結晶、層状物質、そして二次元材料を積み重ねたファンデルワールス・ヘテロ構造では、系をまとめている力の主役が van der Waals 力（ロンドン分散力）です。分散力は離れた領域の電子密度ゆらぎが相関することで生じる非局所的な効果であり、局所的な密度と勾配だけを見る LDA / GGA では原理的に記述できません。実際、後述するように fcc Ar の凝集を担っているのは非局所相関項だけで、グラフェン/h-BN ヘテロ二層では GGA はほとんど結合を与えません。この欠落を、経験的なパラメータを足すのではなく汎関数そのものの中で補うのが vdW-DF 系の非局所相関汎関数です。第一原理計算ソフトウェア Advance/PHASE (ver.5.0) には vdW-DF / vdW-DF2 が自己無撞着に実装されており、さらに Román-Pérez–Soler の波数空間法によって、二重積分の計算コストが GGA と同じ $\mathrm{O}(N\log N)$ スケーリングに抑えられています。本解析では、希ガス固体 fcc Ar の状態方程式で 5 種類の vdW 系汎関数を文献値と突き合わせ、続いてグラフェン/h-BN ヘテロ二層の層間結合で GGA との違いと実用上の計算コストを検証しました。

Keywords: 第一原理計算 (DFT), van der Waals 力, 分散力, vdW-DF, vdW-DF2, 非局所相関汎関数, 状態方程式, 希ガス固体, 層状物質, グラフェン, 六方晶窒化ホウ素 (h-BN), ヘテロ構造, 積層依存性

vdW-DF 法について#

vdW-DF [1] では、交換相関エネルギーを

$E_{xc} = E_{x}^{\rm GGA} + E_{c}^{\rm LDA} + E_{c}^{\rm nl}$

と分割し、最後の非局所 (nl) 相関項

$\displaystyle E_{c}^{\rm nl} = \frac{1}{2}\iint n(\mathbf{r})\,\phi\!\left(q_0(\mathbf{r}),\,q_0(\mathbf{r}'),\,|\mathbf{r}-\mathbf{r}'|\right) n(\mathbf{r}')\, d\mathbf{r}\, d\mathbf{r}'$

が分散力を担います。ここで $n$ は電子密度、 $q_0$ は各点の密度とその勾配から決まる局所的な波数、 $\phi$ は離れた 2 点の密度ゆらぎの相関を与える普遍的なカーネルで、系にも元素にも依存しません。この二重積分を実空間でそのまま評価すると、電荷密度メッシュ点数 $N$ に対して $\mathrm{O}(N^2)$ 回のカーネル評価が必要になり、たとえば 60³ メッシュでは 4.7×10¹⁰ 回に達して実用になりません。Román-Pérez と Soler [2] は、 $q_0$ 依存性を固定 $q$ メッシュ上のスプライン基底 $p_{\alpha}$ で展開することで、上式が

$\displaystyle E_{c}^{\rm nl} = \frac{\Omega}{2}\sum_{\mathbf{G}}\sum_{\alpha\beta} \theta_{\alpha}^{*}(\mathbf{G})\, \phi_{\alpha\beta}(|\mathbf{G}|)\, \theta_{\beta}(\mathbf{G})$

と波数空間の単純な和に書き換えられることを示しました。 $\Omega$ はセル体積、 $\mathbf{G}$ は逆格子ベクトル、 $\alpha,\beta$ は $q$ メッシュの番号で、 $\theta_{\alpha}(\mathbf{G})$ は密度から作った $\theta_{\alpha}(\mathbf{r}) = n(\mathbf{r})\,p_{\alpha}(q_0(\mathbf{r}))$ を FFT したものです。動径平均されたカーネル $\phi_{\alpha\beta}(k)$ は電荷密度に依存しないので計算の初期段階で一度だけ表にしておけばよく、SCF の各ステップでは通常のGGA計算のようなコストでエネルギー・ポテンシャル・応力が得られます。表1はAdvance/PHASE で利用できる主な vdW 系汎関数を示します。

表1. Advance/PHASE で利用できる主な vdW 系汎関数

汎関数	交換	半局所相関
vdW-DF	revPBE	PW92 (LDA)
vdW-DF2	PW86R	PW92 (LDA)
optB88-vdW／optB86b-vdW／optPBE-vdW ／vdW-DF-cx／rev-vdW-DF2	libxc で選択（GUIで自動設定）	PW92 (LDA)

本解析では、2 つの系で目的に応じて汎関数を選びました。fcc Ar では 5 種類の vdW 系汎関数（内蔵の vdW-DF・vdW-DF2 と、libxc [3] 経由の optB88-vdW・optB86b-vdW・rev-vdW-DF2）で計算を行い、文献と比較します。グラフェン/h-BN では、分散力の有無そのものを見るために GGA の PBE と、vdW-DF の交換項と同じ revPBE で計算を行い、vdW-DF・vdW-DF2 と比較します。

計算モデルと計算条件#

図1. 計算モデル。(a) fcc Ar：立方セルで格子定数 $a$ を走査し、状態方程式を求めます。(b) グラフェン/h-BN ヘテロ二層の上面図。最安定の AB 積層（一方の C が B の直上、他方の C が六角形の中心）を示しています。(c) 側面図。層間距離 $d$ を走査し、 $c$ 軸長は 20 Å に固定して真空層を確保しています。

fcc Ar は、分散力だけで凝集している系として知られています。立方セルを用い、各汎関数の平衡体積のまわり ±12 % の体積範囲を 11 点走査して、Birch–Murnaghan 3 次式で平衡格子定数 $a_0$ を求めました。凝集エネルギーは、一辺 20 Å の立方箱に置いた孤立 Ar 原子を別に計算し、 $E_{\rm coh} = E_{\rm atom} - E_{\rm bulk}/4$ として求めています。

グラフェン/h-BN ヘテロ二層は、二次元材料のファンデルワールス・ヘテロ構造の代表例です。六方晶セルを用い、面内格子定数は両者の中間である $a$ = 2.50 Å に固定しました。層間距離 $d$ を 2.80〜5.00 Å で走査し、結合エネルギーは

$E_{\rm bind}(d) = E_{\rm bilayer}(d) - E_{\rm graphene} - E_{\text{h-BN}}$

という定義式から計算されます。

表2. 計算条件

項目	fcc Ar	グラフェン/h-BN
擬ポテンシャル	PAW
交換相関汎関数	vdW-DF / vdW-DF2 / optB88-vdW / optB86b-vdW / rev-vdW-DF2	PBE / revPBE / vdW-DF / vdW-DF2
カットオフ（波動関数 / 電荷密度）	55 / 500 Ry	36 / 324 Ry
k点サンプリング	6×6×6	12×12×1
vdW-DF カーネル表	20 点 $q$ メッシュ、$q_{\rm cut}$ = 5（生成は 1 回だけ）

計算結果と考察#

1. fcc Ar の状態方程式#

図2. (a) fcc Ar の状態方程式。記号が計算点、実線が Birch–Murnaghan 3 次式のフィット、★が平衡点です。破線は実験の平衡格子定数 5.23 Å。(b) 平衡点での凝集エネルギーを、非局所相関項 $\Delta E_c^{\rm nl}$ とそれ以外（半局所交換相関＋静電＋運動エネルギー）に分解したもの。

図2(a) の 5 本の曲線は、そのまま vdW 系汎関数の「性質の違い」を表しています。平衡格子定数は optB88-vdW の 5.224 Å から vdW-DF の 5.492 Å まで 0.27 Å の幅で分かれ、凝集エネルギーも 87〜160 meV/atom と 2 倍近く違います。そして両者の順位は一致しません。格子定数が実験にいちばん近い optB88-vdW・optB86b-vdW は凝集エネルギーを 1.4〜1.6 倍に過大評価し、逆に凝集エネルギーをよく評価できる rev-vdW-DF2 は、格子定数を 2 % 程度大きく見積もります（表3）。

表3. fcc Ar の平衡格子定数・凝集エネルギー（ $a_0$ の小さい順）

汎関数	$a_0$ (Å)	$E_{\rm coh}$ (meV/atom)	文献の $a_0$ (Å) PAW [4] / 全電子 [5]	文献の $E_{\rm coh}$ (meV/atom) PAW [4] / 全電子 [5]
optB88-vdW	5.224	120.0	5.22 / 5.23	143 / 138
vdW-DF2	5.270	132.1	5.26 / 5.28	142 / 124
optB86b-vdW	5.303	138.9	5.30 / 5.32	131 / 137
rev-vdW-DF2	5.351	87.2	5.35 / 5.37	90 / 82
vdW-DF	5.492	159.7	— / 5.50	— / 156
比較の基準
CCSD(T) [4]	5.284	82.8	—
実験 [4]	5.23	88.9	—

平衡格子定数は、文献 [4] の PAW 計算結果と ±0.01 Å の差で一致しており、全電子計算 [5] の結果とも整合しています。序列（optB88-vdW < vdW-DF2 < optB86b-vdW < rev-vdW-DF2 ≪ vdW-DF）も再現しました。凝集エネルギーはいずれも文献の範囲内です。

図2(b) は、この凝集エネルギーの寄与を分解して示したものです。結合を作っているのは非局所相関項 $E_c^{\rm nl}$ だけです。その寄与は凝集エネルギーの 1.4〜2.5 倍にのぼり、半局所交換相関と静電・運動エネルギーの寄与はどの汎関数でも反発です。つまり半局所部分だけを取り出せば fcc Ar は凝集しないことになります。

実用上の結論は「用途によって汎関数を使い分ける」ということになります。吸着エネルギーや分子間結合エネルギーを定量的に出したいなら rev-vdW-DF2、格子定数や吸着距離を合わせたいなら optB86b-vdW か optB88-vdW、というのが文献の共通見解で、本計算はそれを再現しています。

2. グラフェン/h-BN ヘテロ二層の層間結合#

グラフェン/h-BN の層間結合曲線
図3. (a) グラフェン/h-BN ヘテロ二層の層間結合曲線（★は平衡点）。灰色の帯はグラファイト（3.35 Å）および h-BN（3.33 Å）の実験的な層間距離です。(b) vdW-DF の結合エネルギーを、非局所相関項 $\Delta E_c^{\rm nl}$ とそれ以外（半局所交換相関＋静電＋運動エネルギー）に分解したもの。

表4. グラフェン/h-BN ヘテロ二層の平衡層間距離と結合エネルギー（AB 積層）

汎関数	$d_0$ (Å)	$E_{\rm bind}$ (meV / C 原子)
revPBE	—	極小なし
PBE	4.28	−1.7
vdW-DF	3.61	−49.7
vdW-DF2	3.53	−49.2

ここでの GGA の失敗はさらに劇的です。revPBE は走査範囲のほぼ全域で斥力で、極小を持ちません。PBE でも結合エネルギーは −1.7 meV/C 原子と実質ゼロで、しかも極小は 4.28 Å と実験的な層間距離（グラファイト 3.35 Å、h-BN 3.33 Å）から 0.9 Å 以上も外れます。GGA だけで層状物質の $c$ 軸を構造最適化すると層がどこまでも離れていく、という現象の正体がこれです。

vdW-DF では層間距離 3.61 Å、結合エネルギー −49.7 meV/C 原子という、実験値に近い明瞭な結合が得られました。vdW-DF2 は 3.53 Å とやや短く、これは vdW-DF が層間距離をわずかに長く見積もる傾向（PW86R 交換を使う vdW-DF2 で改善される）という、文献で知られた振る舞いと一致しています。

図3(b) の分解も fcc Ar と同じ構図です。平衡点付近（ $d$ = 3.60 Å）で非局所相関 $\Delta E_c^{\rm nl}$ が −83 meV/C 原子の引力を与え、半局所交換相関と静電・運動エネルギーの寄与は +33 meV/C 原子の斥力、差し引き −50 meV/C 原子が結合エネルギーになります。結合エネルギー自体より大きな引力項と斥力項の差し引きで決まっているため、半局所汎関数だけでは層状物質を扱えず、vdW-DF 系でも交換汎関数の選択が結果を大きく左右することになります。

積層の依存性#

次に、vdW-DFを用いて、層間の相対的な平行移動（積層構造）に対するエネルギーの変化を評価しました。

積層依存性
図4. vdW-DF による積層構造の相対エネルギー（ $d$ = 3.40 Å に固定）。C が B の直上に来る AB 積層が最安定です。

図4はvdW-DF による積層構造の相対エネルギーを示します。3 種類の高対称積層を $d$ = 3.40 Å に固定して比較したところ、一方の C 原子が B 原子の直上、他方の C 原子が六角形の中心に位置する AB 積層が最安定となり、C が N の直上に来る配置より 7.6 meV/C 原子、C が B と N の両方に重なる AA 積層より 9.4 meV/C 原子だけ安定でした。この序列は報告されている第一原理計算の結果 [6,7] と一致します。数 meV/原子という微小なエネルギー差が、実験で観測されるモアレ構造や整合／非整合の議論を左右します。この程度の差を安定に取り出せることは、実用上の大きな意味を持ちます。

3. 計算コスト#

vdW-DF の計算コスト
図5. (a) 非局所相関 $E_c^{\rm nl}$ 1 回あたりの評価時間と電荷密度 FFT 格子点数 $N$ の関係（両対数）。実測値は $N\log N$ の線に沿っており、実空間の直接二重和が持つ $N^2$ 依存性から明確に外れています。(b) 同じ系を revPBE (GGA) と vdW-DF で計算したときの SCF 1 反復あたりの実時間。

図5(a)で示すように、グラフェン/h-BNで cutoff_cd を 196〜576 Ry まで変えて FFT 格子点数 $N$ を 10.4 万点から 48.0 万点まで（4.6 倍）増やしたとき、 $E_c^{\rm nl}$ 1 回の評価時間は 0.140 秒から 0.712 秒へ 5.1 倍になりました。6 点の両対数フィットで得られる指数は $t \propto N^{1.10}$ （決定係数 0.97）で、この $N$ の範囲で $N\log N$ が持つ実効指数 1.08 と一致しています。実空間積分の $\mathrm{O}(N^2)$ スケーリングから劇的に改善されており、波数空間実装が意図どおりに働いていることが確認できます。

SCF 全体で見ると、vdW-DF は revPBE の 1.1〜1.5 倍の実時間で完了しました（図5(b)）。カーネル表の生成は 8.4 秒で、これは計算全体を通じて 1 回だけです。分散力を正しく扱うための追加コストが数十パーセントに収まることは、層状物質や吸着系の構造最適化をルーチンとして vdW-DF で行えることを意味します。

まとめ#

第一原理計算ソフトウェア Advance/PHASE の vdW-DF系汎関数（Román-Pérez–Soler の波数空間実装）を用いて、van der Waals 力が支配する 2 つの系を解析しました。fcc Ar では 5 種類の vdW 系汎関数の平衡格子定数が文献値と非常によく一致し、凝集エネルギーの分解から、結合を作っているのは非局所相関項 $E_c^{\rm nl}$ だけで半局所部分はどの汎関数でも反発であることが確認できました。グラフェン/h-BN ヘテロ二層では、revPBE が極小を持たず PBE も実質非結合であるのに対し、vdW-DF は層間距離 3.61 Å・結合エネルギー −49.7 meV/C 原子という実験に近い明瞭な結合を与え、AB 積層の安定性も数 meV/原子の分解能で再現しました。しかもこの精度を、SCF 全体で GGA の 1.1〜1.5 倍という計算コストで得ています。二次元材料のヘテロ積層、有機分子の表面吸着、層状電極材料へのイオン挿入、MOF や分子結晶など、構造も反応性も分散力が決める系に対して、vdW-DF は経験的なパラメータに頼らずに定量的な解析を行う実用的な手段となります。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献#

M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, "Van der Waals Density Functional for General Geometries", Phys. Rev. Lett. 92, 246401 (2004).
G. Román-Pérez and J. M. Soler, "Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes", Phys. Rev. Lett. 103, 096102 (2009).
S. Lehtola, C. Steigemann, M. J. T. Oliveira, and M. A. L. Marques, "Recent developments in Libxc - A comprehensive library of functionals for density functional theory", Software X 7, 1 (2018).
M. Callsen and I. Hamada, "Assessing the accuracy of the van der Waals density functionals for rare-gas and small molecular systems", Phys. Rev. B 91, 195103 (2015).
F. Tran, L. Kalantari, B. Traoré, X. Rocquefelte, and P. Blaha, "Nonlocal van der Waals functionals for solids: Choosing an appropriate one", Phys. Rev. Materials 3, 063602 (2019).
G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, "Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations", Phys. Rev. B 76, 073103 (2007).
B. Sachs, T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein, "Adhesion and electronic structure of graphene on hexagonal boron nitride substrates", Phys. Rev. B 84, 195414 (2011).