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

山

正

憲

(第二回からの続き)

るが間違い)．証明は以下． G• ′－ G• ＝ l1 •b1 ＋ l2 •b2 ＋ l3 •b3 ( li は整.

ら基本格子ベクトル …a1 ， …a2 ， …a3 の三方向に沿う平行六面体数)とし，単位胞内(体積 Qc ＝ …a1 ・ …a2 × …a3 )の積分は，原点か

第一原理擬ポテンシャル法での平面波基底とハミルトニアン
格子周期関数のフーリエ級数展開

なる(16)(39)．

擬ポテンシャルのもと，全領域でスムーズな価電子波動関数を平面波基底で展開する．本章では，波動関数やハミルト
ニアンの平面波基底表現を説明する．関連する重要項目として，格子周期関数の逆格子ベクトル・フーリエ級数展開から説明する．格子周期関数は，式( 14 )の周期系のポテンシャ
ル Veff ( r… )，式( 17 )の固有関数内の U •kn ( r… )，電子密度分布r(r… )など多数あり，各単位胞で同じものが繰り返す．

エ級数展開は
(38)

G•

f

1
Qc

f(r… ) exp [－iG• ・r… ]d…r

(39)

Qc

で与えられる．Qc は単位胞の体積．G• の exp [iG• ・r… ]自体が格子周期関数である(式( 16 )から exp [ iG• ・( r… ＋ R• )]＝ exp[iG• ・r… ])．式(39)で f(r… )から様々の G• の{f(G• )}を求めることをフーリエ変換，逆に式( 38)で{f( G• )}から f( r… )を組み立てることをフーリエ逆変換という(本ノートでは括弧{ }は様々な G• のフーリエ成分やベクトル係数の集合の意味)．
式(38)，(39)の展開には，以下の関係式が必要である．

f

1
Qc

Qc

・r… ]d…r＝dG• G• ′exp [－iG• ・r… ] exp [iG• ′・r… ]dr…
exp [－iG• ・r… ] exp [iG• ′

Qc

＝

1
Qc

－G• )・r… ]d…r
exp [i(G• ′

Qc

1
＝ Qc
Qc

3

1
0

dl1dl2dl3 exp [ilj •bj・lj …aj]

j＝1

＝

0

＝

dl1dl2dl3 exp [i2pl1l1] exp [i2pl2l2] exp [i2pl3l3]

exp [i2pl1l1]
i 2p l1

10

exp [i2pl2l2]
i 2p l2] [exp [i2i2ppl l l ]]10＝dG• ′－G• , 0＝dG• ,G• ′

3 3

3

1
0

(41)

式(15)から …ai・ •bj＝2pdij を使っている．
exp [iG• ・ r… ]は，同じ位相の点 r… が G• に垂直な平面になるとなる．フーリエ係数(フーリエ成分)はf(G• )＝

f
f
∑
f
f
[
][

1
Qc

1

一般に格子周期関数 f(r… ＋R• )＝f(r… )につき，R• の格子系から組み立てた逆格子の逆格子ベクトル G• を用いて，フーリ
f(r… )＝∑f(G• ) exp [iG• ・r… ]内の積分で， l1 …a1， l2 …a2， l3 …a3 として各々 l1， l2， l3 の 0～
1 の積分に変換する．変数変換のヤコビアンが Qc で以下と(40)

積分範囲は単位胞内である(全空間としているテキストもあ波(平面波)で，波長(同じ位相の平面間の間隔)が 2p/|G• |である．例えば G• ＝m1 •b1 の場合，exp [iG• ・r… ]は，式(15)から…a1， …a2， …a3 の平行六面体で形成される単位胞のうち， …a2， …a3

で作る面と面の間を， …a1 に沿って平行に m1 分割する波長の波である．格子周期関数の単位胞内での値の変化を様々な波長と進行方向の平面波の重ね合わせで表わす．変化がスムー
ズでない場合，式( 38 )の展開は， G• ＝ 0 から大きな G• (短波長の G• )まで必要で，G• ＝ m1 •b1＋m2 •b2＋m3 •b3(mi は整数)の大
きな絶対値の mi まで含めた展開になる．
こうして，格子周期関数は，そのフーリエ係数のセット，
例えば{Veff(G• )}，{U •kn(G• )}，{r(G• )}で扱われる．実際の演算は，高速フーリエ変換(fast Fourier transformation; FFT)

国立研究開発法人 産業技術総合研究所 エネルギー・環境領域 電池技術研究部門名誉リサーチャー(〒5638577 池田市緑が丘 18
31 産総研関西センター)
Lecture Notes on FirstPrinciples Methods Using a PlaneWave Basis Set (Part 3); Masanori Kohyama(Research Institute of Electrochemical Energy, Department of Energy and Environment, National Institute of Advanced Industrial Science and Technology, Ikeda,
Osaka)
Keywords: plane-wave basis, Fourier series expansion, Hamiltonian matrix, total energy, separable non-local pseudopotential, Ewald method, fast
Fourier transformation
2022年 2 月21日受理[doi:10.2320/materia.61.782]




講義ノート

により効率化される．FFT の詳細と G• の条件(mi の範囲)は〈f|A|ψ〉＝ff(r… )Aψ(r… )d…r である．

 節)．
後述する(第 6 章

ここで，平面波基底間の規格直交性は

平面波基底での展開

ff1〈 •k＋G• | •k＋G• ′
〉
＝
Q

Q

1
Q

Q

式( 17 )の固有関数形ψ •kn ( r… )＝ e i •k・r… U •kn ( r… )で，式( 18 )の格
子周期関数 U •kn ( r… )のフーリエ級数展開({ U •kn ( G• )}が展開係＝)・r… ]d…r
exp [－i( •k＋G• )・r… ] exp [i( •k＋G• ′－G• )・r… ]d…r＝dG• ,G• ′exp [i(G• ′(44)となり，式( 40 )，( 41 )と同様に証明される．この式の積分数)から
U •kn(G• ) exp [iG• ・r… ]∑•ψ •kn(r… )＝e i •k・r… U •kn(r… )＝e i •k・r…

領域は，結晶 Q 全体である(各単位胞 Qc 内でも直交性が証明

G＝∑U •kn(G• ) exp [i( •k＋G• )・r… ]

(42)

G•

となる．この式は波動関数の平面波基底展開と見なせる．
exp [ i ( •k ＋ G• )・ r… ]が平面波基底，{ U •kn ( G• )}が展開係数で，固有状態を指定する •k を固定し( •k は第一ブリルアンゾーン内)， G• を変えた和である．当然だが，平面波基底自体がe i •k・r… e iG• ・r… ，e iG• ・(r… ＋R• )＝e iG• ・r… でブロッホの定理を満たす．

一方，前述(第 3 章図 1 )の結晶のマクロの周期境界条件(ボルン・フォンカルマンの周期条件(16))で，単位胞が …a1，
…a2， …a3 方向に各々 N1 …a1, N2 …a2, N3 …a3 の周期で繰り返す，体積
Q＝NQc(N＝N1N2N3，N は巨視的な数)の結晶部分で存在確
率 1 に規格化された平面波基底| •k ＋ G• 〉＝ Q－1/2 exp [ i ( •k ＋G• )・r… ]を考えると

ψ •kn(r… )＝∑Q1/2U •kn(G• )| •k＋G• 〉＝∑C n•k＋G• | •k＋G• 〉
G•

(43)

G•

される)．固有関数自体の結晶での規格化は

f

f

1
ψ
•kn(r… )ψ •kn(r… )d…r＝∑∑
Q
G•
G• ′

Q

G• C n•k＋G• ′exp [i(G• ′－G• )・r… ]d…r
C n•k＋

Q＝∑∑C n•k＋G• C n•k＋G• ′dG• ,G• ′＝∑|C n•k＋G• |2＝1
G• G• ′

(45)

G•となり，固有ベクトルの規格化条件である．
固有状態は，固定した •k に対し G• の異なる平面波で展開するので， •k の異なる平面波間の直交性は扱われないが，一般に， •k の異なる平面波間の結晶全体での直交性は，式(41)と同様の積分(ただし，積分領域が N1 …a1, N2 …a2, N3 …a3 での平行六面体の結晶 Q)で証明できる．これは，ボルン・フォン
カ ル マ ン の 周 期 条 件 に 関 わ る •k 点 の 定 義(16) か ら ， N1 …a1,
N2 …a2, N3 …a3 の平行六面体の向かい合う二面間で exp [ i ( •k －)・ r… ]の位相が揃い，式(41)の最終形のように消えるため
•k′
である．

となる．固有ベクトルは{ C n•k＋G• }＝{ Q1/2U •kn ( G• )}である．

平面波基底| •k＋G• 〉は，フーリエ級数展開の G• の平面波とC n•k＋G• の右肩の n はバンド指標で， •k について下から n 番目の固有状態の意味( n 乗の意味ではない)．| •k ＋ G• 〉の表記同様， •k＋G• による波長 2p/| •k＋G• |と進行方向を持つ．価電は，規格化された平面波基底を表し，逆向きの表記〈 •k ＋ G• |は，左から複素共役の基底をかけて体積 Q で積分する際に用いる．一般に，波動関数 f ( r… )，ψ( r… )，演算子 A (ハミル子の波動関数がノードを持たないスムーズなものであれば，式( 43 )は，それほど大きくない •k ＋ G• (それほど短くない波長)の平面波の展開で表現できる．格子周期関数のスムーズさとフーリエ級数展開の G• の範囲の関係と同様である．

トニアンやポテンシャル)で，ブラ“〈”，ケット“〉”を用い

平面波の運動エネルギーはた表現の意味は，| f 〉＝ f ( r… )，〈 f |ψ〉＝ ff( r… )ψ( r… ) d…r ，〈 •k＋G• |－

f

2 2
1
: | •k＋G• 〉＝2m
Q( 2m: ) exp [i( •k＋G• )・r… ]d…r＝2m| •k＋G• |exp [－i( •k＋G• )・r… ] －

Q

2

2

2

2

(46)

で，自由電子のそれに対応する．平面波基底展開の G• は，
G• ＝ 0 から始めて，平面波の運動エネルギーの打ち切りエネ

束の様子から決める(図)．通常，元素種毎の擬ポテンシャ
ル形状が Ecut を決める因子である．一方，逆空間の G• 点のルギー(カットオフエネルギー)Ecut から

密度，つまり，逆空間で同じ半径 Gmax の球内に存在する G•22m| •k＋G• |2Ecut(47)を満たす範囲でとる．平面波基底の波長の下限の指定に相当
2

す る． •k は 第一 ブリ ルア ン ゾー ン内 なの で， 実 質的 に 2m
2|G• |2Ecut で，逆格子空間で G 2max＝Ecut から半径 Gmax の
2m

球内の G• 点を用いる．

点数( NG)は，実空間の単位胞体積に比例する．単位胞体
積とブリルアンゾーン体積(逆空間での G• 点当たりの体積)が互いに反比例するからである．なお，平面波のカットオフ
エネルギー Ecut ( Gmax )と前述のフーリエ変換で用いる G• 点
 節で説明する．
の範囲の関係については，第 6 章

ところで，上記では，平面波基底も固有関数も，体積 Q

この条件から，平面波基底の総数 NG が決まる． NG は計の結晶で規格化されているとした．それで問題ないが，後述

算の負荷を 決める(ハミルトニアンの行列サイズが NG ×の単位胞(スーパーセル) Qc 当たりの全エネルギーの議論でNG )．波動関数がノードを持たずスムーズなら(擬ポテンシ

は，単位胞 Qc 当たりで規格化されているとして扱う方が表

ャルが浅くスムーズなら)Ecut が小さくて済み，NG が抑えら現上便利な場合もあり，適宜，そうした扱いも検討する．

れる．必要な Ecut は，扱う物質について幾つかの値で実際に第一原理計算を行い，Ecut に対する全エネルギー Etot の収

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







平面波展開と対称操作

平面波展開の固有関数について，結晶系の持つ対称要素(対称操 作 ) { S | t• S } の 効 果 を 考 え る ． •k と S •k で の 固 有 ベ ク ト ル
節で説{C n•k＋G• }，{C nS •k＋G• }の互いの関係の問題である．第 3 章

明したように{S|t• S}ψ •kn(r… )＝ψS •kn(r… )(48)

である(位相因子 eiu がかかる不定性はある)．式(43)を左辺に入れると，式(20)から

図8

平面波のカットオフエネルギー Ecut に対する全エネル
ギー Etot の収束の様子(35) ． diamond について，作成条件や作成方法を変えた各種ノルム保存擬ポテンシャル
での比較．Present Method は，文献(35)の TM 型擬ポ
テ ンシ ャルの 結果 ．ノル ムを保 存し ない USPP 法 ，
PAW 法では，もう少し小さい Ecut で収束する．なお，
1Ry13.6 eV である．

{S|t• S}ψ •kn(r… )＝∑C n•k＋G• {S|t• S}| •k＋G• 〉＝Q－1/2∑C n•k＋G• exp [i( •k＋G• )・S－1(r… －t• S)]
G•

G•＝Q－1/2∑C n•k＋G• exp [i(S •k＋SG• )・r… ] exp [－iS •k・t• S] exp [－iSG• ・t• S]
G•

C n•k＋S G• exp [i(S •k＋G• )・r… ] exp [－iG• ・t• S]∑•＝Q－1/2 exp [－iS •k・t• S]－1

G＝exp [－iS •k・t• S]∑exp [－iG• ・t• S]C n•k＋S－1G• |S •k＋G• 〉

(49)

G•

3 行目で，G• の和を SG• の和に替えているが，G• は対称操作〉を考える．同じ •k に対し G• を変の行列表現〈 •k＋G• |H| •k＋G• ′の行列 S で一対一に移りあうので問題ない．式(48)の右辺は(50)
ψS •kn(r… )＝∑C nS •k＋G• |S •k＋G• 〉

えた NG × NG 行列である．運動エネルギー項は式( 44 )，(46)からG•〈 •k＋G• |－である．従って以下のようになる．
C nS •k＋G• ＝exp [－iS •k・t• S] exp [－iG• ・t• S]C n•k＋S－1G•(51)

こうして， •k の固有ベクトル{ C n•k＋G• }を求めれば， S •k の固有ベクトル{C nS •k＋G• }も式(51)から組み立てられる．固有ベク
トルには，位相因子 eiu をかける不定性があるので，式(51)の G• に依らない exp [－ iS •k ・ t• S]は無視しても構わない．電子密度分布や全エネルギー，原子に働く力等の計算(第 6章)では，固有ベクトル成分を含む項の •k での積分が行われる．ブリルアンゾーン全体(あるいは半分)の領域の •k 点の固有ベクトルが必要だが，既約領域の •k 点の固有ベクトルを求めれば，既約領域外の S •k の固有ベクトルも式( 51 )で求ま
り，計算量が節約できる．


平面波基底でのハミルトニアン運動エネルギー項と局所ポテンシャル項

第一原理擬ポテンシャル法では，ハミルトニアン，全エネ
ルギー，原子に働く力等が，平面波展開や逆格子ベクトル・
フーリエ級数展開によりシンプルに表現され，特に逆空間表
現と呼ばれる．第 3 章や文献( 16 )で説明したように周期構造の系を扱うからである．
Kohn Sham 方程式( 4 )のハミルトニアンの平面波基底で




2 2
2
〉＝ | •k＋G• |2dG• ,G• ′: | •k＋G• ′
2m
2m

(52)

で，対角項のみゼロでない．
ポ テ ン シ ャ ル Veff は ， 原 子 か ら の 擬 ポ テ ン シ ャ ル
VPS ( r… )，価電子の電子密度分布からの静電ポテンシャル
VH(r… )，交換相関ポテンシャル mxc(r… )の和である(式(26))．
VPS ( r… )は，原子の擬ポテンシャル V PS
a の総和で(式( 27 ))，V PS
a は非局所項と局所項の和(式( 34 ))であり， VPS ( r… )は以下のようになる．
PS
VPS(r… )＝V PS
local(r… )＋V NL(r… )＝∑∑V alocal(r… －t• a－R• )＋∑∑V aNL(r… －t• a－R• ) (53)
R• a

R• a

R• の和は結晶の全格子点，a の和は単位胞内の内部座標 t• a の

全原子．V aNL(r… －t• a－R• )は t• a＋R• の位置にある非局所項∑|l〉l〈 l |を意味する．周期系全体のポテンシャル Veff のDV NLa, l ( r )内， V PS
local ( r… )と VH ( r… )， mxc ( r… )は，射影演算子を含まない局
所形であるので，まとめて Vlocal と表すと，式(26)は以下のようになる．
PS
Veff(r… )＝V PS
local(r… )＋V NL(r… )＋VH(r… )＋mxc(r… )＝Vlocal(r… )＋V PS
NL(r… )

(54)

まず，Vlocal(r… )の平面波基底間の行列要素を考える．格子
周期関数であるから，式(38)，(39)のフーリエ展開表示で講義ノート

Vlocal(r… )＝∑Vlocal(G• ) exp [iG• ・r… ]

(55)

f

(56)

積分は単位胞(スーパーセル)Qc 内でよい．行列要素は，G•

Vlocal(G• )＝

1
Qc

Qc

Vlocal(r… ) exp [－iG• ・r… ]d…r

〈 •k＋G• |Vlocal(r… )| •k＋G• ′
〉
＝

f

1
Q

Q

)・r… ]d…r
exp [－i( •k＋G• )・r… ]∑Vlocal(G• ″) exp [iG• ″・r… ] exp [i( •k＋G• ′
G• ″

f

1
＝∑Vlocal(G• ″)
Q
G• ″

Q

＋G• ″－G• )・r… ]d…r＝∑Vlocal(G• ″)dG• ″,G• －G• ′exp [i(G• ′G• ″＝Vlocal(G• －G• ′)(57)となる．式(40)を使っている．これは Vlocal(r… )のフーリエ変出される．式(39)から r(r… )のフーリエ成分 r(G• )を計算すれ

換，式(56)の G• －G• ′の項である．

ばよい． mxc ( G• )も同様に，実空間メッシュ点毎の r ( r… )の値

Vlocal(r… )＝V PS
local( r… )＋ VH( r… )＋ mxc( r… )の各項も格子周期関数• )， VH ( G• )， mxc ( G• )の G• －で，各項のフーリエ成分 V PS
local ( Gr( G• )の効率的計算法(高速フーリエ変換の利用)については項を与えればよい．静電ポテンシャル VH ( r… )は，ポアソ
G• ′


 節)．
後述する(第 6 章ン(Poisson)方程式:2V

H(r… )＝－4pe2r(r… )(58)

等から mxc ( r… )を求め，式( 39 )から計算する．なお， r ( r… )や

Vlocal 中 の 残 り の 原 子 の 局 所 擬 ポ テ ン シ ャ ル の 総 和• )は，∑∑V alocal(r… － t• a－ R• )のフーリエ成分 V PS
local ( GV PS
local ( r… )＝R• aから，フーリエ成分の関係式として
VH(G• )＝4pe2r(G• )/|G• |2(59)

式(39)からとなる．これは， :2VH (r… )＝∑VH ( G• ) :2 exp [ iG• ・ r… ]から導
G•

f
∑∑
f

• )＝
V PS
local(G

1
Qc

Qc

＝

1
Qc

Qc R• a• ・r… ]d…r
V PS
local(r… ) exp [－iG
V alocal(r… －t• a－R• ) exp [－iG• ・r… ]d…rf

1
＝ ∑∑exp [－iG• ・(t• a＋R• )]Qc R• aQc

V alocal(r… －t• a－R• ) exp [－iG• ・(r… －t• a－R• )]d…rff

1＝ ∑exp [－iG• ・t• a]∑ V alocal(r… －t• a－R• ) exp [－iG• ・(r… －t• a－R• )]d…r
Qc a
Qc
R•
1
＝ ∑exp [－iG• ・t• a]
Qc a

Q

V alocal(r) exp [－iG• ・r… ]d…r

1
＝ ∑exp [－iG• ・t• a]V alocal(G• )Qc a

最後から二行目で，前の行の単位胞 Qc での積分が， R• の全単位胞の和により結晶全体 NQc＝Q(事実上全空間)での積分に変わり，積分の原点を単位胞内の t• a に変えている． V alocal(r )はクーロン形が遠方まで及ぶのでリーズナブルである．
最終形の a の和は，式( 53 )の説明と同様，一つの単位胞内の全原子である．

(60)

f
f
f
f
f

V alocal(r) exp [i|G• |rv]r 2drdvdq

＝

∞

＝ 2p

0

＝

4p
|G• |

V alocal(r)r 2

∞

0

•

[exp i[|iG|• G|r|rv]] dr

V alocal(r)r sin(|G• |r)dr

1

－1

(61)

式( 60 )の最後の V alocal( G• )の計算を説明する．これは，格となる．積分は G• ベクトル方向を z 軸とした極座標積分 d…r子周期関数のフーリエ変換ではなく，単一の原子の(球対称＝ r 2 sin udrdudq で行う． q の積分で 2p ， u の 0 ～ p の積分の)局所ポテンシャル V alocal ( r )の空間全体でのフーリエ変は ， v ＝ － cosu ， v ＝ － 1 ～ 1 の 積 分で ， dv ＝ sinudu で あ換(39)でる．最終的に r の積分を行う．

f
f
f
f

V alocal(G• )＝ V alocal(r) exp [－iG• ・r… ]d…r＝

V alocal(r) exp [－i|G• |r cosu]r 2 sinudrdudq

原子の局所擬ポテンシャル V alocal( r)は， rc 外で－ e2Za/r ，rc 内で底上げされた形で(前章参照)，rc 外で－e2Za/r になる

適当な解析関数と rc 内の残差で表せる．式( 61 )は，解析関
数部分は一行目に戻って積分を解析式で行い，残差部分は最終形の r のメッシュの数値積分で行う(後述式( 83 ) , ( 84 )参

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




照)．

 節)．
が厳密に処理される(第 6 章

ところで，式( 59)のフーリエ成分 VH (G• )は，その形から• )の G• ＝0
も G• ＝0 の項が発散成分を持つ．式(60)の V PS
local( G項も同様である．これは，V alocal(r)のクーロンポテンシャル
形のために式( 61 )の V alocal ( G• )の G• ＝ 0 項が発散成分を持つ• )の G• ＝ 0 項は，式ためである．これら VH ( G• )と V PS
local ( G)にだけ出現するが，(57)のハミルトニアンの対角項(G• ＝G• ′簡単のためゼロとして扱われる．固有値全体の値のシフトに対応するが問題は生じない(固有値の絶対値が不定性を持つことになる)．しかし，次章で扱う全エネルギーでは，• )の G• ＝ 0 項における発散成分と残留成分
VH ( G• )と V PS
local ( G

平面波基底でのハミルトニアン非局所擬ポテンシャ



ル項式(53)に戻って，非局所擬ポテンシャル V PS
∑V aNL
NL(r… )＝∑R• a(r… －t• a－R• )の平面波基底間の行列要素を考える．式(33)から＋l〈l|で，射影演算子|l〉＝ ∑ |Ylm(u,q)〉
V aNL は∑|l〉DV NL
a,l ( r)m＝－llにより， t• a ＋ R• の位置の a 原子の周りで平面波の各 l 波成分を抜き出して DV NL
a,l (r)を作用させる．式(57), (60)を参考に• ′〈 •k＋G• |V PS〉
NL(r… )| •k＋G

＝
＝
＝

f
∑∑f
∑f

1
Q

Q

)・r… ]dr…
exp [－i( •k＋G• )・r… ]∑∑V aNL(r… －t• a－R• )exp [i( •k＋G• ′1
Q R• a

Q

N
Qa

Qc

R• a

)・r… ]d…r
exp [－i( •k＋G• )・r… ]V aNL(r… －t• a－R• ) exp [i( •k＋G• ′)・r… ]d…r
exp [－i( •k＋G• )・r… ]V aNL(r… －t• a) exp [i( •k＋G• ′

f
f
∑f

1
＝ ∑exp [i(G• ′－G• )・t• a] exp [－i( •k＋G• )・r… ]V aNL(r… ) exp [i( •k＋G• ′)・r… ]d…r
Qc a
Qc1＝ ∑exp [i(G• ′－G• )・t• a] exp [－i( •k＋G• )・r… ]∑|l〉DV NL〈l|exp [i( •k＋G• ′)・r… ]d…r
a, l ( r )
Qc a
Qc
l
1
＝ ∑exp [i(G• ′－G• )・t• a]
Qc a
l

Qc

〈l|exp [i( •k＋G• ′)・r… ]d…r
exp [－i( •k＋G• )・r… ]|l〉DV NL
a,l (r)と展開できる．3 行目から 4 行目の変形は，各単位胞での和が一つの単位胞での各原子の作用で代表できるので，単位胞総数 N を掛け，積分範囲を単位胞内としている．非局所擬ポテンシャルは，原子の半径 rc 内のみで作用する( r ＞ rc で分と関わらない)．方位座標 râ の積分を先に実行する．射影＋l演算子| l〉＝ ∑ | Ylm (u,q )〉の作用は，以下の球面調和関数
m＝－lの規格直交性を用いる(u は 0～p，q は 0～2p の積分)(40)(41)．

DV NLa,l (r)＝0)． 5 行目で実空間積分の原点を a 原子位置に移ff(u, q)sinududq＝dll′dmm′Ylm(u,q)Yl′
m′

している．

〈

＋l

exp [i( •k＋G• )・r… ]＝4p∑ ∑ il jl(| •k＋G• |r)Ylm(râ )Yl m( k ＋ G )l m＝－l(63)
l の和は l＝0, 1, 2…，m の和は l 毎に－l から＋l までである．

〈

jl は球ベッセル関数(実数の関数)，Ylm は球面調和関数，râ やk＋G は，r… や •k＋G• の方位座標(u, q)である．この展開式や)jl(| •k＋G• |r)jl(| •k＋G• ′|r)DV NLYl m( k ＋ G ′a,l (r)で，最終的に式(62)は以下のようになる．
• ′〈 •k＋G• |V PS〉
NL(r… )| •k＋G＝(4p)2
－G• )・t• a]×∑exp [i(G• ′
Qc a

f
rc

2
|r)DV NL
∑ jl(| •k＋G• |r)jl(| •k＋G• ′a,l (r)r dr

l

0

＋l

〈

出ており，電子構造計算関連で頻繁に使われる．

m＝－l

〈

特殊関数は，多くの応用数学・物理数学の本(40)(41)に詳しく＋l

式 ( 64 ) の 方 位 座 標 積 分 の 結 果 は ， ( 4p )2 ∑ Ylm ( k＋G )〈平面波の原点の周りの球関数展開(部分波展開)を利用する．

(65)

〈

式( 62 )の最後の行の実空間積分は，原子毎の半径 rc の球内積分である．射影演算子| l〉の平面波への作用は，次式の(62)× ∑ Ylm(k＋G)Y)lm(k＋G′m＝－l(66)

単位胞内 a 原子の周囲で平面波を部分波展開し，l 波成分に(râ )Yl′m′(k＋G)|l〉DV NL〈 l|＝4p (－i)l′jl′(| •k＋G• |r)Yl′m′a,l (r)l′m′

値的に行う．

∑〈 〈式(62)の最後の被積分関数は式(63)の展開を使うと〈l|exp [i( •k＋G• ′)・r… ]
exp [－i( •k＋G• )・r… ]|l〉DV NL
a,l (r)|r)Yl″m″(râ )Yl″m″(k＋G′)
4p ∑ il″jl″(| •k＋G• ′DV NLa,l を作用させている． r の一次元積分はメッシュ点で数

(64)

l″m″

この式(66)の形(非分離形)では，予め •k＋G• と •k＋G• ′の異なる G• の全組み合わせ( NG × NG 個)で積分を実行するため手間がかかる．そこで， Kleinman Bylander の分離形が用

方位座標 râ ＝( u, q )と動径座標 r に分けて行う( k＋G 等は積いられる(36)．これは，式(33)の射影演算子を用いる非局所〈である．この式の rc 球内極座標積分( d…r ＝ r 2 sinudrdudq )を




講義ノート

形において，同じ DV NL
a,l (r)成分で以下の表現を使う．
PS
PS|DV NL〈DV NLa,l (r)Ra,l (r)Ylm〉
a,l (r)Ra,l (r)Ylm|
PS
NL
PS〈Ra,l (r)Ylm|DV a,l (r)|Ra,l (r)Ylm〉
l m＝－l＋l

V aNL(r… )＝∑ ∑PS
PS|DV NL〈DV NLa,l (r)Ra,l (r)Ylm〉
a,l (r)Ra,l (r)Ylm|
(67)
C a,l
l m＝－l

＋l

＝∑ ∑

NL
PS
は，方位座標積分で
C a,l＝〈RPSa,l ( r ) Ylm | DV a,l ( r )| Ra,l ( r ) Ylm 〉

fの Ylm の規格直交性から C a,l ＝

rc

0

PS
2 2
DV NL
a,l ( r ) Ra,l ( r ) r dr とな.

第一原理擬ポテンシャル法での全エネルギーと原
子に働く力




電子密度分布と全エネルギー

平面波基底で行列表示されたハミルトニアン〈 •k ＋ G• | H | •k
＋ G• ′〉の固有値問題を解いて，固有値 E •kn ，固有ベクトル{ C n•k＋G• }を求める(効率的解法は第 7 章で説明)．出力の固有る(原子球内積分)． RPSa,l ( r )は， a 原子の擬ポテンシャルの l

ベクトルから電子密度分布 r ( r… )を組み立てる．式( 19)のよ成分を構築するときの擬動径波動関数で，RPSa,l (r)Ylm(u,q)のうに占有状態の波動関数のノルムの二乗の重ね合わせであ 節参照)．原子近傍 rc
セットが擬原子軌道である(第 4 章る．式(43)をそのまま代入すると

内では，価電子の擬波動関数は，擬原子波動関数 RPSa,l ( r )Ylm ( u,q )のセット で展開される 挙動をする．そ れへの式

occ

r(r… )＝∑∑w •kin|ψ •kin(r… )|2nPSの展開の形で(67)(分離形)の作用は，|DV NL
a,l ( r ) Ra,l ( r ) Ylm 〉

•ki

occ

n
1
＝∑∑w •kin∑∑C n•ki－G• )・r… ]Q－
exp [i(G• ′＋G• C •ki＋G• ′cある．本来の非分離形の式(33)の作用も同様に|RPSa,l ( r ) Ylm 〉

n •kiG• G• ′の展開に DV NL
a,l ( r )を掛けた形で，分離形で扱っても良いこ
分離形の式( 67 )を使うと，式( 62 )の最終形の原子位置をl

尾に記したように，ここでは基底関数や固有関数の規格化の

f1
C a, lrc

0PS
2
jl(| •k＋G• |r)DV NL
a,l (r)Ra,l (r)r dr×

1 を使って範囲を結晶 Q から単位胞 Qc に変えているため Q－cいる．この場合，占有状態の和を単位胞当たりの電子数で勘定することになり，都合が良い．式( 72 )は， G• の和の二重

f
rc

0

PS
2
|r)DV NLjl(| •k＋G• ′a,l (r)Ra,l (r)r dr

〈

〈

＋l

× ∑ Ylm(k＋G)Y
)l m ( k＋ G ′

(68)

m＝－l

となる．演算子の積分は，左右別個に極座標で r を含めて行
える．r 2 は極座標積分による．結局，式(62)の行列要素は，• ′〉〈 •k＋G• |V PS
NL(r… )| •k＋G＋l1• )Aa,lm( •k＋G• ′)＝∑∑C －a,l ∑ Aa,lm( •k＋G

(69)

m＝－l

〈

a l

1/2Y(k＋G)×
Aa,lm( •k＋G• )＝exp [iG• ・t• a]4pQ－
c
lm

f
rc

0

以外の領域の S •ki の固有ベクトルは式(51)で構築できる．
1 は本来は Q－1 であるが，第 5 章
 節の末式( 72 )中の Q－
c原点とした積分部分は，式(63)～(65)から(4p)2∑(72)となる．既約領域の •ki 点で固有ベクトルを求めた後，それ

とになる．

PS
2
jl(| •k＋G• |r)DV NL
a,l (r)Ra,l (r)r dr

ループを含み，NG×NG のオーダーの手間がかかる(NG は平面波基底数)．この問題は高速フーリエ変換で解決できる(第
 節)．
6 章
 節で説明したように，ハミルトニアンの固有
第 3 章

値・固有ベクトルを一度解いただけではだめで， SCF ルー
プで自己無撞着な{ C n•k＋G• }， r ( r… )を求める．その後，式( 1 )の全エネルギー Etot を計算する．以下，Etot の各項を順に説明する．単位胞(スーパーセル)当たりの全エネルギーとする．

(70)

式( 1 )の第一項の運動エネルギー T は，式( 3 )より求める．式(19)，(72)のように占有準位の重み付きの和としてとなる．a の和は単位胞内の原子のみである．この分離形の場合， NG 個のオーダーの Aa,lm ( •k ＋ G• )だけを予め計算して

occ

n •ki

メモリーすればよい．非分離形の式( 66 )の NG × NG のオー

2

Q

ψ
k• in(r… ) －

2

occ

•kin

• |－
＝∑∑w •kin∑∑C •kni＋
G• C n•ki＋G• 〈′ •ki＋G

ダーに比べ，計算負荷が NG のオーダーに低減化する．
これまでの各項をまとめると，平面波基底でのハミルトニn •kiG• G• ′occ

2
|C n•ki＋G• |2| •ki＋G• |2G• 2m＝∑∑w •kin∑

アンの行列要素は以下となる．
〉＝〈 •k＋G• |H| •k＋G• ′
f ( 2m: )ψ (r… )d…r

T＝∑∑w •kin

n •ki

2 2
〉
: | •ki＋G• ′2m(73)である．式( 52 )を利用している．波動関数は体積 Q で規格2• －G• ′
| •k＋G• |2dG• ,G• ′＋V PS)
local(G2m

化し，Q で積分している(単位胞 Qc で規格化した場合も結果• ′＋〈 •k＋G• |V PS〉＋VH(G• －G• ′)＋mxc(G• －G• ′)NL(r… )| •k＋Gは同じ)．占有状態の和が単位胞当たりの電子数に相当する

右辺の第一項は式( 52)，第二項は式(60)，( 61)，第三項は状態数を勘定に入れれば問題ない．式( 51 )から| C nS •k＋G• |2 ＝
|C n•k＋S－1G• |2 で，|S •k＋G• |2＝| •k＋S－1G• |2 から，式(73)の •ki 点式(69)，( 70)(分離形)，第四項は式( 59)を参照．第五項はの和は，ブリルアンゾーン内の既約領域に限ってよい．

(71)

mxc(r… )のフーリエ成分．上述のように行列の対角項に出てく• ＝ 0)， VH (G• ＝ 0)は，発散項のためゼロとして扱る V PS
local ( G

は，式( 53 )から周期系の擬ポテンシャルの局所項と非局所う．

項で別々に計算する．局所項の寄与 EL は， V PS
local ( r… )のフー式( 1 )の 第二 項 の 原子 か ら のポ テ ン シャ ル エネ ル ギ ー

リエ成分の式(60)，r(r… )のフーリエ成分を使い，単位胞当たりでま て り あ 第61巻 第11号(2022)
Materia Japanまれる．その処理は次節で説明する．

f

EL ＝

Qc

V PS
local(r… )r(r… )d…r

式( 1 )の第四項の交換相関エネルギー Exc は，式( 12 )のLDA の場合として， r (r… )を使って求める交換相関エネルギ
ー密度 exc ( r… )のフーリエ成分を exc ( G• )として，式( 74 )，f• )r(G• ′＝∑∑V PS)
local(GG• G• ′

)・r… ]d…r
exp [i(G• ＋G• ′Qc• )r(－G• )＝Qc∑V PS
local(G(74)

G•• )の G•となる．式( 40)を使っている．前述のように V PS
local ( G

＝0 項には発散成分が含まれる．その処理は次節で説明する．

(76)と同様に次式で与えられる．
Exc＝Qc∑exc(G• )r(－G• )

(77)

G•

以上のように，式( 1 )の Etot は，T，EL，ENL，EH，Exc，

非局所項 V PS
NL ( r… )によるポテンシャルエネルギー ENL は，

EI－I の和で与えられ， EI－I (正イオン間静電相互作用)以外

射影演算子を含むので，占有された波動関数に作用させる形の項の詳細を説明した． T と ENL は，占有状態の固有ベクトルを含み，平面波基底の G• に関する和と •k についての積で，運動エネルギー項式(73)に近い表現となる．
occ

n •ki

分(メッシュ点での重み付きの和)がある． •k 点を抽出する領f

ENL＝∑∑w •kin

Q

PS
ψ
k• in(r… )V NL(r… )ψ •kin(r… )d…r

域は系の対称性を利用して既約領域でよい． EL ， EH ， Exc
は，フーリエ級数展開に関する G• の和である．EL，EH の G•の和は， G• ＝ 0 項が発散成分を持つ(その処理は次節で説明

occn• |V PS• ′＝∑∑w •kin∑∑C n•ki
〉＋G• C •ki＋G• 〈′ •ki＋G
NL(r… )| •ki＋Gn •ki

G• G• ′occ＋lする)．なお， USPP 法や PAW 法でも Etot の各項が基本的

n •ki

m＝－l

に同様の表式で表されるが，特に ENL については，様子の

n
－1
• )
＝∑∑w •kin∑∑C n•ki
＋G• C •ki＋G• ∑′ ∑C a,l ∑ Aa,lm( •ki＋GG• G• ′al)
×Aa,lm( •ki＋G• ′かなり違う複雑な表式となる(28)(30)．
＋l

occ

1
＝∑∑w •kin∑∑C －∑ |∑C n•ki＋G• Aa,lm( •ki＋G• )|2
a, l

n •ki

a

(75)




Ewald 法と発散項の処理

m＝－l G•l

式(67)，(69)の分離形を用いた．Aa,lm( •ki＋G• )は，式(70)参式( 1 )の正イオン間クーロン相互作用 EI－I は， Ewald 法

照．式( 73 )の場合と同様，波動関数の積分は結晶全体であるが，原子の和は単位胞内のみ，占有状態の和は，単位胞当で計算される．発散項が含まれるが，上述の EL，EH 内のフーリエ成分の G• ＝0 項に含まれる発散成分と打ち消しあうのたりの電子数分に取る．また，式( 73 )と同様に， •ki 点の和で，Etot 全体では問題ない(42)(43)．本節ではそれを説明する．

は既約領域内のものに限ってよい．これは，式( 51 )や式(70)の Aa,lm( •ki＋G• )内の Ylm やjl の性質から証明される．

単位胞当たりの正イオン間クーロン相互作用の総和は式( 1 )の第三項の静電エネルギー EH は，EL と同様に
EH ＝

1
2

f
Qc

{

e2
Z 2a
Z aZ a ′＋∑∑EI－I＝ ∑ ∑2 a R• ≠0 |R• | a′≠a R• |R• ＋t• a′－t• a|
e2

}

e2

1
1
＋ ∑ ∑ ZaZa∑＝ ∑Z 2a ∑′• ＋r… aa′|
2 a R• ≠0|R• | 2 a a′≠a
•R |R

VH(r… )r(r… )d…rである．a，a′は単位胞内原子，Za，Za′は正イオンのイオンf＝1) exp [i(G• ＋G• ′)・r… ]d…r∑∑VH(G• )r(G• ′2 G• G• ′Qc＝

1Qc∑VH(G• )r(－G• )＝2pe2Qc∑|r(G• )|2/|G• |2
2 G•
G•

(78)

価， r… aa′＝ t• a′－ t• a である．周期系での同値原子間と同値でな(76)

い原子間の和で，単位胞を超えた距離は R• を加えて表現さ1

1

れる． ∑ |R• | と∑|R• ＋r… | の R• の和は，収束が遅い．次式の
R• ≠0

R•

aa′

となる．最後の形は，式( 59 ) VH ( G• )＝ 4pe2r ( G• )/| G• |2 と実関数 r(r… )のフーリエ成分間の関係式 r(－ G• )＝ r(G• )を使っように Ewald 法で実格子と逆格子の和に変換すれば，比較的短範囲の R• ， G• の和で収束する(変換の証明は，本講義ノている．前述のように VH ( G• )の G• ＝ 0 項には発散成分が含ート最終回末尾の補論 A 参照)．

1∑ • ＝∑| R• ≠0R• ≠0 |R

erfc[|R• |g]p
exp [－|G• |2/4g2] 2g
p
4p
＋ 2 ∑
－ － 2 ＋ lim |G• |－2|R• ||G• |2/4g2
g Qc G• ≠0
p g Qc Qc G• →0(79)

1
erfc[|R• ＋r… aa′|g]p
exp [－|G• |2/4g2]＝∑＋ 2 ∑ exp [－iG• ・r… aa′]| R•
|R• ＋r… aa′||G• |2/4g2
g Qc G• ≠0R• |R＋r… aa′

∑ •

p
4p
－ 2 ＋ lim |G• |－2
g Qc Qc G• →0
g は，収束が速くなるように適当に選ぶパラメータ．erfc は補誤差関数 erfc[x]＝

f

2
p

∞

e－t dt である．もともと式(80)の
2

x

(80)

分，それらの項を R• → 0 ， G• → 0 とした極限の残留項，発散項が各式に追加されている．式( 79 )の第三，四，五項，式( 80 )の第三，四項である(導出は補論 A 参照)．式( 79 )，導出の後，r… aa′→0 の極限として式(79)が導出される．
式(79)では R• ，G• の和が，式(80)では G• の和が，R• ＝0，(80)を式(78)の EI－I に入れると 2∑Z 2a や 2∑ ∑ ZaZa′を掛けG• ＝ 0 の項を除いている．発散成分があるからである．そのた形の和になり，残留項，発散項もまとめて記すと次式になる．




e2

e2

a

≠a
a a′

講義ノート

EI－I＝

(

e2 2
∑Z a
2a

＋
－

e2

)(

erfc[|R• |g]p
exp [－|G• |2/4g2]
＋ 2 ∑
|R• |
|G• |2/4g2
g Qc G• ≠0
R• ≠0

)

∑

(

erfc[|R• ＋r… aa′|g]p
exp [－|G• |2/4g2]＋ 2 ∑ exp [－iG• ・r… aa′]|R• ＋r… aa′||G• |2/4g2
g Qc G• ≠0R•∑ ∑ Z aZ a′ ∑2 a a′≠a)e 2g

pe2N 2c 2pe2N 2c
Z 2a－ 2 ＋
lim |G• |－2pa
2g Q cQc G• →0∑(81)

三行目が残留項と発散項で，∑Z 2a ＋∑ ∑ ZaZa′＝∑∑ZaZa′＝
a

a a′
≠a

f

一行目に G• ＝0 を入れた V alocal(r)d…r の形の極座標積分(dr… ＝a a′

N 2c ，∑Za＝Nc の関係を使っている．Nc は単位胞当たりのイ

r 2sinudrdudq)でa

オン価総数(価電子総数)．正負のイオンを含む古典系では
Nc ＝ 0 なので発散項

2pe2N 2c

Qc

f

1
∑4p
Qc a

lim | G• |－2 は消えるが，第一原理計算の EI－I は正イオンのみで，発散項

Qc

0

1
(V alocal(r)－(－e2Za/r))r 2dr＝ ∑aa
Qc a

(83)

となる． aa は残差の数値積分， 4p は u ， q の積分からのもG• →0

2pe2N 2crcの．長距離形の寄与は，式(61)の一行目に入れて G• ≠0 のま

lim | G• |－2 はま変形するとG• →0

EI－I 内に残る．

f(－e2Za/r)exp [－iG• ・r… ]d…r＝－

一方，前述のように式( 74 )，( 76 )の EL ， EH のフーリエ4p e 2Z a|G• |2(84)

成分の積の和のうち， G• ＝0 項が発散成分を含む．これらのである．これは，クーロンポテンシャルのフーリエ変換とし発散項が式( 81 )の末尾の発散項と打ち消しあうことを示て有名な式である(証明は本講義ノート最終回末尾の補論

1• )r(－G• )}と Qc lim {VH(G• )r
す．具体的には Qc lim {V PSlocal( G
2 G• →0G• →0B)．式(84)を式(60)の最終形に入れ，式(83)も考えると，• )は，以下の残留項と発散項となる．
lim V PS
local(GG• →0(－G• )}を検討する．r(－G• )は小さい G• で r(－G• )＝Nc/Qc＋
b|G• |2＋ 高次項と展開できる(42)．Nc/Qc は価電子の平均濃度，

b は平均からの局所的なズレの効果に関する係数．

11
4 p e 2Z a• )＝ ∑aa－ lim ∑lim V PS
local(GQc a
Qc G• →0 a |G• |2G• →01
Qc lim
2 G• →0

1
4 p e 2N c
＝ ∑aa－
lim |G• |－2
Qc a
Qc G• →0

{ VH ( G• ) r (－ G• )}の方から扱うと，式( 59 ) VH ( G• )＝ 4pe2r(G• )/|G• |2 と r(－G• )＝r(G• )から
1
1
4pe2|r(G• )|2
Qc lim {VH(G• )r(－G• )}＝ Qc lim|G• |2
2
2
G• →0
G• →0

{

＝

Za＝Nc を使っている．式(85)を r(－G• )＝Nc/Qc＋b|G• |2＋∑a

高次項の展開と合わせて代入すると}• )r(－G• )}＝－
Qc lim {V PSlocal(G2pe2N 2c
lim |G• |－2＋4pe2Nc b
Qc G• →0

G• →0

(82)

＋

となる．|r(G• )|2N 2c /Q2c ＋2Nc b|G• |2/Qc を使っている．
• )r(－ G• )}のうち lim V PS• )を考え一方， Qc lim {V PSlocal ( Glocal ( G
G• →0

G• →0

4pe2N2c
lim |G• |－2－4pe2Nc b
Qc G• →0
Nc
∑aa
Qc a

以上から，式( 81 )の EI－I 内発散項

G• →0

• )＝
ると，式( 60 )から lim V PS
local ( G(85)

1∑lim V alocal ( G• )である．
Qc a G• →0(86)

2pe2N 2c
lim | G• |－2 と，
Qc G• →0

式(82)，(86)の EH，EL の G• ＝0 項内の発散成分

2pe2N 2clim
Qc G• →0

4pe2N 2c
lim | G• |－2 の総和がゼロとなり，消える．
Qc G• →0lim V alocal ( G• ) の 発 散 項 は ， 原 子 の 局 所 擬 ポ テ ン シ ャ ル| G• |－2，－

V alocal ( r ) が rc 外 で ク ー ロ ン 形 を 持 つ こ と に 起 因 す る ．

同様に EH ， EL の G• ＝ 0 項内残留成分の和は， b の項が消え

V alocal ( r ) を 長 距 離 形 ( － e2Za / r ) と 短 範 囲 の 残 差 の 寄 与

cて， Q∑
aa のみになる．

G• →0

N

c a

(V alocal(r)－(－e2Za/r))に分けて式(61)の積分 V alocal(G• )を考

これらの残留項，発散項に留意しながら EI－I と Etot をまと め て 記 す ． 式 ( 81 ) の EI－I か ら 発 散 項 を 除 い た も の を• )への寄与は，式(61)のえる．短範囲形(残差)の lim V PS
local(GEEwald とすると，次のようになる．

G• →0

•

•

G| /4g ](e2∑Z )( ∑ erfc|[|R• R| |g]＋gpQ ∑ exp [－|)－e pg∑Z －p2eg NQ|G• | /4geerfc[|R• ＋r… |g]pexp [－|G• | /4g ]＋＋ ∑∑ Z Z ∑exp [－iG• ・r… ]∑()|R• ＋r… ||G• | /4g
gQ
2
2

EEwald＝

2
a

a

2

R• ≠0

2

2

• ≠0
cG

2

2

2

2

2

2
a

a

2

≠a
a a′

a′

R•

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

aa′

c

2

aa′

a

2
c

2

c G• ≠0

aa′

2

2

2

(87)




occ

2
N
• )r(－G• )＋ c∑aa
|C n•ki＋G• |2| •ki＋G• |2＋Qc ∑ V PS
local(G
Qc a
G• 2m
G• ≠0

Etot＝∑∑w •kin∑
n •ki

occ

＋l

1
＋∑∑w •kin∑∑C－∑ |∑C n•ki＋G• Aa,lm( •ki＋G• )|2
a,l

n •ki

a

m＝－l G•l＋2pe2Qc ∑ |r(G• )|2/|G• |2＋Qc∑exc(G• )r(－G• )＋EEwald
G• ≠0

(88)

G•

式(88)は，式(73)～(77)を使っている．以上の EI－I と EH，( 38 )，( 39 )のフーリエ変換の演算は，高速フーリエ変換

EL に 関 す る 発 散 成 分 ， 残 留 成 分 の 問 題 は ， USPP 法 や(FFT)を用いて行う．

PAW 法でも共通であり，同様の処理がなされ，式(83)に基FFT は，離散フーリエ変換である． …a1 ， …a2 ， …a3 で形成さNc
づく Q ∑
aa が付加される．なお， Etot の値は，イオン(裸のれる平行六面体の単位胞内で， …a1， …a2， …a3 方向に各々 M1 ，

擬ポテンシャル)と価電子がばらばらに無限遠に遠ざかった

M2 ， M3 分割で， M ＝ M1M2M3 個の実空間メッシュ点 r… n ( r… n静的状態を原点とし，電気的に中性の系なら必ず負の値であ＝ M1 …a1 ＋ M2 …a2 ＋ M3 …a3 ， ni は整数， 0  ni  Mi － 1 ， i ＝ 1 ～ 3 )る．

で，格子周期関数 f ( r… )の単位胞内データが与えられる．

c a




n

n

1

n

2

3

FFT の順変換式(89)では，M 個の実空間データ{f(r… n)}のセットが，同数の M 個の逆格子点 G• m(G• m＝m1 •b1＋m2 •b2＋m3 •b3，

高速フーリエ変換の活用(その 1)

mi は 整 数 ， 0  mi  Mi － 1 ， i ＝ 1 ～ 3 ) 上 の デ ー タ セ ッ ト{f(G• m)}に変換される．

平面波基底の第一原理擬ポテンシャル法では，上述のように格子周期関数が逆格子ベクトル・フーリエ級数展開で表さ
れ ， ハ ミ ル ト ニ ア ン や Etot の 計 算 が 効 率 化 さ れ る ． 式

f(G• m)＝

これは，式(39)のM 1－ 1 M 2－ 1 M 3－ 1
1∑
∑ ∑ f(r… n) exp
M 1M 2M 3 n 1 ＝ 0 n 2 ＝ 0 n 3 ＝ 0[－i 2pMn m ]exp [－i 2pMn m ]exp [－i 2pMn m ]
1

1

2

1

2

3

2

3

(89)

3

f

1
Qc

Qc

f(r… )exp [－iG• m・r… ]d…r の単位胞内の r… 積分を等間隔メッシュ点 r… n の台形公式で行い，ai・ •bj＝2pdij から

exp [－iG• m・r… n]＝exp[－i 2pMn m ]exp [－i 2pMn m ]exp [－i 2pMn m ]
1

1

2

1

2

3

2

3

(90)

3

を用いている．一方，逆変換(逆 FFT)は次式で計算される．
M 1－ 1 M 2－ 1 M3 － 1

f(r… n)＝∑ ∑ ∑ f(G• m)exp
m 1 ＝ 0 m 2 ＝ 0 m3 ＝ 0[i 2pMn m ]]exp [i 2pMn m ]exp [i 2pMn m ]
1

1

2

1

2

2

3

3

(91)

3

M

ば， G• m＝M1 •b1 の exp [iG• m・ r… ]の波長は， …a1， …a2， …a3 の平行m＝ 1六面体の単位胞の …a2， …a3 で作られる面と面の間を …a1 に沿っ実空間のメッシュ点の離散データ{ f ( r… n )}は，格子周期で上の大きな G• m (メッシュ間隔より短い波長の G• m )は不要で単位胞内のものが繰り返すわけだが，逆空間のデータ{ f ( G• m )}も f ( G• m )＝ f ( G• m ± Mi •bi )で，逆格子点の大きな平行ある．また，上述のように離散メッシュ点上の位相が同じに式( 38 )を M 個の G• m 点の和 ∑ f ( G• m ) exp [ iG• m ・ r… n ]で行うものである．

て M1 分割したメッシュ間隔に相当し， •b1 方向では，これ以六面体(M1 •b1，M2 •b2，M3 •b3 で形成)でサイクリックに同じ値なるので(式( 90 ))，図 9 の原点を中心とする大きな平行六面体内の G• mの展開として扱える．

が繰り返すことが，式( 89 )の形からわかる．± Mi •bi は式

 節で触れたように，電子密度分布 r ( r… )な
一方，第 5 章

( 90 )で m1, m2, m3 に各々± M1, ± M2, ± M3 の効果で， exp
内で 2pi ×整数の効果なので不変で， f ( G• m )＝ f ( G• m ± Mi •bi )ど格子周期関数の単位胞内の変化の様子に対応して，どこまが 証 明 さ れ る ． こ う し て ， FFT の 逆 格 子 空 間 デ ー タ{f(G• m)}は，±Mi •bi 等のシフトで，原点を中心とする平行六面体に移し，その中の G• m と f(G• m)を扱うと考えることがででの細かい r… n メッシュを用いるか(どこまでの大きさの(短波長の)G• m を用いるか)，M1, M2, M3 で設定される．平面波
基底は，式( 47 )のように，平面波のカットオフエネルギー2
Ecut＝ G2max で，最大の G• の大きさ Gmax が決まる．VH, mxc
2m

また，実空間データが離散的なので，実空間メッシュよりを含む Vlocal の平面波基底間の行列要素〈 •k ＋G• |Vlocal(r… )| •k ＋〉は，式(57)のように Vlocal(G• －G• ′)だが，平面波基底の G•
G• ′短い波長の G• ( exp [ iG• ・ r… ]の波長は 2p /| G• |)は意味がないの 範 囲 を 考 え る と ， G• － G• ′は ， 原 点 を 中 心 に 2Gmax とため，大きな平行六面体(M1 •b1，M2 •b2，M3 •b3 で形成)の M＝
M1M2M3 個の G• m 点のフーリエ成分だけを扱えばよい．例え• ), VH(G• ), mxc(G• ))が，この G• －G• ′のする Vlocal(G• )(V PS
local( G

きる(図)．



－2Gmax に渡って変化することになる．フーリエ変換で用意

講義ノート

範囲をカバーせねばならないので，原点を中心に半径 2Gmaxの球内の G• 点が，上記の FFT の M1 •b1，M2 •b2，M3 •b3 で構成
される大きな平行六面体(原点が中心に来るようシフトしたもの)に収容される必要がある(図 9)．こうして，大きな平行六面体の二面間隔が 4Gmax 以上ということが，Ecut からのFFT の M1，M2，M3 の条件である．
 節で説明したように Ecut は Etot の収束の様子か
第 5 章

ら決まる(図 8)．主に擬ポテンシャルの形状が支配因子である．通常，Ecut を与えて Gmax が決まると，FFT の M1，M2，
M3 (実空間メッシュと G• m の範囲)は自動的に上記条件を満たすように与えられる．
 高速 フーリエ変換 の活用(そ の(紙数の関係で ，第 6 章
 原子に働く力は，次回にまわします)2)，第 6 章
(次号へつづく)

文

図9

FFT で用意される G• m 点とフーリエ成分 f ( G• m )の値の周期性．単位胞内 M1 ， M2 ， M3 メッシュの実空間デー
タに対し，フーリエ成分 f ( G• m )は G• m ＝ m1 •b1 ＋ m2 •b2 ＋
m3 •b3(0miMi－1，i＝1～3)で与えられる．その値には f ( G• m )＝ f ( G• m ± Mi •bi )の関係があり， M1 •b1 ， M2 •b2 ，
M3 •b3 で形成される逆格子点の並びの大きな平行六面体
で f ( G• m )の値が周期性を持つ． G• m ， G• m － M1 •b1 ， G• m －M2 •b2 等のフーリエ成分は同じ値で，太線のように平行六面体を原点の周りにシフトできる．平面波のカット
オフエネルギー Ecut からの平面波基底の最大の G• ベク
トルの大きさ Gmax について，ハミルトニアンの行列要の範囲は－2Gmax～
素で必要になるフーリエ成分 G• －G• ′＋ 2Gmax になる． M1 •b1 ， M2 •b2 ， M3 •b3 で形成される平行六面体が半径 2Gmax の球を収容するよう M1 ， M2 ， M3
を決める．

献

(39) 今村 勤物理とフーリエ変換，岩波書店，(1976).
(40) 森口繁一，宇田川 久，一松 信岩波 数学公式～，岩波書店，(1987).
( 41 ) ジョージ・アルフケン，ハンス・ウェーバー基礎物理数学
第 4 版 Vol. 1～4，講談社，(1999).
(42) J. Ihm, A. Zunger and M.L. Cohen: J. Phys. C, 12(1979),
44094422.
(43) M.T. Yin and M.L. Cohen: Phys. Rev. B, 26(1982), 3259
3272.
(
全ての文献番号は，第一回からの通しとなっております)★★★★★★★★★★★★★★★★★★★★★

香山正憲

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

ま て り あ 第61巻 第11号(2022)
Materia Japan