ESopt 使用の手引き
ver 1.2.1

平成２３年３月

大阪大学 大学院基礎工学研究科
物質創成専攻 未来物質領域 新物質創製講座

1

ESopt は東京大学物性研究所において開発が進められた平面波展開法による第一原理電子状態計算プログラム (opt と呼ぶ) を元にして内容を整理したプログラム・
パッケージである。物性研からは ISSP FPSOPT として公開版が提供されているが、
ESopt は機能を限定した上で Fortran 95 による修正と対角化ルーチンの追加を行っていることから、その派生版であると言える。
旧 opt の特徴として、以下の２点が挙げられる。
• ソースコードの可読性の高さ• 波動関数のもつ全自由度を同時に動かす CG 法の実現
前者の点は ESopt においても引き継がれ、Fortran 77 によるメモリ・アロケーショ
ンの困難を回避するために生じた可読性の低下の問題を避ける方向で修正がなされている。また、後者は開発者の荻津格博士 (現米国ローレンス・リバモア研究所) によるものである。比較的多くの CG 法で採用されているバンドごとに最適化を行う
方法と異なり、opt における CG 法では波動関数に関する全自由度の同時最適化が実現されている。この CG 法がもつ本来の形式通りのインプリメンテーションが実現されているのは、稀な例である。
この手引きの目的は、ESopt の基礎としている方程式系と実際のプログラムでの実現例を比較することで、第一原理電子状態計算のプログラムを容易に理解し、プ
ログラム動作を理解しながら計算を実行する機会を、平面波基底展開法による計算の初心者に対して提供することである。従って、基礎としている理論の完全な記述は目指さず、必要最低限のものに限りながら、一方で幾つかの具体例において、プ
ログラム上での実現例を紹介することを目指している。

ver. 1.0.0 原子位置最適化ができる Sopt 機能を含む f90 化第一版。
ver. 1.1.0 内部応力計算と結晶構造最適化ができる Allopt 機能を含む第一版の変更版。
ver. 1.2.0 バンド計算機能を含めた変更版。

1

全エネルギーの表式

この節では、平面波展開法における全エネルギーの表式を与える。但し、実際の第一原理計算における全エネルギーの具体的な計算方法を理解することを目的とするため、擬ポテンシャルとして、ノルム保存型の擬ポテンシャルを適用した場合について記述する。RI , ZI , V̂Ipseudo (r − RI ) を I 番目の原子のもつ原子座標、イオン
価数、擬ポテンシャルを表すものとする。
まず、体積 Ω をもつ結晶の全エネルギーの表式として以下のものを用いる。

Etotal = Ekin + Eel−el + Exc + Eel−ion + Eion−ion .

2(1)ここで、Ekin は電子の運動エネルギー、Eel−el は電子間相互作用のうち電子密度分布
間の古典的な Coulomb 相互作用の部分（以下ハートレー・エネルギーと呼ぶ）、Exc
は交換相関エネルギー、Eel−ion は電子とイオン殻の Coulomb 相互作用、Eion−ion は
イオン殻同士の Coulomb 反撥エネルギーである。ただし、イオン殻は球対称ポテン
シャルを与えるとした。
それぞれの項は、実空間において以下のように表される。

Ekin =
Eel−el =

(

∑ ∫

)

1
dr ϕ∗k,n,σ (r) − ∇2 ϕk,n,σ (r) ,2k,n,σ Ω

(2)n(r)n(r′ )1 ∫∫dr dr′,2 Ω|r − r′ |(3)

dr εxc (n(r))n(r) ,

(4)

∫

Exc =
Eel−ion =

Ω

∑ ∑∫

Ω

k,n,σ I

Eion−ion =dr ϕ∗k,n,σ (r)V̂Ipseudo (r − RI )ϕk,n,σ (r) ,1 ∑ ZI ZJ.
2 I,J |RI − RJ |

(5)
(6)

ここで、一電子密度 n(r) は、

n(r) =∑nσ (r) ,

(7)

σ=↑,↓

nσ (r) =∑|ϕk,n,σ (r)|2 ,(8)

εk,n,σ ≤EFと与えられる。ϕk,n,σ (r) はサンプル k 点である k において、スピン σ をもつ n 番目のバンドにおける Kohn-Sham 軌道の波動関数である。そのエネルギーは εk,n,σ であるとした。EF はこの仮想的一電子系におけるフェルミ・エネルギーである。

2

Dual space formalism

本パッケージでは、周期的境界条件のもとで、波動関数を平面波基底によって展開して与える。Bloch の定理と Fourier 変換を用いることによって、それぞれの波動関数は平面波の線型和によって以下のように表される。

ϕk,n,σ (r) =∑ϕk,n,σ (G) exp (i(k + G) · r) .

(9)

G

ここで、G の和は逆格子の上でとる。
実際の計算では、平面波基底展開特有の周期性を利用した高速化を図るため、分子や原子を計算しようとする場合でも、それを周期的に配置した構造を考えて、全体としての周期性を常に仮定した方法を用いる。この方法に特有の制限に関しても、折に触れて説明することにする。
説明を簡単にするため、周期的構造の基本をなす単位を、単位胞と呼び、N を大きいが有限のある整数として単位胞が N 3 個周期的に配列している構造を考える。つ3

まり、単位胞の単位で周期的構造であると同時に、全体としてそれが周期的境界条
件のもとにある、という意味で２重に周期性が現れる。単位胞の体積を表す Ωcell を導入する。即ち、Ω = N 3 Ωcell であると考える。
ここで、Fourier 変換について纏めておこう。後に、カットオフを導入することで、G の有限和が現れるが、まず、r が連続変数である範囲で変換式を定めておくと便利である。そこで、(9) 式で G の和は無限和であると考えている。電子密度は結晶の周期と同じ周期をもつと考えてよいことから、密度に関する Fourier 変換は、

nσ (G) =1 ∫drnσ (r) exp (−iG · r) ,Ωcell Ωcellnσ (r) =∑nσ (G) exp (iG · r) ,

(10)
(11)

G

で与えられる。
実際の計算では、実空間と逆格子空間（運動量空間）の両方を用いて計算が行われる。これは、交換相関エネルギーのように実空間で計算するのに適している積分と、運動エネルギーなどのように逆格子空間で計算する方が容易な積分の双方を得る必要があるからである。この方法は、“dual space formalism” と呼ばれている。この方法と高速 Fourier 変換 (FFT) を用いることによって、計算量を大きく減らすことができる。
Kohn-Sham 方程式の解は、逆格子空間において表現された波動関数 ϕk,n,σ (G) に関して与えられるエネルギー汎関数の最適化問題として解かれる。そこで、ある波動関数を与えたときに一電子密度を通して得られる各エネルギーを項ごとに与えるために、一電子密度の実空間および波数空間での表現を求めておく必要がある。その具体的な計算の手続きは以下の通りである。
数値的に表現するため、G は逆格子空間において有限な領域の内部でのみ考える
必要がある。より具体的には、展開平面波にあるエネルギーカットオフ Ecut を設ける。|G|2 > Ecut では波動関数の Fourier 成分である ϕk,n,σ (G) の値が 0 であるとするのである。このエネルギーカットオフの設定によって、必要とされる逆格子空間の部分空間が定まる。つまり |G|2 ≤ Ecut となる格子点全てがぎりぎり含まれる有限な大きさの逆格子メッシュを定める。
（図 1）すると対応して、離散 Fourier 変換を考えた際に実空間上に格子構造が与えられる。この格子を実空間メッシュと呼ぶ。メッ
シュ上の点を以下では rj とする。すると、以下のような計算手順が与えられる。

1. 逆格子空間における波動関数 ϕk,n,σ (G) を FFT によって Fourier 逆変換し、(9)式により与えられる実空間メッシュ上の波動関数 ϕk,n,σ (rj ) を得る。
2. ϕk,n,σ (rj ) から、実空間メッシュ上における電子密度 nσ (ri ) を (7) 式によって計算する。
3. 実空間メッシュ上の電子密度 nσ (ri ) から、FFT による Fourier 変換によって電子密度の Fourier 級数を得る。

4

図 1: FFT による Fourier 変換で移り変わる、実空間メッシュと逆格子空間メッシュ√の対応関係の概念図。逆格子空間（運動量空間）では、G ≤ Ecut の範囲にある G
点のみで波動関数が展開される。
最後のステップでは、次式を用いている。

nσ (G) =∑nσ (rj ) exp (−iG · rj ) .

(12)

j

3

各エネルギーの表式

前節の方法で得られた {ϕk,n,σ (G)}, {ϕk,n,σ (rj )}, {nσ (rj )}, {nσ (G)} を用いて以下のように具体的な計算が行われる。以下のエネルギーの表式は、すべて単位胞当たりのエネルギーである。

3.1

運動エネルギー

この項は、逆格子空間で計算される。単位胞当たりの運動エネルギーを Ēkin とすることで、Ēkin =∑ ∑ ∫k,n,σ G,G′

=

Ωcell

dr

|k + G|2
exp (i(G − G′ ) · r) ϕ∗k,n,σ (G′ )ϕk,n,σ (G)2Ωcell ∑ ∑|k + G|2 |ϕk,n,σ (G)|2 .
2 k,n,σ G(13)エネルギーカットオフを導入した段階で、G の和を有限和に書き直すことで、この表式を用いてよい。この計算は、実空間で差分近似などによって評価しても計算量が大きく増加することはないが、指数関数に対しては微分演算の結果が数の積で表現できることから、波数空間で行う方が取り扱いが簡単である。

5

3.2

ハートレー・エネルギー

この項も、逆格子空間で計算される。同じく、単位胞当たりのハートレー・エネ
ルギーを Ēel−el とする。

1
Eel−elN3∫n(r)n(r′ )
1∫
dr
dr′
=
2 Ωcell|r − r′ |
Ω
∫
∫
1
1 ∑
=
nσ (G)nσ (G′ ) exp (iG · r + iG′ · (r − s))drds
2 Ωcell|s| G,G′ ,σ,σ′ΩĒel−el =∫Ωcell ∑1=nσ (G)nσ (−G) ds exp (iG · s) .
2 G,σ,σ′|s|Ω

ここで、G → 0 で生じる発散は、Eion−ion , Eel−ion の各項にも現れ、系が中性の時には完全に打ち消し合うことが分かっている。そのため、各項での発散部分を予め除いておく。Ēel−el では、G = 0 の項が対応するため、この項を除くことにする。さらに、N → ∞ と考えて G ̸= 0 のときに、
∫

ds
Ω

4π
1
exp (iG · s) = 2 ,|s|Gであることを用いて、また n(r) が実であることから nσ (−G) = n∗σ (G) が導かれる∑ことを用いると、n(G) = σ nσ (G) と記載することにして、次式を得る。

Ēel−el =Ωcell ∑′ 4π|n(G)|2
.
2 G
G2

(14)

∑

ここで、 ′G は G = 0 を除くことを示している。
この計算は、実空間で行うと 6 次元空間積分となるため、実空間メッシュの格子2
点数 Nmesh に対して、O(Nmesh) の演算が必要となる。それに対して、逆格子空間
で行った場合には、FFT の計算に必要となる演算量 O(Nmesh log Nmesh ) が全体としての演算量を決める。そこで、この計算は逆格子空間で行った方が容易であると言
える。

3.3

交換相関エネルギー

交換相関エネルギー密度 εxc にはいくつかの表式があるが、本パッケージで提供しているのは、Perdew-Zunger の表式 (PZ81) である。この表式はその後の PerdewWang の表式 (PW92) と比較した場合には、両者とも基本的に Ceperley-Alder の量子モンテカルロ計算の結果を解析的に表現することを基本としているが、PW92 では Vosko らのスピン偏極依存性を取り入れ、モンテカルロ計算の結果をパラメータ6

化する際に個々の計算結果の誤差を考慮して重み付けを導入している点などの差異がある。
∫
Ēxc =
drεxc (n(r)) n(r) .
(15)
Ωcell

交換相関エネルギー密度 εxc は、交換エネルギー密度 εx と相関エネルギー密度 εc に分けることができる。
εxc = εx + εc .
(16)
スピン依存性は本パッケージの現バージョンでは考慮されていない。交換エネルギー密度としては、
(

)1

3
9π 3
εx (rs ) = −4πrs 4()1
3 3,
rs = −
4πn(17)が用いられる。相関エネルギー密度の表式は省略する。

3.4

電子−イオン間相互作用

イオンの荷電は古典的な電荷分布として扱う。それを、擬ポテンシャルとして問題に反映する。本パッケージでは、Eel−ion は擬ポテンシャル V̂Ipseudo を用いて計算される。内殻の電子の効果は擬ポテンシャルに含まれるため、擬ポテンシャル法では価電子のみを顕わに扱うことになる。或いは、内殻電子の電荷分布は固定されているとする近似を採用していることになる。本パッケージで採用している擬ポテンシャ
ルによるエネルギーは、局所ポテンシャルと非局所ポテンシャルの和で表される。
non−local
local
Eel−ion = Eel−ion+ Eel−ion

=

∑∫
I

+

Ωc elldr n(r)VIlocal (|r − RI |)∑ ∑ ∫k,n,σ I,l̸=l0Ωc elldr ϕ∗k,n,σ (r)δ V̂Il (r − RI )ϕk,n,σ (r) .

(18)

ここで、l は軌道角運動量を表す。本パッケージに含まれている擬ポテンシャルでは、l = 0, 1, 2 のどれか一つ（それを l0 とする）の寄与を局所ポテンシャルに含め、その他が非局所ポテンシャルとなることを仮定している。

VIlocal (r) = VIcore + δVIl0 .

(19)

ここで、VIcore は内殻のポテンシャルを表す。
逆格子空間における Eel−ion は次のようになる。

Eel−ion = Ωcell

′ ∑
∑
G

SI (G)n(G)VIlocal (G) + Ωcell n(0)

I

∑
I

7

αI

+ Ωcell∑ ∑ ∑SI (G − G′ )ϕ∗k,n,σ (G)δ V̂Il (k + G, k + G′ )ϕk,n,σ (G′ ) .

k,n,σ G,G′ I,l̸=l0(20)
SI (G) = exp (−iG · RI ) .

(21)

ここで、VIlocal (G) は I 番目の原子における擬ポテンシャルの局所部分の Fourier 成分、SI (G) は構造因子を示す。Eel−el , Eion−ion , Eel−ion の G = 0 での発散は、和をとることによって完全に打ち消し合う。しかしながら、擬ポテンシャル法を用いる
ことによって有限の値が残り、全エネルギーにずれが生じる。この G → 0 の極限
local
での全エネルギーのずれは、G → 0 における Eel−ionと I 番目の原子における生の2
Coulomb ポテンシャル −4πZI /Ωcell G との差によるものである。これは αI として
次のように表される。
(

αI = lim

G→0

SI (G)VIlocal (G) −4πZIΩcell G2

)

.

(22)

エネルギーの非局所部分からの寄与を計算するためには、平面波の数 N に対して、O(N 2 ) の演算量が必要である。しかしながら、Kleinman と Bylander による近似 (K-B 近似) を用いることによって、演算量を O(N ) にすることが可能である。ここでは、この K-B 近似を適用した場合の表式を用いる。

δ V̂ non−local =∑ |δ V̂ l ϕl,m ⟩⟨ϕl,m δ V̂ l |⟨ϕl,m |δ V̂ l |ϕl,m ⟩

l,m

.

(23)

ここで ϕl,m は軌道角運動量 l、磁気量子数 m をもつ擬波動関数である。(23) 式は ϕl,mが完全系を張っていれば正しい式であるが、K-B 近似ではこれを擬イオンに対する固有状態として求められた擬波動関数にとる。一般的な系に対しては、従って近似である。多くの系においては物理量に悪影響を与えないことが知られている。しかし、幾つかの場合にはこの近似による副作用が現れることも知られている。一般的な表式における非局所ポテンシャルは、本来動径座標に関しては局所な表示になっており、その意味で semi-local なポテンシャルと呼ばれている。ところが、K-B 近似を行うと、動径方向に関しても非局所なポテンシャルになってしまう。このため
Wronsky の定理が満たされなくなることが知られている。そこで、K-B 近似したポテンシャルにおいては、ノードをもつ偽の束縛状態が、ノードを持たない真の束縛状態より低いエネルギー固有状態として現れることがあり得る。このような異常な束縛状態は ghost states と呼ばれている。
K-B 近似により、非局所ポテンシャル δ V̂ l (k+G, k+G′ ) は以下のように表される。

1∫dr exp (−i(k + G) · r) δ V̂ l (r)P̂l exp (i(k + G) · r)Ω Ω
4π=
(2l + 1)Pl (cos θk+G,k+G′ )Ω∫δ V̂ l (k + G, k + G′ ) =

×

∞

0

dr r2 jl (|k + G|r)jl (|k + G′ |r)δ V̂ l (r)84πχl (k + G)χlI (k + G′ )(2l + 1)Pl (cos θk+G,k+G′ ) I, (24)ΩXIl(k + G) · (k + G′ )=,(25)|k + G||k + G′ |=cos θk+G,k+G′
χlI (k + G)

=

XIl =

∫ ∞

∫0∞
0

dr r2 jl (|k + G|r)δ V̂ l (r)Rl (r) ,(26)

dr r2 δ V̂ l (r)|Rl (r)|2 .

(27)

ここで、Rl (r) は擬ポテンシャル δ V̂ l (r) によって表される擬波動関数 ϕl,m = Yl,m (θ, ϕ)Rl (r)の動径部分である。

3.5

イオン−イオン間相互作用

この項は、Ewald 和の方法によって計算される。
1 1 ∑ ZI ZJ
Eion−ion =2 Ncell I,J |RI − RJ |

=

1 ∑
Zn Zm
2 n,m,L |dn − dm − L|(

1 ∑erfc(κ|dn − dm − L|) erf(κ|dn − dm − L|)
Zn Zm
=
+
2 n,m,L|dn − dm − L|
|dn − dm − L|

)

.(28)

ここで、L は周期系における並進ベクトルである。原子の座標 RI は単位胞内で原子の位置を指定するベクトル dn によって置き換えられ、RI = dn +L で表される。(28)式の第一項は実空間で計算され、第二項は逆格子空間で計算される。f (r) = erf(κr)の Fourier 係数は次のように計算される。
1 ∫f (G) = lim √dr f (r) exp (−iG · r − µr)
µ→0
Ωcell

4π e−(G/2κ)=.
Ωcell G2
このとき、Fourier 係数 f (G) の G → 0 での振る舞いは、
2

[

(

)

(29)
]

4π
G 2
1
f (G) ≃1−
+ O(G4 )
Ωcell
2κ
G2
]
[
1
1
4π
+ O(G2 )(30)−=
2Ωcell G2κとなる。ここで、第一項は Eel−el , Eel−ion と打ち消しあい、第二項のみが残る。それゆえ、発散項を取り除いた Eion−ion は Ewald パラメータ γewald として以下のように与えられる。
γewald =

∑
n,m

+

[

Zn Zm

(

2
4π ∑′ e−(G/2κ) iG·(dn −dm )1e−2Ωcell GG
4κ∑′ erf(κ|dn − dm − L|)L|dn − dm − L|
9

)

]

2κ
− √ δn,m .
π(31)

単位胞内に含まれる原子数が多くなると、Ewald 和の計算時間は無視できなくなる。そこで、計算量を節約するために、以下の考えに基づいて収束に必要な和の取り方を最適化する。(31) 式において、実空間と逆格子空間のそれぞれにおける和にはカットオフがあり、およそ、1/κ, 2κ の程度に見積もることができる。実空間において、単位胞の一辺が L であるとすると、逆格子ベクトルの長さは、G0 = 2π/L である。両空間での γewald の収束に必要な格子ベクトルの最大値は、カットオフ半径と Lmax = LNRmax ∝ 1/κ と Gmax ∝ G0 NGmax ∝ 2κ の関係がある。よって、両空
(

)3

κL
1
3
3
間においてとられる格子点の総数は、Ntotal ∝ NRmax+ NGmaxと見∝ (κL)
3 +ππ
積もられ、κ ≃ L にとれば最小の計算量で計算できることが分かる。単位胞が直方体の場合には、評価式は κ ≃ (Lx LyπLz )1/3 となる。

4

Hellmann-Feynman Force

この節では、原子に働く力の表式を与える。I 番目の原子に働く力 FI は、全エネ
ルギーを原子座標 RI で微分することで得られる。

FI = −
= −

dE
dRI
∑
i

[

(

)

(

)

d
d
⟨ϕi |
H |ϕi ⟩ +⟨ϕi | H|ϕi ⟩ + ⟨ϕi |H
dRI
dRI

(

d
|ϕi ⟩]
dRI

)]

(32)

ここで、i(= (k, n, σ)) は、Kohn-Sham 波動関数を指定する量子数の組を表すものとする。H を Kohn-Sham Hamiltonian とすると、全エネルギーは今の場合、E =∑i ⟨ϕi |H|ϕi ⟩ で表される。それゆえ、波動関数からの寄与が存在するため、一般には計算量が非常に大きくなる。しかしながら、波動関数 ϕi が H の規格化された固有状態の完全系を与えるとすれば、次のように (32) 式の第二項の和は 0 となる。
∑

)

[(

d
⟨ϕi | H|ϕi ⟩ + ⟨ϕi |HdRI
i∑ d⟨ϕi |ϕi ⟩
=
εi
=0
dRI
i

(

)]

d
|ϕi ⟩
dRI

(33)

それゆえ、原子に働く力の表式は以下のようになる。

FI = −

∑
i

⟨ϕi |

∂E
dH
|ϕi ⟩ = −
dRI
∂RI

(34)

(34) 式は、各波動関数がハミルトニアンの固有状態であるとき、エネルギーの原子座標に関する偏微分が原子に働く物理的な力を与えることを意味している。これはHellmann-Feynman の定理と呼ばれており、この定理を用いることによって、原子に働く物理的な力を簡単に計算することができる。しかしながら、この定理を適用するためには、全エネルギーが基底状態の値に充分収束していなければならない。
また、原子に働く力を精度良く計算することは、全エネルギーを精度良く計算する
10

ことに比べて困難になる。その理由は、全エネルギーの誤差が波動関数の誤差の２次のオーダーであるのに対して、力の誤差は波動関数の誤差の１次のオーダーだからである。
波動関数を平面波基底で展開した場合、原子座標が顕わに含まれないという性質がある。そのため、エネルギーを原子座標で微分したとき、ノルム保存型擬ポテンシャ
ルを用いた場合の力の表式としては、顕わに原子座標依存性をもつ Eel−ion , Eion−ionの２項のみが有限の値をもち、他の項は 0 となる。

FI = −
ここで、

5∂(Eel−ion + Eion−ion ) .
∂RI

(35)

∂Eel−ionは局所部分と非局所部分の２つの項で記述される。
∂RI

Quantum Stress

内部応力の表式は力の表式と同様に変分法で与えられる。体積が Ωcell である単位胞をもつ多電子系を考える。全エネルギーは E({ψk,i (r)}, {RI }) で与えられる。変
分原理は、変量である {ψk,i (r)}, {RI } に関して最適化することで、基底状態エネル
ギーが E({ψk,i (r)}, {RI }) の最小値として与えられることを示している。内部応力を得るには、２階テンソルである回転自由度を落とした対称歪みテンソル ε を与える必要がある。無限小の ε の変化に対するスケーリングから、以下を得る。

RI → R′I = (1 + ε)RI ,ψk,i (r) → ψk,i (r′ ) = det(1 + ε) 2 ψk,i (r).
1基底状態変化は、

Etot → E(ε) = Etot + ∆E(ε).
このエネルギーシフト ∆E(ε) は ε に関して展開される。応力テンソル σ はこの展開の１次の項の展開係数として与えられる。

∆E(ε) = −Tr(σε)Ωcell + O(ε2 ),1 ∂E(ε).
σαβ = −
Ωcell ∂εαβ ε→0もし、歪みと応力が等方的であれば、よく知られた P = −(36)dEという関係式をdΩcell得る。
対称テンソルである歪みテンソル εαβ に関して、一次の範囲で G は (1 − ε)G となり、Ωcell は det(1 + ε)Ωcell となる。Ωcell n(G) と SI (G) は不変であるので、以下の表式を得ることになる。
non−local
el−el
ewald
α
local
xc
kin
.
+ σαβ
+ σαβ
+ σαβ
+ σαβ
+ σαβ
+ σαβ
σαβ = σαβ

11

kin
=2
σαβ

∑

|ψi (k + G)|2 (k + G)α (k + G)β ,

k,G,i

1
el−el
σαβ
=−

′
∑
4π|n(G)|2

2 G

xc
σαβ
= δαβ

(

|G|2

)

2Gα Gβ− δαβ ,
|G|2

1 ∫drn(r)(µxc (n) − εxc (n)),Ωcell Ωcell

local
=−
σαβ

′
∑

n∗ (G)SI (G)G,I
non−local
σαβ
= −2

∑ ∑ ∑
k,i G,G′ I,l∂VIlocal (G)
,
∂εαβ

∗
(G)ψk,i (G′ )SI (G − G′ )ψk,iα
σαβ= δαβ n∗ (0)∑∂δ V̂Il (k + G, k′ + G′ )
∂εαβ

αI ,

I
ewald
σαβ
−

1 ∂γewald.
Ωcell ∂εαβ

ここで原子位置に関して参照座標を qI として導入したが、RI (t) = (1 + ε)qI (t) である。

6

ESopt の構造

現バージョンの ESopt は以下のような構成となっている。なお、将来変更される可能性がある。ESopt.tar.gz を展開すると、ESopt というディレクトリーができる。その下に、INPUT_DATA, mknon, commons, params というサブディレクトリができる。以下それぞれのディレクトリについて説明する。

6.1

ESopt

プログラムファイル (*.f90), makefile と関連した実行ディレクトリ作成スクリプ
トファイルが含まれている。ESopt では入力ファイルを以下に説明する INPUT_DATA
において作成したあとで、このメインディレクトリで make を実行し、実行ディレク
トリを作る。そののち、実行ディレクトリ（その名称は makefile で指定する）において実際の実行を行うことになる。

12

6.2

INPUT DATA計算に必要な入力ファイル群が格納されている。基本的に、ユーザが設定しなけ
ればならないのは、CORD というファイルと、*.CNTL というファイルである。CORD では、コンパイル時に必要なパラメータ群（物質の構造や計算条件）を与え、*.CNTL
は実行時に読まれるパラメータ群（動作モードの切り替えなど）を与える。擬ポテ
ンシャルのデータや、k 点サンプリングのやり方を与えるデータ、バンド計算時に計算される波数を指定するファイルもこのディレクトリにある。

6.3

mknon

計算に必要なパラメータセットを計算する二つのサブプログラムが格納されている。このプログラム群を利用することにより、ユーザは本質的に独立な入力変数のみを設定すれば良いようになっている。一つめのサブプログラム、kbexe では、FFT
mesh 数や、フーリエコンポーネントの数が、カットオフとユニットセルサイズから自動的に計算される。また、電子の状態数は原子の価数の総和から計算される。これらの、システム依存のパラメータのうち、プログラムの実行時に配列サイズ等の決定のため必要なものは、parameters.data というファイルにフォートランの unformatted
形式で出力される。それ以外のものは、input.dat という同じく unformatted 形式
で出力される。また、二つめのサブプログラム mknon では、運動量空間での FFT
メッシュ点（直方体型に分布する）と、カットオフによって球形に切り抜かれた範囲内にあるフーリエコンポーネントとの対応関係のリストを求め、gpt_list というファイルに unformatted 形式で出力する。これらのプログラムは make を実行することにより自動的に実行される。

6.4

commons と params

本パッケージにおいては、サブルーチン間での変数と配列の引き渡しを原則として以下のようにしている。配列に関しては、module 中で public 宣言したグローバル
配列を定義して各サブルーチンで参照する。一方変数の一部は、common 文による参照を行っている。これは、ベースにした Fortran 77 プログラムからの変更が完全ではないためである。mknon で計算されたパラメータにより、グローバル配列のア
ロケーションを行う際に配列サイズが決定されるが、mknon/parameters.data を読み込んで参照するための module が parameters.f90 であり、一方配列参照のためのものが、globalarray.f90 である。commons においては、サブルーチン間で共通に参照する変数を幾つかのヘッダーファイルの形で提供している。一方、params においては、LDA ポテンシャルのパラメタのような物質依存でないパラメータが格納さ
れている。

13

プログラム実行の概要

7

プログラムの実行に関する流れは以下のようになっている。

1. INPUT_DATA/CORD により、計算するシステムに関するパラメータ（原子座標、種類、ユニットセルサイズ）および計算条件（平面波のカットオフ半径、k 点
サンプリングの種類、など）を設定する。
2. makefile の先頭にある DESTDIR=destination の destination の部分が実行ディ
レクトリ名になるため、ここを変更する。なお、実習ではマシンを指定するため、マシン依存性に関する記述は省略。
3. make する。
4. cd destination する。
5. 実行モード切替用入力ファイル、SOPT.CNTL, ALLOPT.CNTL, FILE_IO.CNTL を適当に設定する。
6. ./opt で実行する。または、SGE などのキューイングシステムが動作する場合には、qsub qsub としてジョブをサブミットする。ここで、最初の qsub はジョ
ブサブミット用のコマンドであるが、次の qsub は make したときに destination
内に作られたスクリプトを指定していることになる。

8

入出力について基本的には、INPUT_DATA/CORD, makefile を設定して make し、実行ディレクト
リに移動して*.CNTL を設定すればプログラムの実行準備は終わる。以下に、これらのファイルの設定の仕方を示す。

8.1

CORD

INPUT_DATA/CORD の入力形式は以下の通りである。この例はダイヤモンド構造のSi の場合の例である。
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ <<<<<<<<<<<<<<<<<<< Basic Inputs >>>>>>>>>>>>>>>>>>>>>>>>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#Atom coordinates relative to unit vectors,
Atom,mode.
0.250000000D+00 0.250000000D+00 0.250000000D+00 Si o
0.500000000D+00 0.500000000D+00 0.500000000D+00 Si o
#^
^ ^
^ ^
^ ^^ ^
14

この部分が原子の相対座標（ユニット・ベクトルで与えられる座標に対する位置）を与える。#から#までが入力フィールドであり、この間の行数が原子数を与える。各座
標の相対座標は^から^までの桁で数値を与える。また、原子の種類は大文字、小文字を含めて判別しているため、Si のように記載する。mode は原子の動きについて示しており、o は原子位置の最適化を行い、f は原子を固定することを意味する。各元素の荷電子数や擬ポテンシャルの l0 の情報などは、INPUT_DATA/TABLE_OF_ELEMENTSにある。

================<< Other basic inputs >>=========================________________________________________________________________| unit_vec: Primitive vectors of unit cell in real space [Ang]|| kp_file:
Sampling type of primitive F.B.Z.
||[file name with 6 char.]|| E_cut:
Cutoff energy for plane-wave [Ry]|| band_calc: Choice [yes|no|metal|bands|mband]|| f_crt:Critical force [Ry/a.u.] for judging convergency. || f_tol:Parameter for judging convergency of electronic||
deg. freedom.(Relative )|| calc_mode: Calculation mode.(bas[ic] or ext[entended])|| ideb:ideb = 1; Running verbose mode.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&basic_input
unit_vec = 2.715, 0.000, 2.715,
2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
kp_file= ’DP10GM’,
E_cut= 8.0,
band_calc = ’metal’,
f_crt= 5.0D-03,
f_tol= 1.0D-8,
calc_mode = ’bas’,
ideb
= 0
&end
このセクションは、Fortran の namelist 形式で与えられている。

• unit_vec : 単位胞を与えるユニット・ベクトルを
a1x a1y a1z
a2x a2y a2z
a3x a3y a3z15の順序で与える。このとき、プログラム中での逆格子を与えるベクトルは、

ai = (aix , aiy , aiz ),(37)

a2 × a3
a3 × a1
a1 × a2, b2 =, b3 =,
2πa1 · a2 × a3
2πa1 · a2 × a3
2πa1 · a2 × a3

(38)

として、

b1 =

で与えられる。

• kp_file : k 点のリストを与えるファイルの指定。本パッケージでは結晶系によって k 点リストを取り替えて使っている。’DP10GM’ は O 群の対称性のみを考慮した立方晶系用のもの。
• E_cut はカットオフエネルギー Ecut を Ry 単位で与える。
• band_calc は計算法のスイッチ。’metal’ とした場合にはフェルミ準位近傍でフェルミ関数により準位の占有数をぼかした計算を行うことを意味している。
’bands’ とした場合には、以下の説明にある’bas’ モードを使ってバンド計算をすることができる。
• f_crt : 原子間力についての収束判定条件。
• f_tol : 電子系の CG に関する収束判定条件。
• calc_mode : 計算モード。単位胞固定で原子位置の最適化までを行う場合には’bas’ として使う。単位胞の最適化まで行う場合には’ext’ とすることで圧力一定の最適化を行うことが出来る。
• ideb : 標準出力をコントロールするパラメータ。
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<<<<<<<<<<<<< Mode dependant parameters >>>>>>>>>>>>>>>>>>>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++==============<< Case 1;

band = yes or metal. >>=============______________________________________________________________| pwidth: Width of graduall cutoff for occupation number[Ry].|| incw:
Eigen states is sorted in every incw CG steps.
|| ncband: Number of un-occupied states.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&METALDpwidth=0.005,
16incw=10,
ncband=4&END
Case 1 の部分は、’metal’ のときに必要となるパラメータ。pwidth はフェルミ分布関数の幅を与える。incw は準位の入れ替えを行うステップ回数。ncband は計算す
るフェルミ準位以上にある空の軌道の数。
====<<

Case 2;

cmode = ext. (Variable cell simulation.) >>====___[ Constant-cutoff is mimiced with following parameters. ]___||| vke:
Kinetic energy is scaled in vke where|||k+G|^2 is beyond E_cut.
|| sigma: Above scaling is faded in over sigma.
|| defmax: Maximum deformation of unit-cell allowed in this||simulation.
|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Original
vke = 2.0,
sigma = 0.05,
defmax = 3*0.05&mimicd
vke = 2.0,
sigma = 0.10,
defmax = 3*0.05&end
Case 2 の部分は、圧力一定の分子動力学計算において単位胞が変化してもカットオ
フエネルギーを止めている場合に生じうる G 点数の変化を抑える機能に関する部分。
本パッケージでは使われていない。

8.2

SOPT.CNTL

Sopt は計算上設定されるセル（計算上の単位胞のことで、結晶を計算している場
合でも必ずしも結晶単位胞と一致させる必要はない）を固定したモードでの計算を行う。
Sopt の動作を指定するコントロールファイルである INPUT_DATA/SOPT.CNTL は、
Fortran の namelist 形式を用いて書かれている。その入力形式は、以下の通りである。

17

&soptcntl
itmaxs=1,
f_crt=0.005,
ftol = 1.0d-9,
itmaxd=200,
delta = 0.4,
ideb = 0&end• itmaxs : 構造最適化の最大ステップ数。
• f_crt : 構造最適化における収束判定条件。各原子に働く力の絶対値の最大値
が f_crt 以下になったときに収束したとみなす。単位は A.U.
• ftol : 電子系の収束判定条件。1CG ステップあたりのエネルギーの相対変化を単位とする。
• itmaxd : 電子系最適化の最大 CG ステップ数。
• delta : 原子配置を最適化する時のステップ幅。原子に働く力にかかる比例
定数。

8.3

ALLOPT.CNTL

Allopt では計算上設定されるセルも含めて、構造の最適化を行ったり、ある種の分子動力学法で動かすことができる。
INPUT_DATA/ALLOPT.CNTL の入力形式は以下の通りである。
0. Chose a mode.
mode = ’md|op’: Molecular Dynamics or Optimization.
In case that MD is chosen, you decide whether "Nose’s constanttemparature formula" is used or not.
&modesw
mode=’md’,
nose = ’yes’,&end• Allopt 機能を使う場合には、分子間力と内部応力を用いて mode=’md’ として
分子動力学法として動作させる、または mode=’op’ として構造最適化を行わせることができる。nose=’yes’ とすると能勢の方法による温度コントロールを行う。nose=’no’ のときには温度コントロールなしでの分子動力学を行う。
18

0-a. Mode dependant parameters.
nose:
unit:
temp:
q_inv:

In case that nose = ’yes’, constant temparature simulation
with following parameters is performed.
The unit of temparature.[au|K|eV]
Temperature.
This value corresponds to Q^(-1).
In case that q_inv = 0, then zeta = 0.
See, Nose, Solid State Physics, Vol. 24, pp. 98-106(1989).
Nose, Mol. Phys., Vol. 52, p 255(1984).
q_inv = 0.05??

&noseparam
unit = ’K’,
temp = 100.0,
q_inv = 0.01,&end
p_crt: Convergence parameter to judge convergence of pressure in [GPa].
gamma: Friction ratio. gamma is in [0-1]. If gamma =1, MD is done.
&optparam
f_crt=1d-10,
ftol=1d-8,
p_crt=0.0001,
gamma = 0.5,&end

• unit で温度の単位を原子単位、ケルビンまたはエレクトロン・ボルトで指定できる。数値は temp で指定する。また、q_inv は能勢のパラメータ。
• f_crt は構造最適化時における原子間力の収束判定条件。単位は原子単位。
• ftol は電子状態に関するエネルギー収束判定条件。
• p_crt は構造最適化時における圧力の収束判定条件。単位は [Ry]。
• gamma は原子と単位胞の両方に関する速度に対して与えられる減衰定数。
===========================================================================
1. Define external environment.
-------------------------------------------------------------------------19

1-1. Pressure
iso: A switch to select isotropic or anisotropic pressure.
Case; iso = ’no’, pext_3d_[i|f](1:3) is read.
Case; iso = ’yes’, pext_3d_[i|f] is not read.
pext_ini: initial pressure
pext_fin: final pressure&extpress
iso=’no’,pext_init=0.0,
pext_final=0.0,
pext_3d_i=0.0,0.0,0.0,
pext_3d_f=0.0,0.0,0.0,&end• iso により等方的圧力を指定するか否かを決める。
• pext_init, pext_final または pext_3d_i, pext_3d_f により、分子動力学のプロセスにおいて外部圧力の初期値と終状態での値を指定する。
--------------------------------------------------------------------------1-2. Parameters concernig how the pressure is introduced and how MD runs.
nsec_p:nstp_i:
nstp_g:External pressure is gradually enforced.
Sectioned by nsec_p steps between pext_fin and pext_ini.
nstp_i MD steps are performed at the pext_ini.
nstp_g MD steps are performed at each pressure
except for pext_ini.

&prscntl
nsec_p=10,
nstp_i=10,
nstp_g=1,&end• nsec_p, nstp_i, nstp_g により外部圧力の設定値の変化を指定する。
==========================================================================
3. Parameters for controlling evolution.
20

delt_t:ww:itmaxa:itmaxd:Time interval in femto second.
Fictitious mass ratio [w/mass]
Maximam iterations for MD.
Maximam iterations in diagw.

&mdparams
delt_t=5.0,
ww = 100000.0,
itmaxa=500,
itmaxd=200,&end• delt_t：分子動力学シミュレーションにおける時間刻み。
• ww：単位胞に対する仮想的な質量。
• itmaxa：分子動力学計算のステップ総数。
• itmaxd：CD の最適化における最大 CG ステップ数。
=========================================================================
4. Initial condition.
ekin_init:
Initial kinetic energy is scaled in this value.
init_conf:
Switch how is the initial velocities are given.
Case; init_conf = ’given’, read from iveloc namelist in this file.
Case; init_conf = ’rando’, random vectors, generated with ’iseed’,is used.
ivpoint:
Initial kinetic energy is given at ivpoint step.
&initveloc
ekin_init = 10.0,
ckin_init = 0.0,
init_conf = ’given’,
iseed = 180,
ivpoint = 1,
&end-----------------------------------------------------------------------4-1. Define initial direction to move.
21

This parameters are read only in case init_conf = ’given’.
how:Case; how = ’rel’, final conf - initial conf is used.
Case; how = ’abs’, final conf. is used as a vector.

latv_i(3,3): initial primitive vector in real space [Ang]latv_f(3,3): final primitive vector in real space [Ang]coordi(3,natm): initial coordinate
coordf(3,natm): final coordinate
Example for structural transformation from cubic-diamond to beta-tin
of silicon.
&iveloc
how = ’rel’,latv_i =2.715, 0.000, 2.715,
2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
latv_f =2.715, 0.000, 2.715,
2.715, 2.715, 0.000,
0.000, 2.715, 2.715,
coordi =0.000000000, 0.000000000,
0.500000000, 0.500000000,
coordf =0.000000000, 0.000000000,
0.500000000, 0.500000000,&end

0.000000000,
0.000000000,
0.000000000,
0.000000000,• この部分は、初期速度の与え方に関する部分。ESopt では最後の coordi, coordf
は読み込まれない。

8.4

FILE IO.CNTL

INPUT_DATA/FILE_IO.CNTL の入力形式は以下の通りである。
======================================================================
1. Parameter(s) to controll output.
22

&file_io_cntl
read_file=’no’,
initial_file=1&end

intchg: interval of charge output&outcntl
intchg = 10,
intefl = 20,
intwf = 10,&end

2. Debug Flags&debug_fl
ideb = 0&end• read_file : 計算開始時の原子配置および波動関数として以前の計算結果を利用するかどうかを与える。’cont’, ’wyes’, ’cyes’, ’syes’ または’no’ が指定できる。それぞれ以下のような動作をする。
– ’cont’: STRUCT/strct から構造を、波動関数を実行ディレクトリーにある fort.98 から読み込む。
– ’wyes’: 波動関数を実行ディレクトリーにある fort.98 から読み込む。
– ’cyes’: STRUCT/cell-XXX および STRUCT/latcrd-XXX を読む。XXX は
I3.3 によって与えられる数で、initial_file によって与える。波動関数を実行ディレクトリーにある fort.98 から読み込む。
– ’syes’: STRUCT/cell-XXX および STRUCT/latcrd-XXX を読む。XXX は
I3.3 によって与えられる数で、initial_file によって与える。
– ’no’: 構造は CORD で与えられるものを利用し、波動関数は乱数により発生させたものを使う。
• initial_file : read_file=’cyes’, read_file=’syes’ のときの XXX の値。
read_file がそれ以外のときには無視。

23• intchg : 電荷密度分布を出力するタイムステップの間隔を与える。
• intefl : MONIT/energy を出力する CG ステップの間隔を指定する。IO 時間を短縮したいときには大きめの値を指定する。
• intwf : MONIT/wave000 を出力する間隔を指定していた。現バージョンでは波動関数の出力は実行ディレクトリーに fort.99 として出力される。計算を接続するときにはこのファイルを fort.98 としてコピーして使う。

8.5

kpointkpoint の入力形式は以下の通りである。
number of line (i4) =7^^^^

G X K G L K W X
number of mesh (i4) = 10
1.000000000 0.000000000
1.000000000 0.500000000
number of mesh (i4) = 10
1.000000000 0.500000000
0.750000000 0.375000000
number of mesh (i4) = 100.750000000 0.375000000
0.000000000 0.000000000
number of mesh (i4) = 10
0.000000000 0.000000000
0.500000000 0.500000000
number of mesh (i4) = 100.500000000 0.500000000
0.750000000 0.375000000
number of mesh (i4) = 100.750000000 0.375000000
0.750000000 0.250000000
number of mesh (i4) = 100.750000000 0.250000000
0.500000000 0.000000000

0.000000000
0.500000000
0.500000000
0.375000000
0.375000000
0.000000000
0.000000000
0.500000000
0.500000000
0.375000000
0.375000000
0.500000000
0.500000000
0.500000000

ここで、number of line はブリュアン域において何本の k 線上で計算バンド分散を求めるのかを指定する。この数は 20 以内に制限されている。上記のようにこの数が 7 であるとする。この k 線は連続していて一筆書きできるものとする。一行空
いて、次の行にはこの 7 本の k 線の両端点における名称を書くことができる。空白

24で区切られた文字列を指定することができる。対称性によって決まった名称があればそれを指定してよい。次に、ひとつの k 線ごとに 3 行を用いて、k 線上の分割点数と、両端点を逆格子ベクトル b1 , b2 , b3 の与える空間内の点として与える。

8.6

出力について計算の途中経過には、実行ディレクトリにつくられる MONIT というディレクトリ
にエネルギー等の情報が出力される。全エネルギーは、energy というファイルに出力される。波動関数は、intwf ごとに wave000 という名前で出力される。(ESopt ではこのオプションは停止している。) 構造最適化も終了した時には wave999 という名前で出力される。波動関数は unformatted で出力される。
実空間の電荷密度分布は、CHARGE というディレクトリの下に、chgdnsXXX というファイル名で保存される。
結果は、以下の場所に書き出される。

• Sopt を用いた場合：RESULTS– cord***: 原子座標（実空間での値）– eig***: Kohn-Sham 方程式の固有値– gradient*** : 原子間力
– H***: 格子ベクトルが作る３×３行列の対角項、全エネルギー、
単位胞の体積– Htot***: 全エネルギーを与える以下のエネルギーが順番に出力されている。１．電子系の運動エネルギー、２．ハートレー項、３．交換・相関項、４．擬ポテンシャルの局所成分によるポテンシャルエネルギー、
５．擬ポテンシャルの非局所成分によるポテンシャルエネルギー、６．α
項、７．エバルト項、８．１〜７の合計。
• Allopt を用いた場合：STRUCT– cell***:格子ベクトル– latcrd***:内部座標– press***:外部圧力– realcrd***.xyz:原子座標– result***:全エネルギー（Sopt の Htot***と同じ）、エネルギーとエンタルピー、圧力、原子間力の順に出力される。
– strct***:単位胞のパラメータ, cell***, latcrd, と zeta の値などが出力される。
– strn***:歪みテンソルが出力される。
25

プログラム・ソースコード

9

この節では、各サブルーチンにおいて物理量の計算がどのように実現されているかの例を示すことで、プログラムの動作を理解するための情報を提供する。

9.1

運動エネルギー

運動エネルギーは、etotal.f90 に含まれる関数 ekin(inw) の中で以下のように計算されている。

!**********************************************************************
double precision function ekin(inw)use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: ik,ib,ipw,inw
real(8) frm,z,fd,x!

-----<< definition of fermi distribution function >>----frm(z) = 1d0 + vke_mim/( 1d0 + exp( -1d0*beta_mim*z ) )! ----------<< calc. of KE (mimicing) >>---------!* epw[Ry],gk[a.u.],nadd
ekin = 0d0
do ib = 1,mxband_g
do ik = 1,nk_g
fd = fdist_g(ib,ik)
do ipw = 1,mxpwk(ik)x = 0.5d0*( gk_g(ipw,ik) - ec_mim )ekin = ekin &+ gk_g(ipw,ik)*frm(x) &* ( dreal(wvfn_g(ipw,ib,ik,inw))**2 &+ dimag(wvfn_g(ipw,ib,ik,inw))**2 ) * fd
end doend doend do26

ekin = 2d0*ekin/dfloat(natm_g)/effnk
returnend

ここで、• inw: 波動関数は CG ステップを実現するためにステップの経路上で幾つか保存されている。その波動関数を指定するインデックス。
• ib, ik, ipw: それぞれバンドインデックス、k、及び波数ごとに決められる Ecut
以下のエネルギーをもつ平面波 exp(i(k + G) · r)) を定める G を示す。
• gk_g(ipw,ik): |k + G|2 が格納されている配列。
• wvfn_g(ipw,ib,ik,inw): 波動関数 ϕk,n (G) を格納した配列。
である。よって、変数 ekin にEkin = 2∑|k + G|2 |ϕk,n (G)|2 ,k,n,Gが代入される。これは、ekin は単位体積当たりのエネルギーを、単位を Rydberg 単位 (電子質量 m = 1/2 とする) として、スピン自由度の和からくる因子 2 を含めて与えているためである。さらに、因子 frm(x) によって、単位胞の変形 (ESopt では
圧力一定シミュレーションの機能を止めてあるが）を許した際に、実効的にエネル
ギーカットオフ一定の条件が満たされるようにしている。

9.2

ハートレー・エネルギー

ハートレー・エネルギーは、etotal.f90 に含まれる関数 eee(vcell) の中で以下のように計算されている。

!**********************************************************************
double precision function eee(vcell)use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: igreal(8) :: vcell

27

eee = 0d0
do ig = 2,ng_g
eee = eee + (dreal(chgg_g(ig))**2+dimag(chgg_g(ig))**2)*gg_g(ig)end doeee = eee*4d0*pi/vcell/dfloat(natm_g)
returnend

ここで、• chgg_g(ig): 密度のフーリエ変換 n(G) を格納した配列。
• gg_g(ig): 1/|G|2 が格納されている配列。
である。よって、変数 eee にEel−el =′4π ∑|n(G)|2,Ωcell G |G|2が代入される。これは、eee も単位体積当たりのエネルギーを、単位を Rydberg 単位 (電子質量 e2 = 2 とする) として与えているためである。Ωcell でむしろ割られているのは、本来 Eel−el を与えていた実空間積分が２重積分になっていることから来
ている。

9.3

全エネルギー

全エネルギーは、etotal.f90 に含まれるサブルーチン etotal(vcell,inw,etot)の中で以下のように計算されている。

!**********************************************************************
subroutine etotal (vcell,inw, etot)
use parameters
use globalarray
implicit none
integer :: lmx,inw,i,jreal(8) :: etot(8),vcell,vatm, &ekin,eee,exc,evl,evnl,t1,t2,dsecnd

28

parameter( lmx = 2 )include ’commons/common.h’
t1 = dsecnd(’cpu’)vatm = vcell/dfloat(natm_g)etot(1) = ekin(inw)etot(2) = eee(vcell)etot(3) = exc(vcell)etot(4) = evl()
etot(5) = evnl(inw)etot(8) = etot(1)+etot(2)+etot(3)+etot(4)+etot(5)+etot(6)+etot(7)t2 = dsecnd(’cpu’)
cpu_time(1) = cpu_time(1) + t2-t1
returnend

エネルギーは項ごとに計算され、それぞれ以下の値が入る。

• ’ ekin= ’,etot(1) : 運動エネルギー• ’ eee= ’,etot(2) : ハートレー・エネルギー• ’ exc= ’,etot(3) : 交換相関エネルギー• ’ evl= ’,etot(4) : 擬ポテンシャルの局所項の寄与• ’ evnl= ’,etot(5) : 擬ポテンシャルの非局所項の寄与• ’ ealp= ’,etot(6) : αI の寄与• ’ eewd= ’,etot(7) : イオン-イオン間相互作用 (エバルト和)• ’ etot= ’,etot(8) : 全エネルギー

9.4

Sopt コントロール・ルーチン

main ルーチンはもともと幾つかのモードに対して分岐する形となっていた。
• Allopt : 格子構造も含めた最適化、圧力一定の分子動力学などを行うモード。
• Sopt : 格子を止めたままで、原子位置の最適化を行うモード。
29

• DMD : 減衰効果を含む分子動力学モード。
ESopt は、Sopt, Allopt の一部の機能以外を停止した形で提供している。そこで、
Sopt の本体であるサブルーチン sopt を次に考察する。以下では実際のものから無駄な出力等は省いてある。

!**********************************************************************
subroutine sopt()
use parameters
use globalarray
implicit real*8(a-h,o-z)
include ’commons/lattice.h’include ’commons/common.h’logical(4) :: exist1, exist2
real(8) :: etot(8),delta,cpu0,cpu1
namelist / soptcntl / itmaxs, ftol, f_crt, itmaxd, delta, ideb!==================================================================
inquire(file=’SOPT.CNTL’,exist=exist2)
inquire(file=’../INPUT_DATA/SOPT.CNTL’,exist=exist1)
if(exist2) thenopen(10,file=’SOPT.CNTL’,form=’formatted’, &iostat=ios,status=’old’)read(10,soptcntl,iostat=ios)close(10)else
if(exist1) thenopen(10,file=’../INPUT_DATA/SOPT.CNTL’,form=’formatted’,&iostat=ios,status=’old’)read(10,soptcntl,iostat=ios)close(10)elsestop ’Can not open SOPT.CNTL’
endif
endif
!

----------< generate g-points >--------------call gptgen()!

-----< calc. some components of total energy >----------

30

call vinit()!

-----< calc. ewald sum (independent part with coordinates of atoms ) >----calliniewd()! Main loop starts.
do its = 1, itmaxs! ----------< optimization for electronic state >---------call getll(etot,nloop)
call charge(0)! ----------< calc. Hellmann-Feynman force >---------call force(vcell)
call eforce(natm_g,stoa,grd,frc,fn)
call wfile(etot,its)! ----------< output data >---------call wfile0(000,its)! ----------<

judgement of convergence >-------------

if(ideb.eq.1) write(*,*)’fn, f_crt = ’,fn,f_crt
if(fn.lt.f_crt) exit
! ----------< move all atoms >-------------------------------call mvatm(delta)
end do! Main loop ends.
if(fn.lt.f_crt) thencall wfile0(999,its)
end ifwrite(*,*)’elapse time in sopt = ’, cpu1-cpu0
returnend subroutine sopt

31

これにより、アルゴリズムの流れになっていることが分かる。

1. 初期化部分(a) SOPT.CNTL 読み込み(b) gptgen() : G などのデータを作成。
(c) vinit() : 擬ポテンシャル・データの入力。
(d) iniewd() : エバルト和計算のための初期化。
2. 構造最適化ループ(a) getll(etot,nloop) : 電子状態決定用ルーチン(b) charge(0) :電荷密度分布計算ルーチン(c) force(vcell) : 原子間力計算ルーチン(d) eforce(natm_g,stoa,grd,frc,fn) : 最適化の方法により force を調整するルーチン(e) wfile(etot,its), wfile0(999,its), wfile0(000,its) : 出力用ルーチン群(f) mvatm(delta) : 原子の移動を起こさせるルーチン

9.5

対角化用サブルーチン

サブルーチン diagw において、電子状態の CG 法による最適化が行われる。

subroutine diagw (nloop,vcell,alp,etot)
use parameters
use globalarray
implicit none
include ’commons/common.h’integer :: nloop, inw, itmax, ik, in, ig, ist, loop, its, nwf, nnn, &istc, ireal(8) :: etot(8), vcell, alp, eps, ang, d1, d2, dnf, step, ep, e_p,&fret, fp, gd, avm, dmm, t1, t2, dsecnd,e2,e1
complex(8) :: ggc,gam,s, chg1
parameter (eps = 1d-10)32

parameter (inw = 2)
parameter (e_p = 5d0)
! -----< calc. charge density for G-space >----call schmit(inw)call toreal(inw)chg1 = chgg_g(1)etot(6)=alp*dreal(chg1)! -----< calc. total energy(fp) >----call etotal (vcell,inw,etot)
fp = etot(8)! -----< calc. force >----call wforce(vcell,inw)
call precnd(e_p)
call schmi2(inw,1)
itmax = max(itmaxd*2,npw_g*mxband_g*nk_g*4)! nwf(= npw*nband*nk) is num. of wave func.
! -----< store old data of w.f. >----do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,1)
end doend doend do! -----< iteration start >----ist = 1
istc = 0loop = 0! Main loop starts
do its = 1, itmaxs = (0d0,0d0)33

d1= 0d0
d2= 0d0! -----< summation of some value >----nwf = mxband_g*nk_g
do in = 1,mxband_g
do ik = 1,nk_g
do ig = 1,mxpwk(ik)s = s + dconjg(wvfn_g(ig,in,ik,1))*wvfn_g(ig,in,ik,2)
d1= d1+dreal(wvfn_g(ig,in,ik,1))**2+dimag(wvfn_g(ig,in,ik,1))**2
d2= d2+dreal(wvfn_g(ig,in,ik,2))**2+dimag(wvfn_g(ig,in,ik,2))**2
end doend doend doang = acos(cdabs(s)/sqrt(d1*d2))/pi
dnf = sqrt( d1/nwf )if(ang.lt.0.46d0.and.ang.ge.0.385d0) thencall schmi2(inw,2)elseif(ang.lt.0.385d0)thencall schmi2(inw,1)
endif
call linmin(vcell,etot,fret,ist)if(ist.eq.0) istc=istc+1
step=2d0*dabs(fret-fp)/(dabs(fret)+dabs(fp)+eps)!

----------< judgement for convergency >---------if((step.le.ftol.and.ist.eq.1).or.its.gt.itmaxd.or.istc.gt.3) thenloop = loop + 1
if(loop.gt.2) thenetot(6)=alp*dreal(chg1)
call efile(etot,ang,dnf,its,ftol,step,loop)
nloop = its
write(*,*) ’Count (its,loop)’, its,loop
returnelse
ist = 034

endifendif!

--------------------------------------------------ep=fret-fp
fp = fretcall wforce(vcell,inw)
call precnd(e_p)
call eigenv(avm,dmm,its,nnn)
if(ist.eq.0.or.nnn.eq.1) thencall toreal(inw)chg1 = chgg_g(1)call wforce(vcell,inw)
call precnd(e_p)
endif
if(its.eq.0) ist = 0
if(ist.eq.1) thengd = 0d0
ggc = (0d0,0d0)
do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)gd = gd + dreal(wvfn_g(ig,in,ik,4))**2 &+ dimag(wvfn_g(ig,in,ik,4))**2
ggc = ggc + dconjg(wvfn_g(ig,in,ik,1) &- wvfn_g(ig,in,ik,4))*wvfn_g(ig,in,ik,1)
end doend doend doif(gd.eq.0d0) thenwrite(*,*)’warning. really gd = 0 ?’returnendif
gam = ggc/gd
do in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
end do35

end doend dodo in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,4) &+ gam*wvfn_g(ig,in,ik,5)
end doend doend dodo in = 1, mxband_g
do ik = 1, nk_g
do ig = 1, mxpwk(ik)wvfn_g(ig,in,ik,1)=wvfn_g(ig,in,ik,5)
end doend doend doelseif(ist.eq.0) thendo in = 1,mxband_g
do ik = 1,nk_g
do ig = 1,mxpwk(ik)wvfn_g(ig,in,ik,4)=wvfn_g(ig,in,ik,1)
wvfn_g(ig,in,ik,5)=wvfn_g(ig,in,ik,1)
end doend doend doist = 1endif
nnn = mod(its,intefl)
if(nnn.eq.0) thenif(ideb.eq.1) write(*,*) ’efile in’call efile(etot,ang,dnf,its,ftol,step,loop)
if(ideb.eq.1) write(*,*) ’efile out’
endif
end do! Main loop ends
stop ’Maximum iterations exceeded (diagw).’
end subroutine diagw

36