解析機能

部分電荷密度

部分電荷密度を計算するにはタグPostprocessingの中のタグchargeで指定します。タグchargeの中の変数sw_charge_rspaceとタグpartial_chargeの中の変数sw_partial_chargeをOnにします。

Postprocessing{
  charge{
    sw_charge_rspace = on
    partial_charge {
      sw_partial_charge = on
      Erange_min = -0.45 eV
      Erange_max = 0.45 eV
      Erange_delta = 0.05 eV
      partial_charge_filetype=separate !{individual,separate | integrated}
    }
  }
}

変数Erange_minとErange_maxにエネルギー領域の最大値と最小値を入力します。エネルギーは金属の場合フェルミレベルから測り、絶縁体の場合は価電子帯上端のエネルギーから測ります。変数Erange_deltaに入力した値の間隔のエネルギー窓を先のエネルギー領域に作成します。

出力ファイルoutput000には以下の様にエネルギー窓に関する出力があります（”!pc”で始まる行）。

!pc nEwindows =   20, nvb_windows =   10, ncb_windows =   10 <<m_ESoc_set_nEwindows_pc>>
!pc    iw  if_elec_state                erange(hartree)                        erange(eV)
!pc              (asis)                  (shifted)             (shifted)
!pc     1       1       ( 0.094537  0.096374 ) ( -0.018375 -0.016537 ) ( -0.500000 -0.450000 )
!pc     2       1       ( 0.096374  0.098211 ) ( -0.016537 -0.014700 ) ( -0.450000 -0.400000 )
!pc     3       1       ( 0.098211  0.100049 ) ( -0.014700 -0.012862 ) ( -0.400000 -0.350000 )
!pc     4       1       ( 0.100049  0.101886 ) ( -0.012862 -0.011025 ) ( -0.350000 -0.300000 )
!pc     5       0       ( 0.101886  0.103724 ) ( -0.011025 -0.009187 ) ( -0.300000 -0.250000 )
!pc     6       1       ( 0.103724  0.105561 ) ( -0.009187 -0.007350 ) ( -0.250000 -0.200000 )
!pc     7       1       ( 0.105561  0.107399 ) ( -0.007350 -0.005512 ) ( -0.200000 -0.150000 )
!pc     8       0       ( 0.107399  0.109236 ) ( -0.005512 -0.003675 ) ( -0.150000 -0.100000 )
!pc     9       0       ( 0.109236  0.111074 ) ( -0.003675 -0.001837 ) ( -0.100000 -0.050000 )
!pc    10       1       ( 0.111074  0.112911 ) ( -0.001837  0.000000 ) ( -0.050000  0.000000 )
!pc    11       1       ( 0.112911  0.114749 ) (  0.000000  0.001837 ) (  0.000000  0.050000 )
!pc    12       0       ( 0.114749  0.116586 ) (  0.001837  0.003675 ) (  0.050000  0.100000 )
!pc    13       0       ( 0.116586  0.118424 ) (  0.003675  0.005512 ) (  0.100000  0.150000 )
!pc    14       0       ( 0.118424  0.120261 ) (  0.005512  0.007350 ) (  0.150000  0.200000 )
!pc    15       0       ( 0.120261  0.122099 ) (  0.007350  0.009187 ) (  0.200000  0.250000 )
!pc    16       1       ( 0.122099  0.123936 ) (  0.009187  0.011025 ) (  0.250000  0.300000 )
!pc    17       1       ( 0.123936  0.125773 ) (  0.011025  0.012862 ) (  0.300000  0.350000 )
!pc    18       0       ( 0.125773  0.127611 ) (  0.012862  0.014700 ) (  0.350000  0.400000 )
!pc    19       0       ( 0.127611  0.129448 ) (  0.014700  0.016537 ) (  0.400000  0.450000 )
!pc    20       0       ( 0.129448  0.131286 ) (  0.016537  0.018375 ) (  0.450000  0.500000 )

nEwindowsはエネルギー窓の総数です。nvb_windowsとncb_windowsはそれぞれ価電子状態と伝導電子状態を含むエネルギー窓の数です。この設定例では、Erange_min= -0.45 eV、Erange_max = 0.45 eV、Erange_delta = 0.05 eVなので、エネルギー窓の数は18個になるはずですが、実際には低エネルギー側および高エネルギー側にそれぞれErange_delta幅の袖領域を設けて出力するのでエネルギー窓の総数が20個になっていることに、ご注意下さい。iwはエネルギー窓の番号です。if_elec_stateはそのエネルギー窓に電子状態があるかどうかを示しています。この値が0の時は電子状態がなく、1の時には電子状態が存在します。列asisには原子単位でエネルギー窓の範囲が示されています。二つの列shiftedにはエネルギーの基準から測ったときのエネルギー窓の範囲が原子単位とeV単位で示されています。

変数partial_charge_filetypeにindividualまたはseparateを指定すると、各エネルギー窓ごとに計算された電荷密度が番号付けされたファイルに出力されます。その際の名前の付け方は、スピン分極がない場合であれば、F_CHR = nfchr.cubeに対してnfchr.00xx.cube（xxには上の表のiwの値が入る）というようになります。スピン分極がある場合には、F_CHR = nfchr.cubeに対して、nfchr.up.00xx.cube、nfchr.down.00xx.cubeの二種類のファイルが生成されます。上の表でif_elec_stateが0になっているのは、その範囲に固有値がある状態がないことを示しています。その場合、cubeファイルは生成されません。

integratedを選択すると各電荷密度がひとつのファイルに追記され、各電荷密度データの先頭にはPARTIALCHARGEが記述され、終わりにはENDが記述されます。

BaO/Si(001)界面の部分電荷密度を計算した結果を advanced_BaO_image3 に示します。

$BaO/Si(001)界面構造の部分電荷密度。(a)BaO/Si(001)界面構造のモデル図。(b)フェルミレベル直下(固有エネルギーが-0.05eVから0.0eVまで)の電子状態の部分電荷密度。(c)フェルミレベル直上(固有エネルギーが0.0eVから0.05eVまで)の電子状態の部分電荷密度。電子密度は\mathbf{1} \times \mathbf{10}^{- \mathbf{5}}から\mathbf{1} \times \mathbf{10}^{- \mathbf{3}}までが示されています。青い部分には電子が少なく、赤い分部には電子が多くなっています。$

BaO/Si(001)界面構造の部分電荷密度。(a)BaO/Si(001)界面構造のモデル図。(b)フェルミレベル直下(固有エネルギーが-0.05eVから0.0eVまで)の電子状態の部分電荷密度。(c)フェルミレベル直上(固有エネルギーが0.0eVから0.05eVまで)の電子状態の部分電荷密度。電子密度は1 × 10^− 5から1 × 10^− 3までが示されています。青い部分には電子が少なく、赤い分部には電子が多くなっています。

応用例：STM像の解析

部分電荷密度出力機能を利用すると、STM像を模擬することが可能です。解析したいバイアスポテンシャルに対応したエネルギーウィンドウの部分電荷密度を、表面からある程度離れた平面に投影した像が計算上のSTM像です。以下、サンプルデータ( samples/stm_by_pcharge 以下)を利用して計算方法を具体的に説明します。

サンプルは、Si の(001)面に相当するデータです。通常のPHASE入力に、以下のように部分電荷密度出力の設定を加えています。

postprocessing{
  charge{
    sw_charge_rspace = on
    filetype = cube
    partial_charge{
      sw_partial_charge = on
      partial_charge_filetype = individual
      Erange_min = 0 eV
      Erange_max = 0 eV
      Erange_delta = 1 eV
    }
  }
}

このように設定することによって、フェルミエネルギーからみて-1 eV から0 eVまでのデータと0 eV から1 eVのエネルギーウィンドウの部分電荷密度が出力されます。それぞれ、-1 V（占有状態）および1 V（非占有状態）のバイアスポテンシャルに対応したSTM像が得られます。この入力データを利用して計算を実施すると、nfchr.0001.cube (-1 eV から0 eV の電荷密度データファイル) とnfchr.0002.cube (0 eV から1 eV の電荷密度データファイル) が作成されます。それぞれ、表面から5 Å程度離れた地点でのコンター図を占有状態について advanced_STM_image1 (a)に、非占有状態について advanced_STM_image1 (b)に示します。

Si (100)面のSTM像、(a) 占有状態の像、(b)非占有状態の像。

ベーダ―解析向け電荷密度ファイルの出力 (バージョン2019.02以降)

概要

PHASE/0における実空間における電荷密度出力機能では、通常価電子部分のみが出力されますが、ベーダ―解析では、内殻電子の寄与も考慮することにより、解析精度を上げることが出来ます。そこで、CUBE形式ファイルの出力の際に、内殻電子の寄与を加える機能が搭載されています。

ベーダ―解析について

ベーダ―解析とは、CUBE形式のファイルから各原子に割り当てる電子数を求める解析手法です。PHASE/0にベーダ―解析を行う機能そのものは備わっていませんが、ベーダ―解析はbaderプログラム

http://theory.cm.utexas.edu/henkelman/code/bader/

をダウンロードし、コンパイルすることによって簡単に行うことができます。詳しくはbaderプログラムのドキュメントなどを参照してください。

処理内容

価電子密度及び内殻電子密度を、それぞれ $\rho_{\rm v}(\mathbf{r})$ 及び $\rho_{\rm c}(\mathbf{r})$ とします。CUBEファイルに出力する全電子密度は、両者の合計です。さて、$\rho_{\rm c}(\mathbf{r})$ は、サイトIにある原子種 α の内殻電子密度 $\rho_{\rm c}^{\alpha} (\mathbf{r}-\mathbf{R}_I)$ を足し合わせたものである。

$$\rho_{\rm c} (\mathbf{r}) = \sum_{I} \rho_{\rm c}^{\alpha} (\mathbf{r} - \mathbf{R}_I)$$

ここで、は球対称であり、擬ポテンシャルファイル中に出力されています。これを読みこみ、以下の処理をします。

逆格子空間で処理する場合

$$\rho_{\rm c} (\mathbf{G}) = \sum_{I} \rho_{\rm c}^{\alpha} ( \mathbf{G} ) \exp (-i \mathbf{G} \cdot \mathbf{R}_I )$$

を価電子密度 $\rho_{\rm v} (\mathbf{G})$ に加えた $\rho_{\rm tot} (\mathbf{G})$ を計算し、フーリエ変換により $\rho_{\rm tot} (\mathbf{r})$ を求める。

実空間で処理する場合

まず、にフーリエ変換を施しを得る。次に、実空間のメッシュ点 jで

$$\rho_{\rm c} (\mathbf{r}_j) = \sum_{I} \rho_{\rm c}^{\alpha} ( \mathbf{r}_j - \mathbf{R}_I )$$

を計算し、 $\rho_{\rm v} (\mathbf{r}_j)$ に加えて $\rho_{\rm tot} (\mathbf{r}_j)$ を得る。これを全てのメッシュ点に対して行えばよい。なお、内殻電子密度は原子核近傍に局在しているので、メッシュ点 jは、原子位置 R_I から適当な距離 rcut以内にあるものを選ぶ。

入力

filetype = cube 及びsw_add_corecharge_rspace = on とすると、(擬ポテンシャルから読み込んだ) 内殻電子密度を加えた全電子密度分布が、CUBE ファイルとして出力されます。また、eval_corecharge_on_Gspace =on (off) のとき、逆格子空間 (実空間) で内殻電子密度を価電子密度に加えます。

postprocessing{
  charge{
    sw_charge_rspace = on
    filetype = cube
    sw_add_corecharge_rspace = on ( デフォルト：off )
    eval_corecharge_on_Gspace = off ( デフォルト：off )
  }
}

上述の設定はPAW対応擬ポテンシャルを利用する場合に利用できます。PAW非対応の古い擬ポテンシャルを用いる場合、次のように設定します。

postprocessing{
  charge{
    sw_charge_rspace = on
    filetype = cube
    sw_add_corecharge_rspace = on
    sw_read_corecharge_extra_file = on
  }
}

すなわち、パラメーター sw_read_corecharge_extra_file を on にすることによって別の擬ポテンシャルから内殻の情報を読み込みます。読み込むファイル名は、 file_names.data 内で F_CORE_CHARGE(n) で指定します。

&fnames
...
F_CORE_CHARGE(1) = 'Si_ggapbe_paw_nc_01m.pp'
...
/

なお、 F_POT(n) で指定した擬ポテンシャルファイルが内殻電子密度を含む場合には、 F_CORE_CHAEGE(n) で指定したファイルは無視されます。

出力

出力ファイル名は、価電子密度のCUBEファイル名を元に “_ae” を加えた名称です。例えば、価電子密度CUBEファイル名が “nfchr.data” の場合、全電子密度CUBEファイルは ”nfchr_ae.data” となります。

計算例

GaNの計算例

GaNの例を紹介します。入力ファイルは samples/Bader/GaN 以下にあります。

格子定数や原子位置はあらかじめ最適化したものを採用しました。

bader プログラムを用いる際、以下の評価法1-3を試しました。結果を advanced_bader_table1 に示す。

評価法1: bader nfchr.data ( 価電子密度を利用、従来法 )

評価法2: bader nfchr_ae.data ( 全電荷密度を利用)

評価法3: bader nchr.data -ref nfchr_ae.data

(切断する位置は全電荷密度で決定、積分は価電子密度を利用)

GaNのベーダ―電荷

	Ga	N	評価
Elk
	plot3dのメッシュ：160x160x250
*	1.34,, 1.59	-1.55

PHASE/0	eval_corecharge_on_Gspace = off
	cutoff_cd = 270Ry (CD_FFTメッシュ：32x32x108)
評価法1	1.51	-1.51	OK
*2	3.15	-1.44	NG
3	1.47	-1.47	OK

	cutoff_cd = 2700Ry (CD_FFTメッシュ：108x108x80)
評価法1	1.58	-1.58	OK
*2	0.78	-1.51	NG
3	1.51	-1.51	OK

	cutoff_cd = 6150Ry (CD_FFTメッシュ：160x160x250)
評価法1	1.58	-1.58	OK
*2	1.58	-1.52	OK
3	1.51	-1.51	OK

PHASE/0	eval_corecharge_on_Gspace = on
	cutoff_cd = 270Ry (CD_FFTメッシュ：32x32x108)
評価法1	1.51	-1.51	OK
2	1.35, 1.36	-1.36, -1.35	NG
3	1.28	-1.28	NG

	cutoff_cd = 2700Ry (CD_FFTメッシュ：108x108x80)
評価法1	1.58	-1.58	OK
2	1.22	-1.22	NG
3	1.22	-1.22	NG

	cutoff_cd = 6150Ry (CD_FFTメッシュ：160x160x250)
評価法1	1.58	-1.58	OK
2	1.21	-1.21	NG
3	1.21	-1.21	NG

*電荷密度の積分値が全電荷と一致しない点に注意。

eval_corecharge_on_Gspace = onの場合、全電子密度の積分値がメッシュサイズによらず整数値になるものの、原子近傍に微少な電荷(ごみ)が現れてしまい、bader がこれを検知するため精度が悪いようです。

一方、eval_corecharge_on_Gspace = offの場合、かなりcutoff_cdを上げないと、全電子密度の積分値は、全電荷と一致しません。しかし、評価法3を用いることで、それなりの精度でベーダ―電荷が求められることがわかりました。

4H-SiCの計算例

4H-SiCの例を紹介します。入力ファイルは samples/Bader/4H-SiC 以下にあります。

格子定数や原子位置はあらかじめ最適化したものを採用しました。

GaNの場合と同様、bader プログラムを用いる際、以下の評価法1-3を試しました。結果を advanced_bader_table2 に示します。

評価法1: bader nfchr.data ( 価電子密度を利用、従来法 )

評価法2: bader nfchr_ae.data ( 全電荷密度を利用 )

評価法3: bader nchr.data -ref nfchr_ae.data

(切断する位置は全電荷密度で決定、積分は価電子密度を利用)

4H-SiCのベーダ―電荷

	Si	C	評価
Elk
	plot3dのメッシュ：150×150×480
※	2.66	-2.70, -2.69

PHASE/0	eval_corecharge_on_Gspace = off
	cutoff_cd = 270Ry (CD_FFTメッシュ：32×32×108)
評価法1	4.00	-4.00	×
※2	-4.27, 2.88	-2.66, -2.55	×
3	2.56, 2.58	-2.58, -2.56	○

	cutoff_cd = 2700Ry (CD_FFTメッシュ：100×100×320)
評価法1	4.00	-4.00	×
※2	2.59, 2.67	-2.68, -2.67	○
3	2.65, 2.66	-2.66, -2.65	○

	cutoff_cd = 6150Ry (CD_FFTメッシュ：150×150×480)
評価法1	4.00	-4.00	×
※2	2.66, 2.68	-2.69, -2.68	○
3	2.66, 2.67	-2.67, -2.66	○

PHASE/0	eval_corecharge_on_Gspace = on
	cutoff_cd = 270Ry (CD_FFTメッシュ：32×32×108)
評価法1	4.00	-4.00	×
2	0.94, 1.03	-0.92, -1.03	×
3	0.89, 1.06	-0.94, -1.01	×

	cutoff_cd = 2700Ry (CD_FFTメッシュ：100×100×320)
評価法1	4.00	-4.00	×
2	1.18, 1.24	-1.22, -1.19	×
3	1.18, 1.24	-1.22, -1.19	×

	cutoff_cd = 6150Ry (CD_FFTメッシュ：150×150×480)
評価法1	4.00	-4.00	×
2	1.20, 1.24	1.21, 1.23	×
3	1.20, 1.24	1.21, 1.23	×

※電荷密度の積分値が全電荷と一致しない点に注意。

GaNの場合と同様に、eval_corecharge_on_Gspace = onの場合、全電子密度の積分値がメッシュサイズによらず整数値になるものの、原子近傍に微少な電荷(ごみ)が現れてしまい、bader がこれを検知するため精度が悪いようです。

一方、eval_corecharge_on_Gspace = offの場合、かなりcutoff_cdを上げないと、全電子密度の積分値は、全電荷と一致しません。しかし、評価法3を用いることで、それなりの精度でベーダ―電荷が求められています。また、GaNの場合と異なり、価電子密度のみを用いる評価法1は期待される結果を出力しておらず、注意が必要です。

ストレステンソル

機能の概要

PHASEには、ストレステンソルを計算する機能があります。ストレステンソルを計算することにより、安定な格子定数や弾性定数を計算することができます。

入力パラメータ

ストレステンソルを計算するには、入力パラメータファイルnfinp.dataにおいて、structure_evolutionブロックの下のstressブロックで、ストレステンソル計算を有効にする指定をします。

Si（立方晶）の入力パラメータファイルの例を以下に示します。計算例題は、 samples/elastic/Si/s0 以下にあります。

Control{
    condition = initial
}
accuracy{
    cutoff_wf =  20  rydberg
    cutoff_cd =  80  rydberg
    num_bands =   8
    ksampling{
        mesh{ nx = 10, ny = 10, nz = 10 }
    }
}
structure{
    unit_cell_type = primitive
    unit_cell{
        #units angstrom
        a_vector =   0.000000 2.723515 2.723515
        b_vector =   2.723515 0.000000 2.723515
        c_vector =   2.723515 2.723515 0.000000
    }
    symmetry{
        method = automatic
    }
    atom_list{
        coordinate_system = internal
        atoms{       #tag    rx      ry      rz   element
            0.125   0.125   0.125        Si
           -0.125  -0.125  -0.125        Si
        }
    }
    element_list{
        #tag   element   atomicnumber
                 Si            14
    }
}
structure_evolution{
    stress{
        sw_stress=on
    }
}

SCF 計算と同様に PHASE を実行します。

% mpirun ~//bin/phase

計算が終了したら結果を確認します。

% grep -A3 'Total STRESS TENSOR' output000

STRESS TENSOR
0.0000004032 0.0000000000 0.0000000000
0.0000000000 0.0000004032 0.0000000000
0.0000000000 0.0000000000 0.0000004032

ストレステンソルは

$$\begin{aligned} \begin{pmatrix} X_{x} & X_{y} & X_{z} \\\ Y_{x} & Y_{y} & Y_{z} \\\ Z_{x} & Z_{y} & Z_{z} \\\ \end{pmatrix} \end{aligned}$$

の形式で出力されています。出力されている値の単位は [Hartree/Bohr³] です。上の結果では入力データとしてわずかに格子定数を小さく取ってあるため、正の X_x, Y_y, Z_z が出力されています。

また、釣り合いの位置からの格子変形( ≡ e)、スティフネス定数( ≡ c)を用いると次のようなフックの法則が成り立ちます。

$$\begin{aligned} \left. \ \begin{matrix} X_{x} = c_{11}e_{\text{xx}} + c_{12}e_{\text{yy}} + c_{12}e_{\text{zz}} \\\ Y_{y} = c_{12}e_{\text{xx}} + c_{11}e_{\text{yy}} + c_{12}e_{\text{zz}} \\\ Z_{z} = c_{12}e_{\text{xx}} + c_{12}e_{\text{yy}} + c_{11}e_{\text{zz}} \\\ X_{y}( = Y_{x}) = c_{44}e_{\text{xy}} \\\ Y_{z}( = Z_{y}) = c_{44}e_{\text{yz}} \\\ Z_{x}( = X_{z}) = c_{44}e_{\text{zx}} \\\ \end{matrix} \right\} \end{aligned}$$

弾性定数

ストレステンソルの計算結果から、弾性定数の計算を行う例を紹介します。弾性定数は、歪みのない結晶と歪みのある結晶のストレステンソルを利用すると、上述のフックの法則より計算することができます。ここでは、Si（立方晶）を例に説明します。この例題の入力ファイルは、 samples/elastic/Si の s0 および sxx の下にあります。 s0 が歪みのない結晶の入力データ、 sxx がxx方向に歪みを与える場合の入力データです。

弾性定数を計算する場合、ストレステンソルの絶対値がなるべく小さくなる格子定数を利用することが望ましいです。この例題は、ストレステンソルの各成分の絶対値が小さくなる、以下のような格子定数を採用しています。

unit_cell{
  #units angstrom
  a_vector = 0.000000 2.731958 2.731958
  b_vector = 2.731958 0.000000 2.731958
  c_vector = 2.731958 2.731958 0.000000
}

ディレクトリー samples/elastic/Si/s0 においてストレステンソルを計算すると以下のような結果が出力されます。

% grep -A3 'Total STRESS TENSOR' output000

Total STRESS TENSOR
      0.0000000301        0.0000000000        0.0000000000
      0.0000000000        0.0000000301        0.0000000000
      0.0000000000        0.0000000000        0.0000000301

つぎに、ディレクトリsxxへ移ります。

% cd ../sxx

このディレクトリーに置かれている入力ファイルはs0のものとほぼ同じですが、以下の設定によって11方向(xx方向)に0.001という格子歪みを与えた計算を行うことになります。

structure{
  ...
  ...
  strain{
    sw_strained_cell = on
    e11 = 0.001
  }
}

structureブロックの下にstrainブロックを作成し、変数sw_strained_cellをonとすると格子歪みを与えて計算を行うことができます。さらに変数e11, e12, e13, e22, e23, e33で非ゼロとする歪み成分を指定します。このディレクトリで計算を行うと、以下のようなストレステンソルが得られます。

% grep -A3 'Total STRESS TENSOR' output000

Total STRESS TENSOR
     -0.0000051660        0.0000000000        0.0000000000
      0.0000000000       -0.0000018789        0.0000000000
      0.0000000000        0.0000000000       -0.0000018789

この例題の弾性定数の計算には、ストレステンソルの対角成分を用います。この計算例題では、回転や剪断ひずみを与えていませんので、非対角項は0になっています。

得られたストレステンソルから、スティフネス定数c₁₁, c₁₂を次式から計算します。

$$\begin{aligned} \left. \ \begin{matrix} X_{x} = c_{11}e_{\text{xx}} + c_{12}e_{\text{yy}} + c_{12}e_{\text{zz}} \\\ Y_{y} = c_{12}e_{\text{xx}} + c_{11}e_{\text{yy}} + c_{12}e_{\text{zz}} \\\ Z_{z} = c_{12}e_{\text{xx}} + c_{12}e_{\text{yy}} + c_{11}e_{\text{zz}} \\\ X_{y}( = Y_{x}) = c_{44}e_{\text{xy}} \\\ Y_{z}( = Z_{y}) = c_{44}e_{\text{yz}} \\\ Z_{x}( = X_{z}) = c_{44}e_{\text{zx}} \\\ \end{matrix} \right\} \end{aligned}$$

今のケースではe_xx以外は0なので、歪みのない系とある系のX_x, Y_yの差を歪み量 (0.001)で除すればc₁₁とc₁₂を計算することができます。結果は c₁₁ = 5.196 × 10^− 3 Hartree/Bohr³ = 153 GPa, c₁₂ = 1.909 × 10^− 3 Hartree/Bohr³ = 56.1 GPa となります。

一方、弾性定数(ヤング率( ≡ Y)・ポアソン比( ≡ P)・体積弾性率( ≡ B))はスティフネス定数を用いて次のような式で書き表されます。剛性率は Y/(2 + 2P) と書けます。

$$\begin{aligned} \left. \ \begin{matrix} Y = \frac{c_{11}^{2} + c_{11}c_{12} - 2c_{12}^{2}}{c_{11} + c_{12}} \\\ P = \left| \frac{c_{12}}{c_{11} + c_{12}} \right| \\\ B = \frac{c_{11} + 2c_{12}}{3} \\\ \end{matrix} \right\} \end{aligned}$$

これにスティフネス定数c₁₁, c₁₂を代入すれば Siの弾性定数は Y = 123 GPa, P = 0.268, B = 88.4 GPa と求まります。

より精度の高い弾性定数の計算を行ないたい場合、 cutoff_wf, cutoff_cdを大きめにとり、電子状態を充分に収束させる必要があり、計算時間のかかる計算になります。

ストレステンソルの補正

PHASE/0によるストレステンソルの計算は、精度が低い場合があります。原因は、“格子がひずむことによる平面波数の変化の効果”がとりいれられていないからです。この効果を取り入れることによって、ある程度補正を行うことが可能です。

方法１．

運動エネルギーの計算におけるGベクトルの高周波成分をスメアすることによって、“平面波数一定”の状況を“カットオフエネルギー一定”の状況に近づけることができます。 [Bernasconi95] では、運動エネルギーの高周波成分を以下のように置き換えることが提案されています。

$$G^2 \rightarrow G^2 + A \left[ 1 + \mathrm{erf}\left(\frac{\frac{1}{2}G^2-E_0}{\sigma}\right) \right]$$

PHASE/0では、上式を利用したストレステンソルの計算を行うことができます。以下のような設定を入力パラメーターファイルに記述します。

structure_evolution{
  lattice{ sw_optimize_lattice = on }
  stress{
    sw_smear_KE = on
    a = 15 rydberg
    sigma = 0.1 rydberg
    e0 = 35 rydberg
  }
}

structure_evolutionの下にstressブロックを作成し、設定を行います。sw_smear_KE=onとするとこの機能が有効になります。a, sigma, e0には対応するパラメーターを指定します。

デフォルト値はa=0.375, ecut, sigma = 0.1 Rydberg, e0=ecut-1 Rydbergです。

方法２．（バージョン2019.01以降）

複数のカットオフエネルギーによる計算から誤差を見積もることができます。ターゲットカットオフエネルギーを E_c 、変化量を ΔE_c 全エネルギーの変化量を ΔE_t とすると、ストレスの誤差 σ_e は以下のように見積もることができます。

$$\sigma_{\mathrm e} = - \left( \frac{2 E_{\mathrm c}}{3V} \right) \times \left( \frac{\Delta E_{\mathrm t}}{\Delta E_{\mathrm c}} \right)$$

この補正をPHASE/0に計算させるには、以下のようにstressブロックにおいてsw_stress_correctionをonとします。

structure_evolution{
  stress{
    sw_stress = on
    sw_stress_correction = on
  }
}

ストレステンソルの補正は、カットオフエネルギーを変化させてストレステンソルを求めることによって計算します。どの程度カットオフエネルギーを変化させるかはdelta_ecutによって指定します。カットオフエネルギーをecut-delta_ecutとしたケースとecut+delta_ecutとしたケースのストレステンソル計算が行われ、その後補正が計算されます。なお、補正の計算前に計算が終了した場合継続できないので注意が必要です。補正値は以下のようにoutput000ファイルに記録されます。

!** Pulay stress : -0.000194696412156

補正が計算されたあと、入力パラメーターファイルに格子最適化の設定が行われている場合補正を組み込んだ状態で格子最適化計算が始まります。そうでない場合、以下の要領で補正の値を入力ファイルに書き込み、格子最適化などを行う設定にしたうえで再度計算を実行してください（補正が必要なのは対角要素のみ、また誤差が－0.0001auだったとして）。

structure_evolution{
  lattice{
    sw_optimize_lattice = on
    external_stress{
      s11 = -0.0001
      s22 = -0.0001
      s33 = -0.0001
    }
  }
}

※ バージョン 2019.01未満の場合補正は手動で計算する必要があります

検証

これらの補正を利用し、TiO₂の格子定数を計算した結果を以下の表にまとめました。方法１．のパラメーターはデフォルト値、方法２．のは±5 Rydbergとしました。

\

a (bohr)

c (bohr)

カットオフ36 Rydberg, EV曲線 8.8017 カットオフ36 Rydberg, 補正なし 8.6825 カットオフ36 Rydberg, 方法 1. 8.7593 カットオフ36 Rydberg, 方法 2. 8.8052 カットオフ80 Rydberg, 補正なし 8.7918

5.6355

5.5862

5.6072 5.6200

5.6158

方法１．２．とも改善しています。特に、方法２．を使うとEV曲線からもとめた格子定数とほぼ同じ格子定数が得られています。

仕事関数

機能の概要

PHASEを利用して、仕事関数を評価することが可能です。ここでは、仕事関数を計算する方法を説明します。

第一原理計算の枠組み内における仕事関数とは、真空準位とフェルミエネルギーとの差です。真空準位は、表面のSCF計算を実施し、表面から十分に離れた箇所での局所ポテンシャルを利用して算出することができます。

入力パラメータ

仕事関数を計算するためには、表面のモデルを準備する必要があります。対象としたい系の、対象としたい面

方位をもつ表面モデルを用意します。さらに、入力データのpostprocessingブロックにworkfunc ブロックを作成し、設定を行います。

postprocessing{
  workfunc{
    sw_workfunc = on
    sw_add_xc_to_vloc = off
  }
}

各変数は以下の意味をもちます。

sw_workfunc

仕事関数の計算に必要なデータを出力するためのスイッチです。

出力させたい場合にonとします。

sw_add_xc_to_vloc

局所ポテンシャルを出力する際に、交換相関相互作用を含めるかどうかを指定します。交換相関相互作用は表面から十分離れた場所において0 になると考えられるので、局所ポテンシャルに含めなくても正しい仕事関数が得られることが期待できます。デフォルト値はonですが、offにしておくことによってより少ない真空層で収束した仕事関数を得ることが可能です。

このような設定を行ったら、通常通りPHASE を実行します。計算が収束した後に、必要な局所ポテンシャル

データなどが出力されます。すでに収束した計算に対する継続計算として実行することも可能です。

計算の実行方法

計算が終了した段階では、局所ポテンシャルのデータが逆空間のデータとして保存されます。仕事関数を得るた

めには、逆空間のデータを実空間へ逆フーリエ変換し、表面内で平均を計算しその結果を出力する必要がありま

す。このような処理を行うプログラムがworkfuncです。このプログラムのソースコードはsrc_workfuncディレクトリーにあります。コンパイルするためには、Fortran90コンパイラーが必要です。workfuncをコンパイルするには、たとえば以下のようなコマンドを実行します。

% cd src_workfunc
% export F90=ifort
% make

環境変数F90にFortran90コンパイラーを指定します。環境変数F90のデフォルト値はgfortranです。

以下のように利用します。

% workfunc -z ZAXIS

ZAXIS に、表面に垂直とみなす軸" を指定します。a軸の場合1, b軸の場合2, c軸の場合は3 を指定します。

指定しない場合のデフォルト値は3 です。

計算結果の出力

workfunc の処理が終了すると、nfvlcr.cube とnfvlcr_av.data の2 種類のファイルが生成されます。nfvlcr.cube ファイルは、実空間の局所ポテンシャルデータを持つGaussian Cube 形式のデータファイルです。nfvlcr_av.data には表面に垂直な距離と面内で平均した局所ポテンシャルのデータが記録されています。以下のようなデータ形式となっています。

# Fermi energy (eV) -0.37838
# distance along the z-axis(Angstrom) averaged local potential (eV)
0.104167 -0.218799E+01
0.208333 -0.250195E+01
0.312500 -0.331223E+01
0.416667 -0.427665E+01
0.520833 -0.495695E+01
0.625000 -0.496651E+01
0.729167 -0.425552E+01
.....
.....
.....

ファイルの1 行目にフェルミエネルギーがeV 単位で記録されています。3 行目以降が実際のデータです。1 列目にÅ単位で表面に垂直な距離が、2 列目に対応する局所ポテンシャルの面内平均eV 単位で記録されます。局所ポテンシャルは、表面からある程度離れた地点においてはほぼ一定値となります。この時の値とフェルミエネルギーとの差が仕事関数に相当します。

nfvlcr_av.dataファイルから局所ポテンシャルがフラットになる領域を推定し、フェルミエネルギーとの差を計算することによって仕事関数をもとめるPerl スクリプトがworkfunc.pl です。以下のように利用します。

% workfunc.pl nfvlcr_av.data OPTIONS

実行すると、計算された仕事関数の値が標準出力に出力されます。また、workfunc.eps という、局所ポテンシャルと表面に垂直な距離の関係をグラフ化したEPS ファイルも作成されます。

計算例：アルミニウムの仕事関数

アルミニウムの仕事関数の計算例を紹介します。サンプルデータは、 samples/surface/workfunc/Al です。

利用する系は、Al (111) 7 層の表面モデルです。表面に垂直な軸はc 軸とします。c 軸の長さは、50 Å としました。アルミニウム(111) 面はほとんど再構成しないので、構造最適化は施しませんでした。原点を中心に、反転対称性が存在するようにモデルを作成しました。また、交換相関相互作用は局所ポテンシャルに含めない設定で計算を行いました。アルミニウムの表面モデルを advanced_WF_Almodel に示します。

Al(111)面7層モデル

PHASE によるSCF 計算が終了したのちにworkfunc プログラムによってnfvlcr_av.data ファイルを作成し、さらにworkfunc.pl スクリプトを利用して得られた局所ポテンシャルと表面に垂直な距離の関係を advanced_WF_Alpot に示します。仕事関数は、4.05 eV と計算されました。この値は実測値である4.08 eVと近い結果となっています。

表面に垂直な距離と局所ポテンシャルの関係

原子周囲局所ポテンシャル出力機能（2021.02以降）

機能の概要

本機能を用いることによって原子周囲の局所ポテンシャルの平均値を評価することが可能です。

局所ポテンシャルを V(r), 原子A中心の球対称な重み関数をfとすると、原子Aの周囲の局所ポテンシャル平均値は、

V_A = ∫drV(r)f(|r−R_A|)

であらわされます。関数fとしては、積分半径 r_c 内でのみ有限の値を持つ関数

$$\begin{aligned} f\left(r\right) = \begin{cases} \frac{1}{\Omega} & r \leq r_c \\\ 0 & r>r_c \end{cases} \end{aligned}$$

が考えられます。ここで$\Omega=\frac{4\pi r_c^3}{3}$ です。さて、PHASE/0では、局所ポテンシャルは

V(r) = ∑_Ge^iG ⋅ r

のように表現されるので advanced_lpot_eq2 advanced_lpot_eq3 を advanced_lpot_eq1 に代入するとつぎのような表式が得られます。

$$\begin{aligned} V_A = \sum_G V\left(G\right) e^{i G \cdot r} w, \\\ w = \frac{4\pi}{G^3} \frac{1}{\Omega} \left(G r_c\right)^3 j_1 \left(G r_c \right), \\\ j_1 \left(x\right) = \frac{\sin{x}}{x^2} - \frac{\cos{x}}{x} . \end{aligned}$$

advanced_lpot_eq4 式を用いて原子周囲の局所ポテンシャルの平均値を求めることができます。

入力パラメーターファイル

局所ポテンシャル平均値を計算するには、以下のような入力を作成します。

postprocessing{
  potential_average{
     sw_calc_pot_avg = on          ! { on |off },  default : off
     sw_add_xc_pot = on            ! { on | off }, default : off
     cutoff_radius = 2.0 bohr        !  default : 2.0
  }
}

sw_calc_pot_avg = on を指定すると、原子近傍の局所ポテンシャルの平均値を計算します。sw_add_xc_pot = on の場合、局所ポテンシャルに交換・相関ポテンシャルを含めるようにします。原子を中心とした空間積分の積分半径は、cutoff_radius で指定します。

出力ファイル

局所ポテンシャル平均値は、potential_on_atoms.dataに出力されます。

# Potential on atoms
#   id     pot (eV)
    1     -19.343121
    2     -19.539243
 (後略)

第1列及び第2列は、それぞれ原子のインデックス、局所ポテンシャルの平均値 (単位：eV)に対応します。

計算例： Al-Graphene-Al系

Al-Graphene-Al系を用いた計算例を紹介します。入力ファイルは samples/local_pot_av/Al-Graphene-Al 以下にあります。 Al-Graphene-Al系の原子配置は、 local_pot_av_config に示す通りです。

Al-Graphene-Al系の原子配置。茶及び水色の球は、それぞれ C及びAl原子に対応する。C原子のインデックスは、左端から順につけるものとする。

その他の計算条件は以下のとおりです。

Al-Graphene-Al 系の計算に用いた条件

計算条件	値
平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k点サンプリング	monk (1x36x1, Γ -centered)
交換相関相互作用	GGAPBE, PAW
SCF収束条件 [Ha/atom]	1.0E-8

計算結果を以下に示します。青線及び赤線は、それぞれ積分半径1.0及び2.0 の結果に対応します。また、四角 (丸) は、交換相関ポテンシャルを含まない (含む) 計算に対応します。プロファイルは、積分半径や交換相関ポテンシャルの有無に、あまり依存しないことが分かります。

Al-Graphene-Al 系の局所ポテンシャル平均値のプロファイル。左端のC原子のポテンシャルがエネルギーの原点です。

帯電欠損状態の評価 (バージョン2021.01以降)

概要

PHASE/0 では、電荷非中性状態を計算する際、周期的境界条件によるエネルギー発散を抑えるために一様な背景電荷を考慮して処理しています。この処理により、全エネルギーは周期的境界条件による背景電荷間の相互作用を含んだ値となります。このため、例えば、帯電した原子欠損などを含む半導体の生成エネルギーは、実際の値とは異なるものとなります。これを補正する手法として、近年、FNV 法 [Freysoldt09] が注目されており、対応するプログラム( sxdefectalign , CoFFEE など) がいくつか公開されています。PHASE/0 2021年版以降、これらのプログラムに対するインターフェースが配備され、CUBE ファイル形式で出力することができるようになりました。また、これらのプログラムを利用せずに補正量を評価する仕組みも実装されています。さらに、異なる荷電状態間の生成エネルギーを比較した図を描画するためのPythonスクリプトが付属します。

理論

電荷qをもつ欠陥 (D) の生成エネルギーは

$$E_{\rm form}^{\rm D, q} \left( E_{\rm F} \right) = E_{\rm DFT}^{\rm D, q} - E_{\rm DFT}^{\rm bulk} - \sum_{i} n_{i} \mu_{i} + q \left( E_{\rm VBM} + E_{\rm F} \right) + E_{\rm corr}^{q}$$

で書けます。 $E_{\rm DFT}^{\rm D, q}$ および $E_{\rm DFT}^{\rm bulk}$ はDFT計算における全エネルギー値で、前者は帯電欠陥、後者は非帯電・無欠陥の完全結晶 (Pristine) に対応します。 n_i は完全結晶からの元素iの原子数の増減、 μ_i はその化学ポテンシャルです。 $E_{\rm VBM}$ は完全結晶の価電子帯上端のエネルギー、 $E_{\rm corr}^q$ は帯電量に依存する補正項です。

FNV法

FNV法では、 $E_{\rm corr}^{q}$ を2つの寄与に分けて考えます [Freysoldt09] 。

$$E_{\rm corr}^{q} = E_{\rm PC}^q - q \Delta V_{\rm far}$$

第1項は、孤立した点電荷qの静電エネルギーと、一様な背景電荷を加えて周期的境界条件を課した際の静電エネルギーの差です。

$$E_{\rm PC}^{q} = E_{\rm isolated}^q - E_{\rm periodic}$$

この値は、スーパーセルのサイズ・誘電率から、DFT計算とは無関係に決定されます。一方、第2項は、potential alignment 補正項で、静電ポテンシャルの差

$$\Delta V \left( \mathbf{r} \right) = V_{\rm DFT}^{\rm D, q} \left( \mathbf{r} \right) - V_{\rm DFT}^{\rm bulk} \left( \mathbf{r} \right) - V_{\rm PC}^q \left( \mathbf{r} \right)$$

を、欠陥から十分に離れた位置で評価した値です。ここで、 $V_{\rm PC}^q$ は、静電エネルギーに対応する静電ポテンシャルです。

外部解析プログラム

sxdefectalign、CoFFEEなどの解析プログラムは、以下の入力を要求します。

結果の位置情報
誘電率の値
$V_{\rm DFT}^{D,q} \left( \mathbf{r} \right)$ や $V_{\rm DFT}^{\rm bulk}$ の空間分布。たとえばCUBE形式。

出力されるのは以下の情報です。

$E_{\rm PC}^q$ の値
ΔV(r) の値。各格子ベクトル方向に対するプロット。他の格子ベクトル方向については面平均をする。欠陥からの距離に対して、原子位置近傍で平均した値を出力することもある。

既知の問題として、以下が挙げられます。

sxdefectalign では、欠陥からの距離に対するプロットが正常に出力されないことがある。
CoFFEE では、各格子ベクトル方向に対するプロットも、ユーザーが python スクリプトを、都度書く必要がある。また、欠陥からの距離に対するプロットは得られない。
CoFFEE では、 $E_{\rm PC}^q$ の値は、セルサイズを変えて複数計算し、fitting の結果として得る。このため、ユーザーの手間が多い。

extended FNV 法

Kumagaiら [Kumagai14] のextended FNV法では、 $E_{\rm PC}^q$ は

$$E_{\rm PC}^q = - \frac{q}{2} \left[ \sum_{\mathbf{R}_i \neq 0} \frac{q}{\sqrt{|\varepsilon|}} \frac{{\rm erfc} \left( \sqrt{\mathbf{R}_i \mathbf{\varepsilon}^{-1} \mathbf{R}_i} \right)}{\sqrt{\mathbf{R}_i \mathbf{\varepsilon}^{-1} \mathbf{R}_i}} - \frac{\pi q}{\Omega \gamma^2} + \sum_{\mathbf{G}_i \neq 0} \frac{4\pi q}{\Omega} \frac{\exp \left(-\mathbf{G}_i \varepsilon \mathbf{G}_i / 4\gamma^2\right)}{\mathbf{G}_i \varepsilon \mathbf{G}_i} - \frac{2\gamma q}{\sqrt{\pi |\varepsilon|}} \right]$$

で表されます。ここで、 γ は適当に選んだ収束パラメータです。また、 ε は誘電率テンソル、 Ω は系の体積です。一方、ポテンシャル $V_{\rm PC}^q$ は、

$$V_{\rm PC}^q = \sum_{\mathbf{R}_i \neq 0} \frac{q}{\sqrt{|\varepsilon|}} \frac{{\rm erfc} \left( \sqrt{\mathbf{R}_i \mathbf{\varepsilon}^{-1} \mathbf{R}_i} \right)}{\sqrt{\mathbf{R}_i \mathbf{\varepsilon}^{-1} \mathbf{R}_i}} - \frac{\pi q}{\Omega \gamma^2} + \sum_{\mathbf{G}_i \neq 0} \frac{4\pi q}{\Omega} \frac{\exp \left(-\mathbf{G}_i \varepsilon \mathbf{G}_i / 4\gamma^2\right)}{\mathbf{G}_i \varepsilon \mathbf{G}_i} \exp \left[ i \mathbf{G}_i \left( \mathbf{r}-\mathbf{R}_d\right) \right]$$

で表されます。ここで、 R_d は欠陥の位置です。

$V_{\rm far}$ の評価に関しては、以下の手順をとります。まず、各原子サイトにおける ΔV(r) を計算します。このうち、欠陥からの距離 R_WS より遠方にあるデータを平均し、 $\Delta V_{\rm far}$ とします。なお、 R_WS は、欠陥間の最短距離の半分とします。 eq:defectq5 , eq:defectq6 による計算を用いることによってsxdefectalign, CoFFEEなどを用いずとも(自動的に) $E_{\rm corr}^q$ の評価が可能です。

使い方

外部解析プログラムの使用

sxdefectalign や CoFFEE に受け渡す静電ポテンシャルを CUBE 形式で出力するには、入力パラメーターファイルを以下のように設定します。

postprocessing{
    electrostatic_potential{
       sw_write_electrostatic_pot = on
       unit = Rydberg                     ! { Rydberg | Hartree | eV }
    }
}

単位系は unit で指定し、Rydberg、Hartree、eV の選択が可能です。default 値は　eV です。

結果はelectrostatic_pot.cubeファイルに出力されます。ヘッダー部に単位系が表示され、以降、各軸方向の分割数、メッシュの刻み幅、各メッシュ点における値と続きます。

Calculated by phase
Local+Hartree potential in Rydberg
   215         0.0000         0.0000         0.0000
   180    0.180912    0.000000    0.000000
   180    0.000000    0.180912    0.000000
   180    0.000000    0.000000    0.180912
    33  33.000000     2.396007     2.396052     2.396080
(後略)

帯電欠陥、及び完全結晶の静電ポテンシャルを出力しておくようにします(sxdefectalignを利用する場合には、unitは Rydberg にしておく必要があります。)。以下に、これらのCUBEファイルを用いたsxdefectalign の実行例を示します。

sxdefectalign --charge 3 --eps 12.88
--vdef ../electrostatic_pot.cube --vref ../bulk/electrostatic_pot.cube
--center 0.0,0.0,0.0 --relative
--ecut 270 --qe > Log

系の帯電量を q としたとき、charge には -q を渡します。eps には誘電率の値を指定します。vdef及びvrefには、帯電欠陥及びバルクにおける静電ポテンシャルファイルを、それぞれ指定します。center には欠陥位置の内部座標を与えます。ecut には nfinp.data で指定した　cutoff_cd を Ry 単位で渡します。

sxdefectalignのLog ファイルの末尾には、以下のような出力がされます。

Defect correction (eV): 0.828347 (incl. screening & alignment)

ここで表示されているのは、値 ( 式 eq:defectq3 参照 )であり、potential alignment の寄与は含まれていません。 Logファイル以外には、 vline-eV-an.dat ( n=0, 1, 2 ) 及び vAtoms.dat が出力されます。前者は、式 eq:defectq4 の各結晶軸方向に対する、ポテンシャルの1次元プロファイルを出力したものです。一方、後者は、各原子の欠陥位置からの距離、及び原子上でのポテンシャルの値を出力したものです。ただし、 vAtoms.dat は正常に出力されないことが多々あるようです。

extended FNV 法

まず帯電なし・欠陥なしの bulk の計算を行います。入力パラメーターファイルは以下のように設定します。

postprocessing{
    electrostatic_potential{
       sw_write_electrostatic_pot = on
    }
}

結果はelecpot_bin.data ファイルが出力されます。静電ポテンシャルのG空間での値をバイナリ形式で出力しているため、外部ソフトと連携することはできません。ファイル名は、F_ELECPOT_BINにより変更できます。

&fnames
F_ELECPOT_BIN = './elecpot_bin.data'
/

ついで帯電欠陥の計算を行います。入力パラメーターファイルは以下のように設定します。

postprocessing{
    electrostatic_potential{
       sw_write_electrostatic_pot = on
    }
    charged_defect{
        correction{
            sw_calc_extfnv_correction = on          ! default : off
            dielectric_constant{
               exx = 12.88
               eyy = 12.88
               ezz = 12.88
            }
            position{
               x = 0.0, y = 0.0, z = 0.0
            }
        }
    }
}

charged_defect ブロック内の correction ブロックで、sw_calc_extfnv_correction = on を設定します。また、dielectric_constat ブロック内で誘電率テンソルの値、position ブロック内で欠陥の内部座標を指定します。欠損以外の場合には、内部座標ではなく、原子番号( atom_id )による指定も可能です (例：atom_id= 2 )。さらに、参照する bulk の静電ポテンシャルファイル名をF_ELECPOT_BIN_REF で指定します。

&fnames
F_ELECPOT_BIN_REF = '../bulk/elecpot_bin.data'
F_ELECPOT_BIN = './elecpot_bin.data'
/

結果はdefect_pot_correction.direction_n ( n=1,2,3, )、defect_pot_correction.atomsに出力されます。前者は、式 eq:defectq4 の各結晶軸方向に対する、ポテンシャルの1次元プロファイルを出力したものです。一方、後者は、各原子の欠陥位置からの距離、及び原子上でのポテンシャルの値を出力したものです。

以下に、defect_pot_correction.direction_1 の出力例を示します。

# dist. (Ang), pot_diff,    Vpc,  pot_diff -Vpc (eV)
   0.00000  -0.13205  -0.03907  -0.09299
   0.09573  -0.13743  -0.03902  -0.09841
(後略)

第2列及び第3列は、それぞれ式 eq:defectq4 の $V_{\rm DFT}^{D,q} \left( \mathbf{r} \right) - V_{\rm DFT}^{\rm bulk} \left( \mathbf{r} \right)$ 及び $V_{\rm PC}^q \left( \mathbf{r} \right)$ に対応します。第4列は ΔV(r) に対応します。なお、符号については、電子の電荷が負であることを考慮しています。

以下に、defect_pot_correction.atoms の出力例を示します。

# no.,   dist. (Ang),     pot_diff,            Vpc,              pot_diff -Vpc (eV)
    1    2.1961       -0.0550653932       -0.2936706923        0.2386052991
    2    3.9534       -0.1647261198       -0.1693351641        0.0046090443
(中略)

# Correction energy (eV):        0.7827480523

なお、最終行に表示されている “Correction energy” は、生成エネルギー図作成に必要な $E_{\rm corr}^q$ です。

生成エネルギー図の作成1

欠陥ごとの生成エネルギー図の作成方法を説明します。生成エネルギー図作成にはcalc_defect_formation_energy.pyスクリプトを用います。

平衡状態のとき、化学ポテンシャルは

$$\begin{aligned} \mu_{\rm Ga} + \mu_{\rm As} = \mu_{\rm GaAs}^0, \\\ \mu_{\rm Ga} \leq \mu_{\rm Ga}^0 , \\\ \mu_{\rm As} \leq \mu_{\rm As}^0 \end{aligned}$$

を満たします。ここで、 $\mu_{\rm Ga}^0$ 及び $\mu_{\rm As}^0$ は、それぞれGa及びAs単体の1原子あたりのエネルギーです。また、 $\mu_{\rm GaAs}^0$ は、GaAs 結晶の2原子あたりのエネルギーです。特に、Ga-rich 極限の場合は、

$$\mu_{\rm Ga} = \mu_{\rm Ga}^0, \mu_{\rm As} = \mu_{\rm GaAs}^0 - \mu_{\rm Ga}^0$$

となり、As-rich 極限では

$$\mu_{\rm Ga} = \mu_{\rm GaAs}^0 - \mu_{\rm As}^0, \mu_{\rm As} = \mu_{\rm As}^0$$

となります。これら化学ポテンシャルの情報は、以下で説明するcalc_defect_formation_energy.pyスクリプトの入力に記述する必要があるものです。

例えば、GaAs のスーパーセルの中にGa欠損を導入した計算を行い、帯電量 q =-3, -2, -1, 0 の計算結果が得られているとします。生成エネルギー図作成のために、以下のようなファイルを作成します。ファイル名に制限はありません。以下の例ではtmp1.in とします。

&VBM        #(eV)
6.00887

&band_gap        #(eV)
1.424

&Chemical_potentials #(Ha)
Ga  -138.9838873703      # mu_GaAs_bulk -mu_As
As   -87.9848825305

&Defects        #elements, number ( negative==vacancy, positive==impurity )
Ga  -1

&Host_supercell_energy #(Ha)
-7263.0006368270

&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0 -7123.9292947570

&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0

以下、用語の説明をします。

キーワードの説明

キーワード	単位	意味
&VBM	eV	$E_{\rm VBM}$
&band_gap	eV	バンドギャップ値
&Chemical_potentials	Hartree	化学ポテンシャル μ_i 。系を構成する全元素について記述する。
&Defects		導入した欠陥
&Host_supercell_energy	Hartree	$E_{\rm DFT}^{\rm bulk}$
&Defective_supercell_energy	Hartree	$E_{\rm DFT}^{D, q}.$ 計算した全荷電状態についてqが小さい順に記述する。
&Correction	eV	$E_{\rm corr}^q.$ 計算した全荷電状態についてqが小さい順に記述する。

以下の要領でcalc_defect_formation_energy.pyスクリプトを実行します。

python3  calc_defect_formation_energy.py  input
[-o OUTFILE]
[--emin EMIN] [--emax EMAX] [--de DE]
[--vmin VMIN] [--vmax VMAX]
[--image_format IMAGE_FORMAT]

括弧内は省略可能なオプションで、その意味は以下のとおりです。

キーワードの説明

引数	意味	デフォルト値
-o	出力するファイルの名称	result
--emin	エネルギーEF の最小値	-1.0
--emax	エネルギーEF の最大値	6.0
--de	エネルギーEFの刻み幅	0.01
--vmin	生成エネルギーの表示範囲の最小値	-5.0
--ymax	生成エネルギーの表示範囲の最大値	5.0
--image_format	可視化画像の形式 (png/eps)	png

実行例を以下に示します。

python3 calc_defect_formation_energy.py tmp1.in -o result1

上述のコマンドを実行した結果、生成されるファイルは result1.qdep、result1.min、result1.gnu、result1.png です。以下に、result1.qdep及びresult1.min の一部を示します。前者には、各帯電状態における生成エネルギーの、フェルミエネルギー依存性が出力されています。後者は、各フェルミエネルギー値における最小の生成エネルギー値が出力され、ファイル末尾に Charge Transition Level が追記されています。

result1.qdepファイルの内容

# Formation energy
# Ef (eV)    q=-3       q=-2       q=-1       q=0
  -1.00000   6.48430   4.90715   3.53961   2.37978
(後略)

result1.minファイルの内容

# Formation energy
# Ef (eV)    min
  -1.00000   2.37978
(中略)
   6.00000 -14.51570

# Charge   Transtion level [eV]
#-2/-3             0.57715
#-1/-2             0.36754
# 0/-1             0.15983

result1.gnu は可視化のための gnuplot 用ファイル、result1.png は gnuplot の出力です。

生成エネルギー図の作成2

興味ある全ての欠陥構造について、 defectq_forme_fig1 の作業が終了しているものとします。ここでは、これらをまとめた図の作成を行うため、以下のようなファイルを作成します。ファイル名に制限はありません。以下の例ではgather1.in とします。

#
# title   filename ( excluding ".min" )
#
&List
Vac_Ga    Vacancy_Ga/result1
Vac_As    Vacancy_As/result1
Ga_As     Ga_for_As/resultaa
As_Ga     As_for_Ga/resultaa

&band_gap        #(eV)
1.424

以下、用語の説明をします。

キーワードの説明

ワード	単位	意味
&List		欠陥の名称、及び計算の出力( `defectq_forme_fig1` 指定した outfile 名)
&band_gap	eV	バンドギャップ値

以下の要領でplot_multiple_defect_formation_energy.pyスクリプトを実行します。

python3  plot_multiple_defect_formation_energy.py  input
[-o OUTFILE]
[--emin EMIN] [--emax EMAX]
[--vmin VMIN] [--vmax VMAX]
[--image_format IMAGE_FORMAT]
[--keypos_h KEYPOS_H] [--keypos_v KEYPOS_V]

括弧内は省略可能なオプションで、その意味は以下のとおりです。なお、EMIN、EMAX値は、 defectq_forme_fig1 の指定と揃えた方がよいです。

plot_multiple_defect_formation_energy.pyで省略可能なオプション一覧

引数	意味	デフォルト値
-o	出力するファイルの名称	result_all
--emin	エネルギーEF の最小値	なし
--emax	エネルギーEF の最大値	なし
--vmin	生成エネルギーの表示範囲の最小値	なし
--ymax	生成エネルギーの表示範囲の最大値	なし
--image_format	可視化画像の形式 (png/eps)	png
--keypos_h	凡例の水平位置 (left/center/right)	right
--keypos_v	凡例の垂直位置 (top/center/bottom)	top

以下に実行例を示します。

python3 plot_multiple_defect_formation_energy.py gather1.in

実行するとresult_all.gnu 及び result_all.png が生成されます。前者は gnuplot 用ファイルで、後者はこれを可視化したものです。

例題

概要

GaAs 64 原子をホストとして、帯電欠陥の計算を行った例を紹介します。計算条件は以下の通りです。なお、GaAs の格子定数は、基本格子で最適化を行いました。補正エネルギーは、PHASE/0に実装されている extended FNV 法により評価しました。

帯電欠陥の評価で使用した計算条件

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Monkhorst-Pack (2×2×2)
交換相関相互作用	GGAPBE, PAW
単位胞の1辺 [Å]	11.48882
SCF 収束条件	[Ha/atom] 1.0E-8
力の収束条件	[Ha/bohr] 2.0E-4
擬ポテンシャル	Ga_ggapbe_paw_02.pp, As_ggapbe_paw_02.pp
(補正項の計算で使用する)	誘電率 ※12.88

※ https://www.microwaves101.com/encyclopedias/gallium-arsenide

例題のディレクトリー構成および計算のながれ

本例題は複数の計算を行い、その結果をとりまとめ、スクリプトで処理することによって結果が得られる仕組みになっています。ここではディレクトリーの構成と計算のながれについて説明します。

必要な入力ファイルはサンプルディレクトリーの下の samples/defectq 以下のサブディレクトリーに配置されています。以下のようなディレクトリー構成になっています。

defectq 以下には Preparation ディレクトリーと GaAs_64_lattice_opt ディレクトリーが存在します。

Preparationディレクトリー

単体のAs, Ga, GaAsの格子最適化の入力ファイルが納められたディレクトリーです。bulk_As, bulk_Ga, bulk_GaAsディレクトリーがあり、格子最適化の入力ファイルが納められています。結果得られるエネルギーと eq:defectq8 eq:defectq9 から化学ポテンシャルをもとめ、スクリプト入力ファイルの&Chemical_potentialsに記述します。

GaAs_64_lattice_optディレクトリー

各種生成エネルギー計算の入力ファイルが格納されています。ベースとなる結晶は格子最適化によってもとまった格子定数から作成した64原子系です。GaAs_64_lattice_optディレクトリーは、さらに以下のようなサブディレクトリー群が存在します。
GaAs_64_lattice_opt以下のディレクトリー構成

ディレクトリー名説明

Pristine 欠損のない結晶の入力ファイルが納められたディレクトリー。

Ga_for_As AsをGaで置換した欠陥構造の入力ファイルが納められたディレクトリー。

As_for_Ga GaをAsで置換した欠陥構造の入力ファイルが納められたディレクトリー。

As_interstitial_As4 Asが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。

As_interstitial_Ga4 Gaが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。

Ga_interstitial_As4 Asが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。

Ga_interstitial_Ga4 Gaが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。

それぞれのディレクトリーにはさらにq_ q ディレクトリーが存在します。ここで q は電荷をあらわす数値です。

Pristineの場合はq_0とq_0.2ディレクトリーが存在します。 $E_{\rm VBM} = \frac{E\left( N \right) - E\left( N-\Delta N \right)}{\Delta N}$ (今の場合 ΔN = 0.2 )という関係から $E_{\rm VBM}$ をもとめ、スクリプト入力の&VBMに記述します。またq_0のエネルギーをスクリプト入力の&Host_supercell_energyに記述します。

ほかの欠損に対応するディレクトリーでは、エネルギーの計算結果をスクリプト入力の&Defective_supercell_energyに記述します。また、q_0以外の計算ではextended FNV法による補正エネルギー $E_{\rm corr}^q$ の計算がなされます。結果はdefect_pot_correction.atomsファイルの末尾に記録されるので、その計算値をスクリプト入力の&Correction に記述します (q=0の項には0を記述します)。

各欠損ディレクトリーにおいて calc_defect_formation_energy.py スクリプトの入力を作成し、実行することによって各欠損の生成エネルギー図を作成することができます。またすべての欠陥のデータを集約し defectq_forme_fig2 の手続きをふむことによってすべての欠陥の結果をまとめた生成エネルギー図を作成することができます。

補正エネルギー比較

以下に、Ga欠損 (q=-3) における補正エネルギーについて、sxdefectalign とPHASE/0実装のextentend FNV 法による評価の比較を示します。両者がおおよそ一致していることが分かります。

Ga欠損 (q=-3) における補正エネルギーの内訳 (単位：eV)

	sxdefectalign	PHASE/0
Epc	1.24251	1.24272
dV	-0.00731	-0.02779
Ecorr	1.22058	1.15935

生成エネルギー図の作成 1.

以下では、Ga欠損に対して、2つの極限における入力を作成し可視化します。以降、Ga-rich 極限及びAs-rich極限の入力を示します。

calc_defect_formation_energy.py 用の入力 ( Ga-rich 極限におけるGa欠損; samples/defectq/GaAs_64_lattice_opt/Vacancy_Ga/cond_Ga_rich.in )

&VBM        #(eV)
6.00887

&band_gap        #(eV)
1.424

&Chemical_potentials #(Ha)
Ga -138.9586691142        #mu_Ga
As  -88.0101007866        #mu_GaAs -mu_Ga

&Defects        #elements, number ( negative==vacancy, positive==impurity )
Ga  -1

&Host_supercell_energy #(Ha)
-7263.0006368270

&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0  -7123.9292947570

&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0

calc_defect_formation_energy.py 用の入力 ( As-rich 極限におけるGa欠損; samples/defectq/GaAs_64_lattice_opt/Vacancy_Ga/cond_As_rich.in )

&VBM        #(eV)
6.00887

&band_gap        #(eV)
1.424

&Chemical_potentials #(Ha)
Ga  -138.9838873703      # mu_GaAs_bulk -mu_As
As   -87.9848825305       #mu_As

&Defects        #elements, number ( negative==vacancy, positive==impurity )
Ga  -1

&Host_supercell_energy #(Ha)
-7263.0006368270

&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0  -7123.9292947570

&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0

つぎのコマンドを実行すると、生成エネルギーの図 (result-Ga-rich.png及びresult-As-rich.png ) が出来ます。結果を defectq_fig1, defectq_fig2 に示します。

calc_defect_formation_energy.pyの実行 (Ga欠損)

python3  calc_defect_formation_energy.py
cond_Ga_rich.in  -o result_Ga_rich  --emin -0.5  --emax 2.0

python3  calc_defect_formation_energy.py
cond_As_rich.in  -o result_As_rich  --emin -0.5  --emax 2.0

Ga-rich limit (ファイル名：result-Ga-rich.png )

As-rich limit (ファイル名：result-As-rich.png )

生成エネルギー図の作成 2.

以下で、種々の欠陥構造に対する生成エネルギーをまとめた図を作成します。Ga-rich 及びAs-rich 極限のファイルが、それぞれ、 samples/defectq/GaAs_64_lattice_opt の下の gather_Ga_rich.in 及び gather_As_rich.in にあります。後者は、つぎに紹介する入力のresult-Ga-rich を result-As-rich に置換したものです。なお、欠陥の種類としては、文献 [Broberg18] 記載のものを選択しました。

gather_Ga_rich.inの内容

#
# title   filename ( excluding ".min" )
#

&List
Vac_Ga    Vacancy_Ga/result-Ga-rich
Vac_As    Vacancy_As/result-Ga-rich
Ga_As     Ga_for_As/result-Ga-rich
As_Ga     As_for_Ga/result-Ga-rich
Ga_i_As4  Ga_interstitial_As4/result-Ga-rich
As_i_As4  As_interstitial_As4/result-Ga-rich
Ga_i_Ga4  Ga_interstitial_Ga4/result-Ga-rich
As_i_Ga4  As_interstitial_Ga4/result-Ga-rich

&band_gap        #(eV)
1.424

つぎのコマンドを実行すると、生成エネルギーの図 (results-Ga-rich.png及びresults-As-rich.png) が出来ます。結果を defectq_fig3 および defectq_fig4 に示します。

python3  plot_multiple_defect_formation_energy.py
gather_Ga_rich.in  -o results_Ga_rich  --emin -0.2  --emax 2.0 --vmin -4.0 --vmax 5.0
--keypos_h left --keypos_v bottom
python3  plot_multiple_defect_formation_energy.py
gather_As_rich.in  -o results_As_rich  --emin -0.2  --emax 2.0 --vmin -4.0 --vmax 5.0
--keypos_h left --keypos_v bottom

Ga-rich limit (ファイル名：results-Ga-rich.png)

As-rich limit (ファイル名：results-As-rich.png)

文献 [Broberg18] と同様の生成エネルギー図が得られました (生成エネルギーの絶対値や Charge Transition Level は、完全には一致しません)。

XPS

概要

2020年版までのPHASE/0 では、XPS・XANES のピーク位置が実験値と異なることがありました。原因として、内殻電子の相対論的効果、スピン分裂効果が考えられます。前者については、CIAO を改変し相対論的効果を含む (DIRAC型の) 運動エネルギーを出力可能なようにしました。後者については、PHASE/0 のopencore 法を改変し、その効果を取り込めるようにしました。これによって、SiC、AlN、Li7Ti5O12 結晶の内殻軌道の結合エネルギーについて、従来法に比べて実験値との差が小さくなることなどを確認しています。

理論

相対論的運動エネルギー

Kohn-Sham 方程式は、

(T+V)ψ_i = ε_iψ_i

で表されます。 T 及び V はそれぞれ、運動エネルギー及びポテンシャル演算子、また ε_i 及び ψ_i は固有値・固有波動関数です。ここで、運動エネルギー演算子は、相対論的効果の取り入れ方により、いくつかの形状が提案されています [Lenthe96] 。

$$\begin{aligned} T^{\rm NR} = \frac{1}{2} \mathbf{p}^2, \\\ T^{\rm ZORA} = \mathbf{\sigma} \cdot \mathbf{p} \frac{c^2}{2c^2 - V} \bf{\sigma} \cdot \mathbf{p} , \\\ T^{\rm DIRAC} = \mathbf{\sigma} \cdot \mathbf{p} \frac{c^2}{2c^2+\varepsilon_i - V} \mathbf{\sigma} \cdot \mathbf{p}. \end{aligned}$$

$T^{\rm NR}$ は非相対論、 $T^{\rm ZORA}$ はZORA (Zero Order Regular Approximation)、 $T^{\rm DIRAC}$ はDirac相対論によるものです。ZORAの場合は、scalar relativistic 項とスピン軌道相互作用項に分解できます。

$$\begin{aligned} T^{\rm ZORA} = T_{\rm sc}^{\rm ZORA} + T_{\rm SOC}^{\rm ZORA}, \\\ T_{\rm sc}^{ZORA} = \mathbf{p} \frac{c^2}{2c^2 -V} \mathbf{p}, T_{\rm SOC}^{\rm ZORA} = \frac{c^2}{\left( 2c^2 -V \right)^2 } \sigma \cdot \left( \nabla V \times \mathbf{p} \right) \end{aligned}$$

Dirac の場合、上式にて V を V − ε_i と読み替えればよく、

$$\begin{aligned} T^{\rm DIRAC} = T_{\rm sc}^{\rm DIRAC} + T_{\rm SOC}^{\rm DIRAC}, \\\ T_{\rm sc}^{DIRAC} = \mathbf{p} \frac{c^2}{2c^2 -V - \varepsilon_i} \mathbf{p}, T_{\rm SOC}^{\rm DIRAC} = \frac{c^2}{\left( 2c^2 -V - \varepsilon_i \right)^2 } \sigma \cdot \left( \nabla V \times \mathbf{p} \right) \end{aligned}$$

となります。ただし、擬ポテンシャル作成時はscalar-relativisitic で解く必要があり、スピン軌道相互作用による効果は考慮されません。このため、CIAO が擬ポテンシャルファイルに出力する、内殻電子の運動エネルギーは scalar-relativisitic 項のみで評価します。

さて、内殻電子の運動エネルギーは

$$E_{\rm kin}^{\rm core} = \sum_{i}^{\rm core} \left< \psi_i \left| V_{\rm sc}^{\rm DIRAC} \right| \psi_i \right>$$

で表されます。上記は DIRAC の場合ですが、ZORA の場合も同様です。最後に、CIAOではポテンシャルは球対称として扱うので、波動関数は動径方向成分のみをもちます。このため、運動エネルギー演算子は以下のように展開できます。

$$\begin{aligned} T_{\rm sc}^{\rm ZORA} \psi_i \left( r \right) = - \frac{\partial V}{\partial r} \frac{c^2}{\left( 2c^2 - V\right)^2} \frac{\partial \psi_i}{\partial r} - \frac{c^2}{2c^2 - V} \Delta \psi_i, \\\ T_{\rm sc}^{\rm DIRAC} \psi_i \left( r \right) = - \frac{\partial V}{\partial r} \frac{c^2}{\left( 2c^2 + \varepsilon_i - V\right)^2} \frac{\partial \psi_i}{\partial r} - \frac{c^2}{2c^2 +\varepsilon_i - V} \Delta \psi_i, \\\ \Delta \psi_i \left( r \right) = \frac{1}{r} \frac{\partial^2}{\partial r^2} \left( r\psi_i \right) - \frac{l \left( l+1 \right)}{r^2} \psi_i \end{aligned}$$

内殻電子のスピン

CIAOでは、擬ポテンシャルはスピン中性で作成します。内殻電子についても同様であり、この時の内殻電子密度を $\rho_{\rm core}$ とします。ここで、ある軌道 (n,l) におけるスピン占有数 f_n, l^↑ 及び f_n, l^↓ が変わり、 f_n, l^↑ ≠ f_n, l^↓ になったとします。ただし、 f_n, l^↑ + f_n, l^↓ は擬ポテンシャル作成時と同じであるとします。このとき。各スピンの内殻電子密度は以下のように表されます。

$$\begin{aligned} \rho_{\rm core}^{\uparrow} \left( r \right) = \rho_{\rm core} \left( r \right) + \frac{1}{2} \Delta \rho \left( r \right) \\\ \rho_{\rm core}^{\downarrow} \left( r \right) = \rho_{\rm core} \left( r \right) - \frac{1}{2} \Delta \rho \left( r \right) \end{aligned}$$

ここで、

$$\Delta \rho \left( r \right) = \left( f_{n,l}^{\uparrow} - f_{n,l}^{\downarrow} \right) \left| \phi_{n,l}^{\rm AE} \left( r \right) \right|^2$$

はスピン密度の差で、波動関数 $\phi_{n,l}^{\rm AE}$ はスピン状態により不変と仮定しています。このスピン密度変化により影響を受けるエネルギーは、交換相関相互作用 (XC) エネルギー

$$\widetilde{E}_{\rm XC} \left( \rho \right) - E_{\rm XC}^{\rm PS} \left( \rho \right) + E_{\rm XC}^{\rm AE} \left( \rho \right)$$

です。第1項は非PAW項、第2、3項はPAW由来の項です。なお、スピン密度の和は不変であるので、Hartreeポテンシャルは変更を受けません。よって、XC ポテンシャルのみ Δρ(r) を考慮します。

XPSのピークエネルギーの計算方法

非金属の固体のXPSにおけるピークエネルギーは、以下の結合エネルギー $E_{\rm B}$ に対応します [Ozaki17] 。

$$E_{\rm B} = E_{\rm final} \left( N-1 \right) - E_{\rm initial} \left( N \right) + \mu_0$$

ここで、 N は系の電子数です。 μ₀ は系の化学ポテンシャルで、gapのある系では、価電子帯上端 (VBM) から伝導体下端 (CBM) まで取りえます。特に、真性半導体では、バンドギャップを E_g として、

$$\mu_0 = E_{\rm VBM} + \frac{1}{2} E_g$$

となります。なお、VBMのエネルギー $E_{\rm VBM}$ は

$$E_{\rm VBM} = \frac{E_{\rm initial} \left( N \right) - E_{\rm initial} \left( N - \Delta N \right)}{\Delta N}$$

で計算できます。また、帯電状態を周期的境界条件で計算するので、その補正法の1つとして、式 eq:defectq2 の ( q =1 ) を加えることもできます。

$$E_{\rm B} = E_{\rm final} \left( N - 1 \right) - E_{\rm initial} \left( N \right) + \mu_0 + E_{\rm corr}^1$$

一方、系が金属の場合には、

$$E_{\rm B} = E_{\rm final} \left( N \right) - E_{\rm initial} \left( N \right)$$

となります。

なお、固体ではなく分子の場合は、

$$E_{\rm B} = E_{\rm final} \left( N-1 \right) - E_{\rm initial} \left( N \right) + E_{\rm corr}^1$$

となります。

使い方

相対論的運動エネルギー

相対論的運動エネルギーを出力するためのCIAOの入力は下記の通り。

# PAW
sw_paw     1

# CORE ELECTRON INFO
sw_with_dipole_cor2val  1
method_ekin_core  1              ! default:0 (非相対論), 1:DIRAC, 2: ZORA

method_ekin_core 1 或いは 2 で、内殻電子の相対論的運動エネルギーを出力します。なお、”sw_with_dipole_cor2val 1” は、これまでと同様、内殻電子のXPS・XANES計算には必須です。

“method_ekin_core 1”の場合の3つのスイッチの設定は、以下の入力でも可能です。

sw_write_core_full  1

結果はgncpp2ファイル(擬ポテンシャルファイル)に出力されます。

CORE ENERGY CONTRIB
  0.319190836204E+02     kin              ! 内殻電子の運動エネルギー
 -0.677766126155E+02     ion
  0.697542582590E+01     hartr

内殻電子のスピン

ある原子の内殻電子が励起された状態を考える場合、該当する軌道に正孔を入れた擬ポテンシャルを使う必要があります。例えば、1s軌道励起の場合は 1s 軌道に正孔、2p軌道励起の場合は2p軌道に正孔を入れます。これらの原子は、内殻電子は開殻となるので、スピン自由度をもちます。内殻電子のスピンを扱うには、sw_opencore = on とします。デフォルトでは、式 eq:newxps9 の全項でスピンを考慮しますが、sw_xc_opencore_ae_only = onとすると、第3項 (AE部分)のみ考慮します。価電子が磁気モーメントをもつ磁性材料では、spin_orientation の指定により、内殻電子と (原子近傍の) 価電子の磁気モーメントの向きを、平行・反平行のいずれかに指定できます。ただし、SCF iteration で価電子の磁気モーメントが反転したりすると、内殻電子も追随するため、収束性が悪くなる場合があります。そこで、sw_fix_core_spin_pol = onとすると、初期磁気モーメントの向きに固定することが出来ます。

内殻電子のスピンを考慮するための入力例。

accuracy{
    paw = on
    core_electrons{
       sw_opencore = on                        ! default :off
       sw_xc_opencore_ae_only = on             ! default :off
       spin_orientation = anti_parallel        ! anti_parallel or parallel (default)
       sw_fix_core_spin_pol = on               ! default: off
    }
}
postprocessing{
   corelevels{
     sw_calc_core_energy = on           ! XPS 計算に必要
   }
}

また、s軌道以外から励起する場合には、内殻電子のスピン軌道相互作用により、XPSのピークが2つに分裂します。例えば、2p軌道の場合には、2p1/2と2p3/2 に分裂します。この場合、以下のようにすると、両者の全エネルギー値を計算することが出来ます。

内殻電子準位のスピン軌道分裂を考慮するための入力例。

postprocessing{
   corelevels{
     sw_calc_core_energy = on
     corehole{
         atom_id = 1　         ! 励起する原子
         orbital = 2p　　      ! 励起する軌道
     }
   }
}

(内殻電子のエネルギーを含む) 系の全エネルギーは、core_energy.data に出力されます。sw_calc_core_energy = onとした場合には、core_energy.data 末尾に、スピン軌道分裂したエネルギー値が追記されます。 eq:newxps10 eq:newxps13 eq:newxps15 などの結合ネルギーの計算における $E_{\rm final}$ や $E_{\rm initial}$ はここに記録された値を利用します。

core_energy.dataの出力例

# Etotal (Core+Valence)
  -18558.3245482622

# Etotal (Core+Valence+Soc_corehole)
# J = 3/2:             -18558.3331659845
# J = 1/2:             -18558.3073128174

例題

例題の計算結果を紹介します。XPSの例題は samples/XPS ディレクトリーに配置されています。用いる擬ポテンシャルは newxps_usage1 において説明したsw_write_core_full の値を1として作成した擬ポテンシャルで、それぞれの例題のpseudoサブディレクトリーに配置されています。

4H-SiC結晶 (C 1s)

4H-SiC結晶の例は samples/XPS/4H-SiC 以下に配置されています。そのディレクトリー構成は下記の通りです。

4H-SiC (C 1s) XPS計算のディレクトリー構造

第一階層	第二階層	説明
Final		終状態の入力ファイル
	q_0_no_opencore	中性・内殻スピンなし
	q_0_with_opencore	中性・内殻スピンあり
	q_1_no_opencore	電荷+1・内殻スピンなし
	q_1_with_opencore	電荷+1・内殻スピンあり
Initial		初期状態の入力ファイル
	q_0	中性
	q_0.2	電荷+0.2
	q_0.5	電荷+0.5
	q_1	電荷+1

計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.108 Å、c= 10.170 Å )を 3 × 3 × 3 倍しました。

4H-SiC (C 1s) XPS計算に用いた条件

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Monkhorst-Pack 2×2×2
交換相関相互作用	GGAPBE, PAW
SCF 収束条件 [Ha/atom]	1.0E-8
内殻電子のスピン	sw_xc_opencore_ae_only = on
CIAOの入力	sw_write_core_full 1
(補正項計算で使用する)誘電率	※9.76 (ab面内), 10.32 (c軸)

※https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118313534.app3

計算結果は以下の通りです。

4H-SiC (C 1s) のXPS 計算結果 ( 単位：eV )。式 eq:newxps13 により評価。

	内殻電子のスピン考慮なし	内殻電子のスピン考慮あり
$E_{\rm final} (N-1) - E_{\rm initial} (N)$	279.258	272.309
$E_{\rm corr}^1$	0.237	0.237
$E_{\rm VBM}$ ( ΔN = 0.2)	9.307
E_g	2.435
$E_{\rm B}$	289.975	283.026
実験値	283 +/- 0.8	^a)

a) http://www.xpsfitting.com/2012/01/silicon.html

内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。なお、式 eq:newxps14 を用いた場合、 $E_{\rm B}$ は以下のようになりました。

4H-SiC (C 1s) のXPS 計算結果 ( 単位：eV )。式 eq:newxps14 により評価。

	内殻電子のスピン考慮なし	内殻電子のスピン考慮あり
$E_{\rm B}$	290.815	283.864

w-AlN 結晶 (N 1s)

w-AlN結晶 (N 1s) の例は samples/XPS/w-AlN 以下に配置されています。そのディレクトリー構成は下記の通りです。

w-AlN 結晶 (N 1s) XPS計算のディレクトリー構造

第一階層	第二階層	説明
Final_N1s		終状態の入力ファイル
	q_0_no_opencore	中性・内殻スピンなし
	q_0_with_opencore	中性・内殻スピンあり
	q_1_no_opencore	電荷+1・内殻スピンなし
	q_1_with_opencore	電荷+1・内殻スピンあり
Initial		初期状態の入力ファイル
	q_0	中性
	q_0.2	電荷+0.2
	q_0.5	電荷+0.5
	q_1	電荷+1

計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.140 Å、c= 5.040 Å ) 3 × 3 × 2 倍しました。

w-AlN結晶 (N 1s)

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Monkhorst-Pack 2×2×2
交換相関相互作用	GGAPBE, PAW
SCF 収束条件 [Ha/atom]	1.0E-8
内殻電子のスピン	sw_xc_opencore_ae_only = on
CIAOの入力	sw_write_core_full 1
(補正項計算で使用する)誘電率	※8.23 (ab面内), 9.74 (c軸)

※https://materialsproject.org/materials/mp-661/

w-AlN (N 1s) のXPS 計算結果 ( 単位：eV )。式 eq:newxps13 により評価。

内殻	電子のスピン考慮なし\| 内殻	電子のスピン考慮あり
$E_{\rm final} (N-1) - E_{\rm initial} (N)$	395.772	387.703
$E_{\rm corr}^1$	0.310	0.310
$E_{\rm VBM}$ ( ΔN = 0.2)	6.735
E_g	4.404
$E_{\rm B}$	405.019	396.550
実験値	397.4 :sup:	a)

a) [Mahmood03]

内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。なお、式 eq:newxps14 を用いた場合、 $E_{\rm B}$ は以下のようになりました。

AlN (N 1s) のXPS 計算結果 ( 単位：eV )。式 eq:newxps14 により評価。

	内殻電子のスピン考慮なし	内殻電子のスピン考慮あり
$E_{\rm B}$	406.976	398.889

w-AlN 結晶 (Al 2p3/2)

w-AlN結晶 (Al 2p3/2) の例は samples/XPS/w-AlN 以下に配置されています。そのディレクトリー構成は下記の通りです。

w-AlN 結晶 (Al 2p3/2) XPS計算のディレクトリー構造

第一階層	第二階層	説明
Final_Al2p		終状態の入力ファイル
	q_0_no_opencore	中性・内殻スピンなし
	q_0_with_opencore	中性・内殻スピンあり
	q_1_no_opencore	電荷+1・内殻スピンなし
	q_1_with_opencore	電荷+1・内殻スピンあり
Initial		初期状態の入力ファイル
	q_0	中性
	q_0.2	電荷+0.2
	q_0.5	電荷+0.5
	q_1	電荷+1

計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.140 Å、c= 5.040 Å ) 3 × 3 × 2 倍しました。

w-AlN結晶 (Al 2p3/2)

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Monkhorst-Pack 2×2×2
交換相関相互作用	GGAPBE, PAW
SCF 収束条件 [Ha/atom]	1.0E-8
内殻電子のスピン	sw_xc_opencore_ae_only = on
CIAOの入力	sw_write_core_full 1
(補正項計算で使用する)誘電率	※8.23 (ab面内), 9.74 (c軸)

※https://materialsproject.org/materials/mp-661/

w-AlN (Al 2p3/2) のXPS 計算結果 ( 単位：eV )。式 eq:newxps13 により評価。

内殻	電子のスピン考慮なし\| 内殻	電子のスピン考慮あり
$E_{\rm final} (N-1) - E_{\rm initial} (N)$	64.980	64.212
$E_{\rm corr}^1$	0.327	0.327
$E_{\rm VBM}$ ( ΔN = 0.2)	6.735
E_g	4.404
$E_{\rm B}$	74.244	73.476
実験値	73.3 :sup:`	a)`

a) [Mahmood03]

内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。なお、式 eq:newxps14 を用いた場合、 $E_{\rm B}$ は以下のようになりました。

w-AlN (Al 2p3/2) のXPS 計算結果 ( 単位：eV )。式 eq:newxps14 により評価。

	内殻電子のスピン考慮なし	内殻電子のスピン考慮あり
$E_{\rm B}$	76.221	75.446

fcc Pt結晶 (Pt 4f)

fcc PtN結晶 (Pt 4f) の例は samples/XPS/Pt 以下に配置されています。そのディレクトリー構成は下記の通りです。

fcc Pt結晶 (Pt 4f) XPS計算のディレクトリー構造

第一階層	第二階層	説明
Final		終状態の入力ファイル
	q_0_no_opencore	中性・内殻スピンなし
	q_0_with_opencore	中性・内殻スピンあり
Initial		初期状態の入力ファイル
	q_0	中性

計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.963 Å )を 2 × 2 × 2 倍しました。

fcc-Pt結晶 (Pt 4f) XPS計算に用いた条件

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Monkhorst-Pack 2×2×2
交換相関相互作用	GGAPBE, PAW
SCF 収束条件 [Ha/atom]	1.0E-8
内殻電子のスピン	sw_xc_opencore_ae_only = on
CIAOの入力	sw_write_core_full 1

fcc Pt (Pt 4f 7/2) のXPS 計算結果 ( 単位：eV )。式 eq:newxps14 により評価。

内殻	電子のスピン考慮なし \| 内殻	電子のスピン考慮あり
$E_{\rm B}$	70.673	70.380
実験値	71 ^a)

a) http://techdb.podzone.net/xps/index.cgi?element=Pt

fcc Pt (Pt 4f 5/2) のXPS 計算結果 ( 単位：eV )。式 eq:newxps14 により評価。

内殻	電子のスピン考慮なし \| 内殻	電子のスピン考慮あり
$E_{\rm B}$	74.333	74.040
実験値	71 ^a)

a) http://techdb.podzone.net/xps/index.cgi?element=Pt

なお、CIAOにおいて内殻電子の運動エネルギーを非相対論で評価した擬ポテンシャルを用いた場合、4f 軌道の結合エネルギー ( 4f7/2 と 4f5/2 の重みつき平均 ) は 83.052 eV となり、実験値と離れてしまいます。

O ₂ 分子 (O 1s)

O ₂ 分子 (O 1s XPS計算のディレクトリー構造

第一階層	説明
final_spin_antiparallel	終状態・スピン反平行
final_spin_parallel	終状態・スピン平行
initial	始状態

計算条件は以下のとおりです。

O ₂ 分子 (O 1s) XPS計算に用いた条件

平面波カットオフ [Ry]	30.0
電荷密度カットオフ [Ry]	270.0
k 点サンプリング	Γ 点のみ
交換相関相互作用	GGAPBE, PAW
単位胞の1辺 [Å]	14.0
SCF 収束条件 [Ha/atom]	1.0E-8
内殻電子のスピン	sw_xc_opencore_ae_only = on
CIAOの入力	sw_write_core_full 1
（補正後計算で使用する）誘電率	※1.00

※真空中

O ₂ 分子 (O 1s) XPS計算結果 ( 単位：eV )。式 eq:newxps15 により評価。

	S=1/2	S=3/2
spin orientationの指定	anti parallel	parallel
$E_{\rm final} \left( N-1 \right)- E_{\rm initial} \left( N \right)$	542.000	541.297
$E_{\rm corr}^1$	1.436	1.439
$E_{\rm B}$	543.436	542.736
実験値 ^a)	544.2	543.1

a) https://t-ozaki.issp.u-tokyo.ac.jp/vps_pao_core2019/O/index.html and references therein.

2つのスピン状態とも、実験値に近い値が得られました。

Core Level Shift (CLS)の解析

内殻準位シフト(Core Level Shift, CLS) の計算はエネルギーの相対値が分かればよく、上で説明したような手続きは不要です。ここではCLSを計算する方法について解説します。

内殻準位シフト

内殻準位とは原子の深い電子準位のことであり、化学結合には寄与しないくらい原子に強く局在したものです。例えばシリコン原子の場合、14個の電子は(1s)²(2s)²(2p)⁶(3s)²(3p)² のように5個の準位を占有しますが、この中で1s、2s、2pが内殻準位です。相対論的なシリコン原子の電子準位をCIAOでGGA計算すると

Energy levels [All-electron]
Element ---> Si

symm j Energy (Ha) Energy (eV) nocc focc

1s 1/2 -65.6258330748 -1785.7697073691 2 2.00000 2s 1/2 -5.1250077353 -139.4585506190 2 2.00000 2p 1/2 -3.5260321902 -95.9482139488 2 2.00000 2p 3/2 -3.5022484901 -95.3010265676 4 4.00000 3s 1/2 -0.3967820153 -10.7969875601 2 2.00000 3p 1/2 -0.1503011244 -4.0899015276 2 2.00000 3p 3/2 -0.1491437813 -4.0584086215 4 0.00000

Total number of electrons 14.00000

となり、内殻準位は化学結合の目安である数eVより圧倒的に深いことがわかります。そのため内殻準位の波動関数は隣接した原子と重なりを持たず、エネルギー分散の無い離散準位を生じます。

シリコン原子の内殻準位を観測する手段として 100～130 eVの軟X線領域の単色光を照射し、放出される光電子の運動エネルギーを測定する実験方法があります [Landemark92] 。ちなみに、エネルギー可変の単色光は加速器のアンジュレータから放射されたシンクロトロン放射光を用います。上記のSi原子の電子準位を参考にすれば、この方法では2p準位から放出される光電子に着目していることがわかります。 2p準位はスピン軌道相互作用のために2p_3/2と2p_1/2に0.64 eV程度分裂します。

光電子放出の理論からよく知られているように、照射光のエネルギーをhνとすれば放出される光電子の運動エネルギーE_kinは

E_kin = hν − W − (E_F − E_c)

となります。ここで、Wは仕事関数、E_FはFermi準位、E_cは内殻準位です。一般に結晶表面では電子のしみ出しによる電気二重層が形成されるため、内殻電子の感じるポテンシャルは表面から外側になればなるほど浅くなります。そのため、内殻準位も表面では浅くなります。この関係を模式的に示したのが advanced_XPS_image1 です。その他、化学結合にともない原子のポテンシャルが上下するので、内殻準位はこれに連動した効果も受ける。この2つの効果のため内殻エネルギー準位E_cは原子により異なる値をとります。内殻準位シフトの実験では、表面から十分内部に入ったバルク位置の原子による内殻準位を基準にして表面付近の化学結合が異なる原子の内殻準位のエネルギーの差を測定することにより、表面付近の原子の化学結合状態や構造を推定します。

$光電子放出過程のエネルギープロファイル：図の右側が結晶表面、左側がバルクである。Fermi準位\mathbf{E}_{\mathbf{F}}は一定であるのに対して、内殻準位\mathbf{E}_{\mathbf{c}}は原子により異なり、表面付近ではバルクにくらべて内殻準位が浅くなる。$

光電子放出過程のエネルギープロファイル：図の右側が結晶表面、左側がバルクである。Fermi準位E_Fは一定であるのに対して、内殻準位E_cは原子により異なり、表面付近ではバルクにくらべて内殻準位が浅くなる。

さて、内殻準位シフトを第一原理計算するためには、eq:xps1 によれば個々の原子の内殻準位E_cを計算します。このエネルギー値E_cは内殻電子も扱う全電子計算を行えば直ちに得られますが、価電子のみを扱う擬ポテンシャル法からは得られません。この理由から、擬ポテンシャル法では計算できないと思われますが、この問題を見事に解決したのがSchefflerのグループ [Pehlke93] です。 Schefflerらは advanced_XPS_image2 に示したように、始状態(initial state)と終状態(final state)の違いに着目しています。

$光電子放出過程の始状態と終状態：(I)はバルク位置の原子による光電子放出、(II)は表面付近の原子による光電子放出を表している。始状態(initial state)では、入射光(\mathbf{\text{hν}})と結晶が存在するのみなので(I)と(II)は同じ状態である。終状態(final state)では、放出された光電子の運動エネルギー\mathbf{E}_{\mathbf{\text{kin}}}と内殻正孔(赤丸)1個をもった結晶の全エネルギー\mathbf{E}_{\mathbf{\text{tot}}}は(I)と(II)で異なるが、和\mathbf{E}_{\mathbf{\text{kin}}} + \mathbf{E}_{\mathbf{\text{tot}}}は等しい。$

光電子放出過程の始状態と終状態：(I)はバルク位置の原子による光電子放出、(II)は表面付近の原子による光電子放出を表している。始状態(initial state)では、入射光(hν)と結晶が存在するのみなので(I)と(II)は同じ状態である。終状態(final state)では、放出された光電子の運動エネルギーE_kinと内殻正孔(赤丸)1個をもった結晶の全エネルギーE_totは(I)と(II)で異なるが、和E_kin + E_totは等しい。

advanced_XPS_image2 では、光電子が放出される原子の位置が異なる2つの場合を示しています。(I)は原子がバルク位置にある場合、(II)は原子が表面付近にある場合です。始状態では、考えている系には入射光(hν)と結晶が存在するのみなので(I)と(II)は同じ状態です。一方、終状態では、系には放出された光電子(運動エネルギーE_kin)と内殻正孔(赤丸)が残された結晶(全エネルギーE_tot)が存在します。内殻正孔とは、光電子が放出される前に占有していた内殻準位に残された電子のぬけがらです。 (I)と(II)ではE_kinとE_totはそれぞれ異なるが、始状態が同一ということから、終状態の系の全エネルギーも同じ値にならなければならないという次の重要な関係式が導かれます。

$$E_{\rm{kin}}\left(\rm{I}\right)+E_{\rm{tot}}\left(\rm{I}\right) = E_{\rm{kin}}\left(\rm{II}\right) + E_{\rm{tot}}\left(\rm{II}\right)$$

内殻準位シフトΔE_kinを次式で定義します。

$$\Delta E_{\text{kin}} \equiv E_{\text{kin}}(\rm{II}) - E_{\text{kin}}(\rm{I})$$

ゆえに、eq:xps2 を用いれば

$$\Delta E_{\text{kin}} \equiv E_{\text{tot}}(\rm{I}) - E_{\text{tot}}(\rm{II}) = - \Delta E_{\text{tot}}$$

が得られます。 eq:xps2 の右辺は擬ポテンシャル法でも計算することができます。 Schefflerら [Pehlke93] はこの考えに基づきSi(100)表面の表面内殻準位シフトを擬ポテンシャル法で計算し、Landemarkらの実験結果 [Landemark92] を理論的に説明しています。内殻準位E_cと結晶の全エネルギーE_totとの関係は、eq:xps1 により

ΔE_kin = ΔE_c

が成り立つので、eq:xps4 と比較して

ΔE_c = − ΔE_tot(6)

が導かれます。

実験的に得られる情報を模式的に表した図。

内殻光電子放出のスペクトル

実験的に得られる情報を模式的に表すと advanced_XPS_image3 に示したようなスペクトルです。ここでは、電子間クーロン相互作用によるオージェ過程は考えません。横軸は光電子の運動エネルギーE_kin、縦軸は光電子数(強度)です。また、横軸は内殻準位の結合エネルギーE_bindとみることもできて、この場合は左側ほど結合が強くなっています。結晶バルクからの光電子スペクトルは赤色で示したようにE_kin(I)に唯一のピークをもちます。これに対して、結晶表面付近からの光電子スペクトルは青色で示したようにいくつかのピークE_kin(II)に分かれます。各ピークの強度はそのピークに関係した原子数の比から定まります。電気二重層の効果のためピーク位置は通常右側にずれます。光電子は結晶内部の非弾性散乱でエネルギーを失うので、バルクのピークは結晶表面からせいぜい数nm程度の深さの原子によるものであることを注意してください。計算できるものは、バルク位置からピーク位置のずれと各ピークの相対的な強度です。

内殻正孔を含む原子の擬ポテンシャル

第一原理擬ポテンシャルバンド計算法では、内殻電子状態を凍結させ、価電子状態を自己無撞着に計算します。そのためバンド計算の段階では内殻に正孔を生じさせることができません。そのかわり内殻正孔を含む擬ポテンシャルを作成することは可能です。 Schefflerら [Pehlke93] は、光電子放出時に生成された内殻正孔はまわりの電子によりすみやかに遮蔽(screening)されるものと仮定しました。しかし、格子を変形させるだけの時間はないとしています。また、電子はドープされた不純物から無尽蔵に補給される、すなわちFermi面は不純物準位にピンニングされる、と考えて電気的に中性が保たれているものとしています。これらの仮定により計算されたスペクトルは実験結果とよく一致しています。

以上の考え方に従えば、内殻正孔を含む原子の擬ポテンシャルを次のように作成します。例えば2p準位に内殻正孔を含むシリコン原子の場合、

14個の電子からなる(1s)²(2s)²(2p)⁶(3s)²(3p)²の電子配置の中性シリコン原子において、2p準位の電子1個を3p準位に移動して (1s)²(2s)²(2p)⁵(3s)²(3p)³とした励起状態の中性シリコン原子を考える
その原子の全電子計算を自己無撞着におこなう
価電子として(3s)²(3p)³の5個の電子をはがしてイオン化する
このようにして、内殻正孔を有し5価にイオン化したシリコン原子の擬ポテンシャルができる

とします。これは advanced_XPS_image4 の2つの終状態(final state)のうち右側のscreened hole(遮蔽正孔)に対応します。

スカラー相対論での擬ポテンシャル計算では、内殻準位のスピン軌道分裂は考慮せず縮重度に対して重み平均をとったものを内殻準位とします。内殻準位のシフト量の計算では、この仮定により結果は変わりません。

Si原子の電子配置：左図は基底状態の電子配置、右図は内殻正孔が生じた場合の終状態の電子配置である。

計算の実行方法

内殻正孔を含む原子の擬ポテンシャルを用いて内殻準位シフトをバンド計算します。

バンド計算は次のように行います。

シリコン表面の座標を作成し、通常のシリコン擬ポテンシャルを用いて各原子に力が働かなくなるまで十分に格子緩和させる
その座標を用いて、原子1個を内殻正孔を含む擬ポテンシャルに置き換えて全エネルギーを計算する
すべての原子に対して順番に置き換えて全エネルギーを計算する
バルク位置の原子を決めて、これを基準にして他の配置の全エネルギーの差をとり、運動エネルギーのシフト量を計算する
強度を含めたスペクトルとして表示する

計算例：Si(100)表面

Si(100)表面のモデル

まずはじめに、Si(100) p(2 × 2)再構成表面のモデルを作成します。 Si(100)表面はダングリングボンド数を減らすためダイマー構造をとり、それらのダイマーは列方向に交互に傾くこと(バックリング)で安定化することが知られています。この時、バックリングに連動してダイマーの低い方から高い方に電子が移動します。実験的にはSi(100) c(4 × 2)構造が安定となるが、計算ではバックリングの性質が似通ったp(2 × 2)構造を扱います。この方が計算量を減らせるからです。この仮定はSchefflerら [Pehlke93] と同様です。

$Si(100)表面のモデル：表面平行方向の周期は\mathbf{2} \times \mathbf{2}(左図)、深さ方向は8層(右図)のモデルを採用する。8層めと9層目の中間に反転中心を考えているので、実際の計算では全部で16層64原子のスラブモデルとして扱う。単位胞の等価でない各原子に\mathbf{1}\mathbf{u}、\mathbf{1}\mathbf{d}、\mathbf{2}、\mathbf{3}\mathbf{u}、\mathbf{3}\mathbf{d}、などのラベルをつける。$

Si(100)表面のモデル：表面平行方向の周期は2 × 2(左図)、深さ方向は8層(右図)のモデルを採用する。8層めと9層目の中間に反転中心を考えているので、実際の計算では全部で16層64原子のスラブモデルとして扱う。単位胞の等価でない各原子に1u、1d、2、3u、3d、などのラベルをつける。

出発点の表面モデルは安定構造となってなければならず、十分に格子緩和します。原子に働く力の最大値を5 × 10^− 4程度以下に抑えるのに長時間を要しました。 advanced_XPS_image5 は格子緩和されたされた後の安定構造を表示しています。計算条件は(原子単位)、擬ポテンシャル：Si_ggapbe_nc_01.pp、交換相関ポテンシャル：GGA-PBE、カットオフエネルギー：k_c(wf) = 3.5、 k_c(chg) = 7.0、ユニットセル：a₁ = 14.6816015290、a₂ = 14.6816015290、a₃ = 60.0000000000、 k点：4 × 4 × 1、です。その結果、緩和された座標と力は以下のようになります。

x             y              z             fx          fy          fz

1 11.654665468 7.340800747 19.731672033 -0.000108 0.000000 0.000136 2 10.943522273 0.000000002 18.308554343 -0.000156 0.000000 -0.000260 3 7.408260709 7.340800738 18.308043049 0.000156 0.000000 -0.000288 4 6.696493833 -0.000000008 19.730653256 0.000095 0.000000 0.000136 5 12.644826676 10.798805752 16.929970763 -0.000186 -0.000125 -0.000051 6 12.644826664 3.882795799 16.929970692 -0.000187 0.000125 -0.000052 7 5.707073513 11.223532172 16.928492210 0.000184 0.000135 -0.000066 8 5.707073576 3.458069376 16.928492250 0.000184 -0.000135 -0.000066 9 9.176755229 11.011875138 14.003484695 0.000002 -0.000042 -0.000291

10 9.176755074 3.669726457 14.003484706 0.000002 0.000042 -0.000291 11 1.834472402 11.011012299 14.542706201 -0.000005 -0.000002 0.000195 12 1.834472400 3.670589245 14.542706195 -0.000005 0.000002 0.000196 13 9.235196168 7.340800802 11.460786586 0.000018 0.000000 -0.000178 14 9.118076207 0.000000001 11.459447663 -0.000002 0.000000 -0.000198 15 1.882358440 7.340800774 11.863174990 0.000085 0.000000 0.000607 16 1.785205792 -0.000000009 11.863763387 -0.000075 0.000000 0.000566 17 12.980433730 7.340800764 9.066215802 0.000364 0.000000 -0.000060 18 12.895937822 0.000000005 9.076888750 0.000103 0.000000 -0.000067 19 5.455961640 7.340800779 9.081857152 -0.000145 0.000000 -0.000096 20 5.370066620 -0.000000012 9.069872950 -0.000307 0.000000 0.000025 21 12.895952884 11.013570694 6.478028955 0.000077 -0.000071 0.000099 22 12.895952919 3.668030852 6.478029019 0.000078 0.000071 0.000099 23 5.455699091 11.008682705 6.484530647 -0.000072 0.000047 0.000112 24 5.455699061 3.672918843 6.484530696 -0.000072 -0.000047 0.000112 25 9.172291984 11.011379982 3.920759253 -0.000013 -0.000002 0.000122 26 9.172291976 3.670221565 3.920759374 -0.000013 0.000002 0.000122 27 1.838705787 11.010982242 3.857185244 0.000014 0.000000 -0.000054 28 1.838705789 3.670619311 3.857185195 0.000014 0.000000 -0.000053 29 9.176000956 7.340800765 1.297682500 -0.000017 0.000000 0.002104 30 9.176000956 0.0 1.297682500 0.000472 0.000000 0.002118 31 1.835200191 7.340800765 1.297682500 -0.000774 0.000000 -0.002772 32 1.835200191 0.0 1.297682500 0.000325 0.000000 -0.002761

ここでは、反転対称を考えて半分の32原子の結果を示しています。第8層の4個の原子(29番～32番)は結晶がつぶれないように固定したので力が発生しています。

表面内殻準位シフトの計算

次に、表面内殻準位シフト(Surface Core Level Shift, SCLS)の計算を行います。格子緩和された表面モデルの原子位置を固定したまま、シリコン原子を順に内殻正孔を含む擬ポテンシャルに置き換えて全エネルギーを計算します。

この結果を advanced_XPS_image6 に示します。縦軸は光電子の運動エネルギー(単位はeV)のバルクからのずれです。この場合、バルクを3_u位置の原子にとっています。光電子は結晶内部で非弾性散乱されるので、あまり深い位置の原子からの光電子は実験的に観測されません。観測されるのは表面から数nmといわれています。このため、バルクを3_u位置としたことは正当であると考えられます。その他、3_uと3_dの中間にとる方法も考えられるが、advanced_XPS_image5 を参考にすれば3_uはダイマー列の外側、3_dはダイマー列の内側にあるので 3_dからの放出強度は実験的には抑えられることが考えられるので、ここでは3_uをバルクとしました。

advanced_XPS_image6 によればscreenedとunscreenedの決定的な違いは第1層めのdownのピーク位置です。両者のちがいは電子数がunscreenedの方がscreenedより少ないことです。そのためunscreenedではdown位置のシリコン原子が電子により遮蔽されません。 screenedでは遮蔽されるので、このちがいのためにdownのピーク位置が異なったと考えられます。

$SCLSの運動エネルギー：それぞれ青丸はunscreened、赤丸はscreenedの内殻正孔擬ポテンシャルを用いた結果である。赤丸はSchefflerと同じ方法である。縦軸は光電子の運動エネルギー(単位はeV)のバルクからのずれである。この場合、バルクを\mathbf{3}_{\mathbf{u}}位置の原子にとっている。$

SCLSの運動エネルギー：それぞれ青丸はunscreened、赤丸はscreenedの内殻正孔擬ポテンシャルを用いた結果である。赤丸はSchefflerと同じ方法である。縦軸は光電子の運動エネルギー(単位はeV)のバルクからのずれである。この場合、バルクを3_u位置の原子にとっている。

Schefflerらの論文 [Pehlke93] にある結果と比較するために、advanced_XPS_image7 に表面から3層までの原子によるSCLSを表示します。縦軸は強度であり、各強度は原子数の比をとっています。上から順に、unscreened、screenedの各計算値、実験値を示しています。実験値はLandemarkら [Landemark92] によるものです。

SCLSの強度：縦軸はSCLSの強度であり、各強度は原子数の比をとっている。上から順に、unscreened、screenedの各計算値、実験値を表示した。実験値はLandemarkら [Landemark92] によるものである。実験に合うのはscreenedである。XPS用の内殻正孔擬ポテンシャルはscreenedで作成しなければならない。

advanced_XPS_image7 の計算結果は実験結果 [Landemark92] をよく再現し、また、Schefflerらの論文の結果 [Pehlke93] に一致しています。実験に合うのはscreenedであることがわかります。このため、XPS用の内殻正孔擬ポテンシャルはscreenedで作成しなければなりません。

not latex

参考文献

ELNES / XANES解析機能

はじめに

EELSは電子が試料を透過する際のエネルギー損失を測定し、XAFSは試料にX線を照射した際に現れるスペクトルを測定する実験手法です。ここでは、電子線あるいはX線により、内殻電子が非占有状態に遷移する現象を対象とします。また、広範なエネルギー領域の中で、特に、吸収端付近のスペクトル構造を計算するものとします。この狭いエネルギー領域では、EELS は ELNES、XAFS は XANES と呼ばれています。

ELNESとXANESの二重微分散乱断面積は、双極子近似のもとで、以下のように表されます。

$$\begin{aligned} \frac{\partial^2 \sigma}{\partial E \partial\Omega} &\propto \sum_{i,f} \left| \left< i \left| \mathbf{q} \cdot \mathbf{R} \right| f \right> \right|^2 ({\rm ELNES}), \\\ \frac{\partial^2 \sigma}{\partial E \partial \Omega} &\propto \sum_{i,f} \left| \left< i \left| \mathbf{e} \cdot \mathbf{R} \right| f \right> \right|^2 ({\rm XANES}) \end{aligned}$$

ここで |i> |f> は始状態 (基底状態) と終状態（ある原子の内殻電子が飛び出し、非占有状態に遷移した状態）、qは非弾性散乱時における波数の変化、eは入射X線の偏光ベクトルである。したがって、XANESの計算は、ELNES計算にてqを eと読み替えます。内殻軌道から状態n, kで指定される伝導バンドへの遷移確率は、 [Gao08] に従い、

$$\left< \phi_c \left|r_{\alpha}\right|\psi_{nk} \right> = \left<\phi_c\left|r_{\alpha} \right| \tilde{\psi}_{nk} \right> + \sum_i \left( \left<\phi_c \left| r_{\alpha} \right| \phi_i \right> - \left<\phi_c \left|r_{\alpha}\right|\tilde{\phi}_i\right>\right) \left<\tilde{p}_i | \tilde{\psi}_{nk}\right>$$

で与えられるものとします。ここで右辺第1項は、波動関数の soft part　の寄与、第2項は擬ポテンシャルを用いたことによる補正項です。 |ϕ_i> は原子軌道iの全電子波動関数、ϕ̃_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の間の遷移モーメントは

$$M_{i\rightarrow f} = \Braket{\psi_f | D | \psi_i}, D=\hat{\mathbf{\varepsilon}} \cdot \mathbf{r} + \frac{i}{2} \left(\hat{\mathbf{\varepsilon}} \cdot \mathbf{r}\right) \left( \mathbf{k}\cdot\mathbf{r}\right)$$

で表されます [Gougoussis09] 。ψ_i 及び ψ_f は、それぞれ始状態i及び終状態fの波動関数、$\hat{\mathbf{\varepsilon}}$ 及び k は入射光の偏光ベクトル及び波数ベクトルです。D の第1項が双極子、第2項が四重極子成分に対応します。

ただし、PAW法では、波動関数が規格直交化されないため、は以下のように書き替えられます。

$$\begin{aligned} M_{i \rightarrow f} &= \Braket{\tilde{\psi}_f | D | \tilde{\psi}_i} \\\ &+ \sum_{\mathbf{R}} \sum_{l'm'\tau '} \sum_{nlm} \Braket{\tilde{\psi}_f | \beta_{\mathbf{R},l'm'\tau '}^{v}} \Braket{\phi_{\mathbf{R},l'm'\tau '}^v | D | \phi_{\mathbf{R},nlm}^c} \Braket{\beta_{\mathbf{R},nlm}^c | \tilde{\psi}_i} \\\ &+ \sum_{\mathbf{R}} \sum_{l'm'\tau '} \sum_{nlm} \Braket{\tilde{\psi}_f | \beta_{\mathbf{R},l'm'\tau '}^{v}} \Braket{\tilde{\phi}_{\mathbf{R},l'm'\tau '}^v | D | \tilde{\phi}_{\mathbf{R},nlm}^c} \Braket{\beta_{\mathbf{R},nlm}^c | \tilde{\psi}_i} \end{aligned}$$

ここで、 β_{R, l′m′τ′}^v 及び β_R, nlm^c は、それぞれ価電子及び内殻電子軌道へのプロジェクタです。また、 ϕ_{R, l′m′τ′}^v 及び ϕ̃_{R, l′m′τ′}^v は、それぞれ価電子軌道の全電子及び擬波動関数です。同様に、 ϕ_R, nlm^c 及び ϕ̃_R, nlm^c は、それぞれ内殻電子軌道の全電子及び擬波動関数です。

　なお、XANESの計算は、始状態iとして、原子位置 (R₀) 及び内殻電子軌道 (n,l) を指定して行います。さらに、この軌道は原子近傍に局在していることを考えると、結局、

$$M_{i\rightarrow f} \simeq \sum_{i'm'\tau'} \sum_{m} \Braket{\tilde{\psi}_f | \beta_{\mathbf{R}_0,l'm'\tau'}^v} \Braket{\phi_{\mathbf{R}_0,l'm'\tau'}^v | D | \phi_{\mathbf{R}_0, nlm}^c} \Braket{\beta_{\mathbf{R}_0,nlm}^c | \tilde{\psi}_i}$$

を計算すればよいことがわかります [Gougoussis09] 。

スペクトル強度の計算

スペクトル強度は、

$$\begin{aligned} \varepsilon_2 \left(\omega\right) &= 4 \pi^2 \sum_{f} \left| M_{i\rightarrow f} \right|^2 \delta \left( E_f - E_i - \hbar \omega \right) \\\ &= 4 \pi^2 \sum_f \left[ \left| \Braket{\phi_f | \hat{\mathbf{\varepsilon}} \cdot \mathbf{r} | \psi_i} \right|^2 + \frac{1}{4} \left| \Braket{\psi_f | \left( \hat{\mathbf{\varepsilon}} \cdot \mathbf{r} \right) \left( \mathbf{k} \cdot \mathbf{r} \right) | \psi_i} \right|^2 \right] \delta \left( E_f - E_i - \hbar \omega \right) \end{aligned}$$

で与えられます [Cabaret10] 。ここで、 E_i 及び E_f は、始状態及び終状態の全エネルギー、 ℏω は入射光のエネルギーです。advanced_XANES_quad_eq3 を用いれば、双極子及び四重極子成分は、それぞれ

$$\begin{aligned} \varepsilon_2^{\rm E1-E1} \left(\omega\right) &= 4 \pi^2 \sum_f \left| \sum_{l'm'\tau'} \sum_m \Braket{\tilde{\psi}_f | \beta_{\mathbf{R}_0,l'm'\tau'}} \Braket{\phi_{\mathbf{R}_0,l'm'\tau'}^v | D^{(1)} | \phi_{\mathbf{R}_0,nlm}^c} \Braket{\beta_{\mathbf{R}_0,nlm} | \tilde{\psi}_i} \right|^2 \\\ &\times \delta \left(E_f - E_i - \hbar \omega\right), \\\ \varepsilon_2^{\rm E2-E2} \left(\omega\right) &= 4 \pi^2 \sum_f \left| \sum_{l'm'\tau'} \sum_m \Braket{\tilde{\psi}_f | \beta_{\mathbf{R}_0,l'm'\tau'}} \Braket{\phi_{\mathbf{R}_0,l'm'\tau'}^v | D^{(2)} | \phi_{\mathbf{R}_0,nlm}^c} \Braket{\beta_{\mathbf{R}_0,nlm} | \tilde{\psi}_i} \right|^2 \\\ &\times \delta \left(E_f - E_i - \hbar \omega\right), \\\ D^{(1)} &= \hat{\mathbf{\varepsilon}} \cdot \mathbf{r}, D^{(2)} = \left( \hat{\mathbf{\varepsilon}} \cdot \mathbf{r}\right) \left( \mathbf{k} \cdot \mathbf{r} \right) \end{aligned}$$

となります。なお、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
(後略)

計算例

SrTiO₃結晶の計算例を紹介します。入力ファイルは 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: Mo	nk (2×2×2) XANES: Monk (4×4×4)
汎関数 PAW	, GGAPBE
擬ポテンシャル Sr_ggap	be_paw_02.pp Ti_ggapbe_paw_002_nocorehole.gncpp2 Ti_ggapbe_paw_005_1s_corehole.gncpp2 O_ggapbe_paw_02.pp
スペクトルのブロードニング Gaussian, wid	th = 0.5 eV

また、偏光ベクトル、波数ベクトルの組み合わせとして、以下の2つを試しました。四重極子遷移において、aはt2g, bはeg 軌道成分をプローブします。

\	a	b
偏光ベクトルEの方向 [100] [11 波数ベクトルkの方向 [010] [1-	0] 10]

結果を以下に示します。 SrTiO₃のd軌道は、t2gの方がeg軌道よりもエネルギーが低いため、四重極子遷移で最初にピークが現れるのは t2g 軌道をプローブするa の場合でです ( advanced_XANES_quad_fig1 参照)。advanced_XANES_quad_fig2 は、吸収端近傍の拡大図です。

SrTiO₃ のXANESスペクトル

SrTiO₃ のXANESスペクトル拡大図。

文献

ボルン有効電荷

概要

平面波基底を採用する第一原理計算において、原子の電荷を一意に求めることは簡単なことではありません。ボルン有効電荷はベリー位相理論に基づいて計算されるものであり、電荷を決めるための最良の方法と考えられます。その基本的な考え方は、“原子が微小に変位した場合に発生する分極”を変位量で割ることによって電荷を求めることができる、というものです。分極も変位も3方向あり得るので、ボルン有効電荷は2階のテンソルとなります。なお、ボルン有効電荷の計算はギャップのある系でのみ計算することが可能な点にご注意ください。

ボルン有効電荷の計算の流れ

ボルン有効電荷は、次のような表式によって計算します。

$$Z_{\alpha \beta}^{*} = - \frac{\Omega}{q_e} \frac{\partial P_{\alpha}}{\partial u_{\beta}} = Z_{\rm ion} \delta_{\alpha \beta} + \sum_i \frac{f}{2\pi} a_{i\alpha} \cdot \frac{\partial \phi_i \left( u_{\beta} \right)}{\partial u_{\beta}} .$$

ここでは $Z_{\rm ion}$ イオンの電荷、はi番目の格子ベクトルのα成分、はi番目の逆格子ベクトルに沿ったベリー位相、は変位のβ成分です。ベリー位相の原子変位に対する微分は、差分近似から計算します。このことから、3通りの変位に対して3つの逆格子ベクトルのベリー位相を計算する必要があり、ボルン電荷を計算するためには考慮する原子数の9倍のベリー位相計算が必要になることが分かります（実際には変位しないケースも必要なので、必要なベリー位相計算は原子数×9 + 3 です）。また、ベリー位相の計算はekcalの固定電荷計算によって行うので、さらに固定電荷計算用の電荷を作成するSCF計算をベリー位相計算に先立って行う必要があります。この一連の計算を行うために、PHASE/0に付属するberry.plというPerlスクリプトを利用します。

計算は、PHASE/0に付属するberry.plというPerlスクリプトを利用して行います。berry.plは、binディレクトリーの下にあります。ボルン有効電荷の場合、berry.plは以下のような流れで計算を行います。

対象とする原子すべてのSCF計算を行う。
対象とする原子すべてを3方向変位させた計算を行う。
1., 2.で得られた電荷密度を入力とした固定電荷計算によって、ベリー位相計算を3つの逆格子ベクトルに対して行う。
3. で得られたベリー位相データを結合し、 advanced_eq_born1 を利用してボルン有効電荷を計算する。

入力として必要なのは、SCF計算およびベリー位相計算の入力の元となるテンプレート入力ファイルと、berry.plの動作を制御するコントロールファイルです。

入力データ

概要

ボルン有効電荷の計算は、複数回のPHASEによるSCF計算と、ekcalによるベリー位相計算が必要となります。PHASE/0には、これらの計算を正しい手順で実行するberry.plというPerlスクリプトが付属しています。berry.plを実行するためには、SCF計算およびベリー位相計算用のテンプレートとなる入力ファイルと、berry.plの振る舞いを制御するコントロールファイルを作成する必要があります。

SCF計算およびベリー位相計算用テンプレート入力データの作成

ボルン有効電荷の計算を行うため、“テンプレート入力データ”を作成します。作成するテンプレート入力は、SCF計算用とベリー位相計算用のテンプレート入力です。それぞれについてberry.plを実行するディレクトリーの下にディレクトリーを作成し、その下にSCF計算用の入力データと固定電荷計算用の入力データを作成します。SCF計算はPHASE, 固定電荷計算は2次元版の場合はPHASEもしくはekcal, 3次元版の場合はPHASEプログラムを利用します。作成する入力データは通常の計算と同じですが、以下のような点に留意して入力を作成してください。

座標データは、構造最適化や格子最適化がなされたものを採用してください。
構造最適化などが行われない設定にしてください。
計算は、実行ディレクトリーの一階層下において行われます。すなわち、テンプレート入力が置かれるディレクトリーと同じ階層で行われます。この点に留意してfile_names.dataにおける擬ポテンシャルの指定を行ってください。
通常、固定電荷計算用の入力のfile_names.dataファイルにおいてはファイルポインターF_CHGTを利用してSCF計算によって得られたF_CHGTファイル（電荷密度ファイル）を指定します。ボルン有効電荷計算においては実際に参照するSCF計算はberry.plが動的に決めるので、この部分を正しく指定する必要はありません（指定があっても、berry.plによって書き換えられてしまいます）。
固定電荷計算でPHASEを利用する場合以下のような設定をしてください。この設定を施すことによって、固定電荷計算は“k点を一点ずつ処理する”モードで動作します。ベリー位相計算はこのモードでないと正しく計算できません。

control{
  condition = fixed_charge
  fixed_charge_option{
    kparallel = one_by_one
  }
}

ベリー位相計算のk点サンプリングは、berry.plのコントロールファイルにおけるmesh1, mesh2, mesh3パラメーターによって決まるので、ekcalの入力パラメーターファイルにおけるksamplingブロックの設定は無視されます。

berry.pl用のコントロールファイルの作成

berry.plの振る舞いは、コントロールファイルを介して指定します。ボルン有効電荷計算の場合、典型的には以下のような内容になります（#で始まる文はコメント文）。コントロールファイルのファイル名はberry.pl実行時に指定する仕様となっているので、任意のものを採用してください。

#overall control
property = zeff
cpumax = 1000

#directories under which the template files reside
template_scf = scf
template_berry = berry

#parameters for the berry-phase calculation
atom_list = 1 3
displacement = 0.1
mesh1 = 6 6 15
mesh2 = 6 6 15
mesh3 = 6 6 15

#execution control
np = 4
ndir = 2
ne = 1
nk = 2
ne_b = 2
scf_command = mpiexec –np NP phase ne=NE nk=NK
berry_command = mpiexec –np NP ekcal ne=NE_B

#unit cell info, optional
a_vector = 5.01 0.0 0.0
b_vector = 0.0 5.01 0.0
c_vector = 0.0 0.0 5.01

この例からわかるように、パラメータ1つにつき1行を利用し、“キーワード=値”という形式でパラメータを指定します。ボルン有効電荷計算の場合に考慮する必要のあるキーワードとその説明を以下に記します。

キーワード	説明
property	どのような計算を行うかを指定します。zeff, piezo, strfrcのいずれかです。zeffを指定するとボルン有効電荷用のベリー位相計算が行われます。デフォルト値はzeffなので、ボルン電荷計算の場合未指定でも構いません。
cpumax	計算の最大時間を秒の単位で指定します。ここで指定した時間よりも経過時間が長い場合、計算はすみやかに終了します。0以下の値を指定すると、この条件では計算は終了しません。デフォルト値は-1。
stopcheck	計算停止条件を満たしているかどうかをチェックする間隔を秒の単位で指定します。デフォルト値は10。
length_unit	コントロールファイル中に利用される長さの単位を指定します。bohr, angstrom, nmのいずれかを指定します。デフォルト値はbohr。
template_scf	SCF計算用のテンプレートディレクトリーのディレクトリ名を指定します。デフォルト値はtemplate_scf。
template_berry	ベリー位相計算用のテンプレートディレクトリのディレクトリ名を指定します。デフォルト値はtemplate_berry。
atom_list	変位させる原子のIDを、空白区切りで指定します。
displacement	原子の変位量を指定します。property = zeffの場合に指定します。デフォルト値は0.1 bohr。
mesh1	1番目の逆格子ベクトルに沿ったベリー位相計算のメッシュパラメーターを空白区切りで n1 n2 J のように指定します。property = zeff および piezoの場合必須の指定です。
mesh2	2番目の逆格子ベクトルに沿ったベリー位相計算のメッシュパラメーターを空白区切りで n1 n2 J のように指定します。property = zeff および piezo の場合必須の指定です。
mesh3	3番目の逆格子ベクトルに沿ったベリー位相計算のメッシュパラメーターを空白区切りで n1 n2 J のように指定します。property = zeff および piezo の場合必須の指定です。
np	MPIプロセス数を指定します。デフォルト値は1。
ndir	ディレクトリー並列数を指定します。デフォルト値は1。
ne	バンド並列数を指定します。デフォルト値は1。
nk	k点並列数を指定します。デフォルト値は1。
ng	(三次元版のみ)G点並列数を指定します。デフォルト値は1。
ne_b	ベリー位相計算時におけるバンド並列数を指定します。デフォルト値は1。
ng_b	（三次元版のみ）ベリー位相計算時におけるG並列数を指定します。デフォルト値は1。
scf_command	SCF計算の実行方法を指定します。たとえば、 scf_command = mpirun –np NP phase ne=NE nk=NK (2D版) scf_command = mpirun -np NP phase ne=NE nk=NK ng=NG (3D版) などと指定します。NP, NE, NK, NGは、計算実行時に上述のnp, ne, nk, ngに置き換わります。ただし、ディレクトリ並列の数によっては“あまり”が発生することがあり、その場合はne=NE nk=NKng=NGは省略されて計算が投入されます。デフォルト値はmpirun phaseですが、利用している環境に合わせて適切な指定をする必要があります。
berry_command	ベリー位相計算の実行方法を指定します。たとえば、 berry_command = mpirun –np NP ekcal ne=NE_B (2D版) berry_command = mpirun -np NP phase ne=NE_B ng=NG_B (3D版) などと指定します。NE_Bは上述のne_bに、NG_Bはng_bに置き換わります。デフォルト値はmpirun ekcalです、利用している環境に合わせて適切な指定をする必要があります。
a_vector	a軸の三成分を空白区切りで指定します。必須ではありませんが、指定しておくとベリー位相計算のメッシュパラメーターの参照値を算出します。
b_vector	b軸の三成分を空白区切りで指定します。必須ではありませんが、指定しておくとベリー位相計算のメッシュパラメーターの参照値を算出します。
c_vector	c軸の三成分を空白区切りで指定します。必須ではありませんが、指定しておくとベリー位相計算のメッシュパラメーターの参照値を算出します。

注意点

mesh1, mesh2, mesh3パラメーターについて

ベリー位相の計算においては、対象としたい逆格子ベクトル b_i に垂直な面の面積分と、b_i に沿った線積分が実行されます。この積分のメッシュは、コントロールファイルのmesh1, mesh2, mesh3パラメーターを利用して行います。i番目の逆格子ベクトルに対し、meshi = n1 n2 Jと空白区切りでメッシュを指定します。ここで面積分のメッシュがn1 n2、線積分のメッシュがJです。

b_i 以外の2つの逆格子ベクトルを b_j とすると b_j を b_i に垂直な面に射影したベクトル、${\mathbf b}_j - \left| \mathbf{b}_j \right| \cos \theta_{ji} \frac{\mathbf{b}_i}{\left| \mathbf{b}_i \right|}$ 、が面積分のメッシュの見積もりの基準となるので、その長さから決めます。線積分のパラメータは、b_i の長さをもとに決定します。コントロールファイルにa_vector, b_vector, c_vectorの指定を行っておくと、この長さの計算（bohr^-1単位）とそこから見積もられる参考のメッシュパラメーターが以下のように標準出力に出力されます（あくまで参照値であり、得られる結果の妥当性を保証するものではありません）。

|b_para1|, |b_para2| and |b_perp| (in bohr^-1 units)
for reciprocal vector no. 1 : 0.172224346323159, 0.107572987734313, 0.198867545420854
for reciprocal vector no. 2 : 0.172224346322494, 0.107572987734313, 0.198867545421622
for reciprocal vector no. 3 : 0.198867545420854, 0.198867545421622, 0.107572987734313

reference value for mesh parameters n1, n2 and J
for reciprocal vector no. 1 : 8, 5, 19
for reciprocal vector no. 2 : 8, 5, 19
for reciprocal vector no. 3 : 9, 9, 10

並列計算について

並列計算の設定について注意すべき点を挙げます。berry.plによる計算は、通常のバンド、k点による並列(3D版の場合さらにG点並列)に加え、複数のディレクトリにまたがって並列計算を行う“ディレクトリ並列”によって行われます。したがって、ディレクトリ並列数（パラメーターndir）を2以上にする場合、SCF計算の場合はnp=ne×nkではなくnp=ndir×ne×nk (3D版の場合np=ne x nk x ngではなくnp=ndir x ne x nk x ng)となるように、ベリー位相計算の場合はnp=ne_bではなくnp=ne_b×ndir(3D版の場合np=ne_b x ng_bではなくnp=ne_b x ng_b x ndir)となるように並列数を調整してください。SCF計算とベリー位相計算とでバンド並列数が異なるのは、ベリー位相計算はk点並列に未対応のためです。また、ディレクトリー並列数を2以上にする場合、以下の要領でScaLAPACKを無効にしてください。

wavefunction_solver{
  submat{
    scalapack{
      sw_scalapack = off
    }
  }
}

3D版について

バージョン2019.02から、3D版でもベリー位相計算が実行できるようになりました。3D版でベリー位相計算を行うには、固定電荷計算にekcalではなくphaseを使います。また、以下のような設定を施すようにしてください。

control{
  fixed_charge_option{
    kparallel = one_by_one
  }
}

berry.plの実行

berry.plを引数なしで実行すると、以下のようなメッセージが得られます。

% berry.pl
Usage : berry.pl control [OPTIONS]

第一引数にコントロールファイルのファイル名を指定し、さらに必要に応じてオプションを指定して制御する仕組みになっています。

以下のようなコマンドを実行すると、コントロールファイルの解釈と解析のみ行います。

% berry.pl control –-mode=analyze

以下のようなコマンドを実行すると、コントロールファイルの解釈と解析のあと、計算用のディレクトリを作成します。

% berry.pl control --mode=gendir

以下のようなコマンドを実行すると、コントロールファイルの解釈と解析のあと、計算用のディレクトリを作成し、さらに計算を実行します。

% berry.pl control --mode=exec

--mode オプションのデフォルト値はgendir です。また、

% berry.pl --clean

とすると、berry.plが作成したディレクトリーなどを一括削除することができます。

実行すると、以下のようなログを出力しながらPHASEもしくはekcalによる計算が進行します。

berry.pl : script to calculate the berry phase for the PHASE System
Copyright (c) 2012-2013, IIS, The University of Tokyo
script start time : Tue Mar 21 18:12:13 2017
-- parsing the control file –
…
…
-- generating directories --
number of SCF directories : 7
number of berry directories : 21
-- doing SCF calculations --
running [mpiexec -n 4 $HOME/phase0/binN/phase ne=2 nk=2] under the
following directories
scf_a0
time spent in this calculation : 9 (s), total time : 9 (s)
…
…

まずコントロールファイルの中身が読み込まれ、その内容が出力されます。ついで、計算に必要なディレクトリーが作成されます。その後SCF計算が必要な回数だけ行われ、さらにベリー位相計算が必要な回数だけ行われます。実行中のコマンドと計算中のディレクトリーは、以下のような形式で出力されます。

running [mpiexec -n 4 $HOME/phase0/bin/phase ne=2 nk=2] under the following directories
scf_a0

作成される計算用ディレクトリーは、以下のようなものです。

ディレクトリー名	説明
scf_a0	参照用の、原子を変位させない系の SCF計算が行われるディレクトリー。
scf_aaid_uuid	aid番目の原子を、uidの方向に変位させた系のSCF 計算が行われるディレクトリー。aidは1始まりであり、テンプレート入力データに記述した原子配置の順序に対応する。uidは1, 2, 3のいずれかの値をとり、それぞれx, y, z方向に相当する。
berry_a0_ggid	参照用の、原子を変位させない系における、gid方向の逆格子ベクトルのベリー位相計算が行われるディレクトリー。gidは1, 2, 3のいずれかの値をとり、それぞれ1番目, 2番目, 3 番目の逆格子ベクトルに相当する。
berry_ aaid_uuid_ggid	aid番目の原子を、uidの方向に変位させた系における、gid方向の逆格子ベクトルのベリー位相計算が行われるディレクトリー。aidは1始まりであり、テンプレート入力データに記述した原子配置の順序に対応する。uidは1, 2, 3のいずれかの値をとり、それぞれx, y, z 方向に相当する。gidは1, 2, 3のいずれかの値をとり、それぞれ1番目, 2番目, 3 番目の逆格子ベクトルに相当する。

ボルン有効電荷計算の場合、ベリー位相計算のあとにボルン有効電荷計算用のPHASE計算が行われます。この計算は、zeffディレクトリーにおいて行われます。この計算は非並列で行われ、通常数秒から数十秒程度で終了します。

計算結果

計算が最後まで実行されると、以下のようなログが得られます。

--- Calculated effective charges ---
           [ 2.98266 0.00512 -0.00454 ]
Zeff( 2) = [ 0.00001 3.62664 -0.32666 ]
           [ -0.00010 0.27925 3.42264 ]
...
...
--- Symmetrized effective charges ---
           [ 2.98266 0.00000 -0.00000 ]
Zsym( 2) = [ 0.00000 3.62664 -0.32666 ]
           [ 0.00000 0.27925  3.42264 ]
...
...
--- Effective charges of all atoms ---
           [ 3.46565 -0.27885 -0.28289 ]
Zeff( 1) = [ -0.27885 3.14366  0.16333 ]
           [ 0.24184 -0.13963  3.42264 ]
...
...
--- Averaged effective charges ---
       [ 0.00151 -0.00000 0.00000 ]
Zave = [ 0.00000 0.00151 -0.00000 ]
       [ 0.00000 0.00000 -0.00587 ]
...
...
--- Corrected effective charges ---
           [ 3.46414 -0.27885 -0.28289 ]
Zeff( 1) = [ -0.27885 3.14215  0.16333 ]
           [ 0.24184 -0.13963  3.42850 ]
...
...

--- Calculated effective charges --- 以下に、計算された生のボルン電荷テンソルが、計算対象となった原子数分出力されます。次に、 --- Symmetrized effective charges ---以下に対称化を施したボルン有効電荷テンソルが出力されます。その次に、 ---Effective charges of all atoms ---以下に、あらわに計算しなかった原子のボルン有効電荷も含めた結果が出力されます。あらわに計算しなかった原子のボルン有効電荷テンソルは、対称性を利用して計算されます。 --- Averaged effective charges ---には全原子で平均したボルン有効電荷テンソルが出力されます。全原子で平均したボルン有効電荷はすべての要素がゼロになるはずなので、この総和則を満たすように補正したボルン有効電荷テンソルが---Corrected effective charges ---以下に出力されます。

基本的には、--- Corrected effective charges ---以下の結果が最もよい結果となります。ただし、対称性の与え方が適切でない場合、またはあらわに計算した原子の数が中途半端な場合に「あらわに計算しなかった原子のボルン有効電荷を、対称性を利用して計算する」ことができなくなってしまうことがあります。このような場合は、--- Calculated effective charges ---以下の結果を利用するようにしてください。

例題

AlNのボルン有効電荷を計算する例を紹介します。

入力データ

入力データは、PHASEによるSCF計算用のテンプレート入力、ekcalによるベリー位相計算によるテンプレート入力、そしてberry.pl用のコントロールファイルから成ります。

PHASEによるSCF計算用のテンプレート入力は、 samples/dielectric/lattice/born_revised/AlN/born/template_scf の下にあります。この入力に特殊な設定はありませんが、構造最適化がなされないよう原子のmobile属性値はすべてoffと設定されています。また、対称性を自動的に処理する機能を有効にしています。

ekcalによるベリー位相計算のテンプレート入力は、 samples/dielectric/lattice/born_revised/AlN/born/template_berry の下にあります。この入力も、通常の固定電荷計算の入力と全く同じです。

berry.plのコントロールファイルは、 samples/dielectric/lattice/born_revised/AlN/born/control です。その内容は以下の通り。

property=zeff
atom_list = 1 3
mesh1 = 6 6 15
mesh2 = 6 6 15
mesh3 = 6 6 15
scf_command = mpiexec -n NP $HOME/phase0/bin/phase ne=NE nk=NK
berry_command = mpiexec -n NP $HOME/phase0/bin/ekcal ne=NE_B
np = 1
ndir = 1
ne = 1
nk = 1
ne_b = 1

property = zeffとすることによって、ボルン有効電荷を計算することを指定しています。atom_listに13 とすることによって対象とする原子を指定しています。1番目の原子はAl, 3番目の原子はNです。AlNは計4つの原子から成る結晶ですが、2番目と4番目の原子のボルン有効電荷は対称性から求まるので計算の対象にはしていません。mesh1 mesh2, mesh3によって逆格子ベクトルのメッシュを指定しています。この例では、すべての逆格子ベクトルに対してn1 = n2 = 6, J = 15としています。scf_command, berry_commandによってPHASEおよびekcalの実行方法を指定しています。ここは、計算に利用する環境に合わせて適宜修正する必要があります。np, ndir, ne, nk, ne_bによって並列のさせ方を指定することができます。この例では、すべて1, すなわち非並列で計算を行う指定となっています。

計算の実行

以下の要領で、berry.plを実行することができます。

% berry.pl control --mode=exec

計算時間は利用するCPUやコンパイラーなどに大きく依存します。Intel Core i7-2600@3.40 GHzのCPUを搭載したマシンで実行したところ、ボルン有効電荷が得られるまで1,760秒かかりました。

計算結果

この計算によって得られた各原子のボルン有効電荷テンソルは、下記の通りです。

[ 2.50916 0.00000 -0.00000 ]

Zeff( 1) = [ 0.00000 2.50916 0.00000 ]

[ 0.00000 0.00000 2.64124 ] [ 2.50916 0.00000 -0.00000 ]

Zeff( 2) = [ 0.00000 2.50916 -0.00000 ]

[ 0.00000 0.00000 2.64124 ] [ -2.50916 -0.00000 0.00000 ]

Zeff( 3) = [ -0.00000 -2.50916 -0.00000 ]

[ 0.00000 -0.00000 -2.64124 ] [ -2.50916 -0.00000 0.00000 ]

Zeff( 4) = [ -0.00000 -2.50916 -0.00000 ]

[ 0.00000 -0.00000 -2.64124 ]

1番目と2番目の原子がAl, 3番目と4番目の原子がNです。

berry.plを使わずにボルン有効電荷を計算する方法

ボルン有効電荷の計算は、berry.plを利用せずとも行うことはできます。その手続きを説明します。

SCF計算

ボルン有効電荷を計算するためには、対象としたい原子を変位させたSCF計算を実行する必要があります。また、原子を変位させていないSCF計算も実行しておく必要があります。まずは変位させないSCF計算を通常通りに行います。ついで、対象とする原子ごとにその位置をx, y, z方向に微小量（たとえば0.1 bohr）変位させた計算を行います。変位させる計算を行うためには、入力パラメーターファイルのatom_listブロック以下に定義できるdisplacementブロックを利用すると便利です。

structure{
  …
  …
  atom_list{
    …
    …
    displacement{
      sw_displace_atom = on
      displaced_atom = 1
      ux = 0.1
      uy = 0
      uz = 0
    }
  }
}

sw_displace_atom = onとすると原子を変位させます。displaced_atomによって変位させる原子を指定します。ux, uy, uzによって変位のx,y,zの値を長さの指定します。

以上の入力が作成できたら、通常通りPHASEを必要な回数実行します。

ベリー位相計算

SCF計算のあとに、ベリー位相計算を行います。前節で行ったすべてのSCF計算をもとに、固定電荷計算の入力を作成します。各々のSCF計算について、3方向の逆格子ベクトルのベリー位相計算用入力を作成します。ベリー位相計算を実行するには、入力パラメーターファイルに以下のように最上位にberry_phaseブロックを作成し、設定を行います。

berry_phase{
  sw_berry_phase = on
  g_index = 1
  mesh{ n1 = 3, n2 = 3, J = 5 }
}

sw_berry_phase = onとするとベリー位相計算が有効になります。g_indexによって対象とする逆格子ベクトルを指定します。meshブロックを作成し、n1, n2, Jによってメッシュパラメーターを指定します。このほか、SCF計算に対応するdisplacementブロックを作成しておくことと、file_names.dataにおいて対応するSCF計算のディレクトリーの下の電荷密度データを指すようにすることなどを忘れないようにしてください。

以上の入力が作成できたら、通常通りekcalを必要な回数実行します。

ボルン有効電荷の計算

最後にボルン有効電荷を計算します。以下の手続きを踏みます。

ボルン有効電荷計算用のディレクトリーを作成します（たとえば、zeff）。
ベリー位相データファイルberry.dataをすべて連結し、ファイルの先頭にその数を指定します。たとえば、21のベリー位相計算を行ったならばファイルの先頭に21と記述します。
ボルン有効電荷計算用のディレクトリーに、SCF計算用の入力ファイルと2.で得たberry.dataファイルをコピーします。
コピーした入力ファイルからdisplacementブロックを削除し、以下を挿入します。ボルン有効電荷の計算においてSCF計算を行う必要はないので、controlブロックにおいてmax_iteration = 0とします。postprocessing以下の設定が、ボルン有効電荷を計算するための指示です。

control{
  condition = initial
  max_iteration = 0
}
postprocessing{
  polarization{
    sw_bp_property = on
    property = effective_charge
  }
}

PHASEを非並列で実行します。結果は、outputxxxファイルに、advanced_berry_results 節で説明する形式で出力されます。

有限電場計算機能 (バージョン2019.02以降)

概要

有限電場計算機能とは、一様有限電場下において分極と電場のエネルギーを通常のエネルギーに加えた“エンタルピー”を極小化する計算機能です。この計算機能によって、印加した電場に対する応答としての巨視的な分極を求めることができます。得られた分極から誘電率を計算することや、有限電場下での原子に働く力よりボルン有効電荷を得ることも可能となっています。

計算理論

一様有限電場下の計算は、 [Souza02] において提案されている処方箋に従って行われます。一様電場下におけるエンタルピーを、以下のように定義します。

$$F = E_{\rm KS} - \Omega \mathbf{P}_{\rm mac} \cdot \mathbf{E}$$

$E_{\rm KS}$ が通常のコーン・シャムハミルトニアンから得られる全エネルギー、Ω が単位胞の体積、$\mathbf{P}_{\rm mac}$ が巨視的な電気分極、E が一様有限電場です。$\mathbf{P}_{\rm mac}$ はベリー位相分極理論から計算することができます。E は入力パラメーターです。このエンタルピーを極小化するために必要な微分は下記のように計算できます。

$$\begin{aligned} \ket{g_{n\mathbf{k}}} &= \frac{f}{N_k} \hat{H}_{\rm KS} \ket{u_{n\mathbf{k}}} \\\ &+ \frac{ife}{4\pi} \sum_{i=1}^3 \frac{\mathbf{E}\cdot\mathbf{a}_i}{N_{\perp}^{(i)}} \sum_{m=1}^{M} \left[ \ket{u_{m\mathbf{k}_+^{(i)}}} S_{mn}^{-1} \left( \mathbf{k},\mathbf{k}_+^{(i)} \right) - \ket{u_{m\mathbf{k}_-^{(i)}}} S_{mn}^{-1} \left( \mathbf{k},\mathbf{k}_-^{(i)} \right) \right] \end{aligned}$$

$\ket{u_{n\mathbf{k}}}$ が周期的な波動関数です。k₊, k₋ は、逆格子ベクトルをb_iとするとそれぞれk+b_i/N_{i ,} k-b_i/N_iですS_mn(k,k′) は波動関数の重なりであり、$\Braket{u_{m\mathbf{k}} | u_{n\mathbf{k}'}}$ と表すことができます。

このようにして計算された勾配を利用してエンタルピーを最適化すると、一様電場下での巨視的な分極 $\mathbf{P}_{\rm mac}$ を計算することができます。この量から、さらに様々な誘電応答解析を実施することも可能です。たとえば、電場を印加した状態の分極から差分によって誘電率を得ることが可能です（微弱な電場 δE を与えた際に発生する分極 δ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}$

電場が0なら分極も0, もしくは原子間力も0というケースでは一点の計算で上述の量を計算することができますが、そうでない場合は2点以上の計算の結果から $\frac{\delta P}{\delta E}$ や $\frac{\delta F}{\delta E}$ の値を計算するようにしてください。

計算例

AlP結晶の計算を有限電場法および従来法で行いました。その入力ファイルは samples/FEF/AlP にあります。8原子の系を採用しました。計算結果を以下にまとめます。

計算結果

物性値計算結	果
電子分極 0.00 誘電率（有限電場法） 7.40 誘電率（従来法） 7.40	0050916
ボルン電荷（有限電場法） Al : 2.25, P	: -2.27
ボルン電荷（従来法） Al : 2.227	, P : -2.227

従来法とよく一致する結果が得られました。ボルン電荷の方も従来法とよく一致しています。有限電場法のボルン電荷の和がゼロになっていないのは、数値誤差によるものと考えられます。

注意点

有限電場機能を利用する場合、まずは通常のSCF計算を実施し、収束した波動関数と電荷密度を得ておくことが推奨されます。したがって、計算の手続きは以下のようになります。
1. 対象としたい系のSCF計算を行う。
2. 有限電場計算用のディレクトリーを作成し、そこに1. で利用した入力ファイルと得られた波動関数データと電荷密度データをコピーする（ファイル名はそれぞれzaj.dataとnfchgt.data）
3. 入力パラメーターファイルに有限電場計算用の入力の設定をほどこし、さらに初期波動関数・初期電荷密度はファイルから読み込むように設定する。具体的には下記のような記述を行う。
分極を収束させるためにはSCF計算の収束判定として通常よりも厳しいものを採用する必要があることが多いです。Pelの収束具合を確かめながら必要な収束判定条件を決めていくことが推奨されます。上述のようにinitial_wavefunctionsとinitial_charge_densityをfileに設定しておけばそれまでの計算の結果を読み込んで動作するので、それまで行った計算は無駄にはなりません。
原子間力の計算結果からボルン電荷を得ることができますが、この機能はノルム保存擬ポテンシャルのみ対応しています。ウルトラソフト擬ポテンシャルの場合に得られる原子間力は不正確なので利用しないでください。

not latex

参考文献

陽電子寿命解析

機能の概要

電子の反粒子である陽電子は、電子と同じ質量を持ち、正の電荷を持ちます。陽電子は電子と対消滅し、γ線を放出します。この陽電子消滅の現象を利用して、結晶の品質や、欠陥の研究が可能です。ここで、陽電子消滅実験から有用な情報を引き出すのに、第一原理計算に基づいて陽電子寿命を予測し、実験結果と比較することが重要です。PHASEには、完全結晶における陽電子寿命を予測する機能があります。

陽電子寿命解析、以下の手順で計算を行います。

はじめに通常の電子状態計算（バンド計算）を実行します。PHASEには、擬ポテンシャル・平面波法が実装されているので、この手法に基づいて計算を行います。バンド計算により、価電子の電荷密度 ρ_vが得られます。全電子の電荷密度は次式で与えられます。

ρ_e = ρ_v + ρ_c

ここで、ρ_cは、コア電子の電荷密度です。CIAOで作成し公開されている擬ポテンシャルのデータファイルには、自由な原子におけるコア電子の電荷密度の情報が含まれています。このデータを読み込み advanced_positron_eq1 を評価します。

陽電子波動関数ψ₊は次式で与えられます（原子単位）。

$$\left\lbrack - \frac{1}{2}\Delta - \int_{}^{}d\overrightarrow{r^{'}}\frac{\rho_{e}\left( \overrightarrow{r^{'}} \right) - \rho_{n}\left( \overrightarrow{r^{'}} \right)}{\left| \overrightarrow{r} - \overrightarrow{r^{'}} \right|} + \mu_{c}\left( \rho\left( \overrightarrow{r} \right) \right) \right\rbrack\psi_{+}\left( \overrightarrow{r} \right) = \varepsilon\psi_{+}\left( \text{cr} \right)$$

ここで、μ_c は、電子・陽電子相関に由来するポテンシャルエネルギーであり、ρ_nは、原子核の点電荷を表します。いま、陽電子は固体中に1個しかないと仮定するので、最も安定な固有状態のみを求めればよいことになります。したがって、陽電子の固有状態はブリルアンゾーン中のΓ点に属します。この波動関数を平面波によって展開します。

ψ₊ = ∑_G⃗C_G⃗exp (iG⃗⋅r⃗)

ここで、逆格子周期ベクトルG⃗の和を有限に抑えるため、平面波の運動エネルギーの上限を設定します。

陽電子密度を求めます。

ρ_p(r⃗) = |ψ₊|²

電子及び陽電子の電荷密度を用いて、下記の式を評価し、陽電子寿命を計算します。

$$\frac{1}{\tau} = \pi r_{e}^{2}c\int_{}^{}dr\rho_{e}\left( r\right)\rho_{p}\left( r \right)\Gamma\left( \rho_{e} \right)$$

ここで、r_e は、電子の古典半径、cは光速を表します。Γ は、増大因子であり、電子・陽電子間の相関に由来するものです。PHASEでは、上式の評価において、下記の近似を用いています。

ρ_eΓ(ρ_e) ≅ ρ_vΓ(ρ_v) + ρ_cΓ(ρ_c)

この近似が成り立つためには、価電子とコア電子の分布の重なりが小さいことが条件となります。

本計算では、電子・陽電子間の相関に対して、局所密度近似を用いています。すなわち、相関ポテンシャルと増大因子は、均一な電子ガス中に陽電子が1個ある場合に計算された結果をもとに、電子密度の関数として与えられます。増大因子に関しては、次式が提案されています [Puska95] 。

Γ = 1 + 1.23r_s + 0.9889r_s^3/2 − 1.482r_s² + 0.3956_s^5/2 + r_s³/6

ここで、$\frac{4\pi}{3}r_{s}^{3} = 1/\rho_{e}$です。また、ギャップのある系、すなわち誘電体においては、金属よりも電子のスクリーニング効果が小さいため、次の補正を行う事を推奨します [Puska91] , [Nakamoto08] 。

Γ = 1 + 1.23r_s + 0.9889r_s^3/2 − 1.482r_s² + 0.3956_s^5/2 + (1−1/ε_ele)r_s³/6

ここで、ε_eleは電子系誘電率です。その値が実験により測定されていない場合には、UVSORにより、密度汎関数理論に基づいて評価することができます。計算手法の詳細については、 [Nakamoto08] を参照ください。

入力パラメータ

Si結晶中における陽電子寿命の計算を例とします。計算例題は、 samples/positron/Si です。

Si結晶中における陽電寿命計算の入力パラメータファイルにおいて、陽電子計算に関する部分のみを説明します。

Controlタグで陽電子寿命計算を有効に設定

Control{
   positron = BULK
}

Controlタグ中に　positron=BULK　と宣言すると、通常の電子状態計算(バンド計算）を行った後に、陽電子寿命計算を行います。

accuracy　タグで陽電子寿命計算のオプションを指定

accuracy{
  cutoff_pwf = 50.00 rydberg
  positron_convergence{
    num_extra_bands = 8
    delta_eigenvalue = 1.d-8 rydberg
    succession = 6
    num_max_iteration = 32000
    dtim = 0.01
    epsilon_ele = 12
  }
}

cutoff_pwf = 50.00 rydberg	これは、陽電子の波動関数を展開する際の[ `advanced_positron_eq3` 参照]カットオフエネルギーです。
positron_convergence｛｝	このタグの中で、陽電子波動関数を反復計算によって求める、すなわち `advanced_positron_eq2` を解く、際に、どのように収束解を得ようとするのか、その指定を行います。
num_extra_bands = 8	陽電子の固有状態は、基底状態1個のみを計算すれば充分です。しかし、反復計算で収束解を得るには、それ以外に、基底状態よりもエネルギーの高い状態の波動関数も計算する必要があり、その個数を指定します。なお、求める波動関数は全てブリルアンゾーン中の Γ点に属します。
delta_eigenvalue = 1.d-8 rydberg	5.の説明参照
succession = 6	反復計算において、前回と今回の物理量（7.参照）が4.で与えられた範囲内で一致し、5.で指定された回数だけ連続して、この条件を満たせば、計算は収束したとみなします。
num_max_iteration = 32000	計算は、この数繰り返すと、収束していなくても終了します。
dtim = 0.01	繰り返し計算において、次の波動関数をどれだけ大きく変化させるかの尺度であり、dtimが大きいほど収束が早くなります。しかし、あまり大きいと収束解が得られなくなり、この値が小さいほど、安定に収束解が得られますが、小さくするほど収束が遅くなります。したがって、この値は計算する系により、ユーザーが最適な値を探すことを推奨します。
epsilon_ele = 12	epsilon_eleは、ギャップのある系に対して、LDAの電子系誘電率補正を行う際に用いるtagです。=の後には、誘電率（Siの場合12）を用います。もし、誘電率の補正を行わないのであれば（たとえば、金属の計算を行う場合）、8.の行は削除します。

計算結果の出力

陽電子寿命予測計算を行うと、output000フィルと3個のcube fileが出力されます。

ログ出力ファイル　output000

このファイルの最初の部分は、Siの電子バンドの計算に関するものです。電子のバンド計算が終わり、電子の電荷密度が得られた後、陽電子の計算が行われます。

出力における“--- initial positron energy eigen values ---”からが陽電子計算に関するものです。繰り返し計算により、陽電子の波動関数が決定されます。下記の出力は、繰り返し計算のはじめにおいて、固有値が、14.6379(eV)であることを示しています。extra bandsはそれよりもエネルギーの高い固有値(14.9628460558-15.0292289699）を表します。繰り返し計算の2回目では、固有値が、 0.0021898139(eV)　となっています。

--- initial positron energy eigen values ---
 === positron eigen values ===
     14.6378982055
 -- extra_bands --
     14.9628460558     14.6842242625     14.9879179620     15.2755174303
     14.8070539395     14.6061318397     14.8086346971     15.0292289699
 === positron eigen values ===
      0.0021898139
 -- extra_bands --
      0.0892687578      0.1056325893      0.2037689630      0.2140559068
      0.3115605599      0.3359746459      0.3540270556      0.4738130045

ファイルの下の方には、次のような出力があります。

***************************************************
 positron lifetime(ps)   220.184723312044
 core rate   3.79328791767622      %

これは、陽電子の固有値を求める計算が収束し、陽電子の寿命が220psと計算されたことを示しています。core rateは、全消滅速度に対するコア電子の消滅速度の割合を示します。

Cubeファイル

計算が終了すると、電子の電荷分布、陽電子の電荷分布、電子・陽電子ペアの分布が、ファイルelectron.cube、positron.cube、ep_pair.cubeに出力されます。これらのファイルはGaussian cube 形式であり、可視化できます。図 5‑22にSi結晶における計算結果を示します。価電子は主として、結合領域に存在し、陽電子は、隙間領域に存在することが分かります。陽電子の波動関数が広がり運動エネルギーが低下した方がエネルギー的に有利であることから、一般に陽電子は隙間領域に存在する傾向があります。電子・陽電子対分布を advanced_positron1 (c)に示します。この分布が高いところで、陽電子が大きな確率で消滅することになります。

Si結晶中の価電子分布(a)、陽電子分布(b)、電子・陽電子対分布(c)

使用上の注意

陽電子寿命の計算における注意点です。

擬ポテンシャルの選択

元素によっては、セミコア状態を持つものがあります。ここで、セミコア状態とは、コア電子の内、軌道が空間的に広がり、コア電子の分布と価電子の分布との重なりが無視できない場合のことです。この場合、セミコア電子を価電子として取り扱い、擬ポテンシャルを作成することが望まれます。公開されている擬ポテンシャルには、元素によっては、このようにして作成されたものがあり、その際はこれを利用することを推奨します。もし、そのような擬ポテンシャルが作成されていなければ、CIAOを用いて作成することもできます。

カットオフエネルギーの選択

Si結晶のバンド計算において、電子波動関数、電荷密度、陽電子波動関数に対するカットオフエネルギーを以下のように指定します。

accuracy{
  cutoff_wf =  50.00  rydberg  ! cke_wf
  cutoff_cd =  200.00  rydberg  ! cke_cd
  cutoff_pwf = 50.00 rydberg
}

これらの値を変化させて、計算された寿命が十分収束していることを確かめる必要があります。

出力

陽電子の波動関数は、電子の波動関数と同様、繰り返し計算によって求められます。各繰り返しにおいてoutput000ファイルには、次のような出力があります（ここでは、繰り返しの最後における出力を表示しています）。

=== positron eigen values ===
     -0.5674635596
 -- extra_bands --
     -0.0490686179     -0.0460091253     -0.0446118499     -0.0275856742
     -0.0102856694      0.0069403602      0.0274419414      0.2284487012
lifetime:    220.180365487100        220.179503204077

ここで、計算終了近くで、陽電子の固有値（positron eigen value）が十分収束していることを確かめます。サンプルでの出力（output000)では、

-0.5674635596
-0.5674635638

などとなっており、十分収束していることが分かります。また、上記の出力で220.180365487100, 220.179503204077といった数値がありますが、これは、繰り返し計算において、前回計算された、寿命と今回計算された寿命を示しており、繰り返し計算により、寿命の値が収束に近づいている事を示唆しています。

通常の電子のバンド計算が十分収束しており、かつ使用上の注意の事項を確認できれば、計算は十分考慮されたものであると考えてよいでしょう。

not latex

参考文献

ワニエ関数

機能の概要

電子状態は結晶全体に広がったブロッホ波と呼ばれる波として表されます。ブロッホ波をk空間に関してフーリエ変換することにより得られる局在した関数をワニエ関数と呼びます。ワニエ関数の中心位置は電子分布の平均位置を表すため、その和から結晶の分極を容易に知ることができます。また、ワニエ関数の2乗は電子の分布を表すため、化学結合に関する知見が得られます。一般にワニエ関数は一意に定まりません。ワニエ関数の広がりが最小になるように変換することにより得られる、一意に定まる関数を最大局在ワニエ関数といいます。

最大局在ワニエ関数はワニエ関数の広がりを表す汎関数（広がり汎関数）を最小にするようにブロッホ波をユニタリー変換して求めます。広がり汎関数はユニタリー変換行列の関数で、広がり汎関数をユニタリー変換行列に関して微分して得られる行列はワニエ関数の広がりを狭めるユニタリー変換の方向になっています。この方向にわずかにユニタリー変換していくことで、ワニエ関数の広がりを最小にすることができます。

計算例：Siの最大局在ワニエ関数

最大局在ワニエ関数を計算するには、PostprocessingでWannier 関数を計算することを指定します。

Postprocessing{
    wannier{
        sw_wannier = ON
        eps_grad = 1.d-3
        dt = 1.d-4
        max_iteration = 1000
        filetype = cube
    }
}

汎関数の勾配の大きさがeps_grad 以下になったら計算は終了します。dt は最急降下法の仮想的な時間刻みです。繰り返しがmax_iteration を超えたら計算は停止します。ワニエ関数の出力ファイルはGauusian cube形式に指定します。ワニエ関数のファイルの拡張子がcube になるように、以下のようにfile_names.data にファイルポインタF_WANNIER を記述します。

&fnames
F_INP = './nfinput.data'
F_POT(1) = './Si_ldapw91_nc_01.pp'
F_WANNIER = './nfwannier.cube'
/

この機能はΓ点のみの計算でしか使用できません。

ksampling{
    method = gamma
}

また、並列計算には対応していません。並列計算で収束した結果を得ている場合は、非並列の継続処理としてワニエ関数の計算をしてください。

Si のワニエ関数の計算例題は、 samples/wannier/Si8 です。計算を実行すると、ワニエ関数の出力として、nfwannier.00001.cube といったファイルが16 個出力されます。それらの一つを可視化した結果を advanced_Wan_for_Si_and_GaAs に示します。最大ワニエ関数はSi-Si 結合間に局在していることがわかります。これは、Si 結晶の結合が共有結合であることを示しています。

Si 結晶(a)とGaAs 結晶(b)の最大局在ワニエ関数

期待しない局所極小に収束してしまうことがありますが、sw_random_wannier をON にしてランダムな状態から計算を始めることでこの問題が解決できることがあります。

Postprocessing{
    wannier{
        ...
        sw_randomize = ON
        ...
    }
}

また、収束が不十分な状態で計算が終了しているようでしたら、sw_continue をON にして計算を継続してください。

Postprocessing{
    wannier{
        ...
        sw_continue = ON
        ...
    }
}

Wannier90を用いたワニエ関数解析

機能の概要

Wannier90プログラム(http://www.wannier.org/) と連携することによってワニエ関数を出力することも可能です。Wannier90プログラムと連携することによって、PHASEに組み込まれたワニエ関数解析機能では実施することのできない、以下の計算を行うことが可能です。

Γ点だけでなく、一般のk点サンプリングでワニエ関数を出力することが可能
ワニエ補間により、バンド構造を得ることが可能
ワニエ補間により、フェルミ面や逆空間における等エネルギー面の描画が可能

ここでは、PHASEをWannier90と連携させて上記のような計算を行う方法を説明します。本機能を利用する場合は、上述のウェブサイトからWannier90プログラムを入手し、コンパイル作業を行っておいてください。

計算方法

Wannier90と連携して解析を行うには、以下のような処理を行います。Wannier90の入力ファイル作成方法は、Wannier90のユーザーマニュアルを参照してください。

Wannier90の入力データを作成します。ファイル名はseedname.winとします。seednameは任意の文字列で、系の名前などを指定します。
Wannier90を、以下のように“プリプロセスモード”で実行します。

% wannier90.x -pp seedname

この作業によって、Wannier90が必要とするデータの情報を記録したファイルが作成されます。このファイルの情報をもとにPHASEの入力ファイルを作成します。

PHASEの入力パラメータファイルを作成します。格子定数やk点サンプリングをWannier90の入力に合わせて編集します。その具体例は例題において紹介します。さらに、postprocessingブロックに以下を追記します。

postprocessing{
  wannier{
    seedname = "seedname"
    sw_wannier90 = on
    nb_wan90 = wannier90で利用するバンド数
  }
}

変数seednameに上述のseednameを二重引用符で指定します。sw_wannier90をonとすることによってWannier90が必要とするファイルを出力することができます。nb_wan90には、Wannier90で利用するバンド数を指定します。PHASEのバンド数以下の数値にする必要があります。

PHASEを通常通り実行します。その結果、Wannier90が必要とする行列要素データや固有値データの記録されたファイルが得られます。ただし、Wannier90が利用するデータの出力は非並列計算で行う必要があるので、並列計算をしたい場合はひとまずsw_wannier90をoffとして計算を収束させたあとにsw_wannier90をonとし、非並列の継続計算でWannier90用の出力を行います。
Wannier90を実行します。

% wannier90.x seedname

その結果、Wannier90のログ出力やワニエ関数データのほか、Wannier基底のタイトバインディングハミルトニアンデータなどが得られます。seedname.winの設定によっては、そのほかバンド構造やフェルミ面を構築するために必要なデータ、1次元伝導解析の結果なども得られます。

ワニエ補間

wan_interpプログラムは、ワニエ補間を行ってバンド構造やフェルミ面の計算を行うプログラムです。そのソースコードは、src_wan_interp以下にあります。コンパイルするにはこのディレクトリーへ移り、コンパイル用のシェルスクリプト、make.sh中のFortran90コンパイラーとLAPACKライブラリーの設定を行ったあとmake.shを実行してください。なお、ワニエ補間によるバンド構造やフェルミ面の計算は、Wannier90本体によって実施することも可能です。その方法は、Wannier90のユーザーマニュアルを参照してください。Wannier90によって行う場合は、band.plによるバンド構造図の作成やPHASE-Viewerによるフェルミ面の描画などは行えません。

wan_interpプログラムは、以下のファイルを必要とします。

seedname_hr.datファイル：Wannier90プログラムによって出力される、ワニエ基底のタイトバインディングハミルトニアンデータです。
seedname.nnkpファイル：Wannier90が必要とするデータを作成するために必要な情報が記録されたファイルです。Wannier90を“プリプロセスモード”で実行すると得られるファイルです。
kpoint.dataファイル：計算したいk点の情報が記録されたファイル。データフォーマットは、PHASEの標準のkpoint.dataファイルのフォーマットです。その作成方法は、バンド構造図の場合はバンド計算の項を、フェルミ面の場合はPHASE-Viewerユーザーマニュアルの逆空間ビューアーの項を参照してください。
nfefermi.dataファイル：フェルミエネルギーの値が記録されたファイルです。PHASEから出力されますので、wan_interpプログラムを実行するディレクトリーにコピーしておいてください。
fs.dataファイル（フェルミ面の場合）：PHASE-Viewerがフェルミ面を描画する際に利用する補助データが記録されたファイルです。PHASE-Viewerの手続きによってフェルミ面用のk点データファイルを作成した場合自動的に作成されます。

以上のファイルを揃えたら、つぎの要領で実行します。

% wan_interp seedname

計算例題

具体例として、GaAs結晶のケースを説明します。この例題の入力ファイルは、 samples/wan90/GaAs 以下にあります。

ワニエ関数の出力

まず、Wannier90を以下のように実行します。

% wannier90.x –pp gaa

この操作によって、必要な情報がgaas.nnkpというファイルに出力されます。このファイルの記述をもとにPHASEの入力ファイルを作成します。

まず、gaas.nnkpファイルには次のような書式で単位胞の情報がÅ単位で記述されます。

begin real_lattice
0.0000000 2.8265001 2.8265001
2.8265001 0.0000000 2.8265001
2.8265001 2.8265001 0.0000000
end real_lattice

この指定に従い、PHASEの入力パラメーターファイルの単位胞の指定は以下のようになっています。

structure{
  unit_cell{
    #units angstrom
    a=5.653, b=5.653, c=5.653
    alpha=90, beta=90, gamma=90
  }
  symmetry{
    tspace{
      lattice_system = facecentered
    }
  }
}

この例では、ブラベー格子を利用して単位胞の指定を行っています。GaAsは面心立方格子なので、lattice_systemパラメーターにfacecenteredが指定されています。

また、gaas.nnkpにおいてはk点サンプリングの方法が次のように指定されています。

begin kpoints
64
0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.25000000
0.00000000 0.00000000 0.50000000
0.00000000 0.00000000 0.75000000
0.00000000 0.25000000 0.00000000
0.00000000 0.25000000 0.25000000
0.00000000 0.25000000 0.50000000
…
…
end kpoints

この通りに指定するため、PHASEの入力ファイルにおいては以下のようにk点サンプリングを“直接指定モード”を利用して指定しています。

accuracy{
  ksampling{
    method = directin
    kpoints{
      #tag kx ky kz denom weight
          0 0 0 4 1
          0 0 1 4 1
          0 0 2 4 1
          0 0 3 4 1
          0 1 0 4 1
          0 1 1 4 1
          0 1 2 4 1
            …
            …
    }
  }
}

“method=directin”によってサンプルするk点を直接指定することが可能です。kpointsテーブルにおいて、kx, ky, kz, denom,weightによってk点を指定します。k点の座標は(kx/denom, ky/denom, kz/denom)、その重みがweightです。

Postprocessingブロックでは、Wannier90用の出力を行う設定が施されています。

postprocessing{
  wannier{
    seedname = "gaas"
    sw_wannier90 = ON
    nb_wan90 = 8
  }
}

Wannier90用の出力は非並列計算で行う必要があるので、並列計算をしたい場合はひとまずsw_wannier90をoffとして計算を収束させたあとにsw_wannier90をonとし、継続計算で非並列の計算を行い、Wannier90用の出力を行います。

PHASEを実行し、Wannier90用のデータが得られたら、Wannier90を-ppをつけずに行します。

% wannier90.x gaas

Wannier90による計算が終了すると、gaas_00001.xsf, gaas_00002.xsf,… などのファイルが作成されますが、これはXCrysDenプログラム（http://www.xcrysden.org/\ ）を利用して可視化することのできるワニエ関数データです。XCrysDenプログラムを利用すればワニエ関数を可視化することができますが、PHASE-Viewerによって可視化を行う場合はGaussian cube形式に変換する必要があります。この作業は、PHASEパッケージに含まれるconv.pyというPythonスクリプトを利用することによって実現することができます。advanced_Wan90_GaAs に、本例題によって得られるワニエ関数をXCrysDenで可視化した図を示します。

GaAsのワニエ関数

ワニエ補間

Wannier90には、計算したワニエ関数を基底にタイトバインディングハミルトニアンを構築する機能が備わっており、この機能を利用することによって通常の方法よりも少ない計算量でバンド構造やフェルミ面、等エネルギー面などの描画などを行うことが可能です。この作業は、wan_interpプログラムを利用して行います。この例では、GaAsの伝導体下端付近の等エネルギー面を描画する方法を紹介します。

等エネルギー面の描画に必要なデータを計算するためには、PHASE-Viewerの“逆空間ビューアー”を利用します。その方法はPHASE-Viewerユーザーマニュアルを参照してください。“逆空間ビューアー”によってkpoint.dataファイルが作成されるので、このファイルと、ワニエ関数計算によって得られたgaas.nnkpファイル、gaas_hr.datファイルのほか、PHASEによって得られたnfefermi.dataファイルを一つのディレクトリーにまとめ、そのディレクトリーにおいてwan_interpプログラムを実行します。

% wan_interp gaas

すると、ワニエ補間によってkpoint.dataファイルに記録されたk点の固有値データが計算され、nfenergy.dataというファイルに記録されます。得られたファイルを、“逆空間ビューアー”によってkpoint.dataファイルを作成したディレクトリーにコピーします。以上の作業で等エネルギー面を描画するための準備は完了です。得られた伝導体下端付近の等エネルギー面をPHASE-Viewerで可視化した図 advanced_Wan90_GaAs_isosurface に示します。

GaAs結晶の伝導体下端の等エネルギー面

使用上の注意

ワニエ関数出力は、非並列計算で行う必要があります。並列計算を行う場合はひとまずsw_wannier90をoffとして計算を収束させたあとにsw_wannier90をonとし、非並列の継続計算でWannier90用の出力を行ってください。
本機能で利用できる擬ポテンシャルは、ノルム保存型のみです。

BoltzTraPを利用した解析(バージョン2019.01以降)

BoltzTraP [Madsen06] , [Madsen18] とは、Boltzman理論を利用し、バンド構造から輸送係数などを計算するプログラムです。必要な情報は結晶の対称性とバンドエネルギーです。BoltzTraPを利用すると、温度および化学ポテンシャルの関数として以下のような量を求めることができます。

ゼーベック係数
伝導度（緩和時間で割った値が出力される）
ホール係数
熱伝導率（電子からの寄与分のみ；緩和時間で割った値が出力される）
電子比熱
磁化率

ここでは、PHASE/0が計算するバンド構造からBoltzTraPを利用する方法を説明します。BoltzTraPにはその使い方が大きく異なるバージョン1.とバージョン2.がありますが、PHASE/0はどちらのバージョンにも対応しています。

入力ファイル

BoltzTraPが必要とするファイルを出力するには、入力パラメーターファイルのPostprocessingブロックにおいてboltztrapブロックを作成し、設定を施します。

Postprocessing{
  boltztrap{
    sw_boltztrap = on
    version = 2
    prefix = CoSb3
    header = "CoSb3 calculated by PHASE/0"
  }
}

boltztrapで定義できるタグは下記の通り。

タグ名	説明
sw_boltztrap	Boltz TraPで利用できるデータファイルを出力するかどうかを指定するスイッチ。 onとするとデータファイルを出力する。
version	BoltzTraPのバージョンを指定する整数。1もしくは2を指定する。デフォルト値は1.
prefix	BoltzTraP用データファイルの接頭辞を指定する文字列。デフォルト値は実行ディレクトリー名。
header	BoltzTraP用データファイルのヘッダーを指定する文字列。デフォルト値はprefixと同じ値。

PHASE/0の実行

PHASE/0の実行は通常通り行います。計算機能の制限などは特にありません。構造最適化、格子最適化が収束するとBoltzTraPで必要なデータファイルが出力されます。SCF計算でも固定電荷計算でもBoltzTraPによる計算は可能ですが、BoltzTraPの解析はk点数に対する収束性が悪い場合があり、k点数を変えながら収束を確かめる、という手続きを踏む必要があります。そのため、固定電荷計算を活用することを推奨します。

BoltzTraPの実行

ここでは、BoltzTraPによる解析を実施する方法を簡単に説明します。詳細はBoltzTraPのマニュアル、チュートリアルをご覧ください。なお、ここでの説明はBoltzTraPのインストールは終了していることを前提としています。

バージョン1.

バージョン1.を実行するために必要なファイルはprefix.intrans, prefix.struct, prefix.energyの3つです。この3つのファイルはいずれもPHASE/0から出力されますが、prefix.intransファイルはBoltzTraPの動作を制御する“入力ファイル”なので実行したい解析に応じてユーザーが編集する必要があります。prefix.intransの内容は、典型的には下記のようになっています。なお、行頭の数値は説明用なので、実際には記述しないでください。

1:  GENE # Format of DOS
2:  0 0 0 0.0 # iskip idebug setgap shiftgap
3:  0.17267 0.0005 0.4 58 # Fermilevel (Ry), energygrid, …
4:  CALC # CALC (calculate expansion coeff), NOCALC
5:  5 # lpfac, number of latt-points per k-point
6:  BOLTZ # run mode (only BOLTZ is supported)
7:  .15 # (efcut) energy range of chemical potential
8:  800. 50. # Tmax, temperature grid
9:  -1. # energyrange of bands
10: TETRA
11: 0 0 0 0 0 # For scattering model undocumented
12: 2 # number of
13: 1E20 -1E20 # fixed doping levels in cm-3

この入力の#以降はコメントとなっています。また、ファイルは固定形式なので、記述された数値や文字列を削除してしまうとエラーになります。

重要と思われる（変更する機会が比較的多い）設定値について説明します。2行目後半2つの数値は、バンドギャップのシフトを行うかどうかを指定します。3つめの数値に1を指定するとバンドギャップの補正が行われます。4つめの数値にシフト量をRydberg単位で指定します。7行目は、考慮する化学ポテンシャルの幅をRydberg単位で指定します。8行目は、考慮する温度を指定します。まず最高温度を指定し、次いで温度の刻み幅を指定します。0Kから最高温度まで、温度の刻み幅で分割した温度について計算がなされます。12行目、13行目では考慮するドーピングレベルを指定します。まず12行目でドーピングレベルの数を指定し、さらに13行目でドーピング濃度をcm^-3単位でドーピングレベル数分指定します。なお、ドーピングを考慮するといっても特別な計算が別途行われるわけではありません。ドーピングレベルに応じた化学ポテンシャルによってデータを並べ替え、結果を出力してくれるという計算機能です。

prefix.intransファイルを望みのものに変更したら、以下の要領でBoltzTraPを実行します。PATH_TO_BOLTZTRAPはBoltzTraPをインストールしたディレクトリーに読み替えてください。

%PATH_TO_BOLTZTRAP/src/x_trans BoltzTraP
 (スピン軌道を考慮していない場合)
%PATH_TO_BOLTZTRAP/src/x_trans BoltzTraP -so
 (スピン軌道を考慮している場合)

なお、BoltzTraPの計算結果はk点数をある程度多くしないと収束しません。結果が収束しているかどうか、k点数を変えた計算を行い確認することが推奨されます。BoltzTraPのマニュアルによると、目安として$\frac{16 \times 10^{6}}{V}\ (V\mathrm{は単位胞の体積})$ 程度のk点を取ることが推奨されています。

BoltzTraPによる計算が終了すると、下記のようないろいろなファイルが結果として出力されます。

ファイル名	説明
prefix.outtrace	ログファイル
prefix.trace	各種テンソルのtraceが出力される
prefix.condtens	伝導テンソルの全要素が出力される
prefix.halltens	ホール係数テンソルの全要素が出力される
prefix.trace_fixdoping, prefix.condtens_fixdoping, prefix.halltens_fixdoping	ドーピングを考慮している場合に、prefix.xx の内容を濃度に合わせて加工したデータが出力される

これらのファイルの形式については、BoltzTraPのマニュアルを参照してください。

バージョン2.

BoltzTraPバージョン2はPythonのモジュール群であり、Pythonスクリプトの中で利用できるように設計されていますが、btp2というフロントエンドとなるコマンドも用意されています。btp2コマンドに様々な副コマンドやオプションを渡すことによってBoltzTraPバージョン2.による計算を行うことができます。ここでは、このコマンドの使い方について説明します。

btp2を実行するために必要なファイルはprefix.energyとprefix.structureファイルです。これらのファイルが置かれたディレクトリーをprefixとすると、prefixの親ディレクトリーにおいて以下の要領で“フーリエ補間”を行います。

% btp2 interpolate -m 5 prefix

interpolateは“フーリエ補間を行う”ことをbtp2に伝える副コマンド、-m 5はk点1点あたり5点の補間を行う、というオプションです。これによってinterpolation.bt2というファイルが作成されます。この補間の結果が記録されたファイルを入力として様々な解析が行えるようになっています。たとえば、ある温度範囲で輸送係数の計算を行いたい場合以下のようなコマンドを利用します。

% btp2 integrate interpolation.bt2 300:500:1

integrateは“輸送係数の計算を行う”ことをbtp2に伝える副コマンドです。そのあとに作成した補間の結果ファイルinterpolation.bt2を渡し、さらに温度範囲を最低温度:最高温度:温度の刻み幅という形式指定しています。得られるファイルなどについては、BoltzTraPのバージョン1.と違いはありません。

計算例

体心立方CoSb₃結晶のゼーベック係数の計算例を紹介します。CoSb₃の原子配置は advanced_BT_CoSb3 に示すとおり。

CoSb₃の結晶構造（体心立方構造）

なお、この例題の入力ファイルは samples/BoltzTraP/CoSb3 以下に置かれています。

まずSCFディレクトリーにある入力ファイルを用いて通常のSCF計算を行います。このSCF計算においてはk点サンプリングは4×4×4のMonkhorst-Packメッシュを採用しています。つづいて、fc_5x5x5, fc_10x10x10, fc_20x20x20, fc_30x30x30ディレクトリーにおいてk点メッシュを変更した固定電荷の計算を行います。固定電荷計算の入力ファイルには、以下のようにBoltTraPが読み込める形式の固有値データなどが出力されるような設定が施されています。

postprocessing{
  boltztrap{
    sw_boltztrap = on
    version = 1
    header = "CoSb3 ksamp 10x10x10 bands 96"
  }
}

ekcalもしくはphaseによる固定電荷計算を行うと、以下の3つのファイルが得られるはずです（fc_5x5x5ディレクトリーの場合）

fc_5x5x5.energy fc_5x5x5.struct, fc_5x5x5.intrans

*.energy, *.structファイルは固有エネルギーと座標データが記録されているファイルですが、通常編集する必要はありません。*.intransはBoltTraPの入力ファイルなので目的によっては編集する必要があるかもしれませんが、編集せずともBoltzTraPによる基本的な計算は行うことができるようになっています。

BoltzTraPは以下の要領で実行できます (BoltzTraPの実行ファイルであるx_transにパスが通っていると仮定)。

$ x_trans BoltzTraP

これによって、ディレクトリー名を接頭辞としてもつ多くのファイルが出力されます。ゼーベック係数などの輸送係数の計算結果が出力されるのが*.traceファイルです。たとえば、fc_5x5x5の場合その内容は下記のようになります。

# Ef[Ry] T [K] N DOS(Ef) S s/t R_H kappa0 ...
0.02274  50.0000 34.92246648 0.45246722E+03 -0.16967269E-04 0.94428193E+20 0.15114046E-08 ...
0.02274 100.0000 34.89576852 0.46772674E+03 -0.82371993E-05 0.95266319E+20 0.11470114E-08 ...
0.02274 150.0000 34.86708620 0.47412414E+03 -0.53033332E-05 0.94512268E+20 0.10363499E-08 ...
0.02274 200.0000 34.83810551 0.47509694E+03 -0.38465950E-05 0.94664477E+20 0.93685045E-09 ...
0.02274 250.0000 34.80983692 0.47363201E+03 -0.18452224E-05 0.95321416E+20 0.85951902E-09 ...
...
...

結果は、同じ化学ポテンシャルのデータがセットで出力されます。ゼーベック係数は5カラム目のデータです。300Kにおけるゼーベック係数を、k点サンプリングを変化させながら計算した結果を advanced_BT_CoSb3_S に図示します。この結果から、5×5×5のケースはそのほかの場合と比較し明らかに異なる傾向を示しているため、不十分なk点サンプリングであるといえます。10×10×10, 20×20×20, 30×30×30はおおよそ一致していますが、10×10×10のケースではたとえば化学ポテンシャルが0.25から0.3 Rydbergの領域で異なっているように見えます。20×20×20, 30×30×30はほぼ見分けがつきません。したがって、この例では20×20×20のサンプリングでほぼ収束した結果が得られていると考えられます。

CoSb₃の300Kにおけるゼーベック係数

not latex

参考文献

Energy Density Analysis(バージョン2022.01以降)

概要

Energy Density Analysis法（以降EDA法）とは、平面波基底に基づく第一原理擬ポテンシャル法によって得られるエネルギーを原子ごとに分割する計算手法です [Imamura10] 。

使い方

入力

Energy Density Analysis法を用いる場合、入力パラメーターファイルに以下のような記述を施します。

postprocessing{
  eda{
    sw_eda = on
  }
}

この設定によってSCF計算が収束したあとEDAによるエネルギー分割が行われます。

出力

EDAの結果は、以下のような形式でログファイル(output000ファイル)に記録されます。

atom        density        kinetic        hartree            exc          local       nonlocal          ewald            epc          total
   1  Si        5.1757754      2.3466162      1.7801385     -1.9146166     -2.9861858      0.1762697     -3.5142150      0.0000000     -4.1119930
   2  Si        6.4835609      4.5278124      3.6036561     -2.8902079     -6.0099149      0.1762698     -3.5142143      0.0000000     -4.1065988
   3  Si        5.0811369      2.2787530      1.7134662     -1.8694366     -2.9032419      0.1785505     -3.5142143      0.0000000     -4.1161231
   4  O         5.3960516      6.5742616      5.4242108     -3.2629480    -17.7140909      3.1079633    -10.1029037      0.0000000    -15.9735069
   5  O         5.3919542      6.5371767      5.4250548     -3.2620495    -17.7156335      3.0784449    -10.1029039      0.0000000    -16.0399105
   6  O         4.8630732      5.5310349      4.5725356     -2.8276341    -16.2893323      3.0828468    -10.1029028      0.0000000    -16.0334519
   7  O         4.8607947      5.8746188      4.8724768     -2.9345568    -16.8295440      3.0784455    -10.1029028      0.0000000    -16.0414625
   8  O         5.4033197      6.5416558      5.4209469     -3.2641048    -17.7094035      3.0828468    -10.1029037      0.0000000    -16.0309626
   9  O         5.3443334      6.1435650      5.0619129     -3.1196637    -17.0611132      3.1079635    -10.1029039      0.0000000    -15.9702393
  total        48.0000000     46.3554944     37.8743987    -25.3452180   -115.2184601     19.0696009    -71.1600644      0.0000000   -108.4242486

エネルギー種別ごとに結果が記録されます。一行が1つの原子の結果に対応します。 1カラム目が原子のID, 2カラム目が元素名, 3カラム目が電子密度, 4カラム目が電子の運動エネルギー, 5カラム目がハートリーエネルギー, 6カラム目が交換相関エネルギー, 7カラム目が局所エネルギー, 8カラム目が非局所エネルギー, 9カラム目が原子核間のエネルギー, 10カラム目がpartial core correctionエネルギー, 11カラム目がすべての貢献分の和に対応します。

計算例

EDAを用いた計算例として、Si(001)面にC2H2分子が吸着した系の最適化計算例を紹介します。吸着様式としてはdimerized model (DM, EDA_example_dm_dcm (a)) およびdimer-cleaved model (DCM, EDA_example_dm_dcm (b))を採用しましたた。サンプルの入力ファイルは samples/EDA/Si001_C2H2 の DCM および DM 以下に配置されています。

採用した計算条件は下記の通り。

計算条件

平面波カットオフ [Ry]	30
電荷密度カットオフ [Ry]	300
k 点サンプリング	monk (2×2×2)
交換相関相互作用	GGAPBE
SCF計算収束判定条件	1e-10 hartree
構造最適化収束判定条件	2e-4 hartree/bohr
擬ポテンシャル	Si_ggapbe_paw_nc_01m.pp C_ggapbe_paw_us_01.pp H_ggapbe_paw_nc_01m.pp

Si(001)面にC2H2分子が吸着した系。(a) dimerized model, (b) dimer-cleaved model

構造最適化の結果DMのエネルギーがDCMと比較して1.03 eVほど低いという結果になりました。これは文献値 [Imamura10] の1.04 eVと近い結果です。

EDA法を適用した結果得られた各原子のエネルギーは次表に報告する通り。この結果から、最表面のSi原子の違いが二つの構造のエネルギー差に最も寄与しているなどの情報を読み取ることができます。

DMおよびCDMの原子ごとのエネルギー。原子の並びはz座標の降順。

原子	DM (Hatree)	DCM (Hartree)	DMとDCMの差 (Hartree)
H1	-0.5685776	-0.5637801	-0.0047975
H2	-0.5685893	-0.5637711	-0.0048182
C1	-5.6106157	-5.6098881	-0.0007276
C2	-5.6106366	-5.6099246	-0.000712
Si1	-3.9844977	-3.9695297	-0.014968
Si2	-3.9844253	-3.9695649	-0.0148604
Si3	-3.9372327	-3.9407006	0.0034679
Si4	-3.9371766	-3.9407006	0.003524
Si5	-3.9401731	-3.9412468	0.0010737
Si6	-3.9411641	-3.9368933	-0.0042708
Si7	-3.9425975	-3.9430295	0.000432
Si8	-3.9427341	-3.9419992	-0.0007349
Si9	-4.0002908	-4.0002038	-0.0000870
Si10	-4.0002711	-4.0002176	-0.0000535
H1	-0.5335011	-0.5334118	-0.0000893
H2	-0.5334932	-0.5333809	-0.0001123
H3	-0.5335100	-0.5334081	-0.0001019
H4	-0.5335023	-0.5333772	-0.0001251

not latex

参考文献

Bernasconi95

M. Bernasconi, G.L.Chiarotti, P.Focher,S.Scandolo,E.Tosatti,M.Parrinello Journal of Physics and Chemistry of Solids, 56 501-505 (1995).

Broberg18

D. Broberg, B. Medasani, N. Zimmermann, A. Canning, M. Haranczyk, M. Asta, and G. Hautier,, Comput. Phys. Commun. 226 (2018) 165.

Cabaret10

Cabaret, A. Bordage, A. Juhin, M. Arfaoui, and E. Gaudry, Phys. Chem. Chem. Phys. 12 (2010) 5619.

Freysoldt09

C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. 102 (2009) 016402.

Gao08

-P. Gao, C. J. Pickard, M. C. Payne, J. Zhu, and J. Yuan, Phys. Rev. B 77 (2008) 115122.

Gougoussis09

Gougoussis, M. Calandra, A. P. Seitsonen, and F. Mauri, Phys. Rev. B 80 (2009) 075102.

Imamura10

Yutaka Imamura,Asuka Takahashi, Takeshi Okada, Takahisa Ohno, and Hiromi Nakai "Extension of energy density analysis to periodic-boundary-condition calculations with plane-wave basis functions" Physical Review Letters vol. 89 (2002) pp. 117602 1-4.

Kumagai14

Y. Kumagai and F. Oba, Phys. Rev. B 89 (2014) 195205.

Landemark92

“Core-level spectroscopy of the clean Si(001) surface: Charge transfer within asymmetric dimers of the 2 × 1 and c(4 × 2) reconstructions”, E. Landemark, C.J. Karlsson, Y.-C. Chao, and R.I.G. Uhrberg, Phys. Rev. Lett. 69, 1588 (1992).

Lenthe96

https://www.scm.com/wp-content/uploads/vlenthe.pdf

Madsen06

Georg K.H.Madsen and David J.Singh, “BoltzTraP. A code for calculating band-structure dependent quantities”, Computer Physics Communications 175, 67, (2006)URL: https://www.imc.tuwien.ac.at/forschungsbereich_theoretische_chemie/forschungsgruppen/prof_dr_gkh_madsen_theoretical_materials_chemistry/boltztrap/

Madsen18

Georg K.H.Madsen, JesúsCarrete, and Matthieu J.Verstraete, “BoltzTraP2, a program for interpolating band structures and calculating semi-classical transport coefficients”, Computer Physics Communications 231, 140, (2018) URL: https://www.imc.tuwien.ac.at/index.php?id=21094

Mahmood03

A. Mahmood, R. Machorro, S. Muhl, J. Heiras, F. F. Castillon, M. H. Farıas, and E. Andrade, DIAM RELAT MATER 12 (2003) 1315.

Nakamoto08

A. Nakamoto, M. Saito, T. Yamasaki, M. Okamoto, T. Hamada, and T. Ohno, ``Two-Component Density Functional Calculations on Positron Lifetimes for Band-Gap Crystals'', Jpn. J. Appl. Phys. 47 (2008) p. 2213.

Ozaki17

T. Ozaki and C. Lee, Phys. Rev. B 118 (2017) 026401.

Pehlke93

“Evidence for site-sensitive screening of core holes at the Si and Ge (001) surface”, E. Pehlke and M. Scheffler, Phys. Rev. Lett. 71, 2338 (1993).

Puska91

1. Puska, ``Ab-initio Calculations of Positron Annihilation Rates in Solids'', J. Phys. Condens. Matter 3 (1991) p. 3455.

Puska95

M. J. Puska, A. P. Seitsonen, and R. M. Nieminen, ``Electron-positron Car-Parrinello Methods: Self-consistent Treatment of Charge Densities and Ionic Rel axations'', Phys. Rev. B 52 (1995) p. 10947.

Souza02

I. 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.

ディレクトリー名	説明
Pristine	欠損のない結晶の入力ファイルが納められたディレクトリー。
Ga_for_As	AsをGaで置換した欠陥構造の入力ファイルが納められたディレクトリー。
As_for_Ga	GaをAsで置換した欠陥構造の入力ファイルが納められたディレクトリー。
As_interstitial_As4	Asが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。
As_interstitial_Ga4	Gaが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。
Ga_interstitial_As4	Asが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。
Ga_interstitial_Ga4	Gaが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。

Files

analysis.rst

Latest commit

History

analysis.rst

File metadata and controls

解析機能

部分電荷密度

ベーダ―解析向け電荷密度ファイルの出力 (バージョン2019.02以降)

概要

ベーダ―解析について

処理内容

入力

出力

計算例

GaNの計算例

4H-SiCの計算例

ストレステンソル

機能の概要

入力パラメータ

弾性定数

ストレステンソルの補正

仕事関数

機能の概要

入力パラメータ

計算の実行方法

計算結果の出力

計算例：アルミニウムの仕事関数

原子周囲局所ポテンシャル出力機能（2021.02以降）

機能の概要

入力パラメーターファイル

出力ファイル

計算例： Al-Graphene-Al系

帯電欠損状態の評価 (バージョン2021.01以降)

概要

理論

FNV法

外部解析プログラム

extended FNV 法

使い方

外部解析プログラムの使用

extended FNV 法

生成エネルギー図の作成1

生成エネルギー図の作成2

例題

概要

例題のディレクトリー構成および計算のながれ

補正エネルギー比較

生成エネルギー図の作成 1.

生成エネルギー図の作成 2.

XPS

概要

理論

相対論的運動エネルギー

内殻電子のスピン

XPSのピークエネルギーの計算方法

使い方

相対論的運動エネルギー

内殻電子のスピン

例題

4H-SiC結晶 (C 1s)

w-AlN 結晶 (N 1s)

w-AlN 結晶 (Al 2p3/2)

fcc Pt結晶 (Pt 4f)

O 2 分子 (O 1s)

Core Level Shift (CLS)の解析

内殻準位シフト

symm j Energy (Ha) Energy (eV) nocc focc

Total number of electrons 14.00000

内殻光電子放出のスペクトル

内殻正孔を含む原子の擬ポテンシャル

計算の実行方法

計算例：Si(100)表面

ELNES / XANES解析機能

はじめに

準備：擬ポテンシャル

計算の流れ

基底状態 (A) のSCF 計算

入力

出力

O ₂ 分子 (O 1s)