第一原理電子状態計算ソフトウェアPHASE/0 2023 チュートリアル

はじめに

概要

本稿は第一原理バンド計算ソフトウェアPHASE/0のチュートリアルです。第一原理バンド計算は多大な計算リソースが必要な場合が多いですが，本チュートリアルは一般的なパソコンでも実行できるような規模の問題のみを扱い，その実行方法や結果の見方について解説します。環境としては，Windows上で構築したWSL環境 を想定していますが，一般的なGNU/LinuxやmacOSでも実行できるはずです。お手元のパソコンにおいて実際に例題を実行することを通じて，初学者がPHASE/0の基本的な使い方を自習することができることを目的としています。紹介されている例題の入力ファイルは GitHubで公開されているため 自由にダウンロードすることができるようになっています。

本チュートリアルの構成

シリコン結晶を用いた基本的な計算 では簡単に計算できる例として2原子からなるシリコン結晶の例を扱います。基本となるself-consistent field (SCF) 計算を実行し，それをもとに電子状態密度やバンド構造，誘電関数の計算などを実行します。また，振動解析や \(EV\) 曲線の計算を行う方法も紹介します。磁性を考慮した計算 では磁性の扱い方を紹介します。磁性を有効にする方法や磁性を考慮している場合の結果の解析方法などを説明します。表面の計算 では表面の計算の方法を紹介します。PHASE/0は平面波基底に基づいて計算を行うため，扱えるのは厳密にいうと周期系のみです。しかし表面に垂直な方向に十分な厚さの"真空層"を設けることによって事実上表面と変わらない系を扱うことはできます。ここではその方法について説明します。また，単純な結晶では必要ない「構造最適化」の方法を説明し，さらに原子ごとの電子状態を解析する方法についても説明します。点欠陥の生成エネルギーの計算 では点欠陥のエネルギー計算を扱います。ここでは，第一原理計算において「生成エネルギー」はどのように計算するかを説明します。最後に，付録 においてはPHASE/0を活用するためのちょっとしたコツなどを紹介していきます。

環境設定など

本チュートリアルを実行するにはPHASE/0がお使いのパソコンにインストールされている必要があります。もしインストールされていないのであれば，このサイトの情報 などを参考にWSLとPHASE/0をインストールしてください。本稿では以降ホームディレクトリーに phase0_2023.01 がインストールされていることを仮定して説明を進めます。また，原子配置の可視化に VESTA を用いるので，こちらもインストールしておくことが推奨されます。

PHASE/0の引用情報など

シリコン結晶を用いた基本的な計算

まずはもっとも簡単な例題として，2原子からなるシリコン結晶の計算を行います。最初に基本となるself-consistent field (SCF) 計算を実施します。この計算が終了すると電荷密度が得られるので，VESTAを利用して電荷密度を可視化する作業を行います。さらに，この計算によって得られた電荷密度を入力とし，状態密度・バンド構造・誘電関数の計算を実行し，結果の解析方法を学んでいただく予定です。

SCF計算

すべての計算の基本となる，SCF計算の実行方法を説明します。この計算は，あらかじめ用意したサンプルデータをもとに実行していただきます。留意する必要のある入力パラメーターファイルの設定項目を説明したあとに計算を実行し，さらに結果として得られる電荷密度の可視化の方法を説明します。例題の入力ファイルは wsl/Si2/scf 以下に配置されています。

入力データの確認

nfinp.dataファイル

nfinp.dataファイルは，PHASE/0の計算設定を行う「入力パラメーターファイル」です。この例題のnfinp.dataファイルの内容は下記の通り。

Control{
  condition = initial
}
accuracy{
  cutoff_wf = 20.0 rydberg
  cutoff_cd = 80.0 rydberg
  ksampling{
    mesh{ nx = 4, ny = 4, nz = 4 }
  }
}
structure{
  unit_cell_type = primitive
  unit_cell{
    a_vector = 0.00 5.13 5.13
    b_vector = 5.13 0.00 5.13
    c_vector = 5.13 5.13 0.00
  }
  symmetry{
    method = automatic
    tspace{
      lattice_system = facecentered
    }
  }
  atom_list{
    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
  }
}
postprocessing{
  charge{
    sw_charge_rspace = on
    filetype = cube
  }
}

PHASE/0の入力パラメーターファイルは，このように波カッコ{}で囲った“ブロック”とそこで設定される値から構成されます。以下のような設定が施されています。

  • accuracyブロックにおいて計算精度に関わる設定が施されています。cutoff_wf, cutoff_cdにおいてそれぞれ波動関数および電荷密度のカットオフエネルギー，ksamplingブロックにおいてk点サンプリングの設定が施されています。

  • structureブロックにおいて原子配置が指定されています。

    • unit_cellブロックで単位胞を指定しています。

    • symmetryにおいて対称性が指定されています。method = automaticとすることによってプログラムが自動的に対称性を検出し，それを活用します。この際tspaceブロックのlattice_systemキーワードで結晶がfacecenteredであることを伝えています。

    • atom_listのatomsテーブルにおいて原子配置そのものを指定しています。

  • postprocessingブロックのchargeブロックにおいてsw_charge_rspace = onとしています。こうすることによって実空間にマップした電荷密度をファイルに出力します。また，filetype = cubeとすることによってその形式をGaussian cube形式にしています。

file_names.dataファイル

「ファイル名を指定するためのファイル」がfile_names.dataファイルです。PHASE/0の入出力ファイルはこのファイルを通じてファイル名や場所を指定することができますが，このファイルのみ実行ディレクトリーにこの名前で存在する必要があります。この例題のfile_names.dataファイルの内容は下記の通り。

&fnames
F_POT(1) = '../../pp/Si_ggapbe_paw_nc_01m.pp'
F_CHR = './nfchr.cube'
/

F_POT(1)によって1番目の元素の擬ポテンシャルファイルを指定しています。F_CHRによって実空間の電荷密度を出力するファイルを指定しています。このファイルのデフォルト値はnfchr.dataですが，Gaussian cube形式で出力するので拡張子をcubeに変更しておいた方が便利です。

計算の実行

実際に計算を実行してみましょう。たとえば以下のようなコマンドを用いて計算を起動します。この際 タブ補完 を活用することが推奨されます。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

この例ではk点2並列で計算を実行しました。お使いのコンピューターのCPUに搭載されたコアの数に応じて並列数を決めるようにしてください。

計算結果の解析

計算の結果，次のようなファイル群が出力ファイルとして得られます。

output000

計算のログを記録するファイルです。SCF計算の各iterationにおけるエネルギーなどが記録されます。

nfefn.data

全エネルギーの計算結果が記録されるファイルです。構造最適化や分子動力学シミュレーションを実行している場合，エネルギーの履歴が記録されます。

nfdynm.data

座標データ履歴が記録されるファイルです。原子座標と原子間力が記録されます。

nfchr.cube

実空間の電荷密度データが記録されるファイルです。Gaussian cube形式で記録されます。VESTAなどのアプリケーションによって可視化することができます。

そのほか，継続計算に必要なファイルなどが得られます。また，計算機能によってはほかにも様々な出力ファイルが得られる場合があります。output000ファイルからSCF計算の収束具合を調べてみましょう。

$ grep TH output000
TOTAL ENERGY FOR  1 -TH ITER= -7.843775135159 EDEL = -0.784378D+01 : ...
TOTAL ENERGY FOR  2 -TH ITER= -7.851142435464 EDEL = -0.736730D-02 : ...
TOTAL ENERGY FOR  3 -TH ITER= -7.869822316521 EDEL = -0.186799D-01 : ...
TOTAL ENERGY FOR  4 -TH ITER= -7.873273739596 EDEL = -0.345142D-02 : ...
TOTAL ENERGY FOR  5 -TH ITER= -7.874467564388 EDEL = -0.119382D-02 : ...
...
...
TOTAL ENERGY FOR 12 -TH ITER= -7.875384283153 EDEL = -0.446631D-08 : ...
TOTAL ENERGY FOR 13 -TH ITER= -7.875384283383 EDEL = -0.230161D-09 : ...
TOTAL ENERGY FOR 14 -TH ITER= -7.875384283401 EDEL = -0.179563D-10 : ...

TOTAL ENERGY FOR xx -TH ITER=の後の数値が全エネルギーの計算結果，EDEL = のあとの数値が1ステップ前と現ステップの全エネルギーの差です。この差が閾値（今の場合デフォルト値の1e-9 hartree/atom）よりも小さいステップがある回数（今の場合デフォルト値の2回）続くとSCF計算が収束したとみなされます。構造最適化などを行う場合はこのあと原子配置を更新し次のSCF計算へ移行しますが，今の場合は一点計算なのでこれで計算終了です。

VESTAを用いて電荷密度を可視化してみましょう。VESTAを起動し，FileメニューからOpenを選びます。するとファイル選択ダイアログが得られるので，nfchr.cubeファイルを選択します。すると 図 2.1 で表示するような電荷密度の等値面が描画されます。等値面に採用する値や描画色などの設定は，Propertiesボタンをクリックすると得られる画面のIsosurfaceセクションを選ぶと得られるインターフェースから設定することができます（ 図 2.2 ）。

_images/image2.png

シリコン結晶の電荷密度の等値面

_images/image3.png

等値面設定画面

バンド計算

SCF計算で電荷密度が得られたので，これを入力として“バンド構造”を計算します。すなわち，電荷密度は正しいものとしてこれを固定し，第一ブリユアンゾーン内の対称線上の固有エネルギーの計算を行います。このような計算手法を“固定電荷計算”とよびます。例題の入力ファイルは wsl/Si2/band 以下に配置されています。

入力データの確認

nfinp.dataファイル

バンド計算用の入力ファイルは以下のようなものとなります(SCF計算と共通する部分は省略)。

Control{
  condition = fixed_charge
}
accuracy{
  ...
  ksampling{
    method = file
  }
}
...

  • controlブロックにおいてcondition = fixed_chargeとすることによって固定電荷計算を行う指示をしています。

  • accuracyのksamplingブロックにおいてmethod = fileとすることによってk点生成方法をファイルからの読み込みとしています。

file_names.dataファイル

file_names.dataファイルは入出力ファイルのファイル名を指定するためのファイルです。このファイルにおいて，以下の要領でSCF計算の電荷密度ファイルの位置を指定します。

&fnames
...
F_CHGT = '../scf/nfchgt.data'
...
/

F_CHGTという識別子を用いて，1階層上のscfディレクトリーのnfchgt.dataファイルを指定しています。

kpoint.dataファイル

kpoint.dataファイルにバンド構造を計算したい対称線上のk点の座標が記述されています。その内容は下記の通り。

141 141
0 50 50 100 1
0 49 49 100 1
0 48 48 100 1
...
...

1行目の1カラム目にk点の総数を指定します。2カラム目はこの例題では未使用のデータです。2行目以降がk点の座標で，最初の3カラムがx, y, z座標，4カラム目が規格化のファクター，5カラム目が重みとなります。実際のk点の座標はたとえば1から3カラム目の値をkxkykz 4カラム目の値をdとするとkx/d, ky/d, kz/d となります。また，バンド計算の場合重みは常に1です。

特殊k点の座標を入力とするとkpoint.dataファイルを作成してくれるPerlスクリプトがband_kpoint.plスクリプトです。その入力ファイルがbandkpt.inです。その内容は以下のようになります。

0.02
-1.0 1.0 1.0
1.0 -1.0 1.0
1.0 1.0 -1.0
0 1 1 2 # X
0 0 0 1 # {/Symbol G}
1 1 1 2 # L
...

1行目に対称線上のk点の“間隔”を指定します。続く3行が逆格子の指定ですが，これは比があっていれば問題ありません。以降で対称線上の特殊k点の値を指定します。この例ではまず(0 1/2 1/2)すなわちFCCのX点から始まり (0 0 0)すなわちΓ点に向かい，そこから(1/2 1/2 1/2)すなわちFCCのL点に向かい，... という対称線の指定となります。各特殊k点指定において#のあとに特殊k点をあらわす記号を指定することができます。この情報はk点生成には用いられませんが，後述するバンド図に反映させることができます。このようなファイルを入力とし，band_kpoint.plを以下の要領で実行するとkpoint.dataファイルを作成することができます。

$ ~/phase0_2023.01/bin/band_kpoint.pl bandkpt.in
Distance of 1 = 1
Distance of 2 = 0.866025403784439
Distance of 3 = 0.612372435695794
Distance of 4 = 0.353553390593274
division numbers = 50 43 30 17
i=0 0 -0.01 -0.01
i=1 0.0116279069767442 0.0116279069767442 0.0116279069767442
i=2 0.00416666666666667 -0.00833333333333333 0.00416666666666667
i=3 -0.00735294117647059 -0.0147058823529412 -0.00735294117647059
1 : 0 0.5 0.5
=> 0/100 50/100 50/100
...

計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

SCF計算の場合と同様，k点2並列で実行してみました。

計算結果の解析

バンド計算の結果はnfenergy.dataファイルに記録されます。このファイルから“バンド図”を作成するためのPerlスクリプトがband.plです。以下の要領で実行します。

~/phase0_2023.01/bin/band.pl nfenergy.data bandkpt.in -color -with_fermi

band.plにはいくつかオプションがありますが，そのうち-color (カラーのバンド図を作成する)と -with_fermi (フェルミエネルギーの位置をあらわす線を描画する)を有効にしました。バンド図はband_structure.epsというEPS形式の画像ファイルとして得られます。これを表示するにはevinceコマンドを利用します。

evince band_structure.eps

以上の操作によって 図 2.3 で示すようなバンド構造図が描画されます。

_images/image4.svg

シリコン結晶のバンド構造図

状態密度計算

バンド計算の場合と同様固定電荷計算を用いて状態密度の計算を行います。状態密度の計算はSCF計算収束後のポスト処理として行うことも可能ですが，SCF計算とは異なるk点サンプリングやバンド数を用いて状態密度を求めたい場合は固定電荷計算によって計算することができます。この例題の入力ファイルは wsl/Si2/dos 以下にあります。

入力データの確認

nfinp.dataファイル

状態密度計算用の入力ファイルは以下のようなものとなります(SCF計算・バンド構造計算と共通する部分は省略)。

accuracy{
  ksampling{
    method = mesh
    mesh{ nx = 8, ny = 8, nz = 8 }
  }
  smearing{
    method = tetrahedral
  }
}
postprocessing{
  dos{
    sw_dos = ON
    method = tetrahedral
  }
}

  • accuracyブロックのksamplingブロックにおいてmethodをmeshとし，さらにsmearingブロックにおいてmethodをtetrahedralに設定しています。これらの設定は四面体法による状態密度計算を行うために必要なものです。van-Hove特異点を精度よく再現するために，k点分割数もSCF計算よりも増やしています。

  • postprocessingブロックのdosブロックにおいてsw_dosをonとし，methodをtetrahedralとしています。

計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

SCF計算の場合と同様，k点2並列で実行してみました。

計算結果の解析

状態密度の計算結果はdos.dataファイルに記録されます。その内容は下記のようになります。

No.   E(hr.)  dos(hr.)   E(eV)    dos(eV)   sum
    6 -0.3404 0.00000000 -12.1937 0.0000000000 0.0000
   16 -0.3400 0.00000000 -12.1837 0.0000000000 0.0000
...
...
12176  0.1069 0.00891554  -0.0237 0.0003276399 8.0000
12186  0.1072 0.00297770  -0.0137 0.0001094283 8.0000
12196  0.1076 0.00029899  -0.0037 0.0000109877 8.0000
12206  0.1080 0.00000000   0.0063 0.0000000000 8.0000
12216  0.1083 0.00000000   0.0163 0.0000000000 8.0000
...
...

1行があるエネルギーにおける状態密度に対応します。1カラム目はエネルギーの識別番号，2カラム目および3カラム目がハートリー単位のエネルギーと対応する状態密度，4カラム目と5カラム目がeV単位のエネルギーと対応する状態密度，6カラム目が積算状態密度です。4カラム目はフェルミエネルギー（もしくは価電子帯の上端）が0となるようにシフトされています。4カラム目が0となる準位における積算状態密度は電子数と等しくなります。

dos.dataファイルから“状態密度図”を作成するPerlスクリプトがdos.plです。以下のように実行することができます。

~/phase0_2023.01/bin/dos.pl dos.data -color -with_fermi

dos.plにはいくつかオプションがありますが，そのうち-color (カラーのバンド図を作成する)と -with_fermi (フェルミエネルギーの位置をあらわす線を描画する)を有効にしました。状態密度図はdensity_of_states.epsというEPS形式の画像ファイルとして得られます。これを表示するにはevinceコマンドを利用します。

evince density_of_states.eps

以上の操作によって 図 2.4 で示すような状態密度図が描画されます。

_images/image5.svg

シリコン結晶の状態密度図

誘電関数の計算

固定電荷計算を用いて誘電関数の計算を行います。誘電関数の計算は状態密度の計算と同様指定のk点セットを用いて波動関数をときなおし，バンド間遷移の遷移確率を計算することによって計算します。この例題の入力は wsl/Si2/eps 以下にあります。

入力データの確認

nfinp.dataファイル

固定電荷状態密度計算用の入力に似ていますが，epsilonブロックを作成し，そこで誘電関数計算の設定を行う点が異なります。

epsilon {
  sw_epsilon = on
  photon{
    polar{ux=1.00, uy=0.00, uz=0.00}
    energy{low=0.000, high=2.000, step=0.002}
  }
  fermi_energy{
    read_efermi = off
    efermi = 0.0000
  }
  transition_moment{
    type = ks
    symmetry = on
  }
  BZ_integration {
    method = t
  }
}

  • sw_epsilon = onとすることによって誘電関数計算を有効にします。

  • photonブロックにおいて電磁波状態の情報を設定します。polorブロックで「直線偏光の分極ベクトル」を指定し，energyブロックで電磁波のエネルギー範囲を指定します。

  • transition_momentで遷移確率の計算方法を指定します。この例ではtype = ksとすることによってKageshima-Shiraishi(KS)型遷移モーメント補正を用いることを指定しています。またsymmetry = onとすることによって対称性を活用する指定を行っています。

  • BZ_integrationブロックにおいてブリユアンゾーン内でどのように積分を行うかを指定します。method = tとすると四面体法を利用するようになります。

計算の実行

誘電関数計算はepsmainプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/epsmain ne=1 nk=2

計算結果の解析

誘電関数の計算結果はeps.dataファイルに記録されます。

Dielectric Function Optical Properties
Photon Energy(eV) Real Part Imaginary Part n k abs(in 10**9 m-1) R
0.00000 13.37808 0.00000 3.65761 0.00000 0.00000 0.32558
0.05442 13.38035 0.00000 3.65792 0.00000 0.00000 0.32561
...
...

各行があるエネルギーの計算結果に対応します。各カラムは，それぞれ電磁波のエネルギー, 誘電関数(実部), 誘電関数(虚部), 屈折率(実部), 屈折率(虚部), 吸収係数, 反射スペクトルに対応します。バンド構造や状態密度計算のような作図用のスクリプトはありませんので，gnuplotを用いてグラフを作成してみます。たとえば以下のようなコマンドによって誘電関数の実部と虚部をエネルギーの関数としてプロットすることができます。

$ gnuplot
...
...
Terminal type is now 'wxt'
gnuplot> plot 'eps.data' using 1:2 with lines
gnuplot> replot 'eps.data' using 1:3 with lines axis x1y2
gnuplot> set y2tics
gnuplot> replot
gnuplot> set xrange [0:20]
gnuplot> replot

図 2.5 のような表示が得られるはずです。

_images/image6.png

上述の手続きによって得られる誘電関数の実部と虚部のプロット。

振動解析

PHASE/0には格子振動の基準モードを計算する機能が備わっています。振動解析は安定な結晶から原子位置をわずかにずらし，力計算を行うことによって行います。この際，ずらす対象の原子と方向は検出した対称性に応じてなるべく少なくなるように自動的に決まります。例題の入力ファイルは wsl/Si2/phonon 以下に配置されています。

入力データの確認

nfinp.dataファイル

SCF計算 の入力にPhononブロックを追加し，振動解析の設定を施します。また，通常のSCF計算と違い質量が重要な意味を持つので，元素の質量を正しく設定します。

structure{
  ...
  element_list{
    #units atomic_mass
    #tag element atomicnumber mass
    Si 14 28.0855
  }
}
Phonon{
  sw_phonon = on
  sw_calc_force = on
  sw_vibrational_modes = on
}

  • element_listの下に#units atomic_massを記述することによってブロック内の質量の単位を原子質量単位に変更しています（デフォルトの質量の単位は原子単位；すなわち電子の質量を1とする単位）またmassカラムを追加することによって質量を指定できるようにしています。

  • Phononブロックにおいてsw_phonon, sw_calc_force, sw_vibrational_modesをonとします。このように設定すると振動解析は自動的に行われます。

計算の実行

振動解析はphaseプログラムで行います。たとえば以下のように実行します。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

SCF計算 と違い，対称性に応じて原子を変位させながら原子間力の計算を行うので複数回のSCF計算を行います。この例では2回のSCF計算が実行されます。原子を変位させると対称性がそこなわれるためSCF計算1回あたりの計算時間も長くなる傾向になります。

計算結果の解析

振動解析の計算結果はmode.dataファイルに記録されます。その内容は下記の通り。

--- primitive lattice vectors ---
0.0000000000 5.1300000000 5.1300000000
5.1300000000 0.0000000000 5.1300000000
5.1300000000 5.1300000000 0.0000000000
--- Equilibrium position and mass of each atom---
Natom= 2
1 1.2825000000 1.2825000000 1.2825000000 51196.42133 Si
2 -1.2825000000 -1.2825000000 -1.2825000000 51196.42133 Si
--- Vibrational modes ---
Nmode= 6 Natom= 2
n= 1 T1u IR
hbarW= -0.00000000E+00 Ha = -0.00000000E+00 eV; nu= -0.00000000E+00 cm^-1
1 0.0000000000 0.7071067812 0.0000000000
2 0.0000000000 0.7071067812 0.0000000000
...
...
n= 6 T2g R
hbarW= 0.23347395E-02 Ha = 0.63531493E-01 eV; nu= 0.51241611E+03 cm^-1
1 0.0000000000 0.0000000000 0.7071067812
2 0.0000000000 0.0000000000 -0.7071067812

--- primitive lattice vectors ---以降の三行にセルベクトルが記録され，--- Equilibrium position and mass of each atom---に続く行に原子座標データが記録されます。--- Vibrational modes ---以降が振動解析の結果です。Nmode= 6はモードの数が6つであることを意味します。n= ...から各モードのデータが記録されます。hbarW= ...の行でそのモードの固有エネルギーが，続く行に固有ベクトルが記録されます。ただし最初の3つのモードは系全体が並進するモードに対応するため意味を成しません。

mode.dataファイルから振動モード図を作成するPerlスクリプトがfreq.plです。以下の要領で実行します。

~/phase0_2023.01/bin/freq.pl mode.data

結果得られる振動モード図は 図 2.6 のようなものです。対称性に応じてモードを分類し，その振動数を表示します。

_images/image7.svg

振動モード図

この例題では有効なモードは1つしかないのであまり意味のある図ではありませんが，ある程度の振動モードの数が多い場合に振動モードを分かりやすく分類してくれるため有用な図となります。

EV曲線

複数の格子定数（単位胞体積）で計算を行い，マーナハンの状態方程式に結果をフィットすることによって安定な単位胞の体積，すなわち格子定数を得ることができます。マーナハンの状態方程式とは，以下に示すようなものです。

\[E_{\text{tot}}\left( V \right) = \frac{\text{BV}}{B^{'}(B^{'} - 1)}\left\lbrack B^{'}\left( 1 - \frac{V_{0}}{V} \right) + \left( \frac{V_{0}}{V} \right)^{B^{'}} - 1 \right\rbrack + E_{\text{tot}}\left( V_{0} \right).\]

ここで\(E_{\text{tot}}\left( V \right),\ B,\ B^{'},\ V_{0}\)はそれぞれ体積における全エネルギー，体積弾性率，体積弾性率の圧力微分，安定な体積に対応します。この表式から分かるように，この手続きによって格子定数だけでなく体積弾性率やその圧力微分を得ることもできます。この例題の入力ファイルは wsl/Si2/ev 以下にあります。

入力データの確認

この計算においては格子定数を変化させながら逐次計算を実施していきます。このような計算を実施する場合，入力パラメーターファイルのテンプレートを用意し，シェルスクリプトによって実際に用いる入力パラメーターファイルを作成しながら計算を実行していくと効率よく計算を行うことができます。

入力パラメーターファイル（のテンプレート）

この例題のために用意したテンプレート入力パラメーターファイルは下記の通り（SCF計算 と共通の部分は省略）。

structure{
  unit_cell{
    a = __A__
    b = __A__
    c = __A__
    alpha = 90
    beta = 90
    gamma = 90
  }
}

格子定数の指定が数値ではなく__A__という文字列になっています。この文字列は，シェルスクリプトによって実行時に格子定数の値に置き換わります。入力パラメーターファイルのほかの箇所と被らなければどのような文字列を採用しても問題ありません。

file_names.dataファイル

file_names.dataファイルの内容は下記のようになっています。

&fnames
F_POT(1) = '../../../pp/Si_ggapbe_nc_01.pp'
/

3階層上のディレクトリーの下のppというディレクトリーにある擬ポテンシャルファイルを指しています。ppディレクトリーはサンプルのディレクトリーの2階層上のディレクトリーですが，スクリプト実行時には子ディレクトリーを作成し，そこで計算が行われることを考慮した結果このような指定になっています。

シェルスクリプト

利用するシェルスクリプト(cubic.sh)は以下のようなものです。なお，行頭の数値は説明用に付加したものであり，実際はありません。

 1 #!/bin/sh
 2
 3 inivol=1000
 4 n=21
 5 dv=10
 6 np=2
 7 nk=2
 8 ne=1
 9 rm -f nfefn.data
10 PHASE0="mpiexec -n ${np} $HOME/phase0_2023.01/bin/phase ne=${ne} nk=${nk}"
11 for v in `seq 1 $n`;do
12   vol=$( echo "($v-1)*$dv + $inivol" | bc -l )
13   a=$( echo "e(l($vol)/3)" | bc -l )
14   echo "volume : $vol"
15   mkdir -p vol$vol
16   cp file_names.data vol$vol
17   sed "s/__A__/$a/g" nfinp.data > vol$vol/nfinp.data
18   cd vol$vol
19   ${PHASE0}
20   cd ..
21   echo -n $vol>> nfefn.data ; tail -1 vol$vol/nfefn.data >>nfefn.data
22 done

  • 3行目で1番目の体積の値，4行目で計算する体積の数，5行目で体積の刻み幅を指定しています。この例では，まず1000 bohr3 からはじめ，10 bohr3 きざみで体積を増やしていき，21回の計算を行うことで1200 bohr3 まで体積を大きくして計算を行うことになります。

  • 6行目はプロセス数，7行目はk点並列数，8行目はバンド並列数です。お使いの環境にあわせて編集してください。

  • 10行目でPHASE/0を起動するコマンドを変数PHASE0に格納しています。

  • 11行目からが実際の計算です。for文とseqコマンドを利用し1からnまでループを回しています。

  • 12行目，13行目ではbcコマンドを利用して体積と対応する格子定数を計算しています。

  • 15行目では計算に使うディレクトリーをmkdirコマンドによって作成しています。

  • 16行目でfile_names.dataファイルを計算用ディレクトリーにコピーしています。

  • 17行目では，sedコマンドによって先ほどのnfinp.dataファイル中の__A__という文字列を目的の格子定数に置き換え，計算用ディレクトリーの下のnfinp.dataというファイルに出力しています。

  • 18行目で計算用ディレクトリーに移り，18行目で計算を実行しています。

  • 20行目でもとのディレクトリーに移った後，21行目ではechoコマンドとtailコマンドを利用して結果をスクリプト実行ディレクトリーのnfefn.dataというファイルに積算しています。

このスクリプトは，3行目から8行目を編集することによって他の立方晶の結晶に適用することもできます。

計算の実行

以下のようなコマンドを実行します。

./cubic.sh > log &

この計算は少し時間がかかるので&によってバックグラウンドで実行するコマンドを採用しました。途中経過はlogファイルに記録されます。

結果の解析

このスクリプトを実行すると得られるnfefn.dataファイルの中身は以下のようになっているはずです。

$ cat nfefn.data
1000 1 11 -7.8711979339 0.0000000000
1010 1 11 -7.8720573907 0.0000000000
1020 1 11 -7.8728270483 0.0000000000
1030 1 11 -7.8734755993 0.0000000000
1040 1 12 -7.8740254367 0.0000000000
1050 1 11 -7.8744947325 0.0000000000
1060 1 11 -7.8748582628 0.0000000000
1070 1 12 -7.8751568272 0.0000000000
1080 1 12 -7.8753836059 0.0000000000
1090 1 12 -7.8755212142 0.0000000000
1100 1 12 -7.8755917399 0.0000000000
1110 1 12 -7.8755984193 0.0000000000
1120 1 12 -7.8755265291 0.0000000000
1130 1 12 -7.8753913443 0.0000000000
1140 1 12 -7.8751961234 0.0000000000
1150 1 12 -7.8749652461 0.0000000000
1160 1 12 -7.8746680963 0.0000000000
1170 1 12 -7.8743133309 0.0000000000
1180 1 12 -7.8739200962 0.0000000000
1190 1 12 -7.8734554597 0.0000000000
1200 1 12 -7.8729795665 0.0000000000

通常得られるnfefn.dataファイルの内容と似ていますが，1列目が原子単位の体積である点が異なります。

得られたnfefn.dataファイルとgnuplotを利用して，マーナハンの状態方程式にフィットしてみます。以下の操作を行ってください。

$ gnuplot
...
gnuplot> f(x) = (a*x/(b*(b-1)))*(b*(1-c/x)+(c/x)**b-1)+d
gnuplot> a=0.001
gnuplot> b=2
gnuplot> c=1100
gnuplot> d=-7.8
gnuplot> fit f(x) 'nfefn.data' using 1:4 via a,b,c,d
...
...
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = 0.00073905 +/- 1.202e-06 (0.1626%)
b = 4.22666 +/- 0.08447 (1.999%)
c = 1104.94 +/- 0.09036 (0.008178%)
d = -7.8756 +/- 2.387e-06 (3.031e-05%)
correlation matrix of the fit parameters:
a b c d
a 1.000
b -0.498 1.000
c 0.273 -0.893 1.000
d -0.744 0.287 -0.172 1.000

Final set of parametersに続くデータがフィッティングの結果です。この例では，aが体積弾性率，bが体積弾性率の圧力微分，cが単位胞の体積，dが安定な単位胞体積における全エネルギーに相当します。したがって，まず単位胞の体積は1104.94と求まったことが分かります。三乗根をとると，格子定数はおおよそ10.34 bohr (5.47 Å)となります。また，aは体積弾性率に相当しますが，少し注意が必要です。格子定数はブラベー格子のもので与えていますが，実際の計算は基本格子で行っています。面心立方格子では，ブラベー格子の体積は基本格子の体積の四倍なのでここで得られた体積弾性率の結果も四倍する必要があります。この点に留意すると，体積弾性率は原子単位で0.00073905×4 = 0.0029562, つまりおおよそ87 GPaとなります。得られるEV曲線は，下図のようなものになります。

_images/image8.svg

本例題で得られるEV曲線。

磁性を考慮した計算

前章では非磁性の計算を行いましたが，磁性を考慮した計算を行うことも可能です。ここでは体心立方格子の鉄結晶を例に，磁性を考慮した計算を行う方法を紹介します。例題の入力ファイルは wsl/Fe 以下に配置されています。計算の流れはシリコン結晶の場合と同様で，SCF計算を行い，結果得られた電荷密度を入力としてバンド構造や状態密度の計算を行います。同様の手続きでEV曲線の計算なども行うことができます。

SCF計算

磁性を考慮する場合も基本となるのはSCF計算です。例題の入力ファイルは wsl/Fe/scf の下にあります。

入力データの確認

SCF計算の入力の内容は下記の通り(前章と同じ部分は一部省略)

accuracy{
  ksampling{
    mesh{ nx = 10, ny = 10, nz = 10 }
  }
}
structure{
  unit_cell{
    #units angstrom
    a = 2.845, b = 2.845, c = 2.845
    alpha = 90, beta = 90, gamma = 90
  }
  symmetry{
    method = automatic
    tspace{
      lattice_system = bodycentered
    }
  }
  magnetic_state = ferro
  atom_list{
    atoms{
      #tag rx ry rz element
      0.000 0.000 0.000 Fe
    }
  }
  element_list{
    #tag element atomicnumber zeta
    Fe 26 0.275
  }
}

  • 原子配置の与え方としてBravais格子を与える方法を採用しています。単位胞としては体心立方格子の基本格子ではなく立方格子のそれをあたえ，体心立方格子であることはsymmetryのtspaceブロックにおけるlattice_system = bodycenteredで指定しています。このように指定すると，プログラムは格子ベクトルを自動的に基本格子のそれに変換します。原子座標の指定もそれに合わせ，原点のみの指定となっています。

  • structureブロックのmagnetic_stateにferroを指定しています。このキーワードはややミスリーディングですが，強磁性体に限られるわけではなく，スピンを考慮した計算を行う，という設定です。

  • element_listにおいてzetaの値を指定しています。zetaは“初期スピン分極”の設定です。スピン分極\(\zeta\)とは，アップスピンの電子数を\(n_{\uparrow}\)ダウンスピンの電子数を\(n_{\downarrow}\)とした場合\(\zeta = \frac{n_{\uparrow} - n_{\downarrow}}{n_{\uparrow} + n_{\downarrow}}\)と定義されます。

  • 金属なので，Si結晶の場合と比較するとサンプリングk点を多くしています。

計算の実行

これまでと同じ手続きで実行することができます。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

以下の要因から，原子数自体は鉄結晶よりも多いシリコン結晶の計算より有意に時間のかかる計算となります。

  • スピンを考慮している

  • サンプリングk点の数が多い

  • バンド数（電子数）が多い

  • d電子を含んでいる

結果の解析

スピンを考慮したので，スピン分極の値をもとめることができます。以下のようなコマンドによってアップスピンおよびダウンスピンの電子数の履歴を調べることができます。

$ grep NEW output001
!*--- input-file style = NEW
!NEW total charge (UP, DOWN, SUM) = 5.09024391 (+) 2.90975609 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.08351594 (+) 2.91648406 (=) 8.00000000
...
!NEW total charge (UP, DOWN, SUM) = 5.17484294 (+) 2.82515706 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483127 (+) 2.82516873 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483897 (+) 2.82516103 (=) 8.00000000
!NEW total charge (UP, DOWN, SUM) = 5.17483935 (+) 2.82516065 (=) 8.00000000

最終的にアップスピンの電子数が5.1748935ダウンスピンの電子数が2.82516065となりました。すなわち，スピン分極は0.294程度ということになります。

バンド計算

バンド計算は，Siの場合と同様固定電荷計算を用いて行います。例題の入力ファイルは wsl/Fe/band 以下に配置されています。

入力データの確認

nfinp.dataファイル

Siの場合と同様，controlブロックにおいてcondition = fixed_chargeという設定がほどこされています。また，accuracyのksamplingブロックにおいてmethod = fileとすることによってk点生成方法をファイルからの読み込みとしています。

kpoint.dataファイル

Siの場合と同様kpoint.dataファイルを生成する必要があります。以下のコマンドを実行してください。

$ ~/phase0_2023.01/bin/band_kpoint.pl bandkpt.in
Distance of 1 = 1
Distance of 2 = 0.707106781186548
Distance of 3 = 0.5
Distance of 4 = 0.866025403784439
Distance of 5 = 0.707106781186548
division numbers = 100 70 50 86 70
i=0 -0.005 0.005 0.005
i=1 0.00714285714285714 -0.00714285714285714 0
i=2 0.005 0.005 -0.005
i=3 -0.00290697674418605 -0.00290697674418605 -0.00290697674418605
i=4 0 0 0.00714285714285714
1 : 0 0 0
    => 0/200 0/200 0/200
...

計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

計算結果の解析

バンド構造の計算結果はnfenergy.dataファイルに記録されます。Siの場合と違い，アップスピンとダウンスピン状態の固有値が両方とも記録されます。

band.plをたとえば以下のように実行します。

~/phase0_2023.01/bin/band.pl nfenergy.data bandkpt.in -erange=-10,10 -color -with_fermi

鉄の場合はエネルギーの範囲が広いので，-erangeオプションを用いてフェルミエネルギーを基準に-10 eVから10 eVの状態を対象としました。-color, -with_fermiオプションも有効にしました。

_images/image11.svg

鉄結晶のバンド構造

状態密度計算

Siの場合と同様，固定電荷計算を用いて状態密度の計算を行います。対応する入力ファイルは wsl/Fe/dos 以下に配置されています。

入力データの確認

nfinp.dataファイル

Siの場合と同様，accuracyブロックのksamplingブロックにおいてmethodをmeshとし，さらにsmearingブロックにおいてmethodをtetrahedralに設定しています。また，postprocessingブロックのdosブロックにおいてsw_dosをonとし，methodをtetrahedralとしています。

計算の実行

固定電荷計算はekcalプログラムで行います。以下の要領で実行してみてください。

mpirun -n 2 ~/phase0_2023.01/bin/ekcal ne=1 nk=2

計算結果の解析

スピン分極を考慮している場合，dos.dataファイルにはアップスピンの状態密度とダウンスピンの状態密度が記録されます。そのため，考慮していない場合と比較して1行のカラム数が増えます。具体的には，各行次のような並びでデータが記録されます。

エネルギーの識別子/ハートリー単位のエネルギー/ハートリー単位の場合のアップスピン電子の状態密度/ハートリー単位の場合のダウンスピン電子の状態密度/eV単位のエネルギー/eV単位の場合のアップスピン電子の状態密度/eV単位の場合のダウンスピン電子の状態密度/アップスピン電子の積算状態密度/ダウンスピン電子の積算状態密度/全電子の積算状態密度

eV単位の場合0がフェルミエネルギーになるようにシフトされている点はスピンを考慮していない場合と同じです。

dos.plをたとえば以下のように実行します。

~/phase0_2023.01/bin/dos.pl dos.data -with_fermi -color -erange=-10,10

バンド計算の場合と同様 -erangeオプションを用いてフェルミエネルギーを基準に-10 eVから10 eVの状態を対象としました。-color, -with_fermiオプションもシリコン結晶の場合と同様有効にしました。図 3.2 のような状態密度図が得られるはずです。

_images/image10.svg

鉄結晶の状態密度

表面の計算

PHASE/0は平面波基底を利用するので扱うことが可能なのは周期的な系に限られますが，ある程度の大きさの“真空層”を設けることによって事実上表面と変わらない系を扱うことが可能です。この例題では，Si(100)表面の計算を行います。例題の入力ファイルは wsl/Si001 以下に配置されています。

構造最適化

本例題の表面は次の図で示すようなものです。なお，この系は原点を中心とした反転対称性があるように設定されており，この図では反転対称性から等価な位置の原子は描画されていません。すなわち，実際に計算で扱う系は描画された原子の反転対称位置にも原子が配置されている系です。まずはこの構造の最適化を行います。例題の入力ファイルは wsl/Si001/relax の下にあります。

_images/image12.png

Si(100)面。反転対称位置の原子は描画されていない。

入力データの確認

本例題の入力ファイルの主要部分は下記の通り。

...
accuracy{
  ksampling{
    mesh{
      nx = 2
      ny = 4
      nz = 1
    }
  }
  force_convergence{
    max_force = 2e-4
  }
  ...
}
structure{
  symmetry{
    method = automatic
    sw_inversion = on
  }
  atom_list{
    atoms{
      #tag element rx ry rz weight mobile
      Si 0.315656104534 0.75 0.237152295345 2 on
      Si 0.598029122106 0.75 0.209636130373 2 on
      Si 0.26602495266 0.25 0.185466403501 2 on
      Si 0.73854288575 0.25 0.182492535522 2 on
      Si 0.00536283827178 0.25 0.136873588205 2 on
      Si 0.492058345921 0.25 0.12432672152 2 on
      Si 0.00251446092889 0.75 0.0840435403934 2 on
      Si 0.501404226089 0.75 0.0727860480334 2 on
      Si 0.238553794997 0.75 0.0258423559513 2 off
      Si 0.765555448939 0.75 0.0267457921579 2 off
    }
  }
  unit_cell{
    #units bohr
    a_vector = 14.512 0.0 0.0
    b_vector = 0.0 7.25600000002 0.0
    c_vector = 0.0 0.0 49.6812599214
  }
  ...
  ...
}
structure_evolution{
  method = lbfgs
}

以下の点が特徴的です。

  • accuracyのforce_covergenceブロックにおいて構造最適化の収束判定条件が設定されています(max_force = 2e-4 ; 単位はデフォルト値のhartree/bohr)。ここで指定した値よりも最大原子間力が小さい場合に収束したとみなされます。

  • 原子配置の指定にweight属性を利用しています。この属性値の値が2の場合，原点からみて反転対称の位置に同じ原子が配置されます。

  • c軸の長さが50 bohr(26.5Å程度) 近くとなっています。この設定によって13から14Å程度の真空層を確保しています。

  • c軸方向のk点サンプリングメッシュが1となっています。表面モデルの場合通常表面に垂直な方向のk点メッシュは1とします。違う言い方をすると，1としてよい（表面に垂直な方向の分散は無視できる）程度には厚い真空層を設ける必要があります。

  • 原子配置の指定にmobile属性を利用しています。この例では，最下層の原子のmobileをoffとすることによって構造最適化の対象から外しています。

  • symmetryブロックにおいて，method = automaticだけでなくsw_inversion = onを利用しています。sw_inversionは系に反転対称性がある場合のみ利用できます。

  • structure_evolutionブロックにおいてmethod = lbfgsという記述を行うことによって構造最適化の手法としてlbfgs法を採用しています。

計算の実行

これまでと同じ手続きで計算を実行することができます。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

この例題のk点数は2なので，多くのコアを搭載したマシンでもk点並列数は2以下とする必要があります。

これまでと違い，“構造最適化”が行われるのでSCF計算が収束すると原子配置が計算された原子間力に応じて更新され，次のSCF計算に移行します。構造最適化計算の進行はnfefn.dataファイルの内容をみることによって確認することができます。

$ less nfefn.data

iter_ion, iter_total, etotal, forcmx
1 27 -78.6015420181 0.0074909920
2 34 -78.6020126027 0.0067584284
3 45 -78.6028324863 0.0060223336
4 52 -78.6030036762 0.0045697291
5 59 -78.6033506302 0.0037802003
...

1行が構造最適化のあるステップの結果に対応します。1カラム目が構造最適化のステップ番号，2カラム目がSCF計算の総繰り返し回数，3カラム目が全エネルギー，4カラム目が原子間力の最大値です。4カラム目の数値が収束判定条件よりも小さな値になると収束したとみなされ計算が終了します。

結果の解析

上述のようにnfefn.dataファイルにステップごとのエネルギーや原子間力が記録されるので，このファイルの中身を確認することによってどのように収束していったかを調べることができます。以下のようにgnuplotを用いて履歴のグラフを作成することもできます。

$ gnuplot
...
Terminal type is now 'wxt'
gnuplot> plot 'nfefn.data'using 1:3 with lines title 'energy'
gnuplot> replot 'nfefn.data'using 1:4 with lines axis x1y2 title 'max.  force
gnuplot> set y2tic
gnuplot> replot

上述の操作の結果 図 4.2 のようなプロットが得られます。

_images/image13.png

エネルギーと原子間力の最大値の履歴

原子座標データはnfdynm.dataファイルに記録されます。このファイルをVESTAでそのまま可視化することはできないので，CIFなどのより一般的な形式に変換します。PHASE/0にはconv.pyという座標データを変換するツールが備わっているので( 座標データ変換 )，これを用いてnfdynm.dataファイルをCIFに変換します。結果得られたCIFをVESTAで読み込むことによって構造最適化の座標履歴を確認することができます。

局所状態密度計算と仕事関数

構造最適化 によって得られた安定な表面構造を入力とし，局所状態密度や仕事関数の計算を行います。

入力データの確認

nfinp.dataファイル

nfinp.dataファイルの内容は以下の通り（構造最適化計算と同じである部分は省略）

accuracy{
  ksampling{
    method = mesh
    mesh{
      nx = 4
      ny = 8
      nz = 1
    }
  }
  smearing{
    method = tetrahedral
  }
  ...
}
structure{
  method = file
  file{
    filetype = phase0_output
  }
  atom_list{
    atoms{
      #tag element rx ry rz weight
      Si 0.315656104534 0.75 0.237152295345 2
      ...
      Si 0.765555448939 0.75 0.0267457921579 2
    }
  }
}
postprocessing{
  workfunc{
    sw_workfunc = on
  }
  dos{
    sw_dos = on
  }
  ldos{
    sw_aldos = on
  }
}

  • 状態密度を四面体法で，かつ構造最適化よりも濃い密度のk点メッシュで計算する設定が施されています。すなわちksamplingのmethodはmesh, smearingのmethodはtetrahedral, k点メッシュは4×8×1となっています。

  • 構造最適化計算の最後の結果を用いるため，structureブロックにおいてmethod = fileを指定しています。さらにfileブロックにおいてfiletypeをphase0_outputとしています。このように設定するとnfdynm.dataファイルの最後のコマを入力座標として読み込みます。このような場合でも原子の属性値の定義を行うためにatom_listブロックのatomsテーブル自体は必要です。

  • postprocessingブロックで状態密度計算，原子分割局所状態密度計算，仕事関数計算の設定が施されています。

file_names.dataファイル

file_names.dataファイルの内容は以下のようになっています。

&fnames
F_POT(1) = '../../pp/Si_ggapbe_paw_nc_01m.pp'
F_POS = '../relax/nfdynm.data'
/

F_POSによって構造最適化計算で得たnfdynm.dataファイルを指しています。

計算の実行

これまでと同じ手続きで計算を実行することができます。

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

この例題のk点数は4なので，多くのコアを搭載したマシンでもk点並列数は4以下とする必要があります。

結果の解析

原子分割局所状態密度の結果はdos.dataファイルに全状態密度のあとに記録されます。dos.plスクリプトを用いることによって各原子に割り当てられた状態密度を得ることができます。

~/phase0_2023.01/dos.pl dos.data -mode=atom -color -with_fermi

この操作の結果dos_a001.eps, dos_a002.eps, ....といったEPSファイルが得られます。参考のため，得られる結果の一部の図を紹介します。

1番目の原子と2番目の原子の局所状態密度(dos_a001.epsおよびdos_a002.eps)

_images/image14.svg

局所状態密度

workfuncプログラムを使用することによって仕事関数を得ることができます。まずはworkfuncプログラムをコンパイルしましょう。

$ pushd ~/phase0_2023.01/src_workfunc
$ make F90=ifort
ifort -c -O  m_Const_Parameters.f90
ifort -c -O   m_ArraySize_Parameters.F90
...
$ make install
mv workfunc ../bin/
$ popd

作業中のディレクトリーにすぐに戻ってこられるよう pushd popd コマンドを使ってみました。この例ではworkfuncプログラムをコンパイルする際に F90=ifort としてIntel Fortranコンパイラーを利用するようにしています。デフォルトで用いられるコンパイラーはgfortranですが，gfortranのバージョン10以上を用いる場合はオプションに -fallow-argument-mismatch を加える必要があります。そこで，gfortran 10以上を用いる場合は make コマンドに F90='gfortran -fallow-argument-mismatch' を渡すようにしてください。

workfuncプログラムを実行すると得られるnfvlcr_av.dataファイルから真空域のポテンシャルを見出します。その値とフェルミエネルギーの差が仕事関数に対応します。もしくは，workfunc.plスクリプトを利用します。

$ ~/phase0_2023.01/bin/workfunc
$ ~/phase0_2023.01/bin/workfunc.pl nfvlcr_av.data
estimated work function : 4.74259 eV

図 4.4 はworkfunc.plスクリプトによって得られたポテンシャルとc軸方向の距離の関係です。真空域でのポテンシャルの値とフェルミエネルギーとの差が仕事関数に対応します。

_images/image15.svg

ポテンシャルとc軸方向の距離の関係

点欠陥の生成エネルギーの計算

PHASE/0による構造最適化計算や原子・分子の計算などを組み合わせることによって様々な生成エネルギーの計算を計算することが可能です。ここでは，体心立方バナジウムの点欠陥を例にPHASE/0を利用して生成エネルギーを計算します。サンプルデータは wsl/V 以下にあります。なお，この例題の入力設定は計算負荷を下げるため粗いk点サンプリングを採用し，また電子数の少ない非推奨の擬ポテンシャルを用いています。そのため得られる結果そのものは信頼できるものではありません。計算手順を理解していただくことが目的の例題です。

Self-interstitialの生成エネルギー

母体の結晶と同じ原子種の元素が格子間に入り込む欠陥がself-interstitial欠陥です。入り込む位置の候補として，体心立方格子の場合tetrahedral位置とoctahedral位置が考えられます。tetrahedral位置とは 図 5.1 の原子ABCDが作る四面体の中心位置であり，たとえば(0.5, 0.25, 0.0)の位置です。体心立方格子におけるoctahedral位置とは右の図の原子ABCDEFが作る八面体の中心位置であり，たとえば(0.5,0.5,0.0)の位置です。Interstitialのほか，単純に格子位置に原子が存在しない空孔欠陥などもサンプルデータとしては用意しています。

_images/image17.png

体心立方格子のサイト

入力データの確認

この例題の入力ファイルは wsl/V/self 以下のサブディレクトリーの下に配置されています。各入力パラメーターファイルは，原子座標データの欠陥が導入されていること以外は特別な点はありません。構造最適化は必須であり，すべての原子が最適化の対象となるような設定が施されています。典型的には以下のような内容です。

structure{
  atom_list{
    atoms{
    #units angstrom
    #tag element rx ry rz mobile
      V -0.0283241509 -0.0283241509 0.00 1
      ...
    }
  }
  element_list{
    #tag element atomicnumber mass zeta deviation
      V 23 92860.1055 0.0 1.83
  }
  symmetry{
    method = automatic
  }
  ...
}
accuracy{
  cutoff_wf = 25 Rydberg
  cutoff_cd = 225 Rydberg
  ksampling{
    mesh{
      nx = 2
      ny = 2
      nz = 2
    }
  }
}

また，参照用の欠陥のないスーパーセルの入力ファイルは wsl/V/ref/V 以下にあります。

計算の実行

これまでと同じように計算を実行します。ただ，この例題は計算時間短縮のためk点サンプリングを粗くとっているため1点しかk点がない場合も多いです（対称性に依存します）特に引数をつけずに実行すると可能な場合はk点並列を併用し，そうでない場合はバンド並列のみで計算が実行されるので引数を渡さずに実行すると効率がよいかもしれません。

mpirun -n 2 ~/phase0_2023.01/bin/phase

結果の解析

Self-interstitialの生成エネルギーは以下のように評価することができます。

\[E_{\text{f}} = E_{\text{tot}}\left( \text{Crystal} + \text{point defect} \right) - N_{\text{atom}}E_{\text{tot}}\left( \text{Crystal} \right)\]

ここで\(E_{\text{tot}}\left( \text{Crystal} \right)\)は結晶の原子1個当たりの全エネルギー，\(E_{\text{tot}}\left( \text{Crystal} + \text{point defect} \right)\)は点欠陥を含むスーパーセルの全エネルギーです。結晶のエネルギーは，k点サンプリングが十分なら基本格子のエネルギーから算出することもできますが，欠陥のスーパーセルと同じ大きさの系を採用する方が無難です。本例題ではrefディレクトリーに欠陥のないスーパーセルの入力データが用意されているのでその結果を用います。

全エネルギーの計算結果はnfefn.dataファイルの最後の行に記録されます。たとえば，tetrahedral欠陥の場合エネルギーの計算値は-110.4633242498 hartree, 欠陥のないスーパーセルの全エネルギーは-104.0740790123 hartreeなので，欠陥生成エネルギーは\(- 110.4633242498 - 17 \times \left( \frac{- 104.0740790123}{16} \right) \approx 0.1154\ \text{hartree}\)

となります。eV単位にすると3.14 eV程度となります。

Impurity interstitialの生成エネルギー

母体の結晶とは異なる原子種の元素が格子間に入り込む欠陥がimpurity interstitial欠陥です。この場合も入り込む位置はtetrahedral位置，octahedral位置が考えられます。母体の結晶の原子と入れ替わるsubstitutional欠陥も考えられます。

入力データの確認

この例題の入力ファイルには，原子座標データの欠陥が導入されていること以外は特別な点はありません。構造最適化は必須であり，すべての原子が最適化の対象となるような設定が施されています。

この例題では異種元素はヘリウム原子を想定しています。そこで，生成エネルギーを求めるためにヘリウム原子の計算を行う必要があります。原子の計算は十分大きな単位胞に孤立した原子を1個置くことによって行うことができます。この際，k点サンプリングはΓ点のみとします。カットオフエネルギーは他の計算のそれと同じ値を採用するようにしてください。具体的には下記のような内容になります。

accuracy{
  ksampling{
    method = gamma
  }
  cutoff_wf = 25 Rydberg
  cutoff_cd = 225 Rydberg
  ...
}
structure{
  atom_list{
    atoms{
      #units angstrom
      #tag element rx ry rz mobile
      He 0.0 0.0 0.0 1
    }
  }
  ...
  unit_cell{
    #units angstrom
    a_vector = 15 0.00 0.00
    b_vector = 0.00 15 0.00
    c_vector = 0.00 0.00 15
  }
  symmetry{
    method = automatic
  }
}

計算の実行

計算の実行もself-interstitialの場合と全く同じ要領で行うことができます。

結果の解析

Impurity-interstitialの生成エネルギーは以下のように評価することができます。

\[E_{\text{f}} = E_{\text{tot}}\left( \text{Crystal} + \text{point defect} \right) - \left( N_{\text{atom}} - 1 \right)E_{\text{tot}}\left( \text{Crystal} \right) - E_{\text{tot}}(\text{atom})\]

ここで\(E_{\text{tot}}\left( \text{Crystal} \right)\)は結晶の原子1個当たりの全エネルギー，\(E_{\text{tot}}\left( \text{Crystal} + \text{point defect} \right)\)は点欠陥を含むスーパーセルの全エネルギー，\(E_{\text{tot}}(\text{atom})\)は原子1個の全エネルギーです。

たとえば，tetrahedral位置にヘリウム原子が入り込んだ系の全エネルギーは-106.8410461476 hartree 結晶のスーパーセルの全エネルギーは前節でも言及したように -104.0740790123 hartreeヘリウム原子の全エネルギーは-2.8627340281 hartreeなので欠陥生成エネルギーは

\[- 106.8410461476 - \left( 17 - 1 \right) \times \left( \frac{- 104.0740790123}{16} \right) - ( - 2.8627340281) \approx 0.0958\ \text{hartree}\]

となります。eV単位に変換すると2.61 eV程度です。

付録

PHASE/0による計算のtips

計算に必要なファイル

PHASE/0による計算を実行するには，以下のファイルが必要です。

  • file_names.dataファイル

  • nfinp.data ファイル（入力パラメーターファイル）

  • 元素分の擬ポテンシャルファイル

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

入力パラメーターファイル（nfinp.dataファイル）は，計算対象の系や計算条件などの設定を行うファイルです。その内容は多岐にわたるので詳細な説明はしませんが，以下の設定は重要なのでサンプルの入力を確認する場合は留意してください。

  • accuracyブロック：カットオフエネルギー，k点サンプリング手法，収束判定条件など，計算精度に関わる設定を行います。

  • structureブロック：原子配置や元素の定義を行います。

  • postprocessingブロック：計算が収束したあとに行う後処理の指定を行います（状態密度計算など）

入力パラメーターファイルの以下の点は間違えやすいので，注意が必要です。

  • 単位の間違い。Å単位のつもりがbohr単位だった，など

  • “coordinate_system”の間違い。cartesianのつもりだったのがinternalだった，もしくはその逆。

  • 元素の定義順の間違い。element_listブロックにおける元素の定義順は，file_names.dataにおける擬ポテンシャルファイルの定義順と同じである必要がある。

計算の実行

並列計算は，nfinp.dataファイルやfile_names.dataファイルのあるディレクトリーに移り，mpiexecコマンドやmpirunコマンドを介して実行します。たとえば以下のようなコマンドになります。

mpirun -n NP ~/phase0_2023.01/bin/phase ne=NE nk=NK

NPにはプロセス数（コア数），NEにはバンド並列数，NKにはk点並列数を指定します。NPはお使いの環境に合わせて入力し，NP = NE × NKとなるようバンド，k点の並列数を指定します。また，最後に&と入力して実行するとバックグラウンドで計算を実行することができます。

並列数を決める

最適な並列数を決めるためのポイントはいくつかあります。

  • 総並列数の上限は，利用している環境のコア数

  • 総並列数をNPとするとバンド並列数NEとk点並列数NKとの間にはNP = NE × NKという関係が成立している必要がある

  • 通常バンド並列数がバンド数を超えることはなく，あったとしてもPHASE/0は問題なく動作する

  • k点数を超えたk点並列数を指定すると異常終了する

  • 多くの場合バンド並列よりもk点並列の方が効率がよい

  • ただしk点数がk点並列数で割り切れないと効率が悪い

以上のことから，“事前にk点数が分かっていれば効率の良い並列のさせ方が分かる”ということが分かります。k点数は入力でのメッシュの指定の仕方や系の対称性などに依存するのでその数は自明ではありません。事前に調べるには，以下のような方法を利用することができます。

  • 入力パラメーターファイルのcontrolブロックのconditionの値をpreparationとする

  • PHASE/0を非並列で実行する

  • 計算はすぐに終了する。”kv3”という文字列をoutput000ファイルから検索する。kv3 = のあとに続く数値がk点数（スピンを考慮した計算の場合，その半分）
    $ grep kv3 output000
    !kp kv3 = 10 nspin = 1

  • 本計算を実行する前にconditionパラメーターの値を元に戻すのを忘れないように

SCF計算の収束具合を確かめる

SCF計算の収束具合を確認するには，以下の要領で文字列THをoutput000ファイルから検索し，全エネルギーや1ステップ前のエネルギーとの差を調べます。

$ grep TH output000
TOTAL ENERGY FOR  1 -TH ITER= -7.837000950503 edel = -0.783700D+01 ...
TOTAL ENERGY FOR  2 -TH ITER= -7.839515220565 edel = -0.251427D-02 ...
TOTAL ENERGY FOR  3 -TH ITER= -7.839908027886 edel = -0.392807D-03 ...
TOTAL ENERGY FOR  4 -TH ITER= -7.840075864037 edel = -0.167836D-03 ...
TOTAL ENERGY FOR  5 -TH ITER= -7.840148564528 edel = -0.727005D-04 ...
TOTAL ENERGY FOR  6 -TH ITER= -7.840207027180 edel = -0.584627D-04 ...
TOTAL ENERGY FOR  7 -TH ITER= -7.840209759424 edel = -0.273224D-05 ...
TOTAL ENERGY FOR  8 -TH ITER= -7.840209810085 edel = -0.506609D-07 ...
TOTAL ENERGY FOR  9 -TH ITER= -7.840209813625 edel = -0.353987D-08 ...
TOTAL ENERGY FOR 10 -TH ITER= -7.840209814325 edel = -0.699737D-09 ...
TOTAL ENERGY FOR 11 -TH ITER= -7.840209814372 edel = -0.470983D-10 ...

TOTAL ENERGY FOR N -TH ITER= の後の数値がN回目のSCF計算で得られているエネルギー，そのあとのedel=の後の数値が1ステップ前のエネルギーとの差です。なお，この数値は全セルのエネルギーですが，収束判定は原子1個あたりのエネルギーで行われるので注意が必要です。edelをプロットするためのgnuplotのコマンド例は下記の通りです。

$ grep TH output000>TH
$ gnuplot
gnuplot> set datafile fortran
gnuplot> plot 'TH' using (abs($10)) with lines title 'edel'
gnuplot> set log y
gnuplot> replot
_images/image18.png

SCF iteration回数とedelの関係

計算のログファイルを監視する

計算のログファイル(output000ファイル)は計算の進行があるたびにその内容が更新されます。その様子を監視するには，tail -fコマンドが便利です。

$ tail -f output000
...
...

なお，計算が終了するとなぜ終了したかが表示され，以降更新されなくなります。

構造最適化の履歴

構造最適化の履歴は，nfefn.dataファイルとnfdynm.dataファイルに記録されます。nfefn.dataファイルの中身は，以下のようなものです。

$ less nfefn.data
iter_ion, iter_total, etotal, forcmx
1 21 -78.5809455694 0.0172429910
2 33 -78.5823489899 0.0145375571
3 51 -78.5856703762 0.0061089674
4 58 -78.5860239085 0.0051604695
5 70 -78.5872784872 0.0046626573
6 79 -78.5874820406 0.0045542295
....
....

1列目は原子配置の更新回数，2列目はSCF計算の総回数，3列目は全エネルギー，4列目は原子間力の最大値です。4列目の数値が閾値 (デフォルト値は1e-3) を下回ると収束したとみなされ構造最適化計算は終了します。

nfefn.dataファイルの内容は単純なので，プロットツールなどを利用してその内容を簡単にプロットすることができます。以下はgnuplotのコマンド例です。

$ gnuplot
gnuplot> plot 'nfefn.data' using 1:3 with lines title 'energy'
gnuplot> replot 'nfefn.data' using 1:4 with lines axis x1y2 title 'max. force'
gnuplot> set log y2
gnuplot> set y2tics
gnuplot> replot
_images/image19.png

構造最適化のステップ数とエネルギーおよび原子間力最大値の関係のプロット

nfdynm.dataファイルには座標データの履歴が記録されます。付属の変換スクリプトconv.py( 6.2 章 )を利用して様々な形式に変換することができます。conv.pyは対話的に利用することができるので，実行し，指示に従い入力していけば座標データファイルの変換を行うことができます。

SCF計算がなかなか収束しない場合の対処

  • ミクシング比を小さくする
    charge_mixing{
    rmx = 0.1 ! デフォルト値はスピンを考慮しない場合0.4, する場合0.1
    }
  • バンド数を増やす（SCF計算の場合）
    accuracy{
    num_bands = ...
    }
  • “num_extra_bands”を増やす（固定電荷計算の場合）
    accuracy{
    ek_convergence{
    num_extra_bands = 20
    }
    }
    num_extra_bandsのデフォルト値はバンド数に依存；output000ファイルのdefault value for num_extra_bandsに続く数値
  • RMMに切り替えるタイミングを厳しくする（もしくは無効にする）
    wavefunction_solver{
    rmm{
    edelta_change_to_rmm = 0 ! default : 1e-3 hartree/atm
    }
    }

単位

PHASE/0においては，原則としてハートリー原子単位系が用いられます。ハートリー原子単位系は電子状態計算の理論の記述には便利ですが，それ以外の用途ではあまり利用されないため変換が必要な場合が多いです。ハートリー原子単位系と別の単位系との変換係数をいくつか記述します。

エネルギー

1 hartree = 2 rydberg = 27.21139615 eV = 4.359745836 \(\times 10^{- 18}\) J

長さ

1 bohr = 0.5291772480 Å= 0.5291772480 \(\times 10^{- 10}\)m

質量

1 au mass = 電子の質量 = 9.1094\(\times 10^{- 31}\) kg

力

1 hartree/bohr = 51.42208259 eV/Å= 8.238725025 \(\times 10^{- 8}\) N

時間

1 au time = 2.418884327 \(\times 10^{- 2}\) fs = 2.418884327 \(\times 10^{- 17}\) s

ストレス

1 au stress = 2.903628623 \(\times 10^{9}\) atm = 2.942101703 \(\times 10^{13}\) Pa

座標データ変換

PHASE/0は独自の形式(F_DYNM形式)で原子座標の履歴などを出力します。これをVESTAなどの一般的な原子配置可視化ソフトで可視化するためには，CIF形式などの形式に座標を変換する必要があります。その変換のためのPythonスクリプトconv.pyはPHASE/0インストールディレクトリーの下のbinディレクトリーに含まれています。

conv.pyは対話的に利用します。たとえば，nfdynm.dataファイルをCIFに変換する手続きは下記の通りです。

conv.pyの利用例

画面に現れる文字列

説明

$ conv.py

atomic configuration converter

utility.

Copyright (C) PHASE System Consortium

select the type of the input atomic coordinate file

  1. phase_input

  2. phase_output

  3. VASP_input

  4. VASP_output

  5. OpenMX_input

  6. OpenMX_output

  7. XSF

  8. xyz

  9. cube

  10. cif

  11. dmol

  12. LAMMPS_output

  1. Exit

Please enter a selection (0/1/2/3/4/5/6/7/8/9/10/11/x) [0]:

変換元のファイル形式

を指定する。nfdynm.dataの場合ph

ase_outputなので1を指定し，Enter

Please enter the name of the

input atomic coordinate file, or

type x to exit. [nfdynm.data]:

nfdynm.data

ファイルのファイル名を指定。nf

dynm.dataでよいならそのままEnter

Please enter the frame no.

(enter a negative value in order

to output all frames when

possible), or type x to exit.

[-1]:

フレームを

選択する。負の値の場合「全フレー

ム」を選択することに相当する。ま

た，カンマ区切りで三つの整数を指

定することによって，

始フレーム, 終フレーム，間隔を

指定することができる。この指定方

法の場合負の値は最終フレームを意

味する。フレームの数値は0始まり

select the type of the output

atomic coordinate file

  1. phase_input

  2. phase_output

  3. VASP_input

  4. OpenMX_input

  5. XSF

  6. xyz

  7. cube

  8. cif

  9. dmol

  10. LAMMPS_input

  1. Exit

Please enter a selection (0/1/2/3/4/5/6/7/8/9/x) [1]:

変換先の形式を指

定する。CIFの場合7と入力しEnter

Please enter the name the output

atomic coordinate file, or type

x to exit. [coord.cif]:

出力ファイル名を指定する。CI

Fの場合，デフォルト値はcoord.cif

以上の操作によって，nfdynm.dataファイルからCIFを得ることができます。そのほかのファイル形式についても同様に変換することができます。

conv.py起動時に，以下のオプションを指定することができます。

オプション

説明

--pack

単位胞の中に原子を押し込めます

--na=NA

a軸をNA倍にしたスーパーセルを作成し，変換します

--nb=NB

b軸をNB倍にしたスーパーセルを作成し，変換します

--nc=NC

c軸をNC倍にしたスーパーセルを作成し，変換します

Linux tips

タブ補完

Linuxのコマンドラインにおいてはコマンド名やファイル名を入力しますが、これらは「タブ補完」によって効率よく入力することができます。たとえば以下のコマンドの場合

mpirun -n 2 ~/phase0_2023.01/bin/phase ne=1 nk=2

~/phaseまで入力してタブキーを押下すると phase0_2023.01 まで補完されます。長いコマンドやファイル名を入力する際に間違いなく簡単に入力することがででるので使ってみてください。

Linuxコマンド

本チュートリアルで用いたLinuxコマンドについて説明します。

本チュートリアルで用いたLinuxコマンド

コマンド

説明

cd

引数で指定したディレクトリーを変更する。引数な

しで実行すると，ホームディレクトリーに移る。

. はカレントディレクトリー，

~ はホームディレクトリー，

/ はルートディレクトリーを表します。

例：

cd abc (abcというディレクトリーに移る)

cd (ホームディレクトリーに移る)

ls

ディレクトリーの下のファイルを表示する。

例：

$ ls

nfinp.data file_names.data

less

ファイルの中身を参照する（編集は不可）

例：nfinp.dataファイルの中身を参照する

$ less nfinp.data

control{

condition = initial

}

accuracy{

...

...

jで下に，kで上にスクロールできる。qで抜けることができる。

pwd

現在のディレクトリーを出力する。

例：

$ pwd

/home/phase0

cp

一つ目の引数のファイルを，二つ目の引数のファイルにコピーする。

コピー先がディレクトリーの場合はそのディレクトリーの下に

コピー元と同じファイル名でコピーされる。

-Rオプションによってディレクトリーごと

コピーすることもできる。コピー先のファイル

が存在する場合上書きされるので，注意が必要。

mv

一つ目の引数のファイルを二つ目の引数のファイルにリネームする。

元のファイルは保持されないので注意が必要。

rm

引数に指定したファイルを削除する。やり直せないので注意が必要。

echo

引数の文字列をそのまま標準出力に出力する。

例：

$ echo "abc"

abc

command > file

コマンドの結果を，ファイルに出力する。ファイルは上書きされる。

例：lsの結果をls.txtに出力する

$ ls > ls.txt

$ less ls.txt

nfinp.data file_names.data

command >> file

コマンドの結果を，ファイルの末尾に追加する。

command1 | command2

command1の結果を，command2の入力とする。

例：実行中のmpiexecコマンドを調べる

$ ps aux |grep mpiexec

... mpiexec -n 4 phase ...

... grep mpiexec ...

tail

ファイルの末尾を表示する。

-fオプションをつけるとファイル更新を監視する

ようになり，更新される度にその内容が

表示される。-nオプション(nは

数値)を指定すると，末尾n行を出力する。

grep

ファイルの中身を検査し，引数に与えた文字列と一致する行を出力する。

例：output000からエネルギーの履歴を抽出する

$ grep TH output000

TOTAL ENERGY FOR 1 -TH ITER= -198.325431056805 ...

TOTAL ENERGY FOR 2 -TH ITER= -198.341016261819 ...

TOTAL ENERGY FOR 3 -TH ITER= -198.379630813150 ...

...

...

mpiexec

MPIジョブを開始する。

例：4プロセス，バンド・k点2並列ずつでPHASE/0を実行する

mpiexec -np 4 ~/phase0_2023.01/bin/phase ne=2 nk=2

gedit

geditエディターを起動する。

例：nfinp.dataを編集する。

gedit nfinp.data

sed

引数で与えるファイルの内容を指定のルールに

従って編集し，その結果を標準出力に出力する。

例：fileというファイルにあるabcという文字列をすべてdefに置換する

$ less file

abcdefghijklmnop

$ sed "s/abc/def/g" file

defdefghijklmnop

seq

連続した数値を出力する。

例：1から10を出力する

$ seq 1 10

1

2

3

4

5

6

7

8

9

10

bc

標準入力から得た数式を評価し，結果を出力する。

例：(1+2)×3 を計算する

$ echo "(1+2)*3" | bc –l

9

pushd/popd

pushdコマンドによってcdコマンドと同様引数に

指定したディレクトリーに移ることができる。

その後popdコマンドを実行するとpushdコマンド

を発行したディレクトリーに移ることができる。