PDielecの使用法メモ の履歴(No.7)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
PDielecはDFT計算結果等から赤外吸収スペクトル等を計算するpythonプログラムのパッケージです。単体で赤外吸収を計算出来る訳ではなく、DFT計算等がまず必要です。PDielecは微小な試料が媒質に埋まっている時の赤外吸収の周波数変化や強度変化を計算することができます。私は実験ではラマンを主に使っていて赤外はやってないのですが、ある事情で使い始めています。私の普段使っているQuantum EspressoのPHonon(ph.x)コード出力にも対応しています。しかし実際使ってみると、マニュアルには書かれていない、いくつか注意すべきことがあるようですので、それらをメモにしておきます。PDielecについての論文がここに出ています。なお、PDielecについては３年前にサマープログラムで研究室に滞在していたChris Gregsonさんに教えてもらいました。
パッケージはgithubにあります。https://github.com/JohnKendrick/PDielec
pythonの環境が必要です。私の場合はM3 MacBookAirを使ってますが、以前Anacondaインストールして問題があったので、今回はminiconda3をインストールしました。condaが使えるとterminalからいくつかコマンド打ち込むことでインストールできます。ドキュメントのインストールでAnacondaのところを参照してください。pythonベースなので、WindowsやLinuxでも動きます。
インストールが終わったら、ターミナルからpdguiとタイプすると、GUIのウィンドウが出てきます。そこから計算が実行できて、結果のスペクトルのプロット、振動モードの表示など全てできます。ただ、前記の論文では端末での使用例が示されていて、GUI版はその後で開発されたようです。
上記の論文の時点ではVASP, CASTEP, GULPのみの対応でしたが、現在はそれら以外にAbinit, Crystal, Quantum Espressoの出力にも対応しています。このメモはQuantum Espressoが対象です。
Quantum EspressoのPHononコードの出力ファイルの１つにダイナミカル・マトリックスや各振動モードの原子変位を記録したdynファイルがありますが、PDielecではこのファイルのみが必要となります。PDielecではこのファイルを拡張子がdynGを持っているとして認識しているので、dynファイルの拡張子をdynGに変えておく必要があります（またはPHononの入力ファイルでdynファイル名を.dynGとして指定しておきます）。pdguiのメインタブでQuantum Espressoを選んでおいて、dynGファイルをファイルダイアログで選択します（拡張子.dynG以外は選択できません）。
既に計算してあったdynファイルをいくつか読ませたところ、読み取り時にエラーが出ました。そこで色々と試してみたのですが、Quantum Espressoの出力ファイルを使うためにはいくつかお約束があるらしいことが分かりました。
PHonon計算する前に、普通PWscfで構造を最適化しますが、その入力ファイルで空間群（spacegroup =）やibrav = 0以外を利用するとどうも良くないようです。
実際ibrav = 0とそれ以外の時によって出来るdynファイルの違いを見てみたのですが、ファイルの最初の方のBasic vectorsのところが空間群を使った場合にはなくて、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も同じでさらに悪いかも…外部モニターを繋いでみたら最下部まで見られました。ウィンドウ上下のリサイズはできないようになってます。
私の作っているxtl2pw.pyではセルの対称性からibravを定義するためにibrav = 0にはなりません。それでpw.xで計算するとPDielecの計算には使えない出力ファイルになってしまいます。ibrav = 0に対応させるにはpw.xの入力ファイルを少し編集する必要があります。
まずはibrav = 0とします。そうすると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コードにあります。
3行目の３つ目のコラムを0にして、４つ目のコラム（= celldm(1)）はそのままで、それ以降の値をゼロにします。
そして4行目に以下を挿入。
Basis vectors
そしてその次にpw.xの出力の構造最適化されたCELL_PARAMETERS部分をコピーします。pw.xの出力ファイルを"Final"で検索するとFinal Enthalpyが検索されますが、それの7行目くらい後にcellのマトリックスがあります。その3x3マトリックス部分をコピーして、Basis vectors次の行に挿入します。それでファイル自体をセーブして、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のところで正しく読めていることを確認します。