9. PHASE/0プログラム、ツールの実行方法

9.1. プログラムphase

  • プログラムphaseの実行

PHASEはSCF計算、分子動力学法計算を行います。また収束した電荷密度分布から状態密度やバンド分散を計算することができます。

入力パラメータファイル,擬ポテンシャルファイルを実行ディレクトリに置きます。file_names.dataを使用する場合には、それも同じディレクトリに置いてください。

1プロセッサ(1コア)の逐次計算を行う場合には、次のようにプログラムphase を実行します。”../../phase0_2020.01/bin/”は、PHASEがインストールされているディレクトリです。

% ../../phase0_2020.01/bin/phase

並列計算を行う場合には、お使いの計算機の利用するMPIライブラリの実行コマンドを使用します。詳細はお使いの計算機システムのマニュアルを参照ください。一般的なコマンドはmpirunです。

% mpirun -np NP ../../phase0_2020.01/bin/phase ne=NE nk=NK (2次元版)

% mpirun -np NP ../../phase0_2020.01/bin/phase.3d ne=NE nk=NK ng=NG (3次元版)

ここで,NP はMPI プロセス数,NE はバンド並列数,NK はk点並列数です。3次元版の場合のNGはG点並列数です。

9.1.1. プログラムphase2次元版の並列計算オプション

9.1.1.1. バンド並列、k点並列

並列計算(バンド並列、k点並列)では、バンド並列数NE、k点並列数NKを指定します。 NP = NE×NK という関係が成立している必要があります。

% mpirun -np NP ../../phase0_2020.01/bin/phase ne=NE nk=NK

通常,バンド並列よりもk点並列の方が効率が良いです。 したがって,可能な場合はk 点並列数を大きくすると良いと考えられます。 ただし,Brillouin領域内にサンプルするk点数は系が大きくなるほど少なくても充分になり,そのk点数が必ずしも利用できるプロセス数で割り切れるわけではない(MPIプロセス数がサンプリングk点よりも大きくなる)という点に注意が必要です。 k点数よりもNK の値が大きいとエラーになります。 また,k点数がNK で割り切れない場合は理想的な並列効率が得られません。 そこで,必要に応じてバンド並列も組み合わせて計算を実行してください。

ne, nk という引数は省略することも可能です。その場合のデフォルト値は下記の通り。

  • バージョン2020.01未満:NE = NP, NK = 1

  • バージョン2020.01以上:対称性を考慮した上で得られるk点数と総並列数NPが割り切れる最大の整数値がNK, NEはNP/NK.

9.1.1.2. レプリカ並列

NEB 法,拘束条件付きダイナミクス,メタダイナミクスなどの機能によっては“レプリカ並列” が利用できる場合があります。 レプリカ並列を実行するには以下のコマンドを利用します。

% mpirun -np NP ../../phase0_2020.01/bin/phase nr=NR ne=NE nk=NK

NR はレプリカ並列数です。NP = NR×NE×NK という関係が成立している必要があります。 レプリカ並列の効率はk 点並列よりも更に良いですが,k 点並列と同じ注意が必要です。 また,“もっとも収束の遅いレプリカ”が律速となるので,実効的には必ずしも効率的とは限りません。

9.1.2. プログラムphase 3次元並列(G点並列)版

PHASE/0 はバンドとk 点の2 軸並列に対応していますが,平面波のG 成分の並列化にも対応しています。 3軸並列版も,2軸並列版と同様インストーラーによってコンパイルすることができます。

% ./install_3d.sh

実行は、次の例のように行います。

% mpirun -np NP ../../phase0_2020.01/bin/phase.3d ng=NG ne=NE nk=NK

ここで,NPは総MPIプロセス数,NG はG 点並列数,NE はバンド並列数,NKはk点並列数を意味します。 NG とNE とNK の積は,総MPI プロセスの数に等しい必要があります。 2次元版の場合と同様,NEB法,拘束条件付きダイナミクス,メタダイナミクスなどの機能を利用する場合はnr=NRによってレプリカ並列をすることも可能です。

ne, nk, ngが省略された場合のデフォルト値は下記の通り。

  • バージョン2020.01未満:NE = NP, NK = 1, NG=1

  • バージョン2020.01以上:対称性を考慮した上で得られるk点数と総並列数NPが割り切れる最大の整数値がNK, NGとNEはNE*NG=NP/NKを満たし,かつNE:NGが1:2に最も近くなる取り方

9.1.3. 3次元版のオプションについて

3次元版は,以下のような設定を施すことによって,おもに低並列時に高速になる場合があります。

FFTを非並列で処理する

3次元版はFFTを並列で処理します。これをあえて非並列で処理することによって,G点並列数が少ない場合に高速になる場合があります。このオプションは以下の要領で有効にすることができます。

control{
  sw_serial_fft = on
}

特にng=1とする場合は高速化が期待できるオプションです。

電荷密度の処理用にコミュニケーターを追加する

3次元版は,電荷密度を波動関数と同じコミュニケーターで扱います。以下の設定を施すことによって電荷密度に専用のコミュニケーターを割り当てることができます。

control{
  sw_communicator_for_chg = on
}

G並列数が少ない場合に,おもに交換相関相互作用の処理が高速化されます。

9.2. プログラムekcal

状態密度計算、バンド計算において、k点の個数が多い場合に使うプログラムとしてekcalがあります。

SCF計算の計算結果の電荷密度を入力として計算できます。

SCF計算の計算結果の電荷密度ファイルnfchgt.dataを実行ディレクトリにコピーします。または、入出力ファイル設定ファイルfile_names.dataにおいて、F_CHGにSCF計算で得られた電荷密度ファイル指定します。

バンド構造計算においては、サンプリングk点の設定ファイルkpoint.dataを用意します。

次のようにプログラムekcal を実行します。”phase0_2020.01/bin/”は、PHASEがインストールされているディレクトリです。

% ../../phase0_2020.01/bin/ekcal

ekcalプログラムは2次元版にしか用意されていませんが、入力パラメーターファイルに以下のような設定を加えることによってphaseをekcalと同じように動作させることができます。

control{
  fixed_charge_option{
    kparallel = one_by_one
  }
}

ONE_BY_ONEモードで動作する場合,あるk点iterationにおける初期波動関数は1つ前のk点iterationの波動関数が採用されます。この際,以下のように変数sw_modified_kpoint_incrementの値をonとすることによってk点の更新方法が変更され,k点並列時により近いk点の波動関数が初期波動関数として採用されるようになります。

control{
  fixed_charge_option{
    kparallel = one_by_one
    sw_modified_kpoint_increment = on
  }
}

9.3. プログラムepsmain

電子系の誘電関数の計算に利用するプログラムがepsmainです。その動作はekcalとほぼ同じですが、入力パラメーターの設定に応じて誘電関数計算用の処理が行われる点が異なります。3次元版のバイナリー名はepsmain.3dです。

9.4. プログラムtdlrmain

線形応答時間依存密度汎関数法による励起スペクトル計算に利用するプログラムがtdlrmainです。3次元版は用意されていません。

9.5. 状態密度図作成ツール dos.pl

状態密度図の作成

PHASEあるいはEKCALによって状態密度データを出力させることが出来ます。 それについては、たとえば3.10.1をご覧ください。 その状態密度データdos.dataを可視化するプログラムがdos.pl exampleのdos.dataをworkにコピーします。

$ cp ../example/dos.data

間違いなくdos.dataがコピーされていることをlsで確認してください。 を用いて,この状態密度データdos.dataを可視化してみましょう。次のようにコマンドを入力して下さい。

$ dos.pl dos.data -erange

こうすると,EPSファイルdensity_of_states.epsが生成されます。 UNIX環境で,これを見るにはghostviewやgvなどを必要とします。

$ ghostview
または $ gv
../_images/ch8_10_image1.svg

図 9.1 バルクSiの状態密度図

dos.plを実行するときに状態密度データdos.dataの後に付加した-erangeは表示するエネルギーの範囲を 制御するオプション,-colorはカラー出力を行うためのオプションです。

9.5.1. dos.plのオプション

なにも付加せずにdos.plを実行すると利用方法が表示されます。

$ dos.pl

Version: 3.00
Usage: dos.pl DosD
ata -erange=Emin,Emax -einc=dE -dosrange=DOSmin,DOSmax -dosinc=dDOS
-title=STRING -with_fer
mi -width=SIZE -font=SIZE -color -mode={total|layer|atom|projected}
-epsf={yes|no} -data={yes|no}

DosDataに状態密度データが記録されたファイル(通常dos.data)を指定します。 その後に作図を制御するオプションを指定します。

-erange=Emin,Emax

表示 するエネルギーの範囲をeV単位で指定する。

たとえ ば,-10 eVから5 eVまで表示したい場合は, -erange=-10,5 とします。指定をしないと,データの 最小値・最大値から自動的に決定されます。

-einc=dE

目盛りの間隔を指定する。 た とえば,2eV間隔に目盛りをふりたいなら,   -einc=2 とします。

-dosrange=DOSmin,DOSmax

表示する状態密度の範囲を変える。 たとえば,0 states/ eVから12 states/eVまで表示したい場合は, -dosrange=0,12 とします。

-dosinc=dDOS

縦軸(状態密度)の目盛りの間隔を指定する。 たとえば, 2 states/eV間隔に目盛りをふりたいなら, -dosinc=2 とします。

-title=STRING

グ ラフにタイトルを付けたいときに設定する。 たとえ ば,タイトルを「Total DOS」とするなら, -titile="Total DOS" とします。

-with_fermi

デフォルトでは描か ないフェルミレベルまたは価電子帯上端のエ ネルギーレベルを描く。金属ではフェルミレ ベルを表示し,絶縁体・半導体であれば価電 子帯上端のエネルギーレベルを表示します。

-width=SIZE

図の幅のデフォルト 値は1であるが,その値を変更したい場合は このオプション を使う。たとえば,0.8に変更したい場合は -width=0.8 とします。

-font=SIZE

フォントのサイズを変更したいと きには,これを設定する。既定値は14です。 たと えば,フォントサイズを28にしたいならば, -font=28 とします。

-color

グラフをカラー表示します。

-mode={total|layer|atom}

totalを 指定すると,全状態密度図が作成されます。 layerを指定すると, 層分割の局所状態密度の図が作成されます。 atomを指定すると,原 子分割の局所状態密度の図が作成されます。 projectedを指定すると,原子軌 道分割の局所状態密度の図が作成されます。 既定値はtotalです。

-epsf={yes|no}

ポストスクリプトファイルを作成しない ときには,noを指定します。既定値はyesで 指定がなければ, ポストスクリプトファイルが作成されます。

-data={yes|no}

層分割や原子分割の状態密度データからe psファイルを作成するのではなく個別のファ イルに出力するときにはyesを指定します。

9.6. 改良版状態密度図作成ツールdos.py (バージョン2020.01以降)

9.6.1. 概要

PHASE/0には状態密度可視化スクリプトdos.plが付属します。このスクリプトを用いると,全状態密度・原子分割局所状態密度・層分割局所状態密度・射影状態密度の状態密度図をEPS形式で得ることができます。dos.pyスクリプトはその改良版です。スクリプト名の拡張子から示唆されるように,Pythonスクリプトです。下記のような機能が搭載されています。

  • dos.plが持っている全機能

  • EPS以外の画像ファイルの対応

  • 状態密度を加算する機能;たとえば,原子分割局在状態密度において指定の原子群の状態密度を足し上げて状態密度図を作成したりデータファイルを出力したりすることができる機能

  • 層分割局所状態密度のヒートマップ作成機能

  • 動作モードの追加:dos.plが提供するバッチモードのほか,対話モードとGUIモード

層分割状態密度のヒートマップは,x座標が層の座標,y座標がエネルギー,z座標が状態密度という三次元データをヒートマッププロットすることによって得られます。このような可視化を行うことによって,層ごとに状態密度がどのように変化するかを視覚的に明らかにすることができるようになります。

状態密度を加算する機能などを利用する場合,すべてのオプションを引数で渡すのは煩雑な場合がありえます。そこで,dos.pyは対話的にも利用できるようになっています。たとえば加算機能を利用する場合,加算対象となりえる状態密度データがリストアップされるので,所望の状態密度データをそこから選択します。

dos.pyは簡易的なGUIも提供します。利用するGUIのフレームワークはtkinter (https://docs.python.org/ja/3/library/tkinter.html) です。tkinterは通常Pythonに標準的に備わっているので,利用の際特別な準備は必要ありません。

状態密度プロットの際,dos.plはgnuplotを用います。これに対し,dos.pyはmatplotlib (https://matplotlib.org/) をプロットの際に用いる仕様となっているため,動作にはPythonのほかmatplotlibが必要となります。また,numpyも必要です。 Matplotlibやnumpyはpip(https://www.python.jp/install/windows/pip.html)などの仕組みを用いてインストールすることが可能です。 また,両方ともAnaconda (https://www.anaconda.com/) のようなPython distributionにはプリインストールされています。 なお,Python3以降が必須であり,Python2系列では動作しません。

9.6.2. 使い方

9.6.2.1. バッチモード

バッチモードにおいては,-fもしくは--fileオプションによって状態密度データファイルを指定し,さらに様々なオプションを指定します。利用可能なオプションは下記の通り。なお,オプションは原則としてハイフン(-)が一つの場合は空白の後に,ハイフンが二つ続く場合(--)は=の後に値を指定します。 たとえば-m total, --mode=total など。 また,真偽値を指定するタイプのオプションの場合そのオプションがあるかどうかで真偽が判定されるため,値は指定しません。 オプションは1. スクリプト全般の振る舞いを制御するオプション,2. 状態密度を加算する場合に加算対象を選択するオプション,3. 描画に関するオプションの三種類があります。

スクリプト全般の振る舞いを制御するオプション

スクリプトの全体的な振る舞いを以下のオプションによって制御することができます。

オプション

説明

--version

バージョンを表示する。

-h, --help

ヘルプを表示する。

-i, --interactive

対話モードで 実行する場合に指定するオプション。

-g, --gui

簡易GUIモードで 実行する場合に指定するオプション。

-f FILE --file=FILE

状態密度データファイル を指定する。デフォルト値はdos.data.

--output000=OUTPUT000

output000ファイルを指 定する。無指定の場合,作業ディレク トリーにおいてタイムスタンプが最も 若いoutput000ファイルが採用される。

-m MODE --mode=MODE

状態 密度データの種類を指定する。total, atom, layer, projectedのいずれか。

-a ACTION --action=ACTION

スクリ プトの振る舞いを指定する。analyze, split, sumのいずれか。

analyze の場合は状態密度データを解析し, 結果を出力する。splitの場合は局所状 態密度などを分割し,画像ファイルな どを作成する。sumの場合フィルターオ プションの指定に応じて状態密度を加 算して画像ファイルなどを作成する。 カンマ区切りで複数指定してもよい。

--heatmap

層分割局所状態 密度の場合に,ヒートマップを作成し たい場合このオプションを指定する。

-o OUTPUT,

--output_action=OUTPUT_ACTION

出力の振る舞いを指定する。genfig, storedata, bothのいずれか。genfigの場合 画像ファイルが出力される。storedat aの場合テキストファイルに加工した状 態密度データが出力される。bothの場 合両方行われる。デフォルト値はboth.

加算対象の状態密度データを選択するためのオプション

状態密度を足し上げてその結果を画像ファイルに出力したりテキストファイルに出力したりすることができます。この時に加算対象とする状態密度データをオプションによって選択します。カンマ(,)とハイフン(-)を用いて複数の整数値を選択することができます。カンマ区切りで値を一つずつ,ハイフン区切りで連続する値を選択することができます。たとえば1,3,4,8-11 などとすると1,3,4,8,9,10,11と展開されます。ただし--elemidは文字列の指定なので,カンマ区切りのみ利用可能です。

オプション

説明

--dosid=DOSID

dosidによっ て加算対象の状態密度を選択する。

dosidは,得に射 影状態密度の場合は分かりづらいので後述のatomid, lid, mid, tidを利用してもよい。

--atomid=ATOMID

加 算対象とする原子のIDを指定する。原子分割局所状 態密度および射影状態密度の場合のみ意味をもつ。

--layerid=LAYERID

加算対象とする層のIDを指定 する。層分割局所状態密度の場合のみ意味を持つ。

--elemid=ELEMID

加算対象とする元素を指定する。

--lid=LAYERID

加算対象とする方位量子数 を指定する。射影状態密度の場合のみ意味を持つ。

--mid=MID

加算対象とする磁気量子数 を指定する。射影状態密度の場合のみ意味を持つ。

--tid=TID

加算対象とする主量子数 を指定する。射影状態密度の場合のみ意味を持つ。

描画オプション

状態密度図を作成する際にどのように描画するかをオプションによって制御することができます。

オプション

説明

-e ERANGE, --erange=ERANGE

態密度図描画の際のエネルギーの範囲を 指定する。emin,emaxという形式で指定す る。eminが下限値,emaxが上限値。emin, emaxは片方の み指定することもできる。すなわちemin, もしくは,emax. この場合指定されなかった方の値はmatp lotlibのデフォルト値が割り当てられる。

--einc=EINC

エネルギーの軸の目盛り値を指定する。

-d DRANGE,

--drange=DRANGE

状態密度描画の際の状態密度の範囲を指 定する。指定の形式はエネルギーと同じ。

--dinc=DINC

状態密度の軸の目盛り値を指定する。

--lrange=LRANGE

層分割局所状態密度の ヒートマップ作成の場合に,層の範囲を指 定する。指定の形式はエネルギーと同じ。

--linc=LINC

層分割局所状態密度のヒートマップ作成 の場合に,層の軸の目盛り値を指定する。

--with_fermi

フェルミエネルギーを表す点線を描画し たい場合にこのオプションを有効にする。

--title

状態密度図にタイトルを表示 したい場合このオプションを有効にする。

--cmap=CMAP

層分割局所状態密度のヒートマップ作成の 場合に,カラーマップの種類を指定する。

可能な選択肢が※のウェブサイトに 掲載されている。デフォルト値はviridis

--imgtype=IMGTYPE

画像ファ イルの種類を指定する。以下のいずれか。

eps, ps, png, jpg, pdf, svg

https://matplotlib.org/3.1.1/gallery/color/colormap_reference.html

9.6.2.2. 対話モード

dos.pyを-i もしくは--interactiveをつけて実行すると対話モードで利用することができます。対話モードでは選択肢がいろいろと提示されるので,所望の振る舞いに応じて選択します。対話モードでは,おおよそ以下のように動作します。

  • 状態密度の種類を選ぶ。状態密度の種類とは,total, atom, layer, projectedのいずれか。

  • atom, layer, projectedの場合何を実行するかを選ぶ。atom, projectedの場合splitもしくはsum, layerの場合はこれにheatmapが加わる。

  • sumの場合,加算対象のIDを選ぶ。

  • 描画オプションの選択:erange, drange, with_fermiを入力する。

  • 画像ファイルの種類:eps, ps, png, jpg, pdf, svgのいずれかを選ぶ

9.6.2.3. GUIモード

dos.pyに-g もしくは--guiオプションをつけて実行すると 図 9.2 で示すようなGUIが得られます。

../_images/ch8_10_image2.png

図 9.2 dos.pyが提供するGUI.

図のアルファベットを用いてこのGUIについて説明します。

  1. プロット表示域。

  2. 状態密度の種類などを選択するタブ域。totalは全状態密度,atomは原子分割局所状態密度,atom (sum)は原子分割局所状態密度加算,layerは層分割局所状態密度,layer (sum)は層分割局所状態密度加算,layer (heatmap)は層分割局所状態密度のヒートマップ,projectedは射影状態密度,projected (sum)は射影状態密度加算,figure optionsは描画オプション設定域。

  3. 状態密度選択域。atom ID, layer ID, pdos IDなどのチェックボックスを選択するとそれに応じて描画が更新される。加算モードの場合複数選択することが可能で,選択状態のすべての状態密度が加算され描画される。

  4. プロット操作域。プロットの描画のされ方などを変えることができる。フロッピーディスクのアイコンをクリックするとファイル選択ダイアログが得られ,画像ファイルにエクスポートすることができる。

9.6.3. 出力

出力としては,分割もしくは加工した状態密度データが記録されたテキストファイルと状態密度図のデータが得られます。二つのファイルのファイル名は拡張子以外共通であり,拡張子は前者がdata, 後者は--imgtypeの指定に応じたそれになります。拡張子を除いた分は以下に示すように場合によって異なります。

  • --action=sumでなく,--heatmapもつけていない場合

  • 全状態密度:dos_total

  • 局所状態密度:dos_atomatomid

  • 層分割局所状態密度:dos_layerlayerid

  • 射影状態密度:dos_atomatomid_llid_mmid_ttid

  • --heatmapを付けている場合

  • 層分割局所状態密度:layer_dos_heatmap

  • --action=sumの場合

  • 局所状態密度:dos_summed_atom

  • 層分割局所状態密度:dos_summed_layer

  • 射影状態密度:dos_summed_projected

9.6.4. 利用例

使用例と結果得られる状態密度図を紹介します。利用する状態密度データはPHASE/0のサンプルにあるBaO-Si界面とBaTiOの状態密度データです。

../_images/ch8_10_BaOSi_BTO.svg

図 9.3 (a) BaO-Si界面と(b) BaTiO の原子配置。

9.6.4.1. BaO-Si界面

dos.py -f dos.data --imgtype=svg

全状態密度の状態密度図が得られる。画像ファイルはsvg形式。

../_images/ch8_10_image5.png

図 9.4 全状態密度。

dos.py -f dos.data -m atom --imgtype=svg

原子に分割した状態密度図が得られる。原子数分の画像ファイルが得られる。一番目の原子(最下層のBa)と二番目の原子(最下層のO)の状態密度を図示する。

../_images/ch8_10_aldos_Ba_O.svg

図 9.5 1番目の原子と2番目の原子の局所状態密度。いずれも 図 9.3 (a)の最下層の原子で,元素は1番目がBa, 2番目がOである。

dos.py -f dos.data -m layer --imgtype=svg

層に分割した状態密度図が得られる。層の数分の画像ファイルが得られる。最下層と中央付近の層の状態密度を図示する。

../_images/ch8_10_layerdos_Ba_Si.svg

図 9.6 1層目( 図 9.3 の最下層)と65層目( 図 9.3 (a)の中央付近)の状態密度

dos.py -f dos.data -m atom -a sum --elemid=Si --imgtype=svg

Si原子の状態密度をすべて足しこんだ状態密度図が得られる。

../_images/ch8_10_image10.png

図 9.7 Si原子の状態密度をすべて足しこんだ状態密度。

dos.py -f dos.data -m layer --heatmap --imgtype=png --drange=,2

層分割状態密度のヒートマップが得られる。色の違いを際立たせるため,状態密度の上限を2としている。

../_images/ch8_10_image12.png

図 9.8 層分割状態密度のヒートマップ。横軸が最低面を原点とした場合の層の原点からの距離,縦軸がエネルギーであり,状態密度の値の違いは描画色の違いによって表されている。

9.6.4.2. BaTiO結晶

dos.py -f dos.data -m projected --imgtype=svg

射影状態密度の図が軌道ごとに出力される。Tiのd軌道の状態密度を図示する。

../_images/ch8_10_BTO_Ti3d.svg

図 9.9 Ti d軌道の射影状態密度。

dos.py -f dos.data -a sum --lid=2 -m projected --imgtype=svg

Tiのd軌道の状態密度をすべて足し上げた状態密度図が得られる。

../_images/ch8_10_image15.png

図 9.10 Ti d軌道の射影状態密度。

9.7. k点ファイル生成ツール band_kpoint.pl

バンド構造図を描くには,対称線に沿った\(\mathbf{k}\)点の列を生成し,その各\(\mathbf{k}\)点での固有エネルギーを ekcalで計算します。ekcalは\(\mathbf{k}\)点のデータが書き込まれたファイルkpoint.dataを 読み込み各\(\mathbf{k}\)点での固有エネルギーを計算します。その\(\mathbf{k}\)点のファイルの生成を支援するプログラムが band_kpoint.plです。band_kpoint.plの入力ファイルの記述形式は以下の様になっています。

dkv
b1x b2x b3x
b1y b2y b3y
b1z b2z b3z
n1 n2 n3 nd # Symbol
...

dkvが\(\mathbf{k}\)点の間隔,b1x,b1y,b1zは逆格子ベクトル\(\mathbf{b}_{\mathbf{1}}\)のx,y,z成分。 逆格子ベクトル\(\mathbf{b}_{\mathbf{2}}\),\(\mathbf{b}_{\mathbf{3}}\)についても同様です。 五行目以降に特殊\(\mathbf{k}\)点とそのシンボルの指定をします。シンボルの指定は必須ではありませんが,指定がある場合,バンド構造図作成の際に利用されます。 整数\(n_{1},n_{2},n_{3},n_{d}\)を用いて\(\mathbf{k}\)ベクトルを

\[\mathbf{k =}\frac{n_{\mathbf{1}}}{n_{\mathbf{d}}}\mathbf{b}_{\mathbf{1}}\mathbf{+}\frac{n_{\mathbf{2}}}{n_{\mathbf{d}}}\mathbf{b}_{\mathbf{2}}\mathbf{+}\frac{n_{\mathbf{3}}}{n_{\mathbf{d}}}\mathbf{b}_{\mathbf{3}}\]

のように指定します。シンボルは#の後に書いてください。面心立方格子の場合の例を示します。

0.02                        <---- k点の間隔
-1.0  1.0  1.0
1.0 -1.0  1.0               <---- 逆格子ベクトル
1.0  1.0 -1.0
0 1 1 2 # X                 <---- n1 n2 n3 nd # Symbol
0 0 0 1 # {/Symbol G}
1 1 1 2 # L
5 2 5 8 # U
1 0 1 2 # X

これと同じものがディレクトリexampleにあるので,それをコピーしてband_kpoint.plを実行してみましょう。

$ cd PHASE_INST_DIR/samples/tools/work
$ cp ../example/bandkpt_fcc_xglux.in .
$ band_kpoint.pl bandkpt_fcc_xglux.in > output

こうするとkpoint.dataが生成されます。これがバンド構造計算用の\(\mathbf{k}\)点のファイルです。 この\(\mathbf{k}\)点のファイルを入力に加えて,ekcalで\(\mathbf{k}\)点での固有エネルギーを計算してください。

9.8. バンド構造図作成ツール band.pl

9.8.1. band.plの実行

band.pl PHASE/0のekcalの出力nfenergy.dataとband_kpoint.plの入力ファイルがband.plの 入力になります。 前節の入力例で生成したkpoint.dataを入力とし、ekcalで固有エネルギー計算を行い、 結果得られた固有エネルギーファイルnfenergy.dataが ディレクトリexampleにあります。 このファイルを使ってバンド構造図を描いてみましょう。exampleにあるnfenergy.dataと bandkpt_fcc_xglux.inをworkにコピーし、それらを入力としてband.plを実行します。
$ cp ../example/nfenergy.data
$ cp
$ band.pl nfenergy.data bandkpt_fcc_xglux.in

こうすると,EPSファイルband_structure.epsが生成されます。 このファイルをご覧になるには,ghostviewやgvなどのソフトウェアが必要です。

$ ghostview
または $ gv
$ gv band_structure.eps
../_images/ch8_10_image17.svg

図 9.11 バルクSiのバンド構造図

band.plを実行するときにいくつかのオプションを付加することができます(次節)。

9.8.2. band.plのオプション

なにも付加せずにband.plを実行すると利用方法が表示されます。

$ band.pl

 Usage: band.pl EnergyDataFile KpointFile -erange=Emin,Emax
 -einc=dE -ptype={solid_circles|lines} -with_fermi
 -width=SIZE -color

KpointFileの後が作図を制御するオプションです。

-erange=Emin,Emax

表示するエネルギーの範囲をeV単位で指定する。

たとえば,-10 eVから5 eVまで表示したい場合は, -erange=-10,5 とします。

-einc=dE

目盛りの間隔を指定する。 たとえば2eV間隔に目盛りをふりたいなら, -einc=2 とします。

-ptype=TYPE

描画種を選択する。 -ptype=solid_circles  : 黒く塗りつぶされた円で表示する。 -ptype=lines : 直線でつなぐ(デフォルト)。

-with_fermi

デフォルトでは描かないフェルミレベルまた は価電子帯上端のエネルギーレベルを描く。金属で はフェルミレベルを表示し,絶縁体・半導体であれ ば価電子帯上端のエネルギーレベルを表示します。

-width=SIZE

図の幅のデフォルト 値は0.5であるが,その値を変更したい場合はこのオ プションを使う。たとえば,0.3に変更したい場合は -width=0.3 とします。

-color

グラフをカラー表示する。

9.9. 原子構造の拡張trajectory形式への変換ツール dynm2tr2.pl

Perlスクリプトdynm2tr2.plは、構造最適化、分子動力学法計算のデータ(nfdym.data)を拡張trajectory形式に変換します。

ツールdynm2tr2.plを以下のように実行します。

$ dynm2tr2.pl nfdynm.data

このようにすると,dynm.tr2というファイルとgrid.mol2というファイルが生成されます。前者は原子の座標などが記述されたファイルであり, 後者は対応するセルの情報などが記述されたファイルです。

FCCのプリミティブセルにSiが二原子入った非平衡状態を初期構造とし,構造最適化した結果を拡張trajectory形式に変化し、可視化した例です。

../_images/ch8_10_image18.svg

図 9.12 バルクSiの構造最適化過程の可視化例

図 9.12 の矢印は原子に作用する力を表しています。 力が極大になったあとは,原子座標の更新が進む毎に原子に作用する力が小さくなり,原子構造が最適化されていく様子が分かります。 図 9.12 ではプリミティブセルで表示されますが, 以下のようなcontrol.inpというファイルを作成すれば,原点の移動やセルの変更ができます。

origin  1.2825 1.2825 1.2825
vector1 10.26  0     0
vector2  0    10.26  0
vector3  0     0    10.26

このcontrol.inpを使用してdynm2tr2.plでdynm.tr2を作成すると原点が(1.2825,1.2825,1.2825) bohrに移り, セルのベクトルが(10.26,0,0),(0,10.26,0),(0,0,10.26) bohrになります。 以下のようにして,dynm.tr2を作成します。

$ dynm2tr2.pl nfdynm.data control.inp

ブラベーセルで構造最適化過程のstep 10を図示したのが,図 9.13 です。

../_images/ch8_10_image19.svg

図 9.13 ブラベーセルで表したバルクSiの構造最適化過程(step 10)

9.10. 振動数レベル図作成ツール freq.pl

PHASEの振動解析機能を使用すると,結晶の基準振動モードの振動数と固有ベクトルが得られます。 振動解析の結果は、ファイルmode.dataに出力されます。Perlスクリプトfreq.plmode.dataから振動数のデータを取り出し 振動数レベル図を作成します。freq.pl実行すると,EPS形式の画像ファイルfreq.epsが出力されます。されます。

$ freq.pl [options] mode.data

バルクSiの振動数解析結果の振動数のレベル図を 図 9.14 に示します。

../_images/ch8_10_image20.svg

図 9.14 バルクSiの振動数レベル図

振動数レベルを表す横棒は既約表現ごとに列にまとめて分類され,その各列の上には既約表現の名称と活性を表す記号(IR,R,IR&R,NON)が 表示されます。IRは赤外活性を表し,Rはラマン活性を表します。IR&Rは赤外活性とラマン活性があることを示します。 NONはサイレントモードであることを示しています。作成された振動数レベル図では,横線の右側には振動数がcm-1単位で表示されます。 既約表現ごとに振動数の低い順に番号付けされ,横線の左側に表示されます。

9.10.1. freq.plのオプション

なにも付加せずにfreq.plを実行すると利用方法が表示されます。

$ freq.pl

*** A visualization program for vibrational freqencies ***
Usage: freq.pl [-width=W] [-height=H] [-nrep=N] {-solid|-mol|-ignored_modes=LIST} mode.data

freq.plのオプションです。

-width=W

図の幅のデフォルト値は1であるが,そ の値を変更したい場合はこのオプションを使う。

たとえば,0.3に変更したい場合は -width=0.3 とします。

-height=H

図の幅のデフォルト値は1であるが,そ の値を変更したい場合はこのオプションを使う。 たとえば,2.5に変更したい場合は -height=2.5 とします。

-nrep=N

一つの図に表示する既約表 現の数。振動モードの既約表現かこの数よりも多 いときには,複数のEPSファイルが作成されます。

-solid

固体の場合に並進を非表示にするオプション。 これはデフォルトで設定されています。

-mol

分子 の場合に回転と並進を非表示にするオプション。

-ignored_modes=LIST

LISTのところにコンマで区切っ て並べた番号のモードは表示されなくなります。 たとえば, -ignored_modes=1,2,3 と すると1,2,3番のモードは表示されなくなります。

9.11. 基準振動の軌跡の拡張trajectory形式ファイル変換ツール animate.pl

Perlスクリプトanimate.plmode.dataに出力されている振動モードの固有ベクトルのデータを読み込み,基準振動の軌跡を拡張trajectory形式ファイルに変換します。

control.inpというファイルを用意すると、原点の移動とセルベクトルの変更ができます。

control.inpの例です。

origin  1.27189 1.27189 1.27189
vector1 10.17512 0 0
vector2 0 10.17512 0
vector3 0 0 10.17512

この例では、ブラベーセルで表示するために,原点を(1.27189,1.27189,1.27189) bohrに移し,セルベクトルを(10.17512,0,0),(0,10.17512,0),(0,0,10.17512) bohr に変更します。

animate.plを以下のように実行します。

$ animate.pl mode.data control.inp

各振動モードごとの拡張trajectory形式ファイルmode_1.tr2,mode_2.tr2,...,mode_6.tr2というファイルとgrid.mol2というファイルが作成されます。 拡張trajectory形式のファイルは振動モードの数だけ出力されます。

バルクSiの振動解析の6番目の基準振動の固有ベクトルmode_6.tr2を可視化した図を 図 9.15 に示します。

../_images/ch8_10_image21.png

図 9.15 バルクSiの基準振動の固有ベクトル

9.12. 座標データ変換ツールconv.py

conv.pyというPythonスクリプトを使って,座標データを変換することができます。binディレクトリーの下にあります。conv.pyは対話的に利用します。たとえば,nfdynm.dataファイルをCIFに変換する手続きは下記の通りです。

画面に現れる文字列

説明

$ conv.py

atomic configuration converter utility.

Copyright (C) the RISS project, The University of Tokyo

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.da taファイルのファイル名を指定。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

9. LAMMPS_input x. 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倍にしたスーパーセルを作成し,変換します

PHASE/0の入力を変換する場合,ホームディレクトリーに.piouというファイルが存在し,その中身が下記のようになっている必要があります。

pp.default_pp_dir=PATH_TO_PPDIR

ここでPATH_TO_PPDIRは擬ポテンシャルファイルが納められたディレクトリーです。

9.13. 入力ファイル正誤チェックツールinpcheck.py

inpcheck.pyとは,PHASE/0の入力の正誤チェックを行うツールです。binディレクトリーの下にあります。

このスクリプトを利用するには,まずホームディレクトリーに.piouというファイルを作成し,その中身を以下のようにする必要があります。

pp.default_pp_dir=PATH_TO_PPDIR

ここでPATH_TO_PPDIRは擬ポテンシャルファイルが納められたディレクトリーです。

以下のように計算の実行ディレクトリーにおいて実行すると診断をはじめ,結果が標準エラー出力に出力されます。たとえば,以下のような出力が得られます。

$ inpcheck.py

input data validator utility for PHASE
Copyright (C) the RISS project, The University of Tokyo
INFO: -- running the input validator --
INFO: specfile : phase.spec
INFO: checking directory : /home/jkoga/phase0_2018.01/samples/Si8
INFO: validating input...
INFO:
INFO: checking if required/recommended entries exist...
INFO:
INFO: found [accuracy.matrix_diagon.cutoff_wf] at line 31, a recommended
entry when [accuracy.initial_wavefunctions] is 'matrix_diagon'
...
...
INFO: checking whether each entires are valid...
INFO:
INFO: line 2: entry [control.cpumax] matches the specification.
INFO: line 3: entry [control.condition] matches the
...
...
WARNING:
WARNING: line 95: could not find [postprocessing.dos.nwd_window_width] in specification.
WARNING: candidate entry: [postprocessing.dos.nwd_dos_window_width]
WARNING:
INFO: line 98: entry [postprocessing.charge.sw_charge_rspace] matches the specification.
INFO: line 99: entry [postprocessing.charge.filetype] matches the specification.
INFO: line 100: entry [postprocessing.charge.title] matches the specification.
INFO:
INFO: input validation done.
INFO:
WARNING: found 2 warnings in /home/jkoga/phase0_2018.01/samples/Si8
WARNING: check the log for details.

inpcheck.pyに渡せるオプションは下記の通り。

オプション

説明

-l LOGLEVEL, --loglevel=LOGLEVEL

ログレベルを指定する。

0 : 最小出力 1 : デフォルト出力 2 : デバッグ出力

--prop=PROPFILE

プロパティ―ファイルをユーザ ー指定のものにする。(デフォルト 値はpiou/data/props.properties)

--ppdir=PPDIR

擬ポテンシャル格 納ディレクトリーをデフォルト以外 のものにする(デフォルト値は$HOM E/.piouファイルに記述されている)

-s SPECFILE, --specfile=SPECFILE

入力仕様定義ファイルを ユーザー指定のものにする(デフォ ルト値はpiou/config/phase.spec)

-r, --recursive

起動したディレ クトリーの下のすべてのサブディレ クトリ―を再帰的にたどり,PHASE/0 計算フォルダーと判定された場合に 入力ファイル正誤チェックを行う。