PDielecの使用法メモ の履歴(No.15)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
粉末誘電体の赤外吸収スペクトルは粒径、粒子形状、使っている媒質によって大きく変化します。ピーク数の増加、ピークの移動、強度の変化が生じて解釈を難しくしています。PDielecは粉末誘電体の赤外吸収スペクトルを計算してくれるpythonプログラムのパッケージです（現在は結晶も対応してます）。赤外吸収スペクトルの解釈に非常に有用なプログラムだと思います。ただPDielec単体で赤外吸収スペクトル等を全て計算してくれる訳ではなく、別途DFT等プログラムでの振動計算をしておいて、その結果をPDielecが利用します。私は実験ではラマンを主に使っていて赤外はやってないのですが、ある事情で使ってみました。私の普段使っているQuantum EspressoのPHonon(ph.x)コード出力にも対応しています。しかし実際使ってみると、マニュアルには書かれていない、いくつか注意すべきことがあるようですので、それらをメモにしておきます。PDielecについての論文はここに出ています。なお、PDielecについては３年前にサマープログラムで研究室に滞在していたChris Gregsonさんに教えてもらいました（使い始めたのは最近）。
また、以下はQuantum Espressoの出力でPDielecを動かすことがメインで、PDielec自体の使い方はあまり書いてません（まだよく理解していないので）。ただ、GUIになっていて分かりやすいので、試しているうちにすぐ一応は使えるようになると思います。
PDielecのパッケージはgithubにあります。condaなど使う場合はこのパッケージは必要ないのですが、examplesなどもあるのでダウンロードしておく。https://github.com/JohnKendrick/PDielec
PDielecはpythonのプログラムなので、pythonの環境が必要です。私の場合はM3 MacBookAirを使ってますが、以前Anacondaをインストールして問題があったので、今回はminiconda3をインストールしました。condaが使えると、terminalからいくつかコマンド打ち込むことで簡単にインストールできます。ドキュメントのインストールでAnacondaのところを参照してください。pythonベースなので、（同様にpython環境を整えれば）WindowsやLinuxでももちろん動きます。
インストールが終わったら（condaの場合、conda activate pdielecをターミナルから実行後）、ターミナルにpdguiとタイプすると、GUIのPDielecのウィンドウが出てきます。そこから色々な計算が実行できて、結果のスペクトルのプロット、振動モードの表示などができます。結果をExcelファイルに書き込んでくれるので、計算結果をさらに利用できるようです。
PDielecを起動したい時はターミナルを開いて、conda activate pdielecを実行して、pdguiをタイプする。
前記の論文ではターミナルでの使用例が示されていてますが、GUI版はその後で開発されたようです。また、その論文では単結晶については全く触れられていませんが、現在のバージョンでは計算できるようになってます。上記の論文の時点ではVASP, CASTEP, GULPのみの対応でしたが、現在はそれら以外にAbinit, Crystal, Quantum Espressoの出力にも対応しています。このメモはQuantum Espressoが対象です。
Quantum EspressoのPHononコード(ph.x)の出力ファイルの１つにダイナミカル・マトリックスや各振動モードの原子変位を記録したdynファイルがありますが、PDielecではこのファイルのみが必要となります。PDielecではこのファイルを拡張子がdynGを持っているとして認識するので、dynファイルの拡張子をdynGに変えておく必要があります（またはph.xの入力ファイルでdynファイル名を*.dynGとして指定しておきます）。pdguiのメインタブでQuantum Espressoを選んでおいて、dynGファイルをファイルダイアログで選択します（拡張子.dynG以外は選択できません）。
既に計算してあったdynファイルをいくつか読ませたところ、読み取り時にエラーが出ました。そこで色々と試してみたのですが、Quantum Espressoの出力ファイルを使うためにはいくつかお約束があるらしいことが分かりました。
PHonon計算する前に、普通PWscfで構造を最適化しますが、その入力ファイルで空間群（spacegroup =）やibrav = 0以外を利用するとどうも良くないようです。実際ibrav = 0とそれ以外の時によって出来るdynファイルの違いを見てみたのですが、ファイルの最初の方のBasic vectorsのところが空間群やibrav = 0以外を使った場合にはなくて、ibarv = 0で指定した時には存在しました。どうもこのBasic vectors部分がないとエラーがでるようです。実際、空間群を使った出力にBasic vectorsのところ（それ以外にも少し変更必要で最後参照）を追加したところ、エラーが出なくなりました。なので、再計算するのが（時間がすごくかかるなどで）難しい場合、Basic vectorsを追加することでエラーを回避できます。そのためにはPWscfの出力からBasic vectors部分を取り出す必要があります。また、３行目もちょっと変更が必要です。このやり方については最後にまとめています。またpythonで自動的にやってくれるスクリプトを一応作ってみました。
PHonon計算の方ですが、入力ファイルでepsil = .true.で計算する必要があります。epsilは誘電率を計算させるフラグです。epsil = .false.でも読み込み時にはエラーは出ないようですが、誘電率テンソルや有効電荷のリストがdynGにないためにPDielecの計算ができないようです（まあそれはそうだと思いますが）。
一応上記２点に対応するとdynGの読み込みは成功して、粉末赤外吸収スペクトルの計算ができるようになりました。まだまだPDielecの使用方法を理解してないのですが、振動モードをアニメーションで見たり、それをgifやmp4で出力できました。ラマンモードは赤外スペクトルの方からは当然除外されてますが、振動モードを見ることはできます。
なお、M3 MacBookAirだと、pdguiのウィンドウ下部が画面からはみ出していて、リサイズ出来ないので、特にSettingタブの時にモードが多いと一番下が見えないことが生じました。また、スペクトル表示のところの下側のボタンが見えません。古いIntel MacBookAirも同じでした（さらにひどい）。仕方ないので、外部モニターを繋いでみたら最下部まで見られました。ウィンドウ上下のリサイズはできないようになっています。pythonなので対応するコード部分を修正すればいいと思いますが、ウィンドウサイズを設定しているコードの場所がまだ特定できていません…
さて実際に計算する場合にはibrav = 0で計算するか、またはibrav = 0以外で計算した結果をPDielecで読めるように変換する必要があります。まずは変換する方法から。後の方でibrav = 0で計算する方法を２つ紹介してます。
3行目の３つ目のコラムを0にして、４つ目のコラム（= celldm(1) = alat）はそのままで、それ以降の値をゼロにします。
そして4行目に以下を挿入。
Basis vectors
そしてその次にpw.xの出力の構造最適化されたCELL_PARAMETERS部分をコピーします。pw.xの出力ファイルを"Final"で検索するとFinal Enthalpyが検索されますが、それの7行目くらい後にcellのマトリックスがあります。その3x3マトリックス部分をコピーして、Basis vectors次の行に挿入します。その時にalatもメモしておいて、3行目の４つ目のコラムをその値で置き換えます。それでファイル自体をセーブして、PDielecからオープンします。
以下はdiasporeでの計算例。最初は空間群を使った入力ファイルから出来たdynGのファイル最初の部分。このファイルはPDielecで読み取り時にエラーが出ます。
Dynamical matrix file
phonon of diaspore pbesol
3 16 8 17.7689343 0.3023712 0.4634348 0.0000000 0.0000000 0.0000000
1 'Al ' 24592.588567557767
2 'O ' 14582.196445495396
3 'H ' 918.73579705352574
以下略
次は手直しして、Basis vectorsなどを追加したもの。
Dynamical matrix file
phonon of diaspore pbesol
3 16 0 17.7689343 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
Basis vectors
0.998080364 0.000000000 0.000000000
0.000000000 0.301790745 0.000000000
0.000000000 0.000000000 0.462545137
1 'Al ' 24592.588567557767
2 'O ' 14582.196445495396
3 'H ' 918.73579705352574
以下略
この編集したファイルはエラーを出さずに読み込めました。MainタブのUnit cellのところで正しく読めていることを確認します。なお上記リストの格子定数部分（最初の17.7689343はa軸に対応）は原子単位(au)で、PDielecの出力はオングストローム単位で表示されているので値が同じにはなりません。またその次の２つはb, cなのですが、b/a, c/aで示されています。最後の３つは軸間の角度で、こちらはcos(angle)となってます。
上記マニュアルでやっていたことを自動化したpythonコード(dynG.py)を使ってdynGを作る方法です。この場合ph計算では普通に*.dynを出力させます。以下のリンクでdynG.pyをダウンロードしてください。
download
最初の行のpython3の場所を正しく設定します。なおalatとマトリックス部分をpwの出力から取るので、そのファイルも必要です。なので以下のように２つのファイルを引数として取ります（順番重要）。
% dynG.py test-pw.out test-ph.dyn
これをターミナルで実行すると、test-ph.dynGが作成されます。これをPDielecから読みます。
ただまだ問題があって、面心構造などで格子定数がおかしくなることを確認しています。１つの回避策は構造をVestaからxtlで出力して、それをxtl2pw.pyで出力ファイルに変換してQuantum Espressoの計算を進めることです。
私の作っているxtl2pw.pyはxtlファイルからpw.x用の入力ファイルを作成しますが、セルの対称性からibravを適当に設定する仕様のためにibrav = 0にはなりません。それでそのままpw.xで計算するとPDielecの計算には使えない出力ファイルとなってしまいます。ibrav = 0に対応させるにはpw.xの入力ファイルを少し編集する必要があります。
まずは入力ファイルを開いて、ibrav = 0とします。そうなるとcelldmについてはcelldm(1)以外は不要なので、celldm(1)以外がある場合は削除します。また、Atomic Speciesの前にCELL_PARAMETERS (alatと3x3のcellのマトリックス3行を追加する必要があります。この部分はpw.xを実行した時に出力で出てくる部分と同じフォーマットですので、一度xtl2pw.pyで作った入力ファイルでpw.xを実行して、構造最適化ルーチンが１度終わったら、計算をkillして、出力部分を見ると以下のような部分があるはずです（具体的な数値はもちろん異なりますが）。この部分をコピーして、入力ファイルのAtomic Speciesの前に挿入してやります。
CELL_PARAMETERS (alat= 20.00842021) 1.004344960 0.000000000 -0.000388049 0.000000000 1.343906566 0.000000000 -0.233130619 0.000000000 0.929018807
ibrav = 0の入力ファイルを作るもう１つの方法は、まずはvesta2pw.pyかxtl2pw.pyを作ってvc-relaxの入力ファイルを作成します。構造最適化が最低１回は回った後で計算を止めて、それらの出力を使って、pwout2in.pyで入力ファイルを再作成します。pwout2in.pyで作られるpw.xの入力ファイルは常にibrav = 0になるので、これを使ってpw.xの計算を行います。それらのコードはQE用pythonコードにあります。