平面波基底の第一原理計算法(第四回)
香

正

憲

NG の列ベクトル(式( 43 ))で， Hc 計算は N 2G オーダーの演(第三回からの続き)算量である．式(52), (69)から，Hc 計算のうち運動エネル

高速フーリエ変換の活用(その 2)

山ギー項や非局所擬ポテンシャル項の部分は演算量を NG オー

高速フーリエ変換(fast Fourier transformation; FFT)の活

ダーに減らせるが，式( 57 )の局所ポテンシャルの Vlocal × c

用として，固有ベクトルを用いた r(r… )や Hci の演算が独特の演算は G• 毎に∑Vlocal ( G• － G• ′) C n•k＋G• ′で， N 2G オーダーであの工夫で飛躍的に効率化できる(44)(45) ．電子密度分布 r ( r… )の二重ループの演は，式( 72)の固有ベクトル成分の G• ， G• ′る．

算量が N 2G オーダーで膨大である(NG は平面波基底数，第 5

Vlocal(r… )c •kn(r… )＝Vlocal(r… )U •kn(r… ) exp [i •k･r… ]〉＝∑Vlocal(G• ) exp [iG• ･r… ]∑C n•k＋G• ′| •k＋G• ′occ

 節参照)．一方，式(19)から r(r… )＝∑∑w •ki n|U •ki n(r… )|2 で
章

G• ′

ここで，式(17)，(43)，(55)から

G•

n •ki

ある．格子周期関数 U •ki n ( r… )のフーリエ変換 U •ki n ( G• )は，式
1/2
( 43 )の議論から固有ベクトルを用いて{ U •ki n ( G• )}＝{ Q－cG•

1/2 exp [(iG• ＋G• ′)･r… ]＝∑∑Vlocal(G• )C n•k＋G• ′Q－
c
G• G• ′

数の規格化を単位胞 Qc で行うとする)．固有ベクトルからFFT の G• m 点で{U •ki n(G• m)}を組み立て，逆 FFT で単位胞内×exp [i •k･r… ]1/2 exp [iG• ･r… ]
)C n•k＋G• ′Q－＝∑∑Vlocal(G• －G• ′c

実空間メッシュ点{r… l}の{U •ki n(r… l)}が一度に与えられる．これG• G• ′

を|U •ki n(r… l)|2 で式(19)に代入すればよい．式(72)の N 2G オー

章の図 9 から M  16NG として NG が 1000 ， 10000 の時， M
log2 M/N 2G の比は0.22，0.028で，大きく低減化できる．
 節の対称性の議論から S •ki 点の |US •ki n(r… )|2
一方，第 3 章

は •ki 点の|U •ki n(r… )|2 の分布を実空間の対称操作{S|t• S}で移したもので，また|U－ •ki n(r… )|2＝|U •ki n(r… )|2 である．従って，ブリルアンゾーンの既約領域内の •ki 点のみで{U •ki n(r… l)}を求め，G• ′×exp [iG• ′･r… ] exp [i •k･r… ] 節の末尾のように基底関数，固有関C n•ki＋G• }である(第 5 章ダーの演算を回避し， FFT の演算量 M log2 M オーダー( Mは FFT の G• m 点総数)で r ( r… )が得られることになる．第 6

G• ′

1 /2
＝∑Vlocal(G• ) exp [iG• ･r… ]∑C n•k＋G• ′Q－
c×exp [i •k･r… ](92)

1/2 を用いている(第 5 章と表現できる．Q－1/2 の代わりに Q－c 節の議論参照)．一行目と最後の行の比較から格子周期関
数 Vlocal ( r… ) U •kn ( r… ) Q1c /2 のフーリエ変換が G• 毎の∑Vlocal ( G• －G• ′)C n•k＋G• ′と言える．固有ベクトル{C n•k＋G• }が与えられると，G• ′1/ 2C nFFT の G• m 点での{ U •kn ( G• m )}＝{ Q－•k＋G• m }から逆 FFT でc
実空間の{U •kn(r… l)}が求まる．同様に{Vlocal(G• m)}から逆 FFT
で実空間の{Vlocal(r… l)}が求まる．両者から実空間メッシュ点∑∑w •k n | U •k n ( r… l )|2 を計算した後，この和の全体を{ S | t• S }でnでの積{Vlocal(r… l)U •kn(r… l)Q1c /2}のデータが決まる．これを FFT
すれば，全 G• m の∑Vlocal(G• m－G• ′)C n•k＋G• ′が一度に求まること

対称操作した分布の重ね合わせで最終的な{ r ( r… l )}が求まになる．これが目的とする量で，convolution 積分の FFT にる．用いる対称操作は，既約領域をブリルアンゾーンの半分よる効率化である． N 2G オーダーの演算が M log2 M オーダに広げるものだけで良い．

ーになる．

occ

•ki

i

i

次に，第 7 章の大規模系の基底状態計算の高速化技法では，波動関数 c の最適化のために Hc (行列×ベクトル)の演算を繰り返す．式( 71 )の H は NG × NG 行列， c はサイズ

G• ′




原子に働く力

第一原理計算では，原子に働く力は Hellmann Feynman

国立研究開発法人 産業技術総合研究所 エネルギー・環境領域 電池技術研究部門名誉リサーチャー(〒5638577 池田市緑が丘 18
31 産総研関西センター)
Lecture Notes on FirstPrinciples Methods Using a PlaneWave Basis Set (Part 4); Masanori Kohyama(Research Institute of Electrochemical Energy, Department of Energy and Environment, National Institute of Advanced Industrial Science and Technology, Ikeda,
Osaka)
Keywords: HellmannFeynman force, firstprinciples molecular dynamics, CarParrinello method, GramSchmidt orthogonalization, conjugategradient method, subspace diagonalization, block Davidson method, RMMDIIS method, charge mixing
2022年 5 月11日受理[doi:10.2320/materia.61.878]




講義ノート

力(46) と呼ばれ，全エネルギー Etot の原子の位置ベクトルに×exp [－iG• ･r… aa ′]よる微分に負符号をつけたものである．ここで，Etot を原子
位置{R• a}と占有された波動関数のセット{ci}の汎関数とみな×∑R•

せば，式( 1 )～( 3 )から

f (－2m : )c (r… )d…r)er(r… )r(r… ′＋ V (r… )r(r… )d…r＋
f
f|r… －r… ′| dr… d…r ′＋E [{c }]＋E
2focc

Etot[{R• a}, {ci}]＝2∑ c
i (r… )

2

2

i

＋

i

xc

i

I－ I

(93)
occ

r(r… )＝2∑ci(r… )ci(r… )

(94)

i

1
R• ＋r… aa ′erfc (|R• ＋r… aa ′|g)|R• ＋r… aa ′| |R• ＋r… aa ′|2{2g
1
exp [－|R• ＋r… aa ′|2g 2]
p |R• ＋r… aa ′|

}

(96)

r… aa ′＝ t• a′－ t• a である．導出は， r… a ′a の場合や－ R• ，－ G• の和に

2I

exp [－|G• |2/4g2]－ e 2 Za ∑ Za ′≠a
|G• |2/4g2
a′

留意し，derfc(x )/dx＝－2/ p e－x 2 を用いている．Ewald 法と同様，パラメータ g を適当に選べば比較的短範囲の R• ，G•の和で収束する．
非局所項も含めて F• 1a を考える．これは式( 74 )，( 75 )のとなる( occ は占有状態の和， 2 はスピンの和)． VI ( r… )は，

EL ， ENL 内の t• a への直接依存項の微分で，式( 60 )，( 69 )，原子からのポテンシャル(擬ポテンシャル)で，簡単のため局(70)から

所項のみとしている．式( 94 )を踏まえながら，具体的に式( 93 )の原子の位置ベクトルによる微分を考える．原子位置F• 1a＝－に直接依存する項は VI(r… )と EI－I で，さらに ci(r… ), ci(r… ),

f
∑
f [

n •ki

&ci(r… )
&R• a

occ
i

－

G•

2 2
: ＋VI(r… )＋VH(r… )＋mxc(r… )2m]×ci(r… )d…r

f [

2 2
ci(r… ) －: ＋VI(r… )＋VH(r… )＋mxc(r… )
2m

occ

－2∑
i

]

&ci(r… )
×
d…r
&R• a

f{

－2∑ Eii&c &ci(r… )i (r… )
ci(r… )＋ci(r… )
&R• a
&R• a

(97)

となる．一項目が局所擬ポテンシャル，二項目が非局所擬ポ
テ ン シ ャ ル (分 離 形 )に よ る 寄 与 ． 自 己 無 撞 着 な { C n•ki＋G• } ,r(G• )をそのまま入れれば計算できる．
式(97)の右辺第二項の •ki 点の和は，実質上，既約領域(第 節)の •ki 点に絞ることができる．ただし，既約領域の3 章

•ki 点で演算を実行して求めた和(ベクトル)を，対称操作の行証明される．用いる点対称操作 S は既約領域をブリルアン

} d…r

f

&VI(r… )&EI－I＝－＝F• 1a＋F• 2ar(r… )d…r－
&R• a
&R• a

ゾーンの半分に広げるもののみで良い．
原子に働く力は，Etot の原子座標への直接依存項の微分の(95)r(r… )の原子位置依存性は式(94)より ci(r… )，ci(r… )の依存性で扱っている．VI(r… )と EI－I 以外の寄与は，式(10)のように

dEtot / dci， dEtot / dci から二行目～五行目が出てくる(スピみ で ， ci ( r… ) を 通 じ た 依 存 性 は 消 え る ． こ の こ と を 特 に
Hellmann Feynman の定理と呼ぶ．応力(ストレス)テンソ

ルの計算も同様の考え方で行われる．応力はセル当たりの
Etot の歪テンソル ea b による微分でン自由度 2 に注意)．ci は KohnSham 方程式の固有関数で
Hci ＝ Ei ci ， ciも ciH ＝ Ei ciを 満 た し， 下 か ら 二 行 目 の• a ci(r… )＋ c
• a} d…r が出てくる．
{ &c Ei fi ( r… )/ &Ri ( r… ) &ci ( r… )/ &R

規格化条件 fci(r… )ci(r… )d…r＝1 から，この項は最終行で消える．これが消えるので，r(r… )や ci(r… )の原子位置依存性は原子に働く力には入らないのである．例えば式(93)の f
VI(r… )

r(r… )d…r について，VI(r… )の部分の微分だけでよい．
最終 的に， 原子 に働く 力は ，原子 からの ポテ ンシ ャルVI(r… )を原子位置で微分した項 F• 1a，正イオン間静電相互作用s a b＝

で，各単位胞の同じ原子の一斉の変位についての微分である．
F• 2a の方は， Ewald 法の式( 87 )の r… aa′依存項( t• a 依存項)の微分で次式になる．
&EEwalde2 p＝－ 2 Za ∑ Za′∑ iG•
g Qc a ′≠a G• ≠0&t• a

ま て り あ 第61巻 第12号(2022)
Materia Japan(98)

である．Etot の各項の歪への直接依存成分の微分だけを考え固有ベクトル成分 C n•k＋G• の直接依存性はなく，exp [iG• ･t• a]等• )なの歪依存性も G• とt• a で逆なので消える． r(G• )や V PS
local ( G
どの歪依存性は Qc や| G• |を使った表式まで書き下して導出される．詳細は文献(4)，(21)等を参照．

.

に非局所項がある場合は表式が少し変わるが，意味は同じである．具体的には，周期セルの内部座標 t• a についての微分

1 &Etot
Qc &ea b 節)で は ，る ． EI－I 項以 外 の各 項 の逆 空間 表 現(第 6 章22
Qc，|G• |，|G• | ，|k• ＋G• | などが歪への直接依存項である．

EI－I を原子位置で微分した項 F• 2a の二項のみである． VI ( r… )F• 2a＝－

G•〈&VI(r… )&EI－Ir(r… )d…r－
&R• a
&R• a

occ

m＝－l

列 S で回転して平均したベクトルが力への寄与になる．このことは，式(51)や式(70)の Aa, lm( •ki＋G• )内の Ylm(ki＋G )や jl(| •ki＋ G• |r )の性質( •ki を S •ki に替えたものとの関係)から

f

＝－

l

iG• C n•ki＋G• Aa, lm( •ki＋G• ))＋c.c.}

×(∑

&Etot
&VI(r… )&EI－I＝－r(r… )d…r－
&R• a
&R• a
&R• a

－2

＋l

occ

n
1
• ))
－∑∑w •ki n∑C－a, l ∑ {(∑C •ki＋G• Aa, lm( •ki＋G

r(r… )の原子位置依存性を考えると次式のようになる．
F• a＝－&(EL＋ENL)＝ ∑ iG• exp [－iG• ･t• a]V alocal(G• )r(－G• )
&t• a
G• ≠0




大規模電子構造計算の高速化技法

CarParrinello 法と直接最小化法

微小時間刻みの原子の動きのたび毎に，その時の原子配列に対するスーパーセルの電子構造計算を行い，断熱近似(第 節)のもとでの電子系の最安定状態(基底状態)を求
2 章

め，原子に働く力を正確に与えて原子の運動を扱う計算を第

一 原 理 分 子 動 力 学 法 ( first principles molecular dynamics;{C n•k＋G• (t )}(式(43))の運動方程式として，古典分子動力学法

FPMD)と呼ぶ．これは，平面波基底の第一原理擬ポテンシ

での Verlet アルゴリズム(48)に類似した手続きで，微小時間〈 •k＋ G• |H | •k ＋
刻み Dt で時間発展させる． Hci の演算は，∑

ャル法(NCPP 法，USPP 法，PAW 法)において，大規模スーパーセルの基底状態計算を飛躍的に高速化する計算技術に基づくものであり， Car Parrinello ( CP )法(47) をきっかけに急速に発展した．構造緩和で安定構造を求める場合や大規模系の電子構造を求めるだけの場合も，この技術が使われる．
平面波基底の第一原理計算では，平面波基底数 NG は原子
当たり 100 個のレベルで， 100 原子を超えるスーパーセルで 節で論じたように NG はは数万のオーダーである(第 5 章単位胞体積原子数に比例)．式(71)のハミルトニアンは数G• ′ 節の FFT による効率化法を用い〉C n•k＋G• ′( t )で，第 6 章G• ′る． li は期待値〈 ci | H | ci 〉である(この部分は別の扱い方もある(48))．

Dt ステップ毎に全ての状態(ベクトル){C n•k＋G• (t )}を変化させ た 後 ， 各 状 態 の 規 格 直 交 化 を Gram Schmidt の 直 交 化法(49) で行う．これは，同じ •k の状態間で下の準位から順に次式の処理を行う．
|ci〉cj′，ci′＝ci″/|ci″|ci″＝ci－∑〈cj′

万×数万のオーダーとなる．通常の(数値計算ライブラリに

(101)

j＜ i

ある)エルミート行列対角化の固有値・固有状態計算は，NG

同じ •k の状態の間で，ベクトル{C i•k＋G• (t )}から下の準位のベ個全ての固有値を求めるアルゴリズムで，巨大行列への適用

クトル{C j•k＋G• (t )}との重なりを差し引く処理である．プライは困難である．実際は，全ての固有状態は必要でなく，半導ムのついた波動関数 c j′
は，既に下の準位の波動関数と規格体や絶縁体なら占有準位のみ，金属的な系でも占有準位と少直交化済みである．異なる •k の状態間は自動的に直交していし上までの固有状態だけで良い．そこで，行列対角化でなく，る．さらに Dt ステップ毎にハミルトニアン内の VH， mxc を

input し た 初 期 ベ ク ト ル を 繰 り 返 し 法 ( 反 復 法 ， iterative

新しい{ C n•k＋G• ( t )}による r ( r… )を用いて更新する．再び式

method)で徐々に固有状態に収束させる方法が適すると考え(100)の Hci 計算を行い，{C n•k＋G• (t )}を変化させる．

ci ( t )の運動({ C n•k＋G• ( t )}の運動)は(仮想的)運動エネルギられる(次節で論じる)．
一方，KohnSham 方程式のハミルトニアンを扱うのでなー ∑1/2 me|c· i|2 を徐々に減らし，全状態で c· i＝ 0 になるよ

く，前段階の式( 1 )の全エネルギー Etot を占有波動関数のうアニールすれば，固定した原子配列についての基底状態のセット{ci}の汎関数 Etot[{ci}]として直接に最小化すること{ ci }({ C n•k＋G• })に行きつく．図のように{ ci }の配位座標空が考えられる( ci は c •kn と同じ意味)．適当に作った{ ci }を間におけるEtot[{ci}]の最小点を効率的に見つける問題(49)とinput し，規格直交条件の下で Etot [{ ci }]を下げるよう{ ci }考えれば，CP 法は，図10(a)にように{ci(t )}の運動とアニを徐々に変化させ，Etot を最小化する{ci}を求めれば，密度ールでそれを実現する．

汎関数理論の基底状態が(KohnSham 方程式を介さず)求まるはずである．

i

当初の CP 法は，原子系の運動と電子系(波動関数)の運動
を同じ力学系として(同じラグランジアンの式に表して)解く

この方策を最初に実行したのが CP 法である(47) ．第 2 章

手法として提案されたが，固定した原子配列に対して，

 節，式( 8 )～( 10 )の Kohn Sham 方程式の導出過程を思Etot[{ci}]を規格直交条件のもと最小化する{ci}のセット(断い出そう．式( 9 )の Lagrange 未定係数法から，状態間の規格直交化は別過程で付加するとして－dQtot[{ci}]/dci＝－Hci＋li ci＝－(H－li)ci＝zi(99)

が全占有状態でゼロになればよい．式(10)を用いている．li＝〈ci|H|ci〉として， ci が固有状態でなければ Hci≠lici，zi  0 である．この zi (平面波基底でサイズ NG のベクトル)は規格直交条件の下で Etot を下げる方向の ci についての勾配( gradient )で， ci をその方向に変化させれば( zi の成分を少し混ぜれば) Etot が下がるので， ci に働く「力」と見なせる．そこで，適当に作った初期の波動関数セット{ ci }を式(99)の勾配(力)に従って徐々に変化させる．変化過程の{ci}からの r ( r… )を Etot や H に入れることで Etot や H ， li ， zi も同時に変化する．zi＝0 になるまで繰り返せば，Etot[{ci}]を
規格直交条件下で最小化する{ci}に行き着く．
CP 法では，これを「仮想的力学系の運動」として実行する．具体的には，仮想時間依存性を入れた波動関数 ci(t )について，式( 99 )の勾配(力)を用いた二階の運動方程式としてme cä i＝－(H－li)ci＝zi(100)が組み立てられる(47)．cä i は時間での二階微分で加速度の意
味，me は電子質量(これは厳密な値でなくパラメータとして扱う)．この式を平面波基底の展開係数(ベクトル)




図10

波動関数の配位座標空間で，全エネルギー Etot[{ci}]を最小にする占有波動関数のセット{ ci}(基底状態の波動関数)を，繰り返し法による波動関数の逐次更新で効率的に求める方法(直接最小化法)の概念図．曲線または
折れ曲がった実線が{ci}の変化の軌跡を示す．(a) CarParrinello 法，(b) 最急降下法，(c) 共役勾配法．

講義ノート

熱近似での基底状態)を効率的に求める手法と見なせる．そ＋1[{ ci }]に cmを代入した最小化条件から決まる(式( 88)のi

うすると，一階の運動方程式m＋1の寄与で置き換え， u での微分をEtot で cm
i の寄与を ci

me c· i＝－(H－li)ci(102)

計算)．この共役勾配法のステップを各状態で数回ずつ行から， ci を勾配方向に直にステップで変化させる最急降下

い ， 全 •k 点 の 全 状 態 で 巡 回 す る 過 程 を 繰 り 返 せ ば ，(steepest descent; SD)法(図10(b))も考えられる．単純な系

Etot[{ci}]を最小化する{ci}のセットが，CP 法や最急降下法

なら，より簡単に Etot[{ci}]を最小化する{ci}のセットに到よりも効率的に決まる可能性がある(図(a))．

達できる可能性がある．
一方，{ci}の配位座標空間における Etot[{ci}]の最小化問

これら式( 100 )，式( 102 )～( 105 )の各方法の演算は，全て，サイズ NG の列ベクトルと NG × NG 行列〈 •k ＋ G• | H | •k ＋

題は，多変数関数の最小化の逐次最適化( sequential optimi-〉を用いた演算である．行列対角化によらず， input したG• ′zation)法として，共役勾配(conjugate gradient; CG)法(49)の初期ベクトル{C n•k＋G• }の逐次更新(繰り返し法(反復法))によ適用が考えられる(図 10 ( c ))．適当な初期ベクトル c0i からる最適化である．初期ベクトルは乱数でランダムに与えたり，始めて m ステップ目の波動関数を cm
i とし，ステップを進NG を小さくして試行的にハミルトニアンの対角化から作っめることで，全状態について式( 99 )の勾配 zi がゼロになる{ci}を求める．m ステップ目の c mi についての勾配は式(99)たりする．特徴的なことは，{ C n•k＋G• }の更新のたび毎に式(101)による直交化処理を行い，行列〈 •k＋ G• |H | •k＋ G• ′
〉内の

m
m
m
m
から z m〈 cm)，cm
i ＝－( H － li ) ci だが( li ＝i | H | ci 〉
i に与え

波動関数依存項(r(r… )依存項)VH，mxc の更新も行う．その時mる変化の方向は，勾配 zm
i ではなく共役勾配ベクトル qi

点の{C n•k＋G• }から直に r(r… ), VH, mxc を与えることで，SCF を－1
で，前回ステップの勾配ベクトル z m，共役勾配ベクトルi

保ったまま，互いに直交した{ci}のセットが収束に向かう．

－1
も用いて，
qm
i

これは，KohnSham 方程式の固有値問題を行列対角化法で

m
m－1
qm
i ＝ z i ＋ hm q i

－1
－1
m
/〈z m
|z m
〉
hm＝〈z mi |z i 〉
ii(103)(m＞1)，hm＝0

(m＝1)
(104)

m
m
で作られる(逐次補正)．qm
i と下の準位 {cj }及び ci との直交化処理を式( 101 )で行った後，次のステップの波動関数＋1mmを次式のように cm
cm
i と qi の混合で与える(ci と下の準i

位 cm
j との直交化処理は既に行っている)．
＋1
m
＝acm
cm
i ＋bqii解き，r(r… )の更新を SCF ループで分けて行う従来法(第 3 章 節)と大きく異なる．その意味で，本節の諸法は Etot の直接最小化法(direct minimization)(49)と呼ばれる．


大規模行列固有状態計算の高速化技法

前節の一連の手法は， Si など半導体で威力を発揮し，多くの成果が得られた(図11(a))．しかし，金属の取り扱いで(105)

はうまくいかない(50)．前節の手法では，占有状態の数をは混合比の a ， b ( a ＝ cos u ， b ＝ sin u ， u で指定)は直接に Etot

じめから固定して{ci}を扱うが，金属的な系の場合，収束過程で占有状態と非占有状態が入れ替わったり， •k 点によって占有状態の数が変わったりする．そうした変化が{ci}と r の逐次の更新時に取り入れられると，式( 99 )の Hci がそのス
テップで急変し，{ci}の最適化の安定な運行が阻害されるためである(50)．
そこで，固有状態計算と SCF ループでの r の更新を分け 節)に立ち返る． input のて行う従来のアプローチ(第 3 章r ( r… )に固定したハミルトニアン行列 Hfix (式( 71 )の VH ， mxc
を固定)の固有状態計算(第 3 章の図 4 の二重四角の部分)を，行列の対角化でなく，input の波動関数(ベクトル)の繰り返し法(反復法)による最適化で行い，次に SCF ループでの r，Hfix の更新を行い，交互に繰り返すわけである．初期
図11(a)直接最小化法での全エネルギーの収束性の比較(50)．
SiC( 110 )表面スラブ( 10 原子＋真空層)において，適当な input 波動関数から始めて，更新 step での全エネル
ギーの低下・収束の様子(縦軸の単位は Ry13.6 eV)．
SD ， CG は 最 急 降 下 法 ， 共 役 勾 配 法 の 結 果 ， ALL 
BAND CG は，波動関数の更新を一つずつ逐次に行うのではなく，全状態を一斉に更新する共役勾配法の結果を示す．
( b )金属的な系も扱える高速計算技法の比較． BKL 法( Bylander ら の 共 役 勾 配 法 ) ， BD 法 ( block Davidson法)を Kerker の chargemixing 法(式(112))と組み合わ
せて Al/SiC 界面スーパーセル(15原子＋真空層)に適用した結果(50)．SCF ループでの input と output の電子密度分布の差の平均(メッシュ点での差の二乗平均の平方根)の変化を示す．電子密度分布の収束は，全エネルギーの収束と同様である． charge mixing のパラメータ
G 20 をうまく選べば BKL 法が効率的であった(図中のG 20の値の単位は Ry，A は1.0に固定)．

ま て り あ 第61巻 第12号(2022)
Materia Japanの input の r(r… )は，例えば自由原子の価電子密度分布の重ね合わせ等で設定する．
Bylander らは，式(103)～(105)の共役勾配法を Hfix の期

待値最小化の見地で適用し， SCF ループと組み合わせる手法を提案した(51)．勾配ベクトルは式(99)と同様，－(Hfix－li)ci＝zi(106)で与えられる( li ＝〈 ci | Hfix | ci 〉)．これは Etot でなく，期待値の和 ∑〈の規格直交条件付きの最小化のための
i ci|Hfix|ci〉－∑i li
勾配で，式( 9 )の代わりに Qtot[{ci}]＝∑〈
i ci|Hfix|ci〉〈( ci|ci〉－1)についての－dQtot[{ci}]/dciから導出される．
 節の FFT による効
Hfixci の演算は，上記のように第 6 章率化法を用いる．
期待値最小化なので，非占有状態も扱える．金属を扱う場



合，{ci}のセットに占有状態に加え非占有状態も一定数含めて，勾配ベクトル zi に掛ける処理が行われる(49)(52)(53) ．前る． GramSchmidt の直交化処理(式(101))を行いながら，節や本節の各手法では，この preconditioning 処理が行われ式( 103)～(105)の共役勾配法の手順で(式( 105)の混合比はた勾配ベクトルが使用される．

期待値最小化で決定)各状態 ci を各々数回更新する．全てのさて，上記の block Davidson 法は，アルゴリズムが簡便状態でこの過程を行った後，同じ •k の状態の間で subspaceで多用される．この手法の弱点は，状態数 M が大きい場

対角化を行う．{ci}を同じ •k の状態のセット(非占有含めて合，{ci}と{zi}による subspace の 2M × 2M のハミルトニア

M 個)として，その subspace で M×M のハミルトニアン行

ンの行列サイズが大きく，通常法で対角化するのに計算時間列がかかることである．対角化を多数回繰り返すので影響は大[Hfix]ij＝〈ci|Hfix|cj 〉(107)

きい．そこで， •k 点毎の全状態{ci}を低い順にいくつかのグを組み立てる．比較的小さいサイズなので，通常法による対ループに分け， subspace 対角化をグループ毎に行い，行列角化が可能で，固有ベクトル{ C ij }を求め，波動関数を次式

サイズを小さくする方策がある(54) ．状態数 M を n 分割しで更新する(i＝1～M)．

て，各グループの M / n 個の波動関数に関して，対応する
M＝∑ C ij cj
ci′

(108)

j＝ 1

subspace 対角化は，変分自由度を大きく下げた対角化なの{zi}も含めて 2M/n×2M/n 次元の subspace 対角化をグループ毎に行う．全ての波動関数は，GramSchmidt の直交化処理で(グループを超えて)下の準位と直交化させるので，下ので，本来の大規模ハミルトニアンの固有状態が得られるわけ

準位に落ち込む危険は避けられる．適宜，全状態での sub-ではないが，より固有状態に近い，直交化した波動関数セッ

space 対角化も行う．原子数が大きいスーパーセルでは，Mトが得られる．全 •k 点で行い，エネルギー期待値の低い順にも原子数の数倍以上になるので，重要である．並列計算機で並べることで，フェルミ準位や各状態の占有数が更新され，グループ毎に別のノードで扱えば効率的である．

output の r( r… )が求まり， charge mixing で次回 SCF ループ

本節で扱っている，固定した r ( r… )の Hfix の固有状態計算の input の r(r… )が与えられる(例えば第 3 章の式(25)，またを繰り返し法による波動関数の更新で行い， SCF ループのは後述の式( 112 )を用いる)．それにより Hfix 内の VH ， mxc

charge mixing で の r ( r… ) の 更 新 と 組 み 合 わ せ る 手 法 は ，を更新し，再び共役勾配法のループに入る．

{ ci }と r ( r… )の両者が全体として効率的に収束すればよく，

一方，r の固定した大規模ハミルトニアン行列 Hfix の繰り必ずしも r の更新の都度，Hfix の固有状態計算を厳密に行う返し法による固有値解法として， block Davidson 法が知ら

必 要 は な い ． そ の 意 味 で charge mixing 法 も 重 要 で ， 式

れている(52) ．これは，上記の •k 点毎の subspace 対角化を(25)の単純な方法に替わる方法として Kerker 法(20)がある．

{ ci }と式( 106 )の勾配{ zi }を合わせた subspace で(予め式( 101 )で直交化処理をした後)2M × 2M のハミルトニアン行＋1• )＋(G• )＝r N
rN
in
in(G

A|G• |2• )－r N
• ))
(r N
out(G
in(G
|G• |2＋G 20(112)

列(要素が〈ci|Hfix|cj〉，〈ci|Hfix|zj〉，〈zi|Hfix|cj〉，〈zi|Hfix|

これは，フーリエ成分での mixing で， N 回目の SCF ルーzj〉)を組み立て，対角化で得た固有ベクトル{C ij}で波動関数N
プでの input の r N
in と output の r out の差を用いて次回の N ＋
N＋1
1 回目の input の r in を組み立てる． mixing 係数 A | G• |2 /

を
M

2M

i
i
c′
i ＝ ∑ C j c j ＋ ∑ C j zj － M
j＝ 1

(109)

j＝ M ＋ 1(|G• |2＋G 20)は，振動の原因となる長波長成分(小さな G• の成• )－ r N• )に対して小さく，大きな G• では A に分)の r Nout ( Gin ( G

ことで，それらが混合した，より低い固有値の(真の固有状

漸近し，振動を抑える働きをする．パラメータ G 20が大きい
ほど，波長による係数 A|G• |2/(|G• |2＋G 20)の変化が顕著にな態により近い)波動関数が得られる．この過程を数回行い，る(G 20＝0 なら A|G• |2/(|G• |2＋G 20)＝A，式(25)と同じ)．ス

出力の r ( r… )から charge mixing で入力の r ( r… )を組み立て，

ーパーセルの大きさや形状，電子構造に応じて G 20 と A のVH，mxc を更新した Hfix で次回の subspace 対角化を行う．

最適値をテストして調整する．

のように更新する(i＝1～M)．勾配{zi}を subspace に加える

ところで，式( 99 )，式( 106 )で与えられる勾配 zi には，

図 11 ( b )に Bylander らの共役勾配法， block Davidson 法

preconditioning の問題がある(49)(52)(53)．式(106)で，現在のを式(112)の Kerker 法と組み合わせたテスト結果を示す(50)．

更新中の波動関数 c i と正しい固有関数 c0i の間に c 0i ＝ ci ＋

dci の関係があるとすると
zi＝－(Hfix－li)ci＝－(Hfix－li)(c 0i －dci)＝(Hfix－li)dci



残差最小化に基づく高速化技法

 節，
 節で論じた各種の基底状態計算の高速化技
上述の

(110)

法は，対象のサイズや電子構造(金属か非金属か)に応じて使となる．正しい固有状態なので Hfixc 0i ＝ lic 0i の関係を使っ

い分けられる(図 11 )．一方，クラスター計算機など並列計ている( li ≠ l0i は差が小さいとして無視する)．効率的に c i

算機の普及で，並列計算に適した手法が有利である．並列計を変化させ c 0i に近づけるには，勾配ベクトル zi よりも式算機では，各ノード( CPU )で別々の演算を行い，適宜，デ(110)の dci を扱う方が良いわけで，ータのノード間通信で集計する．{ci}のセットの最適化におdci＝(Hfix－li)－1zi(111)

いて，各波動関数 ci を別々のノードで扱い，頻繁なノードとなる．これは，勾配ベクトルに逆行列(Hfix－li)－1 を作用間通信なしに独立に固有状態に収束させられれば理想的な並させたものである．簡便のため，逆行列を対角項の逆数で近列計算になる． GramSchmidt の直交化処理(式(101))は，似し，対角項のみの preconditioning 行列[ K ]G• , G• を組み立

ノード間の大きな波動関数データの通信が必要なので，並列




講義ノート

計算効率を落とす．この直交化処理が少なくて済む方法が望ましい．
RMMDIIS(residual minimization/direct inversion in theiterative subspace)法(53)(55)(56)は，前節で議論した繰り返し法による固有状態計算法の一種であるが，提案者にちなんで
Pulay 法と呼ばれ，直交化処理が減らせると考えられる．前節までの方法は，全エネルギーやエネルギー期待値の最小化の見地から波動関数の更新，最適化を行うが， RMM DIIS
法では，残差( residual )の最小化の見地から波動関数の更新，最適化を行う．残差は，最適化途上の波動関数 ci について
Ri＝(Hfix－li)ci(113)である． li は期待値〈 ci | Hfix | ci 〉． Ri は ci と同じくサイズ
NG の列ベクトルで，固有状態に収束すれば Hfixci ＝ lici でゼロとなるので残差と言う．式(106)の勾配と符号が違うだ
けだが意味が異なる．本手法では，初期に input した波動関
数 c1i から更新した m ステップ目の c mi までの過去の m 個の波動関数{ c ji }( j ＝ 1 ～ m )の線形結合で次のステップの
＋1
を
cm
i
m

＋1
＝∑ a j c ji
cš m
i

(114)

＋1
＋1
＋1
＝cš m
＋hR̃ m
cm
i
i
i

(115)

j＝ 1

図 12

期待値最小化の方法(上)と残差最小化の方法(下)の概念の比較．期待値最小化法では，波動関数 ci の期待値
E i を下げる方向の勾配に従い，逐次に波動関数を変化させる．残差最小化法(RMMDIIS 法)では，過去数回の波動関数の線形結合で次回の波動関数を組み立てる．その線形係数は，過去の波動関数の残差 R i の線形結合の大きさが最小になるように決める．

で構築する．式( 114 )の線形係数{ a j }は，過去の波動関数{c ji}の式(113)による残差{R ji}の同じ線形係数{a j}の予測残いないと，複数の状態が同じ低い固有状態に近づく問題が起差R
m

R＝∑ a j R ji

(116)

j＝ 1

きる．本手法では残差を減らす方向なので必ずしもそうならない． Gram Schmidt の直交化や subspace 対角化は必要だの内積〈 R | R 〉が，式( 114 )の波動関数の規格化条件(ノルムが，頻度を減らすことができる．ノード間通信が減らせるのが 1)のもとで最小になる条件で決める．係数{a j}は，m×mで，並列計算に有利である．

行列[R ]lj，[S ]lj を用いた以下の固有値方程式の解となる．
m

m

j＝ 1

j＝ 1

∑ [R ]lj a j＝e∑[S ]lj a j(117)[R ]lj＝〈R li|R ji 〉，[S ]lj＝〈c li|c ji〉(118)

式(118)の各要素はサイズ NG のベクトルの内積である．

残差最小化法は， SCF ループの charge mixing 法としても有効である(53)(57) ．この場合の残差は， N 回目の SCF ループの input，output の差でN
RN(r… n)＝r Nout(r… n)－r in(r… n)(119)である．{r… n}は FFT の実空間メッシュ点．まず，式(112)の＋1
は式( 114 )の一方，式( 115 )は trial step と呼ばれ， R̃ mi

Kerker 法(20) に よる SCF ル ープを 数回 繰り返 し， input の＋1
を式( 113 )に入れて求めた残差で， preconditioning 処
cš mir jin(r… n)と式(119)の残差 R j(r… n)のデータをメモリーする．今，理した後，式(110)，(111)の dci と同じ意味でパラメータ hj ＝ N から j ＝ N ＋ l までをメモリーした後， N ＋ l ＋ 1 回目ので加える(定義上， h は負)．この step は線形従属を防ぐ意

input の r を以下の線形結合で与える．

＋1
味で不可欠．こうして，式(115)の c m
を加えて拡張したi

N ＋l

N ＋ l＋ 1(r… n)＝ ∑ a j r jin(r… n)
r insubspace { c ji } ( j ＝ 1 ～ m ＋ 1 ) で 次 の RMM DIIS 過 程 ( 式(114))を行う．m が進むほど，[R ]lj と[S ]lj の行列サイズが大きくなるが， m は数ステップで打ち切るので，通常の行列解法で解ける．全ての{ci}で独立にこの過程を行った後，N ＋l

係数{a j}は， ∑ a j＝1 の条件のもと，同じ係数を用いた残差
j＝ Nの線形結合からの予想残差

全 状 態 で Gram Schmidt の 直 交 化 と subspace 対 角 化 を 行

N＋ l

R(r… n)＝ ∑ a j R j(r… n)

い，出力の r ( r… )の計算， charge mixing による次回入力のr(r… )の構築，Hfix の更新を行う．以上の過程を繰り返す．
本手法の特徴は，次回の波動関数を，残差を最小化するように，過去のステップの残差情報から線形性を想定して組み立てる点にある(図)．各波動関数は，近くの固有状態(残差ゼロの状態)に近づいていく．従来手法では，全エネルギ
ーやエネルギー期待値を下げる方向の更新なので，波動関数が Gram Schmidt の直交化で下の準位と厳密に直交化してま て り あ 第61巻 第12号(2022)
Materia Japan

(120)

j＝ N

(121)

j＝ N

の大きさ(∑|R(r… n)|2，メッシュ点データの二乗和)が最小になr… nる条件から決める．式(117)に類似した方程式から係数{aj}が決まり，式( 120)で次回 input の r を構築する．その後，再び数回にわたり式(112)で更新を行い，式(120)，(121)から次回 input の r を更新する．この過程を繰り返す．
以上，本章では，大規模系の基底状態計算の高速化技法と

して， Etot の直接最小化法，大規模行列の繰り返し法(反復が，現在の汎用コードの主要な構成要素であり， 1980 年代法)による固有状態計算法，効率的な SCF ループの charge

から 90 年代にかけて，急速に進歩した．三つの柱のうち，

mixing 法等を論じた．これらは，最適化問題や大規模固有どれが欠けても今日の第一原理計算の隆盛はなかったろう．

値問題のクリロフ部分空間法など，広い分野の数値計算法に今世紀に入り，各種汎用コードの普及が飛躍的に進展し，物含まれるもので，平面波基底法以外も含め，第一原理計算と質・材料研究で不可欠のツールになっている．

密接に関係している(3)(4)．

ここで，第一原理計算の活用に関して，幾つか留意事項を記しておきたい．本講義では，基底状態の電子構造，全エネ

. ま

と

め

ルギー，原子に働く力を求める方法を中心に論じたが，各種
物性や諸性質―弾性的・機械的性質，格子振動や熱的性質，

本講義では，平面波基底の第一原理計算法として，まず，

相安定性，光学的性質や内殻励起スペクトル，電子伝導物

密度汎関数理論やブロッホの定理など第一原理計算の基礎を性，誘電関数や誘電率，磁気的性質，電気化学的性質，化学論じた(第 2，3 章第一回)．密度汎関数理論により電子間反応や拡散など―を平面波基底の第一原理計算の枠内で扱う多体相互作用を取り入れた電子構造や全エネルギーの定量的方法論や関連コードも開発されている．文献(4)や関連文献計算が可能である．全エネルギーを最小化する基底状態の電を参照されたい．

子密度分布を求める問題となり，KohnSham 方程式を自己材料科学・工学では，結晶自体よりも，表面・界面，格子

無撞着に解くことで達成される．結晶やスーパーセルのよう欠陥，粒界，不純物や溶質，異相界面や析出など，局所的にな周期系の KohnSham 方程式では，ブロッホの定理からブ

特別な原子配列や組成を持つ部分の電子状態や原子間結合がリルアンゾーン内の •k 点ごとに固有値，固有関数を求める．

鍵を握る場合が多い．平面波基底の第一原理計算では，大き次に，平面波基底を用いる第一原理擬ポテンシャル法の原理なスーパーセルの計算でそれらを含む系の安定原子配列や電と擬ポテンシャルの構築法を論じた(第 4 章第二回)．固子状態が高精度に得られる．しかし，局所的な原子間結合の体中の価電子の波動関数は，ポテンシャルが比較的平坦な原様子を深く探るには，局所状態密度(local density of states;

子間領域でスムーズであるが，原子核近傍では内殻軌道とのLDOS )や結合次数( bond order )の分析など，原子軌道基底

直交化のため振動する(ノードを持つ)．そこで，原子間領域を用いる局所解析法(58) を平面波基底の方法で求めた構造にで正確な挙動をし，原子核近傍(半径 rc 内)でノードを持た適用することが有益で，汎用コードも開発されている．平面

ない価電子波動関数を出力するような原子ごとの擬ポテンシ

波基底の手法の枠内では，局所エネルギー・局所応力の方法

ャルを構築し，価電子とイオン(擬ポテンシャル)の系に密度も有益である(59)．

汎関数理論を適用する．全体でスムーズな価電子波動関数は一方，第一原理計算の結果の意味を理解し，正しく解釈す

平面波基底で効率的に展開される．擬ポテンシャルは，自由るためには，様々な物質系での原子間結合の支配因子や価電原子の全電子軌道計算から rrc で正しい散乱の性質を再現

子挙動の知識を持っていることが望ましい．筆者の経験でするように構築される．次に，ハミルトニアンや全エネルギ

は，文献( 60 )～( 63 )などの書籍が有益であった．これらのーの具体的表式(逆空間表現)と演算方法を論じた(第 5 ， 6

書籍にあるように，原子間結合の支配因子や価電子挙動は，

章第三回)．格子周期関数であるポテンシャルや電子密度

単純金属，遷移金属，半導体，典型元素の化合物，遷移金属

分布は逆格子ベクトル・フーリエ級数展開で表現され，ハミ化合物等の各々の系で大きく異なる特徴を持つ．強結合近似

ルトニアンや全エネルギーは，波動関数の平面波基底展開の法(第 4 章(1)節の LCAO 法の簡略版)や経験的擬ポテンシャ

ベクトル係数やポテンシャル項・電子密度分布のフーリエ係

ル法(実験的パラメータを用いる簡便法)を用いたモデル解析

数，それらの積や総和でシンプルに表される．擬ポテンシャ

で，本質が捉えられてきた．もちろん，第一原理計算でさら

ルの非局所項は原子周りの軌道角運動量成分毎の射影演算子に豊かな解明が進展している．

を通じて扱われる．無限に繰り返す系の静電相互作用の発散

最後に，筆者は， 1980 年代に研究生活をスタートし， 90

成分は厳密に打ち消しあう．高速フーリエ変換の活用で演算

年代から独力での第一原理計算プログラムの開発に着手しが大幅に効率化される．最後に，第一原理分子動力学法を可た．汎用コードなど望むべくもない時代で，論文情報のみで能とする基底状態計算の高速化技法を論じた(第 7 章今徒手空拳であったが，次第に仲間も広がった．本講義で扱っ回)．これは，大規模系の基底状態計算において， Kohn た手法や計算技術の多くは，実地にコードを開発し，有用性

Sham 方程式の固有状態計算をハミルトニアンの行列対角化を確認してきたもので，また，新しい手法や計算技術への展法で行うのではなく， input した個々の波動関数を，規格直開も進めてきた(59)．CP 法以降，世界的に手法やコードの開交化条件のもと，勾配や残差に従い，繰り返し法(反復法)で発ラッシュの時代となり，その中で研究できたことは excit-

逐次に変化させ，subspace 対角化も駆使して，全エネルギー

ing な体験であった．同時期には，スパコンや並列計算機のやエネルギー期待値，或いは残差の最小化の見地で効率的に進歩，電子顕微鏡や STM 等観察技術の革新もあり，第一原

最適化し，SCF の基底状態に到達させるアプローチである．

理計算の精度と有効性が実証され，計算科学と実験・観察の 密度汎関数理論やブロッホの定理のよう
以上のように，◯連携が広く展開されるようになった．一方，機械学習等 in- 擬ポテンシャル法のような平面波基底にな基礎「理論」，◯

formatics と第一原理計算を組み合わせた材料探索や材料設

 CP 法に端をよる電子構造計算の効率的「手法」，そして◯計が材料開発で注目されている．汎用コードの普及もあり，発する大規模系の基底状態の高速「計算技術」の三つの柱

現在では，実験家も含めて，材料研究者が第一原理計算を実




講義ノート

行したり，第一原理計算との共同を進める機会も多い．そう

一方，前記の数学公式において a ＝ r ， b ＝ 0 ， x ＝ t とするした際に本講義ノートが役に立てれば幸いである．さらにま∞ － r 2t 2
∞ － r 2t 2とf
dt ＝ p / 2r で， 1 / r ＝ 2 / p f
dt から次式
0 e
0 e

た，若い方が従来の理論や手法，計算技術のエッセンスを習が成り立つ．

得し理解することで，新たな手法・計算技術の開発など，計算材料科学の発展に繋げていただくことを期待したい．

1
2＝∑• －r… | R•
p
R• |R

∞

2
p

g

∑

＝

補論 A

f
∑
f
∑
f
∑
[
R•

 節の Ewald 法の式( 79 )，( 80 )の導出を説明す
第 6 章

＝る．まず，次の格子周期関数を考える．
2∑ e－|R• －r… |2t 2
p R•

F(r… , t )＝

(A1)

＝

周期系の格子ベクトルR• の総和で，格子周期関数である．この関数の逆格子ベクトル・フーリエ級数展開∑ F(G• )e iG• ･r… を

e－|R• －r… | t dt

e－|R• －r… | t dt＋∑

f
f∑
f

＝

f∑

2
p Qc

Qc R•

2
＝
p Qc

Q

Qc R•

2p
eiG• ･r…
Qc G•

g

2p
eiG• ･r…
Qc G•

2e－|G• | /4t|G• |2R•

2 2e－|R• －r… | t e－iG• ･(r… －R• )d…r
2 2

f
f
f

e－|r… | t e－iG• ･r… d…r＝
2 2

Q

(A2)

e－r t e－i|G• |r cos u r2 sin u drd u dq
22

ei|G• |rv

e－r t r2
2 2

i|G• |r0∞

2 2
e－r t r 20p

lim eiG• ･r…

c G• →0

＝

v＝1

t＝ 0

erfc [|R• －r… |g]|R• －r… |

(A6)

e－r t
sin|G• |r
2t 2
2 2

－

|G• |
2t 2

exp [－|G• |2/4g 2]p＝ 2 lim|G• |2/4g 2
g Qc G• →0

2

4p
1
p
－
lim
Qc G• →0 |G• |2 g 2Qc

2

(A7)

1
p
exp [－|G• |2/4g 2]＝ 2 ∑ e iG• ･r…
• －r… | g Qc G• ≠0|G• |2/4g 2
R• |R

＋∑
R•

erfc [|R• －r… |g]
p
4p
－ 2 ＋
lim |G• |－2 (A8)|R• －r… |g Qc Qc G• →0

この式の r… を－r… aa ′に替えれば，式(80)になる(証明終)．

r＝∞次に式(79)を証明するため，式(A8)において r… ＝0 を考え
r＝ 0∞

e－r t cos|G• |r dr
2 20

る．式(A8)の左辺の R• の和と右辺の第二項の R• の和におけ}る R• ＝0 の項，つまり，1/|r… |と erfc[|r… |g]/|r… |を取り出し，右辺に移し，| r… |を r に変えると( erfc ( r g )－ 1 )/ r である． r∞

2p p －|G• |2/4t 2
2 2ee－r t cos|G• |r dr＝ 2t 2t
0→0 の極限を考えると
1

(A3)

p

R• ≠0

erfc[|R• |g]erfc(rg)－1
p
＋lim
－ 2
•
|
|
g
r
Qc
R
→0
r
R• ≠0

＋∑

v＝－ cos u，dv＝sin udu の変数変換を行っている．六，七
行目は部分積分で，左の項はゼロになる．最後から二行目は，
22
p /2a e－b /4a の数学公式(40)を用いている．(A2)，(A3)から
2p
2
2
F(G• )＝ 3 e－|G• | /4t
Qc t

(A4)

である．(A1)と組み合わせて以下になる．
2
2p1
F(r… , t )＝∑ e －|R• －r… |2t 2＝ ∑ 3 e－|G• |2/4t 2 eiG• ･r…
p R•
Qc G• texp[－|G• |2/4g 2]|G• |2/4g 2
G• ≠0∑ |R• |＝g 2 Qc ∑となる．極座標の積分は u＝0～p，q＝0～2 p，r＝0～∞で，ま て り あ 第61巻 第12号(2022)
Materia Japan(|G• |1/4g －1)∑

2 sin|G• |rdr|G• |rp 3 /2＝ 3 e－|G• |2/4t 2
t

f

R•

となる．式(A6)に入れると次式になる．

dr

v＝－1

0

∞ － a 2x 2
cos bx dx ＝
0 e•

] ＋∑ erfc|[|R• －R－r… |r… |g]

この式( A6 )の最終形の第一項の G• での和は， G• ＝ 0 の項

2 2

＋

2p
＝ 2
t

t＝ g

が発散成分を持つので，これだけ取り出して検討する．小さい G• での展開 exp [－|G• |2/4g 2]1－|G• |2/4g 2 を使うとr e－r t sin|G• |r dr

4p
＝
|G• |

2

になる．

∞

4p
|G• |

2

える．四行目の第一項で，t の積分を実行．t＝0 の項はゼロ

g2 Q

f
f
f [ ]
f
f
{[
]
f
f
∞

＝

2 2

二行目で t による積分の範囲を g で分ける．三行目では，二2 2
e－r t ei|G• |rvr 2drdv＝2 pe－|R• －r… | t dt

g

∞ －t 2分を補誤差関数 erfc[x ]＝ 2/ p f
dt を使った表現に替x eの単位胞の積分の和を結晶全体(全空間)での積分に替えている．この最終形の積分は，G• 方向を z 軸とする極座標積分で

＝2 p

f

(A5)のフーリエ級数展開に替え，第二項の t＝g から∞の積

2 2e－|r… | t e－iG• ･r… d…r

＝2 p

∞

行目の第一項の t ＝ 0 から g までの積分の被積分関数を式

式( 16)の G• ･R• ＝ 2p ×整数を使っている．最後の形は，各 R•

f2
p

1 －|G• |2/4t 2
erfc [|R• －r… |g]edt＋∑
3|R• －r… |
•
t
R
0

＋∑
e－|R• －r… | t e－iG• ･r… d…rR•0

考えると，式(39)より以下のようになる．
2
F(r… , t )e－iG• ･r… d…r＝
p
Qc
Qc

2 2

p
exp [－|G• |2/4g 2]＝ 2 ∑eiG• ･r…
|G• |2/4g 2
g Qc G•G•1
F(G• )＝
Qc

2 2

0

＋

4p
lim |G• |－2
Qc G• →0となる．右辺第三項は，以下のように変形される．
∞ －t 2∞ －t 2erfc(rg)－1
2
dt－f
dt
f
rg e
0 e
＝
lim
r
r
r→0
p r→ 0

lim

＝－
(A5)

(A9)

＝－

2

fr0g e－t dt
2

lim

r

p r→0
2

lim
p r→0

2 2
ge－r g

1

＝－

2g

p

(A10)




∞ －t 2
∞ －t 2
erfc [ x ]＝ 2 / p f
dt ， erfc [ 0 ]＝ 2 / p f
dt ＝ 1 を
x e
0 e

文

使用し，最後の変形は分子と分母を r で微分している．式献(A9)，(A10)から

p1

exp[－|G• |2/4g 2]|G• |2/4g 2
G• ≠0∑ |R• |＝g 2 Qc ∑R• ≠0erfc[|R• |g] 2gp4p－－ 2 ＋ lim |G• |－2
•
|
|
g
Q
Qc G• →0
R
p
c
R• ≠0

＋∑

(A11)
となる．この式は式(79)にほかならない(証明終)．
補論 B
 節の式( 84 )のクーロンポテンシャルのフーリエ
第 6 章

変換の導出を説明する．一般の波数ベクトルを …q としてf(－e 2 Za/r)exp[－i …q･r… ]d…r＝－

4p e 2 Z a| …q|2(B1)の証明である．格子周期関数ではないので，積分は全空間である．クーロン形の積分は，長距離で問題が生じる可能性があるので，k を正の定数として，e－kr をかけた関数(－e 2 Za/r ) e－kr を扱い，後で k → 0 を実行する． …q を z 軸とする極座

標積分で

( )
f
f
f
f
f
f
f [
f
{[
]
}
f
{ f
－

e 2 Z a － kreexp[－i …q･r… ]d…rr＝－e 2 Zae － krexp[－i| …q|r cos u]r2 sin u drd u d qr＝－2pe 2 Za

re－kr exp[i| …q|rv]drdv

∞

re－kr

＝－2pe 2 Za0

－ 4 pe 2 Z a＝| …q|
4 pe 2 Z a
＝－
| …q|

exp[i| …q|rv]
i| …q|r(44) W. E. Pickett: Comp. Phys. Rep., 9(1989), 115197.
(45) J. L. Martins and M. L. Cohen: Phys. Rev. B, 37(1988), 6134
6138.
(46) R. P. Feynman: Phys. Rev., 56(1939), 340343.
(47) R. Car and M. Parrinello: Phys. Rev. Lett., 55(1985), 2471
2474.
( 48 ) 日本機械学会編原子・分子モデルを用いる数値シミュレー
ション，コロナ社，(1996).
(49) M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D.
Joannopoulos: Rev. Mod. Phys., 64(1992), 10451097.
(50) M. Kohyama: Model. Simul. Mater. Sci. Eng., 4(1996), 397
408.
(51) D. M. Bylander, L. Kleinman and S. Lee: Phys. Rev. B, 42(1990), 13941403.
(52) D. Singh: Phys. Rev. B, 40(1989), 5428
5431.
(53) G. Kresse and J. Furthm äuller: Phys. Rev. B, 54(1996), 11169
11186.
(54) N. Kosugi: J. Comp. Phys., 55(1984), 426
436.
(55) P. Pulay: Chem. Phys. Lett., 73(1980), 393398.
(56) 香山正憲，田中真悟，岡崎一行材料，52(2003), 260265.
(57) P. Pulay: J. Comput. Chem., 3(1982), 556
560.
(58) S. Maintz, V.L. Deringer, A. L. Tchougreeff and R.
Dronskowski: J. Comput. Chem., 37(2016), 10301035.
(59) M. Kohyama, S. Tanaka and Y. Shiihara: Mater. Trans., 62(2021), 1
15.
( 60 ) W. A. ハリソン固体の電子構造と物性―化学結合の物理―(上，下)，現代工学社，(1984).
(61) D. G. ペティフォー分子・固体の結合と構造，技報堂出版，(1997).
( 62 ) 津田惟雄，那須奎一郎，藤森 淳，白鳥紀一電気伝導性酸化物(改訂版)，裳華房，(1993).
(63) 小口多美夫遷移金属のバンド理論，内田老鶴圃，(2012).

v＝1

]

dr
(全ての図番号および文献番号等は，第一回からの通しとなってお

v＝－1

∞

ります)

e－kr sin| …q|r dr(完)

0e－kr cos| …q|r－| …q|∞★★★★★★★★★★★★★★★★★★★★★

0∞e－kr cos| …q|r
－k
dr
| …q|
0

＝－

4 pe 2 Z a 1
－k
| …q|
| …q|

∞

e－kr cos| …q|r
dr
| …q|
0

}

(B2)

極座標積分は， u ＝ 0 ～ p ， q ＝ 0 ～ 2 p ， r ＝ 0 ～∞で， v ＝－cos u，dv＝sin u d u の変数変換を行っている．六行目の r＝∞の項はゼロになる．式(B2)の最終形に対して，k→0 で第二項が消え，式(B1)と同じく－4pe 2 Za/| …q|2 となる(証明
終)．

香山正憲

1985年 東京大学大学院工学系研究科(金属材料学専攻)博士課程中退
同年 工業技術院 大阪工業技術試験所(現産業技術総合研究所関西センター)入所
2004年 産業技術総合研究所 ユビキタスエネルギー研究部門 グループ長
2015年 産業技術総合研究所 エネルギー・環境領域
電池技術研究部門 首席研究員
2021年～現職
専門分野計算材料科学，粒界・界面・ナノ構造の材
料科学，新規手法・コードの開発
Modelling and Simulation in Materials Science andEngineering 誌 Editorial Board membermkohyama＠aist.go.jp★★★★★★★★★★★★★★★★★★★★★

計算材料科学に導いていただいた堂山昌男先生，山本良一
先生，多くの議論をしていただいた共同研究者の皆様にあらためて感謝申し上げます．本講義ノートの内容について丁寧な助言をいただいた圓谷貴夫博士にお礼申し上げます．




講義ノート