Quantum-ESPRESSOとVestaの連携(2015/01/07作成)

  • xtl2pw.pyのupdate(2016 01/01):xtlファイルで単位格子外の原子が含まれていても対応するようにしました。psuedopotentialを原子毎に予め指定しておくことができるようにもしました。できたファイルに手直しをあまりせずに動かせるようになってます。
  • xtl2pw.pyのupdate(2015 11月):単斜晶系で通常セッティング(betaが非90度)に対応しました。また、C面心格子を単純格子に変換したxtlファイルで、ibravが正しくならないバグをとりました。
  • (2017/06/07)京大の吉川様の指摘により、xtl2pw.pyにバグがあることが分かりました。直したものと入れ替えました。

前書き

  • Quantum-ESPRESSO(以下QE)は擬ポテンシャルと平面波を使う第一原理計算のプログラムパッケージである。我々のところで最近よく使っている。一方、Vestaは結晶構造や電子密度を表示するプログラムである。QEにおける計算の入力ファイルには結晶構造データが必要であるが、QEでは空間群を利用してないので、単純格子だと単位格子中の全原子の座標を記述する必要がある。私は最初Vesta上の結晶構造から、Fractional Coordinatesをexportして、これをカット&ペーストして、QEの入力ファイルを作っていた。大した手間でもなかったが、Pythonでこの過程を自動化するスクリプトを作ったので公開する。また、QEのpw.xで構造最適化した出力ファイル(.out)から、Vestaで構造を表示できるように.xtlに変換するスクリプトも公開する。これらによって、QEとVestaが相互に連携できるようになり、QEの入力ファイル作成が簡単にでき、QEの計算後に最適化した構造をVestaでチェックすることができる。スクリプトはMacOS XとLinux上でテストした。

QEのpwscf用入力ファイルを作成する

  • pwscf(実行ファイル名pw.x)は1点電子状態計算(scf)、構造最適化計算(vc-relax, relax)などができるQEの最も基本的なプログラムである。これ用の入力ファイルを自動作成する。そのために必要なのはVestaからExportされたxtl(fractional coordinate)ファイルである。以下にQE入力用のxtlファイルの作り方を示す。
    • まずVesta上で正常に結晶構造が表示されていることを確認。第一原理計算では部分的に占有された系の計算はそのままではできないので,占有率は全て1であること。
    • もし空間群の取り方が特殊な場合は標準の設定に直しておく。これはVestaのEditメニューから,Edit Dataで,Unit Cell..のところでできる。Settingを1にする。Applyボタンをクリック,OKボタンでダイアログを抜ける。
    • Vesta上で単位格子内の原子が全て表示されるように、ObjectsメニューのBoundary...で、x(min)=0 ~ x(max)=1.0のように設定する。FileメニューのExport Data...で、File typeとして、Fractional Coordinatesを選択。適当な名前で出力する。
    • なお、複合格子の場合は単純格子に変換することで、原子数を減らせ、計算が速くなる。単純格子への変更はVestaの機能を使えばできる。
  • このxtlファイルを使って、pw.xの入力ファイルを作るPythonコードがこちらで、以下のようにターミナルから使う。なおコードは実行可能にしないといけないので、
    > chmod 744 xtl2pw.py
    をターミナル上で実行しておく。自分以外のユーザーが実行する場合があるなら755にする。
    > ./xtl2pw.py test.xtl vc-relax 1
    このコードには3つ引数があり、最初は.xtlファイル名(.xtlが必須)、2番目は計算方法(scf, relax, vc-relaxのいずれか)、3番目が原子位置に制約をつけるかどうかのフラグで1が制約を自動でつける、0が制約をつけない。これを実行すると、test.inファイルが作成される。引数を忘れた時のために、引数なく実行した時には例が出るようにしている。
    • なお、scfは与えられた構造での1点電子状態の計算を行い、relaxは原子位置の最適化(格子はそのまま)、vc-relaxは原子位置と格子の最適化計算を行う。
    • 制約とは構造最適化時に原子位置を固定するフラグのことで、座標x,y,zの後に0または1をつけて指定できる(x,y,zに応じて3つ、空白でわける)。何もないと1と解釈される。0は原子座標固定で、1または未記入はその原子座標を最適化する。このコードでは0.00, 0.50など特殊位置と考えられる座標の場合に自動的に0フラグを書き込むことを行う。偶然そういう位置になった場合もあるので、その場合は1に変えないといけない。
  • 生成した*.inファイルはあくまでもテンプレートで、エディタで色々と直す必要がある。特に擬ポテンシャルについては私の使う適当なファイル名が入っており、実際に使うものとは異なっているはずなので、それらを変更する。使うファイル名がほぼ決まっている場合は、Pythonコードの最初の方にあるpsを実情に応じて変更すればよい。実際には元素毎にファイル名が異なることも多いので,元素毎にファイル名を指定できるようにps[i]を定義しているので,そこでファイル名を書き込んでおけばよい。また、擬ポテンシャルのディレクトリパス、設定圧力、エネルギー等のカットオフ、K-pointなども実行したい計算によって異なるはずなので、それらも変更する。最適化計算の方法も現在bfgsをデフォルトで使っているが、必要に応じで変更する。K-pointは与えられた格子常数から適当に設定しているので、必要に応じて変える。また格子(ibrav)についても単純格子を仮定しているので、そうでない場合はibravの指定を変更するか、単純格子に変換する。
  • 注意点:
    • 以前、単斜晶系の場合にgammaを非90度角に取る必要があって、事前にVestaで変換する必要があり面倒でしたが、最新のスクリプトでは不要です。ibrav=-12とすることで結晶学の標準(betaを非90度にとる)に対応できました。
    • 複合格子の場合、必要ならば単純格子に変換する。原子数が多い場合には原子数を減らして、計算速度を上げられるので変換した方がよい。これもVesta上のEdit DataのUnit cellでcustamizeすればよい。supercellにする場合も同様。supercellでは問題が起きることがある。その場合、一度supercellとしてxtlで出力しておいて、そのxtlを再度読み込んで作業するとよい。
    • 例えばC面心格子の場合、単純格子への変換は以下のようにする。
      (-0.5 0.5 0) (0.5) 
      (0.5 0.5 0) (0.5) 
      ( 0  0  1) (0.0) 
    • ibravで複合格子を指定した時は、複合格子分は原子座標として与える必要はない。しかしC面心の場合では私はうまくできなかった。上記のように単純格子に直して対応している。または単純格子として全て処理する(計算量は増える)。
    • pwscfの入力ファイルのキーワードやオプションの詳細はDocディレクトリー内のINPUT_PW.htmlをウェブブラウザーで見ればよい。
    • 環境に合わせてコードをカスタマイズしていただければ、編集なしか、最小の編集で実行可能なpw.xの入力ファイルができます。

pw.xによる構造最適化結果のVestaでの表示

  • こちらは逆にVestaで読めるxtlファイルをpw.xの出力から作るスクリプトです。まあ、最適化された原子座標を元のxtlファイルをコピー&ペーストすれば済むような話ですが、格子常数を変換する部分が結構面倒です(特にmonoclic, triclinic)。このスクリプトを使うとそこを自動でやってくれます。出力ファイルの拡張子が.outであることを仮定しています。また、構造最適化が正常に終わっている必要があります(.outファイル中に「Final」があるかどうかで判断している)。
  • pw.xの出力ファイルから構造部分を取り出してxtlファイルを作る。Pythonコードはこちら。以下のようにターミナルから使う。実行可能にしておく。
    > ./pwout2xtl.py test.out
    このコードには1つ引数があり、pw.x計算の出力ファイル名(.outが必須)を与える。これを実行すると、test.xtlファイルが作成される。このファイルはVestaで読むことができる。既存の同名xtlファイルは警告なしに上書きされるので、注意が必要。

pw.x用のユーティリティーコード (以下Vestaとは無関係)

pw.xの構造最適化結果から新しい.inファイルを作成

  • pw.x用の入力ファイルを構造最適化した結果(セルと座標)で更新する。Pythonコードはこちら。条件を変えて計算するような場合に、前の計算での原子位置と格子を新しい入力ファイルに取り込んで計算を行う。前の計算の*.inファイルと*.outファイルが必要。
  • 以下のようにターミナルから使う。実行可能にしておく。
    > ./pwout2in.py test
    このコードには1つ引数があり、pw.x計算のベースファイル名(.outは不要)を与える。この場合、test.outおよびtest.inが存在している必要がある。これを実行すると、test.inの原子位置と格子の情報がtest.outの構造最適化の結果で置き代えられたtest.new.inが作られる。それをエディタで条件などを変えて、pw.xで実行することができる。
  • 構造最適化が途中で終わっているファイルの最後の構造から続けて計算したい場合は、その最後の構造部分のまえにFinalの文字を追記して、実行すれば入力ファイルが更新される。
  • 場合によっては生じるファイル名を直接指定したい場合があります。コードの最初のファイル名を指定しているところを変更すれば可能です。

.inファイルのパラメーター値の変更

  • pw.x用の入力ファイル(.in)中の1つのパラメーター(press, conv_thrなど)の数値を変更するだけのコード。条件を変えて一連の計算を実行したい場合などに使用する。sedで出来るでしょうが… Pythonコードはこちら。。
  • 以下のようにターミナルから使う。実行可能にしておく。
    > ./pwinupdate.py test.0kbar.in test.10kbar.in press 10.0
    このコードには4つ引数があり、最初は既存の.inファイル名、2番目は新たに作成する.inファイル名、3番目はパラメーター名、4番目はその値。上記例では圧力(press)を10 kbarに設定している(Pwscfで使っている圧力の単位はkbar)。
  • 1つのパラメータしか1度に変更できない。複数変更したい場合は繰り返し実行する。 普通よく変更するパラメータとしては、press(圧力), ecutwfc(平面波のエネルギーcutoff), ecutrho(電子密度のcutoff), conv_thr(scfの収束条件), press_conv_thr(圧力の収束条件)などがある。上記のプログラム等をshellか、Pythonのcodeで使ってやると、圧力を変えた一連の計算などが自動で行えるようになる。

計算中の計算状況チェック方法

  • topコマンドで、実行中か調べる。さらにgrepを使えば、最適化計算の出力ファイル中の体積や原子間力の変化を見て、収束状況がわかる。出力ファイル中にFinalがあると、最適化計算がすでに終了したか、終了前の最後のscf計算をおこなっているところだと分かる。
    > grep Ang test.out
    > grep "Total force" test.out
    > grep Final test.out

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2016-01-01 (金) 14:38:18 (693d)