構造最適化
ストレステンソルを利用したユニットセル最適化機能
PHASE/0には,ストレステンソルを利用して単位胞を最適化する機能が備わっています。ここでは,この機能の使い方の説明を行います。
入力パラメータ
まずは,通常通り入力パラメータファイルを記述します。セルを変形させたあとに座標の緩和を行いたい場合は通常通り原子座標最適化用のパラメータを設定すればセルの変形→力が収束していない場合は原子座標の最適化,という動作をするようになります。さらに,単位胞最適化用の,以下のような設定を加えます。
structure_evolution{
lattice{
sw_optimize_lattice = on
}
}
変数sw_optimize_latticeをonとすると本機能を利用することができます。latticeブロックには,以下の変数を定義することが可能です。
sw_optimize_lattice |
単位胞最適化機能を有効に する場合onとします。 デフォルト値はoff です。なお,このスイッチがonの場合 はsw_stressは自動的にonになります。 |
sw_uniform |
単位胞を一様に変化させたい 場合にonとします。デフォルト値はoff です。このパラメータがonの場合,スト レステンソルの対角要素の平均値によっ て体積を変化させるように動作します。 |
sw_rebuild_pws |
単位胞を変 形させた後に平面波基底を作り直すかど うかを指定します。デフォルト値はon, つまり格子が変 形する度に平面波を作り直します。Off とすることによって電子状態計算の収束 性を向上させることができますが,格子 が変形しても同じ平面波セットを利用し ている,ということは厳密にはカットオ フエネルギーが微妙に変化している,と いうことに相当する点に注意が必要です 。また,このパラメータをoffとすると 継続計算ができなくなってしまいます。 |
method |
最適化の手法を選択します。bfgs, quench, sdのいずれかを 指定します。デフォルト値はbfgsです。 |
delta_stress |
me thodがquenchかsdの場合の更新の刻み幅 を指定します。デフォルト値は1です。 |
max_stress |
収束判定に利用する ,ストレステンソルの最大値を圧力の単 位で指定します。デフォルト値は1.e-6 ha rtree/bohr3です。sw_unifo rmがonの場合はストレステンソルの対角 要素の平均が収束判定に採用されます。 |
sw_optimize_coordinates_once |
原子配置の最適化は1回目の格子の 更新時のみ行いたい場合にonとします。 |
sw_fix_length_a |
onとするとa軸 の長さを固定して格子を最適化します。 |
sw_fix_length_b |
onとするとb軸 の長さを固定して格子を最適化します。 |
sw_fix_length_c |
onとするとc軸 の長さを固定して格子を最適化します。 |
sw_fix_angle_alpha |
onとすると格子 定数αを固定して格子を最適化します。 |
sw_fix_angle_beta |
onとすると格子 定数βを固定して格子を最適化します。 |
sw_fix_angle_gamma |
onとすると格子 定数γを固定して格子を最適化します。 |
バージョン2019.01未満の場合の注意点
ストレステンソルのカットオフエネルギーに対する収束性はかなり悪い場合があります。ストレスミニマムと全エネルギーのミニマムが一致しない場合,おそらくカットオフが不十分であることが要因と思われます。このようなケースに遭遇したら,ストレステンソルとカットオフエネルギーの関係をしらべていただくことを推奨します。
ストレステンソル補正(バージョン2019.01以降)
ストレステンソルの計算精度は低い場合が多いので,補正する設定を施すことが推奨されます。
structure_evolution{
stress{
sw_stress_correction = on
}
}
電荷密度および波動関数の再利用(バージョン2019.02以降)
格子最適化を行う場合,格子が変形するたびに電荷密度や波動関数のメッシュは変化します。そのため,2019.02未満のバージョンにおいては波動関数や電荷密度は格子が変形するたびに作り直す振る舞いがデフォルトの振る舞いでしたが,バージョン2019.02以降より新しいメッシュに補間して古い電荷密度および波動関数データを割り当てる動作がデフォルトの振る舞いとなりました。この変更により,SCF計算の収束に要する回数を大幅に減らすことができるようになりました。この動作を抑制したい場合は以下のような設定を行います。
structure_evolution {
lattice {
sw_optimize_lattice = on
sw_interpolate_charge = off
sw_interpolate_wfs = off
}
}
計算結果の出力
結果はoutput000ファイル,nfefn.dataファイル,nfdynm.dataファイルに記録されます。
output000ファイルには,ストレステンソルが記録されます。以下のコマンドによってその情報を抽出することができます。
% grep –A3 ‘STRESS TENSOR$’ output000
STRESS TENSOR
0.0002326236 0.0000000000 0.0000000000
0.0000000000 0.0002326236 0.0000000000
0.0000000000 0.0000000000 0.0002142790
--
STRESS TENSOR
0.0002272841 0.0000000000 0.0000000000
0.0000000000 0.0002272841 0.0000000000
0.0000000000 0.0000000000 0.0002077216
...
...
通常の計算の場合ストレステンソルが1組出力されるのみですが,本機能を利用している場合はストレステンソルの履歴が出力されます。
nfefn.dataファイルには,通常通り全エネルギーや原子に働く力の最大値のほか,ストレステンソルの最大値(sw_uniformがonの場合は対角要素の平均値)が記録されます。たとえば,以下のような出力が得られます。
iter_unitcell, iter_ion, iter_total, etotal, forcmx, stressmx
1 1 18 -181.4043211413 0.0020128619
1 2 27 -181.4043355689 0.0015666906
1 3 36 -181.4043464493 0.0011267018
1 4 44 -181.4043509953 0.0008837770
1 5 53 -181.4043582176 0.0000137026 0.0002326236
2 1 73 -181.4044226903 0.0000645338 0.0002272841
nfdynm.dataファイルも通常のものとほぼ同じですが,通常の計算の場合は一度しか出力されないヘッダーが,セルベクトルが変形される度に出力されます。
#
# a_vector = 8.6795114819 0.0000000000 0.0000000000
# b_vector = 0.0000000000 8.6795114819 0.0000000000
# c_vector = 0.0000000000 0.0000000000 5.5916992108
# ntyp = 2 natm = 6
# (natm->type) 2 2 1 1 1 1
# (speciesname) 1 : O
# (speciesname) 2 : Ti
#
cps and forc at (iter_ion, iter_total = 1 18 )
1 0.000000000 0.000000000 0.000000000 0.000000 0.000000 0.000000
2 4.339755741 4.339755741 2.795849605 0.000000 0.000000 0.000000
3 2.643779197 2.643779197 0.000000000 -0.001423 -0.001423 0.000000
4 6.983534938 1.695976544 2.795849605 -0.001423 0.001423 0.000000
……
……
#
# a_vector = 8.7672856463 0.0000000000 0.0000000000
# b_vector = 0.0000000000 8.7672856463 0.0000000000
# c_vector = 0.0000000000 0.0000000000 5.6429940606
# ntyp = 2 natm = 6
# (natm->type) 2 2 1 1 1 1
# (speciesname) 1 : O
# (speciesname) 2 : T
#
cps and forc at (iter_ion, iter_total = 1 111 )
1 0.000000000 0.000000000 0.000000000 0.000000 0.000000 0.000000
2 4.383642823 4.383642823 2.821497030 0.000000 0.000000 0.000000
3 2.663907294 2.663907294 0.000000000 0.001773 0.001773 0.000000
4 7.047550117 1.719735530 2.821497030 0.001773 -0.001773 0.000000
5 1.719735530 7.047550117 2.821497030 -0.001773 0.001773 0.000000
6 -2.663907294 -2.663907294 0.000000000 -0.001773 -0.001773 0.000000
……
……
計算例:ルチル型TiO2 (ストレス補正なし)
ルチル型TiO 2 の格子最適化を行った例を紹介します。入力データは samples/unitcel_optimization/TiO2
以下にあります。
入力パラメータファイルには,以下のような設定を施しました。
カットオフエネルギーは80 Rydberg
擬ポテンシャルはポータルサイトにおいて公開されているTi_ggapbe_us_02.ppとO_ggapbe_us_02.pp
原子座標の最適化を施す設定;手法はBFGS法,収束判定となる力の最大値は2e-4
初期原子配置および格子定数は,無機材料データベースAtomWork(http://crystdb.nims.go.jp/)に登録されていたルチル型TiO2のデータを採用
波動関数ソルバー,電荷密度ミキサーは指定せず,デフォルト設定を採用。
採用したカットオフエネルギーは80 Rydbergと比較的大きなものですが,後述のようにTiO2の場合はこれくらい必要であると考えられます。
nfefn.dataファイルの内容は,以下のようになりました。
iter_unitcell, iter_ion, iter_total, etotal, forcmx, stressmx
1 1 18 -181.4043211413 0.0020128619
1 2 27 -181.4043355689 0.0015666906
1 3 36 -181.4043464493 0.0011267018
1 4 44 -181.4043509953 0.0008837770
1 5 53 -181.4043582176 0.0000137026 0.0002326236
2 1 73 -181.4044226903 0.0000645338 0.0002272841
3 1 92 -181.4044839579 0.0001241955 0.0002222588
4 1 111 -181.4056948858 0.0025074070 0.0002222588
4 2 120 -181.4057176163 0.0020195652 0.0002222588
4 3 130 -181.4057600852 0.0000156213 0.0000444895
……
……
9 1 248 -181.4058191217 0.0001647915 0.0000332105
10 1 268 -181.4058328662 0.0000709369 0.0000119789
11 1 287 -181.4058349707 0.0000268520 0.0000015502
12 1 306 -181.4058351835 0.0000244918 0.0000006790
まずは,原子座標の最適化が5回実施されています。その間ストレステンソルは未計算なので,6列目は空欄になっています。5回目で原子に働く力の最大値が閾値より小さくなったので,セルを変形させたのちに計算が進行しています。この際に,単位胞最適化の更新回数を表す1列目の数値が2になっていることがわかります。また,6列目にストレステンソルの最大値が記録されています。2回目と3回目の更新時はセルを変形させても原子に働く力の最大値は閾値以下だったので原子座標の最適化は実施されませんでしたが,4回目セルベクトル更新時にはそうではなかったので原子座標の最適化が行われています。このようにセルの最適化と必要に応じた原子座標の最適化が行われつつ計算が進行し,セルの更新回数が12回となったところでストレステンソルの最大値が閾値以下となったので計算は収束したとみなされ終了しています。単位胞最適化収束の履歴を,図にまとめました。

単位胞最適化の履歴。赤線:全エネルギー,緑線:ストレステンソルの最大成分。
安定な格子定数は,nfdynm.dataファイルに記録された最後のセルベクトル更新の情報からもとめることができます。この例の場合,a=8.7934 bohr, c=5.6164 bohrと得られました。
ストレステンソルとカットオフエネルギー
ストレステンソルは,全エネルギーや原子間力と比較してカットオフエネルギーに対して収束しづらい傾向があります。例として,ルチル型TiO2の,実測値の格子定数で計算したストレステンソルとカットオフエネルギーの関係を図にプロットしました。

ルチル型TiO2の場合の,ストレステンソルとカットオフエネルギーの関係
図からわかるように,カットオフ50 Rydberg程度の場合ストレステンソルの符号が間違ってしまっています。このケースでは,ある程度収束したストレステンソルを得るためには,最低でも80 Rydberg程度以上のカットオフエネルギーが必要であることが示唆されます。
計算例:ルチル型TiO2 (ストレス補正あり,バージョン2019.01以降)
5.6.1.6 章 で行った最適化を,ストレス補正を有効にして実行してみます。
入力データは samples/unitcell_optimization/TiO2_with_correction
にあります。
入力パラメータファイルには,以下のような設定を施しました。 計算条件は下記の通り。
カットオフエネルギーは36 Rydberg
ストレス補正は有効 (structure_evolutionブロックの下のstressブロックにおいてsw_stress_correction = onと設定)
擬ポテンシャルはポータルサイトにおいて公開されているTi_ggapbe_us_02.ppとO_ggapbe_us_02.pp
原子座標の最適化を施す設定;手法はBFGS法,収束判定となる力の最大値は2e-4
初期原子配置および格子定数は,無機材料データベースAtomWork(http://crystdb.nims.go.jp/)に登録されていたルチル型TiO2のデータを採用
波動関数ソルバー,電荷密度ミキサーは指定せず,デフォルト設定を採用。
結果得られる単位胞最適化収束の履歴を次の図に示します。比較的低いカットオフ(36 Rydberg)を採用しているにも関わらずスムーズな最適化が行えました。結果得られる格子定数は,a = 8.806 bohr, c = 5.619 bohrとなりました。
単位胞最適化の履歴。赤線:全エネルギー,緑線:ストレステンソルの最大成分。
ノンコリニア系の計算、スピン軌道相互作用計算
ノンコリニア系の計算
機能の概要
通常の計算では、局所的な磁気モーメントの大きさは、アップ・ダウンスピン電子密度の差で表現され、その向きについては規定されません。向きまで表現するためには、1つの波動関数がアップ及びダウンスピン成分をもつようにする、すなわち2成分スピノールで表す必要があります。
これに対応して電荷密度は
のように、スピン指標に関して2x2行列となります。また、局所的な電荷密度及び磁気モーメントは、
となります。
入力パラメータ
ノンコリニア系の計算、すなわち2成分スピノールでの計算を行うには、structureブロックに “magnetic_state = noncollinear “ と記します。また、各原子種の局所磁気モーメントの方向の初期値を、”mdx mdy mdz” に指定します。特に指定しない場合は、z 方向を向くと判断します。局所磁気モーメントの方向の初期値は、”theta phi” (単位:degree) でも指定可能です。
structure{
...
magnetic_state = noncollinear
element_list{
#units atomic_mass
#tag element atomicnumber zeta deviation mdx mdy mdz
O 8 0.166666 1.5 0.0 0.0 1.0
}
...
}
計算結果の出力
標準出力には、以下のような磁気モーメントに関する情報が出力されます。Tot, Mx, My, Mzは、それぞれ、局所的な電荷密度及び磁気モーメントを単位胞内で足し合わせた値です。これ以外は通常の計算による出力と共通です。
!OLD Chg \*\* Tot: 40.00000000 Mx: 0.00073742 My: 0.00000000 Mz: 16.17742289
!NEW Chg \*\* Tot: 40.00000000 Mx: 0.00075559 My: 0.00000000 Mz:
スピン軌道相互作用計算
機能の概要
2成分スピノールが重要になるのは、スピン軌道相互作用を考慮する場合です。スピン軌道相互作用は、
で表されます。ここで、は原子核周りの球対称なポテンシャルです。このHamiltonianは、波動関数のアップ及びダウンスピン成分の間に相互作用を働かせるため、5.8.1 で説明した2成分スピノールが必要になります。
入力パラメータ
スピン軌道相互作用を利用するには、accuracyブロック内のspin_orbitブロックに “mode = pawpot “ と記します。これ以外は、ノンコリニア系の計算と同様です。
accuracy{
…
spinorbit{
mode = pawpot
}
…
}
計算結果の出力
スピン軌道相互作用によるエネルギーは、標準出力にESpinOrb_old, new で表示されます。
なお、これらの値は、コンパイル時にCPPFLAGに -DUSE_ESPINORB をつけた場合にのみ計算されます。
TOTAL ENERGY FOR 53 -TH ITER= -41.454944288742 edel = -0.170628D-08 : SOLVER = SUBMAT + RMM3
KI= 13.204535394898 HA= 32.283599969986 XC= -6.801519951682 LO=
NL= 7.597059454569 EW= 5.402894293900 PC= 0.000000000000 EN=
PHYSICALLY CORRECT ENERGY = -41.454944288742
EOHXC_PAW= -0.4786729 HA_PAW= 0.0218350
XC_PAW_AE= -15.6619733 XC_PAW_PS= -5.6000049
!XC_PAW_AE-XC_PAW_PS= -10.0619684
ESpinOrb_old= -0.0000308 ESpinOrb_now= -0.0000293
計算例:O2分子、Pt2分子
計算例は、 samples/spin_orbit
以下のサブディレクトリーにあります。これらの例では、分子をx軸方向に配置し、磁気モーメントの向き(
theta )による全エネルギーの違いを計算します。
スピン軌道相互作用計算の計算例の分子配置
O2分子
samples/spin_orbit/O2/Theta_0
samples/spin_orbit/O2/Theta_90
前者の方が、1原子あたり 0.108 meV 安定です。
Pt2分子
samples/spin_orbit/Pt2/Theta_0
samples/spin_orbit/Pt2/Theta_90
後者の方が、1原子あたり 17.215 meV 安定です。
Si 結晶のバンド
スピン軌道相互作用を考慮したバンド計算は、通常のコリニア計算と同じ手順で行います。まず、scf 計算を行い (phase)、得られた電荷密度のもとでバンド計算 (ekcal) を行います。
例として、Si diamond
結晶のバンド計算を示します。scf及びband計算ともに、以下のキーワードを用います
( samples/spin_orbit/Si_band/soc_on
フォルダ参照 )。
accuracy{
paw = on
spinorbit{
mode = pawpot
}
}
structure{
magnetic_state = noncollinear
}
図 5.92 左図は上記入力により得られたバンド構造で、右図はスピン軌道相互作用を考慮しないバンド構造です。スピン軌道相互作用によりgamma 点の縮退が解けていることが分かります。

スピン軌道相互作用を考慮したバンド構造(左)としないバンド構造(右)
使用上の注意
スピン軌道相互作用を用いてエネルギー比較を行う場合、対称操作をEのみにされることをお勧めします。
smearing においてmethod = parabolic のみが使用可能です。
ELNES / XANES解析機能
はじめに
EELSは電子が試料を透過する際のエネルギー損失を測定し、XAFSは試料にX線を照射した際に現れるスペクトルを測定する実験手法です。ここでは、電子線あるいはX線により、内殻電子が非占有状態に遷移する現象を対象とします。また、広範なエネルギー領域の中で、特に、吸収端付近のスペクトル構造を計算するものとします。この狭いエネルギー領域では、EELS は ELNES、XAFS は XANES と呼ばれています。
ELNESとXANESの二重微分散乱断面積は、双極子近似のもとで、以下のように表されます。
ここで \(|i>\) \(|f>\) は始状態 (基底状態) と終状態(ある原子の内殻電子が飛び出し、非占有状態に遷移した状態)、qは非弾性散乱時における波数の変化、eは入射X線の偏光ベクトルである。したがって、XANESの計算は、ELNES計算にてqを eと読み替えます。内殻軌道から状態n, kで指定される伝導バンドへの遷移確率は、 [Gao08] に従い、
で与えられるものとします。ここで 右辺第1項は、波動関数の soft part の寄与、第2項は擬ポテンシャルを用いたことによる補正項です。 \(|\phi_i>\) は原子軌道iの全電子波動関数、\(\tilde{\phi}_i\) は原子軌道iの擬波動関数である。
また、吸収端のエネルギーは、基底状態と内殻電子励起状態の全エネルギー差で与えられるものとします。
準備:擬ポテンシャル
公開されている擬ポテンシャルには、内殻軌道に関する情報が含まれていません。このため、ELNES・XANES の計算を行うためには、擬ポテンシャルをあらためて作成する必要があります。作成すべき擬ポテンシャルは、
基底状態の擬ポテンシャル
内殻に空孔を入れた状態の擬ポテンシャル
です。例えば、LiF 結晶の Li 原子のK 端のスペクトルを計算する場合には、
Li原子の基底状態の擬ポテンシャル
Li原子の1s軌道に空孔を入れた擬ポテンシャル
F原子の基底状態の擬ポテンシャル
を作成します。なお、いずれの擬ポテンシャル作成においても、CIAOの入力ファイルに、
sw_write_core_to_valence 1
を記入する必要があります。なお、上記は”sw_with_dipole_cor2val 1”でも結構です。キーワードについてはCIAO のマニュアルを参照してください。
計算の流れ
計算は以下の順で行います。
基底状態 (A) のSCF計算
内殻電子励起状態 (B) のSCF計算
内殻電子励起状態 (B) のスペクトル計算
①と②で得られた(A)と(B)の全エネルギー差をスペクトルの吸収端エネルギーとし、③の計算に反映します。
なお、実行コマンドは、①及び②のSCF計算は phase、③のスペクトル計算は epsmain を使用します。
基底状態 (A) のSCF 計算
入力
吸収端のエネルギーを計算するために、入力パラメータに以下のキーワードを設定します。
accuracy{
paw = on
}
Postprocessing{
CoreLevels{
sw_calc_core_energy = on
}
}
なお、file_names.data の F_POT に、5.8.2で作成した擬ポテンシャルファイル名を指定します。通常の、内殻軌道情報を含まない擬ポテンシャルを指定した場合には、該当原子の内殻電子の寄与を考慮せずにエネルギー計算が行われます。
&fnames
F_INP = './nfinp1.data'
F_POT(1) = '../pp/Li_ggapbe_paw_02_no_corehole.pp'
F_POT(2) = '../pp/F_ggapbe_paw_02_no_corehole.pp'
/
出力
core_energy.data の最終行に、内殻電子と価電子の寄与を合わせた全エネルギー値が出力されます。core_energy.data以外のファイルに出力したい場合には、file_names.data でF_CORE_ENERGY_OUTにファイル名を指定します。
# Etotal (Core+Valence)
-3437.92292869603
内殻電子励起状態 (B) のSCF 計算
入力
基底状態の計算と同じキーワードを用います。その他の注意事項として、内殻電子を励起する原子には、別の原子種を割り当てる必要があります。また、励起原子は単位胞に 1 原子のみとする。以下の例では、励起原子は Li2 です。
structure{
atom_list{
atoms{
#tag element rx ry rz mobile
Li2 0.0000000000 0.0000000000 0.0000000000
F1 0.2500000000 0.0000000000 0.0000000000
…
(中略)
…
Li1 0.7500000000 0.5000000000 0.7500000000
F1 1.0000000000 0.5000000000 0.7500000000
}
}
element_list{
#units atomic_mass
#tag element atomicnumber zeta deviation
Li1 3 0.00 1.5
Li2 3 0.00 1.5
F1 9 0.00 1.5
}
}
Li2 には、内殻電子を励起させて作成した擬ポテンシャルを割り当てる。内殻電子励起状態 (B) のSCF 計算の file_names.dataの記述例です。
&fnames
F_INP = './nfinp1.data'
F_POT(1) = '../pp/Li_ggapbe_paw_02_no_corehole.pp'
F_POT(2) = '../pp/Li_ggapbe_paw_02_1s_corehole.pp'
F_POT(3) = '../pp/F_ggapbe_paw_02_no_corehole.pp'
/
出力
基底状態の計算と同様に、core_energy.data の最終行に、内殻電子と価電子の寄与を合わせた全エネルギー値が出力されます。
内殻電子励起状態 (B) のスペクトル計算
入力
スペクトル計算の入力では、control ブロックにてcondition = fixed_charge を指定します。
control{
condition = fixed_charge
use_additional_projector=on
}
また、structure ブロックは、基底状態の計算と同様に記述します。
epsilon ブロックでは、スペクトルに関するパラメータ指定を行う。sw_corelevel_spectrum は、内殻軌道からの励起スペクトル計算を行うためのスイッチで、onと指定します。probe ブロック内のatom_id には、励起させる原子の番号を指定します。上記のLi2 は1番目の原子であるため、atom_id = 1 としています。また、orbital には、擬ポテンシャル作成時に空孔を入れた軌道を指定します。この例では orbital に 1s を指定しており、Li K端におけるスペクトルが計算されます。
fermi_energy ブロックでread_efermi = on とし、efermiには、内殻電子励起状態 (B)のSCF計算の nfefermi.data 記載の値を指定します。
energy ブロックには、吸収端を基準にして観察したいエネルギー領域を指定する。low, high, step は、それぞれ、エネルギー領域の下端、上端、及び間隔に対応する。単位は hartree です。
XANESの計算の入力パラメータは、EELS の際は photonを eels と書き換えます。なお、内部の処理は共通であるため、どちらのキーワードを用いても結果は変わりません。
epsilon {
sw_epsilon = on
sw_corelevel_spectrum = on
probe{
atom_id = 1
orbital = 1s
}
fermi_energy{
read_efermi = on
efermi = 0.21930399
}
photon{
energy{
low = -0.10, high = 2.0, step = 0.002
}
}
transition_moment{
type = ks ! {l
symmetry =on
}
BZ_integration {
method = g !{parabolic(p)
width = 1.0 eV
}
}
最後に、file_names.data では、基底状態 (A)及び内殻電子励起状態 (B) の全エネルギー値出力ファイル名を、F_CORE_ENERGY_INITIAL 及びF_CORE_ENERGY_FINAL に指定します。また、F_CHGT 及び F_CNTN_BIN_PAWに、内殻電子励起状態 (B)の電荷密度ファイルを指定ます。
&fnames
F_INP = './nfinp1.data'
F_POT(1) = '../pp/Li_ggapbe_paw_02_no_corehole.pp'
F_POT(2) = '../pp/Li_ggapbe_paw_02_1s_corehole.pp'
F_POT(3) = '../pp/F_ggapbe_paw_02_no_corehole.pp'
F_CHGT = '../scf_excited/nfchgt.data'
F_CNTN_BIN_PAW = '../scf_excited/continue_bin_paw.data'
F_CORE_ENERGY_INITIAL = '../scf_ground/core_energy.data'
F_CORE_ENERGY_FINAL = '../scf_excited/core_energy.data'
/
出力
スペクトルの出力ファイル名はeps.data である。以下に出力例を示す。3行目のエネルギー値は、吸収端エネルギーにe_lowを加えた値です。spectrum はスペクトル強度、すなわち誘電関数の虚部 に対応します。
# Spectrum data
# Energy[eV] Spectrum
0.6018296869E+02 0.7167942086E-06
0.6023739146E+02 0.1031555487E-05
0.6029181423E+02 0.1478755242E-05
…
計算事例
LiF 結晶の Li 原子K端のスペクトル
例題として、LiF 結晶の Li 原子K端の XANES
スペクトルを示します。計算事例は、 samples/XANES/LiF
です。計算は、
基底状態 (A) のSCF計算 ( フォルダ名:scf_ground )
内殻電子励起状態 (B) のSCF計算 ( フォルダ名:scf_excited )
内殻電子励起状態 (B) のスペクトル計算 ( フォルダ名:eps_excited )
の順で行います。
スペクトル計算はk 点数を増やして行うことを推奨します。

(左) LiF結晶。(右)スペクトル
四重極子成分 (バージョン2019.02以降)
概要
XANESのスペクトル計算では、始状態iと終状態fの間の遷移モーメントを計算します。一般に、四極子成分は弱いことが多いですが、双極子極子遷移が禁止されているエネルギー領域に明瞭なピークを生じることがあります。この四重極成分の計算をPHASE/0で実行する方法を説明します。
計算理論
遷移モーメントの計算
四重極まで含めた、始状態iと終状態fの間の遷移モーメントは
で表されます [Gougoussis09] 。\(\psi_i\) 及び \(\psi_f\) は、それぞれ始状態i及び終状態fの波動関数、\(\hat{\mathbf{\varepsilon}}\) 及び \(k\) は入射光の偏光ベクトル及び波数ベクトルです。\(D\) の第1項が双極子、第2項が四重極子成分に対応します。
ただし、PAW法では、波動関数が規格直交化されないため、は以下のように書き替えられます。
ここで、 \(\beta_{\mathbf{R},l'm'\tau'}^v\) 及び \(\beta_{\mathbf{R},nlm}^c\) は、それぞれ価電子及び内殻電子軌道へのプロジェクタです。また、 \(\phi_{\mathbf{R},l'm'\tau'}^v\) 及び \(\tilde{\phi}_{\mathbf{R},l'm'\tau'}^v\) は、それぞれ価電子軌道の全電子及び擬波動関数です。同様に、 \(\phi_{\mathbf{R},nlm}^c\) 及び \(\tilde{\phi}_{\mathbf{R},nlm}^c\) は、それぞれ内殻電子軌道の全電子及び擬波動関数です。
なお、XANESの計算は、始状態iとして、原子位置 (R0) 及び内殻電子軌道 (n,l) を指定して行います。さらに、この軌道は原子近傍に局在していることを考えると、結局、
を計算すればよいことがわかります [Gougoussis09] 。
スペクトル強度の計算
スペクトル強度は、
で与えられます [Cabaret10] 。ここで、 \(E_i\) 及び \(E_f\) は、始状態及び終状態の全エネルギー、 \(\hbar \omega\) は入射光のエネルギーです。(5.76) を用いれば、双極子及び四重極子成分は、それぞれ
となります。なお、E1及びE2は電気双極子、四重極子成分という意味です。
入力
内殻電子励起前後の SCF 計算は済んでいるものとします。四重極子計算機能を利用するには、以下のような入力を作成します。影を付けた箇所が、本機能特有の項目です。sw_local_approx_trans_moment=on とすると、双極子成分に加えて四重極子成分を計算します。wave_vector では、入射光の波数ベクトルの方向を指定します。
四重極子成分を計算するための入力例
epsilon{
sw_epsilon = on
sw_corelevel_spectrum = on
sw_local_approx_trans_moment = on
probe{
atom_id = 2
orbital = 1s
}
photon{
polarization{
ux = 1, uy = 0, uz = 0
}
wave_vector{
kx = 0, ky = 0, kz = 1
}
energy_range{
e_low= -0.10, e_high=2.0, e_step=0.002
}
}
BZ_integration {
method = g !{parabolic(p)
width = 0.2 eV
}
fermi_energy{
read_efermi = on
efermi = 0.17409189
}
}
なお、fermi__energy ブロックは無くてもかまいませんが、このときは、必ず、内殻電子励起状態のSCF計算時に生成されたフェルミ準位ファイル名を file_names.data 内に明記する必要があります。
fermi_energy ブロックを用いない場合の file_names.dataの例
&fnames
(中略)
F_EFERMI = '../scf-final/nfefermi.data'
(中略)
/
出力
双極子及び四重極子成分のスペクトルは、それぞれ、eps_E1_E1.data及びeps_E2_E2.data に出力されます。
eps_E1_E1.data の出力例
# Core level spectrum ( dipole-dipole )
# atom: 2, orbital: n = 1, l = 0
#
# Energy [eV] spectrum
#
(中略)
0.4968388438E+04 0.8211847200E-07
0.4968442860E+04 0.9903378493E-07
0.4968497283E+04 0.1147357290E-06
(後略)
eps_E2_E2.data の出力例
# Core level spectrum ( quadrupole-quadrupole )
# atom: 2, orbital: n = 1, l = 0
#
# Energy [eV] spectrum
#
(中略)
0.5008878978E+04 0.1286045986E-09
0.5008933400E+04 0.1328488558E-09
0.5008987823E+04 0.1279594782E-09
(後略)
計算例
SrTiO3結晶の計算例を紹介します。入力ファイルは samples/XANES/SrTiO3
にあります。
まず、単位胞 (5原子) の格子定数を、k点: Monk(8×8×8) で最適化しa=3.923Åを得ましたた。これを各結晶軸方向に3倍したスーパーセル (135原子、a=11.769Å) を作成し、1つのTi原子の擬ポテンシャルを1sに内殻正孔を持つものに置き換えました。計算条件は以下の通。
cutoff_wf |
25 Ry |
cutoff_cd |
225 Ry |
k点サンプリング |
SCF: Monk (2×2×2) XANES: Monk (4×4×4) |
汎関数 |
PAW, GGAPBE |
擬ポテンシャル |
Sr_ggapbe_paw_02.pp Ti_ggapbe_paw_002_nocorehole.gncpp2 Ti_ggapbe_paw_005_1s_corehole.gncpp2 O_ggapbe_paw_02.pp |
スペクトルのブロードニング |
Gaussian, width = 0.5 eV |
また、偏光ベクトル、波数ベクトルの組み合わせとして、以下の2つを試しました。四重極子遷移において、aはt2g, bはeg 軌道成分をプローブします。
a |
b |
|
---|---|---|
偏光ベクトルEの方向 |
[100] |
[110] |
波数ベクトルkの方向 |
[010] |
[1-10] |
結果を以下に示します。 SrTiO3のd軌道は、t2gの方がeg軌道よりもエネルギーが低いため、四重極子遷移で最初にピークが現れるのは t2g 軌道をプローブするa の場合でです ( 図 5.94 参照)。図 5.95 は、吸収端近傍の拡大図です。
SrTiO3 のXANESスペクトル
SrTiO3 のXANESスペクトル拡大図。
文献
- Gao08
-P. Gao, C. J. Pickard, M. C. Payne, J. Zhu, and J. Yuan, Phys. Rev. B 77 (2008) 115122.
- Gougoussis09(1,2)
Gougoussis, M. Calandra, A. P. Seitsonen, and F. Mauri, Phys. Rev. B 80 (2009) 075102.
- Cabaret10
Cabaret, A. Bordage, A. Juhin, M. Arfaoui, and E. Gaudry, Phys. Chem. Chem. Phys. 12 (2010) 5619.
化学ポテンシャル一定のシミュレーション
概要
PHASE/0による通常のシミュレーションでは,電子数一定の計算を行います。これに対し,化学ポテンシャル(フェルミエネルギー)を一定とし,構造最適化や分子動力学シミュレーション,NEB計算などを行う [Bonnet12] ことができる機能(constant-mu法)もPHASE/0には備わっています。この場合,電子数はシミュレーション中変化します。
入力パラメーター
化学ポテンシャル一定のシミュレーションを行うには,accuracyブロックの下にfcpブロックを作成し,次の例のように設定を行います。
accuracy{
...
fcp{
sw_fcp = ON
mu = -0.1
relax_crit = 1.0d-5
}
...
}
fcpブロックにおいて定義可能なパラメータは下記の通りです。
sw_fcp |
化学ポテンシャル一定の計算を行うかどう かを指定するスイッチです。 行う場合にonとします。 |
---|---|
mu |
ターゲットと なる化学ポテンシャルの値をエネルギーの単位で指定 します。事前に通常の計算で状態密度を計算しておき ,フェルミエネルギー近辺のエネルギー固有値の分布 や電子数などを調べておくと設定しやすいでしょう。 |
temperature |
化学ポテンシャルを 制御するための"粒子浴"の温度を指定します。分子動 力学シミュレーションの場合に意味のある設定です。 |
qmass |
化学ポテンシャルを 制御するための"粒子浴"の質量を指定します。分子動 力学シミュレーションの場合に意味のある設定です。 |
mass |
化学ポテンシャル を制御するための"電子の質量"を指定します。分子動 力学シミュレーションの場合に意味のある設定です。 |
relax_crit |
構造最適化の際の収束判定条件を指定し ます。得られる化学ポテンシャルの値とターゲットの 値の差の絶対値がここで指定する値よりも小さくなっ た場合に収束したと見なされます(これ以外に,通常 の原子間力に関する収束判定条件も考慮されます)。 |
tot_charge_first |
NEB計算 の際に,始点のレプリカに与える電荷を指定します。 |
tot_charge_last |
NEB計算 の際に,終点のレプリカに与える電荷を指定します。 |
なお,本計算機能を使って(ジョブ1とする)得られた(中性でない)電荷を外部電荷としてあたえて通常の計算を行う(ジョブ2とする)場合,(ジョブ1において)与えたmuの値とは異なるフェルミエネルギーが(ジョブ2において)得られます。逆に,外部電荷を与えて通常の計算を行って(ジョブ3とする)得られたフェルミエネルギーをmuに指定して(ジョブ4によって)最適化を行うと,得られる電荷はもとの(ジョブ3における)外部電荷とは異なる電荷となります。このようにつじつまの合わない結果が得られるのは,中性でない電荷を与える場合は初期電荷の与え方が中性の場合と異なるので,エネルギーの原点が変化するためです。つじつまの合った結果を得るためには,外部電荷を与えた計算について以下のような設定を加え,初期電荷の計算方法を合わせるようにしてください。
accuracy{
...
sw_add_qex_to_initial_charge = off
...
}
計算の実行
通常の構造最適化もしくは分子動力学シミュレーションの設定に加え,化学ポテンシャル一定のシミュレーションの設定を施したら,通常通りPHASE/0を実行すれば計算を行うことができます。
計算中に電荷がどのように変化したかは,以下の要領で調べることができます。
% grep 'Total Charge' output000
FCP : Total Charge = 31.98942095
FCP : Total Charge = 31.99795170
FCP : Total Charge = 32.01254363
FCP : Total Charge = 32.02605805
FCP : Total Charge = 32.03237025
FCP : Total Charge = 32.02886985
FCP : Total Charge = 32.01742419
FCP : Total Charge = 32.00372842
FCP : Total Charge = 31.99503316
...
...
また,化学ポテンシャル(フェルミエネルギー)がどのように変化したかは,以下の要領で調べることができます。
% grep 'Fermi Energy' output002
FCP : Fermi Energy = 0.24621583
FCP : Fermi Energy = 0.24679898
FCP : Fermi Energy = 0.24700415
FCP : Fermi Energy = 0.24674191
FCP : Fermi Energy = 0.24618985
...
...
例題
非常に単純な例題として,シリコン8原子の系の分子動力学シミュレーションを取り上げます。この例題の入力ファイルは, samples/dynamics/FCP/NVT_nose_hoover
にあります。
この例題の入力ファイルは,基本的には samples/dynamics/molecular_dynamics/NVT
以下にあるものと同等ですが,以下のように“化学ポテンシャル一定の分子動力学シミュレーション”を行うための設定が施されています。
accuracy{
...
fcp{
sw_fcpopt = ON
mu = 9.0e-3
mass = 1000.0d0
qmass = 4000
}
...
}
この入力ファイルを利用して,通常通りPHASE/0を実行すれば化学ポテンシャル一定の第一原理分子動力学シミュレーションを行うことが可能です。その結果,たとえば 図 5.96 に示すように電荷がシミュレーション中時々刻々と変化します。
分子動力学シミュレーション中の電荷の変化。
使用上の注意
本計算機能は,すべての機能と組み合わせて利用することが可能ですが,本計算機能を利用すると電荷中性ではない計算を行うことになる点には注意が必要です。有効遮蔽体法(ESM法。5.3.6 章 参照)で境界条件をpe1(両側のESMが金属)とすると,隣り合う単位胞の相互作用の影響を取り除くこともできます。
参考文献
- Bonnet12
Bonnet, T. Morishita, O. Sugino, and M. Otani, “First-Principles Molecular Dynamics at a Constant Electrode Potential”, Physical Review Letters 109 266101 (2012).
有限電場計算機能 (バージョン2019.02以降)
概要
有限電場計算機能とは,一様有限電場下において分極と電場のエネルギーを通常のエネルギーに加えた“エンタルピー”を極小化する計算機能です。この計算機能によって,印可した電場に対する応答としての巨視的な分極を求めることができます。得られた分極から誘電率を計算することや,有限電場下での原子に働く力よりボルン有効電荷を得ることも可能となっています。
計算理論
一様有限電場下の計算は, [Souza02] において提案されている処方箋に従って行われます。一様電場下におけるエンタルピーを,以下のように定義します。
\(E_{\rm KS}\) が通常のコーン・シャムハミルトニアンから得られる全エネルギー,\(\Omega\) が単位胞の体積,\(\mathbf{P}_{\rm mac}\) が巨視的な電気分極,\(\mathbf{E}\) が一様有限電場です。\(\mathbf{P}_{\rm mac}\) はベリー位相分極理論から計算することができます。\(\mathbf{E}\) は入力パラメーターです。このエンタルピーを極小化するために必要な微分は下記のように計算できます。
\(\ket{u_{n\mathbf{k}}}\) が周期的な波動関数です。\(\mathbf{k}_+, \mathbf{k}_-\) は,逆格子ベクトルをbiとするとそれぞれk+bi/Ni , k-bi/Niです\(S_{mn}\left( \mathbf{k},\mathbf{k}' \right)\) は波動関数の重なりであり,\(\Braket{u_{m\mathbf{k}} | u_{n\mathbf{k}'}}\) と表すことができます。
このようにして計算された勾配を利用してエンタルピーを最適化すると,一様電場下での巨視的な分極 \(\mathbf{P}_{\rm mac}\) を計算することができます。この量から,さらに様々な誘電応答解析を実施することも可能です。たとえば,電場を印可した状態の分極から差分によって誘電率を得ることが可能です(微弱な電場 \(\delta E\) を与えた際に発生する分極 \(\delta P\) から,誘電率は \(\varepsilon = 1+4\pi \frac{\delta P}{\delta E}\) と計算することができる)。あるいは,有限電場下での原子に働く力から \(eZ = \frac{\delta F}{\delta E}\) という関係を用いてボルン電荷を得ることも可能です。
入力
一様有限電場下の計算をPHASE/0で実行するには,controlブロックの下でsw_fef = onとすることによって一様有限電場機能を利用することを指定し,さらにaccuracyブロックの下のfefブロックの下のelectric_filedブロックにおいて電場を入力します。具体的には以下のように指定します。
control{
sw_fef = on
}
...
...
accuracy{
fef{
electric_field{
ex = 1e-5
ey = 1e-5
ez = 1e-5
}
}
}
...
...
変数ex, ey, ezによって電場のx, y, z成分を指定します。また,有限電場法を利用する場合はk点サンプリング手法に制限があります。以下のようにkshiftパラメーターをすべて0と設定してください(methodは無指定もしくはmonkとする)
ksampling{
kshift{ k1=0, k2=0, k3=0 }
}
計算ソルバーについて
有限電場法においては利用可能な計算ソルバーに本質的な制限があります。これは,有限電場法は非周期的な項 ( \(\mathbf{P}_{\rm mac} \cdot \mathbf{E}\) )を追加するため,波動関数がハミルトニアンの固有状態ではなくなってしまうことに起因します。自由にユニタリ変換することができないため,たとえば部分空間対角化は利用できません(したがって部分空間対角化が必須であるDavidson系ソルバーはすべて利用不可です)。検証の結果「部分空間対角化を伴わないCG法」が最も収束の速い手法であることが分かったので,有限電場計算機能利用時はこのソルバーがデフォルトになります。
出力
PHASEを実行すると,ログファイル(output000ファイル)には以下の要領で各SCFイテレーションにおける電気分極の情報が出力されます。
BP = 0.006502276 0.006502276 0.006502276
Pel = 0.000078647 0.000078647 0.000078647
Pion = 0.000000000 0.000000000 0.000000000
Pmac = 0.000078647 0.000078647 0.000078647
BP = の後にベリー位相が,Pelの後に分極の電子からの寄与が,Pionのあとに分極のイオンからの寄与が,Pmacに電子とイオンの分極の合計が出力されます。単位はすべて原子単位です。得られた分極から,誘電率やボルン電荷を以下の要領で計算することができます。
誘電率 \(\varepsilon = 1+4\pi \frac{\delta P}{\delta E}\)
ボルン電荷 \(eZ = \frac{\delta F}{\delta E}\)
計算例
AlP結晶の計算を有限電場法および従来法で行いました。その入力ファイルは samples/FEF/AlP
にあります。8原子の系を採用しました。計算結果を以下にまとめます。
計算結果
物性値 |
計算結果 |
---|---|
電子分極 |
0.000050916 |
誘電率(有限電場法) |
7.40 |
誘電率(従来法) |
7.40 |
ボルン電荷(有限電場法) |
Al : 2.25, P : -2.27 |
ボルン電荷(従来法) |
Al : 2.227, P : -2.227 |
従来法とよく一致する結果が得られました。ボルン電荷の方も従来法とよく一致しています。有限電場法のボルン電荷の和がゼロになっていないのは,数値誤差によるものと考えられます。
注意点
有限電場機能を利用する場合,まずは通常のSCF計算を実施し,収束した波動関数と電荷密度を得ておくことが推奨されます。したがって,計算の手続きは以下のようになります。
対象としたい系のSCF計算を行う。
有限電場計算用のディレクトリーを作成し,そこに1. で利用した入力ファイルと得られた波動関数データと電荷密度データをコピーする(ファイル名はそれぞれzaj.dataとnfchgt.data)
入力パラメーターファイルに有限電場計算用の入力の設定をほどこし,さらに初期波動関数・初期電荷密度はファイルから読み込むように設定する。具体的には下記のような記述を行う。
分極を収束させるためにはSCF計算の収束判定として通常よりも厳しいものを採用する必要があることが多いです。Pelの収束具合を確かめながら必要な収束判定条件を決めていくことが推奨されます。上述のようにinitial_wavefunctionsとinitial_charge_densityをfileに設定しておけばそれまでの計算の結果を読み込んで動作するので,それまで行った計算は無駄にはなりません。
原子間力の計算結果からボルン電荷を得ることができますが,この機能はノルム保存擬ポテンシャルのみ対応しています。ウルトラソフト擬ポテンシャルの場合に得られる原子間力は不正確なので利用しないでください。
参考文献
- Souza02
Souza, J. Íñiguez and D. Vanderbilt, “First-Principles Approach to Insulators in Finite Electric Fields” Physical Review Letters vol. 89 (2002) pp. 117602 1-4.