Quantum-Espressoの空間群入力 の変更点

#author("2024-03-08T07:17:01+09:00","","")
*Quantum-Espressoのpw.x用入力ファイルで結晶構造に空間群を使う方法 (2018/03/18 created)(2022/06/30 updated)
#author("2024-03-08T07:27:54+09:00;2024-03-08T07:17:01+09:00","","")
*Quantum-Espressoのpw.x用入力ファイルで結晶構造に空間群を使う方法 (2018/03/18 created)(2024/03/08 updated)
　VESTAの出力する.vestaファイルからQuantum Espressoのpw.xで使う入力ファイルを作成するpythonプログラム(vesta2pw.py)を作って使っています。これを使うと空間群で結晶構造を指定できます。
**前書き
　pw.xの入力ファイルでは結晶構造を指定する必要があります。ずっと以前はpw.xは空間群が使えず、格子タイプ（ibrav）を指定して、座標データを入力する必要がありました。しばらく前から空間群が利用できるようになってます。空間群を利用する利点はいくつもあって、１つは入力ファイルの構造データ部分が簡単になります。構造最適化も恩恵を受けます。対称性を正しく生かすことができます（ibravで入力すると点群が正しくならないことがしばしばあります）。ただしpw.xは標準設定の空間群しか受け付けないなど使いづらい面もあります。空間群が使えるようになっても構造データ（cifファイルなど）から入力データの構造部分を作るのは面倒です。そこでPythonでVESTAのファイルから構造部分を読んで、それを使って入力データのテンプレートを作るプログラムを作ってます。QE7.0で試しています。
　pw.xの入力ファイルでは結晶構造を指定する必要があります。ずっと以前はpw.xは空間群が使えず、格子タイプ（ibrav）を指定して、座標データを入力する必要がありました。しばらく前から空間群が利用できるようになってます。空間群を利用する利点はいくつもあって、１つは入力ファイルの構造データ部分が簡単になります。構造最適化も恩恵を受けます。対称性を正しく生かすことができます（ibravで入力すると点群が正しくならないことがたまにあります）。ただしpw.xは標準設定の空間群しか受け付けないなど使いづらい面もあります。空間群が使えるようになっても構造データ（cifファイルなど）から入力データの構造部分を手入力するのは面倒です。そこでPythonでVESTAのファイルから構造部分を読んで、それを使って入力データのテンプレートを作るプログラムを作ってます。QE7.0-7.3で試しています。
**使用方法
　Pythonコードの[[[ダウンロード:https://mkanzaki.sakura.ne.jp/coes/vesta2pw.py]]]
　Pythonコードの[[[ダウンロード:https://mkanzaki.sakura.ne.jp/codes/vesta2pw.py]]]
+このコードはVestaの*.vestaファイルを構造データとして使うので、まずは結晶構造データ(cifなど)をCrystal Structure Open Databaseなどから入手します。それをVESTAで読み込みます。まず構造がちゃんと表示されているか見てください。
+次にVESTAのEditメニューからEdit Data--> Unit Cell...を選びます。そこで表示されるSpace groupとSettingに注目します。もしCustomだとうまく変換できないので、別のcifを探してください。Settingで軸の取替えなど(b,a-c)になっていて、Setting No.が1以外になってない場合はSettingで1を選んでください。なお、Origin choiceが2の場合はそのままで構いません（1に直しても構いませんが）。
+逆にRhombohedralでSetting No.が1で(Hexagonal axes)となっている場合は、Setting No.を2にしてください（なぜかhexagonalのままだとうまく計算できないし、Rhombohedralの方が原子数少なくて計算速いので）。
+Setting No.を変えた場合はOKボタンをクリックします。
+次にUtilitiesメニューからStandarization of Crystal Dataをクリックします。
+最後にFileメニューからSaveまたはSave As...で.vestaファイルをセーブします。
　次にvesta2pw.pyと同じディレクトリに*.vestaファイルを入れて、Terminalでそのディレクトリに移動して、以下を実行します。計算としてscf, relax, vc-relaxの３つから選択できます。
 % ./vesta2pw.py test.vesta scf
 % ./vesta2pw.py test.vesta relax
 % ./vesta2pw.py test.vesta vc-relax
　もしPythonのあるパスがvesta2pw.pyの最初の行と異なる時は正しく設定してください。
　なおPATHを切っておけばvesta2pw.pyはどこにあっても構いません。その場合は./は不要です。vesta2pw.pyに実行権限がついてない場合はchmodで権限を755に変える。
これで上の例だとtest.inファイルが同じディレクトリに出来ているはずです。Setting No.が１でない場合やrhombohedralでhexagonal cellを選んだ時には途中でメッセージを出してストップします。VESTA上でSetting No.を直して再度実行してください。
　なおPATHを切っておけばvesta2pw.pyはどこにあっても構いません。その場合は./は不要です。vesta2pw.pyが実行できない場合、chmodで権限を755に変えます。
　これで上の例だとtest.inファイルが同じディレクトリに出来ているはずです。Setting No.が１でない場合やrhombohedralでhexagonal cellを選んだ時には途中でメッセージを出してストップします。VESTA上でSetting No.を直して再度実行してください。
+*.inをエディタで開いて必要な変更を行った後に、pw.xを実行する。
**コードの注意点
　得られる*.inはあくまでもtemplateです。必要な計算精度等に応じて直す必要があります。擬ポテンシャル名も実際に使っているものに直さないとそのままではpw.xが動きません。普段使っている擬ポテンシャル名やその他の設定についてPythonコードの対応するところを書き換えてください。
　固溶体や占有率が100%でない構造はこのコードでは対象外です。端成分の構造がない場合で固溶体の構造を使わざるを得ない場合は、１つのサイトに１つの原子が100%占めるようにVESTA上で編集してください。
　pw.xは現時点では標準設定の空間群しか扱えません。なので軸を取り替えたりした空間群でそのまま使うと計算がおかしくなります。鉱物分野の典型例としては、forsterite（brigmaniteも）があり、歴史的経緯からほとんどの場合Pbnm空間群で記載されています。しかしこれは標準の取り方ではなく、Pnmaが標準設定です。このような場合はVESTAで標準設定に変換しておくことがpw.xを使う上で必須となります。また、構造データによっては非標準なデータを入れている場合があるので、VESTAの標準化機能を使っておくと多くの場合それらが直されますので、何か計算がおかしくなる場合は標準化機能を試してみてください（最初から使っておいた方がいいです）。
　pw.xはWyckoff symbolを読むことができるのですが、どうも完全ではないので、多くの場合はWyckoff symbolではなく単にx,y,zの座標値を与える方がいいようです。ただし例外的にtrigonalとhexagonalは逆に座標値ではおかしくなるので、その場合だけWyckoff symbolを使う様にしてます。
　pw.xでは空間群によって格子原点として２つ可能な場合がある時に、origin_choice(=1 or 2)でそれらが選べる様になっているのでそれも利用しています。なのでorigin choiceが2でも問題なく計算できるはずです。
　rhombohedral cellの場合は本来hexagonal cellとrhombohedral cellのどちらも使えますが、pw.xではどうもhexagonal cellだとうまくいきません。この場合はSettingでrhombohedral axesを選んでください。rhombohedral axesの方が単位格子中の原子数が少ないので、計算的にも有利になります。
　pw.xの出力は残念ながら元の空間群で使った座標に変換してくれません。また単純格子Pではない場合は、単純格子に変換して計算して、そのまま出力されます。しかし、VESTA上のEditメニューのEdit DataのUnit Cell...のとこのTransform...を使って元の格子に変換することが可能です。
　上記はQuantum Espressoのバージョンアップで変わる可能性があります。私が今回試しているのはversion 7.0です。
　得られる*.inはあくまでもtemplateです。必要な計算精度等に応じて直す必要があります。K-mesh（一番最後）も格子定数から適当に計算してます。擬ポテンシャル名も実際に使っているものに直さないとそのままではpw.xが動きません。普段使っている擬ポテンシャル名やその他の設定についてPythonコードの対応するところを書き換えてください。
　固溶体や占有率が100%でない構造はこのコードでは対象外です。端成分の構造がない場合で固溶体の構造を使わざるを得ない場合は、１つのサイトに１つの原子が100%占めるようにVESTA上で編集してください。この時にネットチャージがゼロになるようにします。
　pw.xは現時点では標準設定の空間群しか扱えません。なので軸を取り替えたりした空間群でそのまま使うと計算がおかしくなります。鉱物分野の典型例としては、forsterite（bridgmaniteも）があり、歴史的経緯からほとんどの場合Pbnm空間群で記載されています。しかしこれは標準の取り方ではなく、Pnmaが標準設定です。このような場合はVESTAで標準設定に変換しておくことがpw.xを使う上で必須となります。また、構造データによっては非標準なデータを入れている場合があるので、VESTAの標準化機能を使っておくと多くの場合それらが直されますので、何か計算がおかしくなる場合は標準化機能を試してみてください（最初から使っておいた方がいいです）。
　pw.xはWyckoff symbolを読むことができるのですが、どうも完全ではないようなので、多くの場合はWyckoff symbolではなく単にx,y,zの座標値を与える方がいいようです（手動でやる必要はありません）。ただし例外的にtrigonalとhexagonalは逆に座標値ではおかしくなるので、その場合だけWyckoff symbolを使う様にしてます。
　pw.xでは空間群によって格子原点として２つ可能な場合がある時に、origin_choice(=1 or 2)でそれらが選べる様になっているので、それも利用しています。なのでorigin choiceが2でも問題なく計算できるはずです。
　rhombohedral cellの場合はhexagonal cellとrhombohedral cellのどちらも使えますが、pw.xでは空間群を使うとhexagonal cellだとうまくいきません。この場合はSettingでrhombohedral axesを選んでください。rhombohedralの方が単位格子中の原子数が少ないので、計算時間的にも有利になります。
　pw.xの出力は残念ながら元の空間群で使った座標に変換までしてくれません。また単純格子Pではない場合は、単純格子に変換して計算して、その座標がそのまま出力されます。しかし、VESTA上のEditメニューのEdit DataのUnit Cell...のとこのTransform...を使って元の格子に変換することが可能です。
　上記はQuantum Espressoのバージョンアップで変わる可能性があります。私が今回試しているのはversion 7.0-7.3です。
**作成した*.inファイルがうまく動かない時
　構造データとInternational Tables Aを見ながら直すのがベストですが、International Tables Aは全ての方が持っている訳ではないので…大学だと図書室にあったり、ネットで読めるように契約している場合もあると思います。
　空間群を使うのを諦めて、ibravを使うのが最後のオプションです。もちろん欠陥の入った構造などは対称性がなくなっているので、空間群がそもそも使えません。ibravを使うためのPythonコードも作ってます[[Quantum-ESPRESSOとVestaの連携]]。内容は古いのですが、プログラム自体は最近も更新してます。
**計算結果の構造をVESTAで見る
　vc-relax, relax計算を行った場合には最適化された構造がどうなっているかVESTAで見たいと思います。pw.xの出力ファイルからVESTAで読める結晶構造を抽出するプログラムも作成しています。これはpw.xの出力ファイルから最適化された時はその構造を、最適化が終わってない場合には最後の構造をxtlファイルに変換するプログラムです。xtlファイルはVESTAで読めます。なお、pw.x出力のフォーマットがバージョンアップで時々変わるので、プログラムでエラーが出ることがありますが、最近のバージョン6.*,7.0では動くはずです。Macの場合は、最後から２、３行目の行頭の#を取るとVESTAが自動的に起動して構造を表示します。このプログラムもダウンロードできます。[[pwout2xtl.py:https://www.misasa.okayama-u.ac.jp/~masami/pwout2xtl.py]] 
　vc-relax, relax計算を行った場合には最適化された構造がどうなっているかVESTAで見たいと思います。pw.xの出力ファイルからVESTAで読める結晶構造を抽出するプログラムも作成しています。これはpw.xの出力ファイルから最適化された時はその構造を、最適化が終わってない場合には最後の構造をxtlファイルに変換するプログラムです。xtlファイルはVESTAで読めます。なお、pw.x出力のフォーマットがバージョンアップで時々変わるので、プログラムでエラーが出ることがありますが、最近のバージョン6.*,7.0-7.3では動くはずです。最後から数行目の行頭の対応するOSの#を取ると、vesta2pw.py実行時の最後にVESTAが自動的に起動して構造を表示します。このプログラムもダウンロードできます。[[pwout2xtl.py:https://mkanzaki.sakura.ne.jp/codes/pwout2xtl.py]] 
　使い方は簡単で、出力ファイル名を引数に取るだけです。*.xtlファイルが出来るので、それをVESTAで読んでください。
 > ./pwout2xtl.py test.out
　前にも書いたように出力ファイルの構造は使った空間群に対応する様には出力してくれないので、複合格子だと単純格子に変換されて計算されていて、それがそのまま出力に出てきます。元の格子に戻したい場合には、VESTAの機能を使うと戻すことが可能です。
　例えばC面心格子からpw.xで計算した後の出力は単純格子になっているはずです。それを元のC面心格子に戻すには、VESTAでEditメニューからEdit Data -> Unit Cell...を選択。Transformボタンをクリックして、出てきた画面で、Rotation matrix Pを
 1 -1  0
 1  1  0
 0  0  1
に変えて、OKボタンをクリック。ダイアログが出るが、そのままOKをクリックして、Applyボタンをクリックすると、C面心格子での構造が表示されているはず。ただVESTA上では構造はP1のままで、C面心とも認識していません。必要ならさらに使ったcif（または.vesta）のコピーを作って、そこの座標の値をこの最適化された座標を拾って、手で入れ替えを行います。こちらには単位格子内全ての座標が出てますので、その中の必要な座標だけをコピーします（多分似た値になっているのですぐ分かります）。これで最適化された結晶構造のcifまたは.vestaが得られます。pw.xでここまでやってくれるといいのですが…上記の変換マトリックスは変換前後の座標系を比べてやるとすぐ導出できるはずです。間違うと右手系から左手系に変換する等の注意が出ます。