Quantum-ESPRESSOで振動モードを計算する時のメモ (2013/10/26 created) (更新 2017/07/09)

更新情報

  • (2017/07/09) 振動方向をVESTAで表示させるコードのVibratz版を作成中。Vibratzは力定数を入力すると、振動モードを計算してくれる有料のアプリケーション。振動モードの表示が一応できるが、簡易な表示なので、VESTAでも表示できるようにしている。
  • (2017/06/24) 振動方向をVESTAで表示されるために、Vectorを手入力するのは大変なので、PythonでVESTAファイルを作成するコードを作った。下にその情報を追加した。
  • (2017/06/07) 以下の方法で計算したstishoviteの振動モードの計算例を追加した。

計算の流れ

  • ラマンテンソルを計算しないなら使う擬ポテンシャルに制約はないが、ラマンテンソルを計算したい場合はノルム保存(NC)擬ポテンシャルのみが使える。まずpw.xで構造最適化計算(構造がすでに最適化されているならscf計算)を行い、続いてph.xで格子振動モード等の計算を行う。ph.xの入力ファイルでは、prefixで指定する名前はpw.xの入力ファイルと同じにする。この時にqをどの範囲で計算するかを指定する。さらに必要に応じて、dynmat.x, matdyn.xなどを使う。
  • ph.xの計算を並列で行う時は、事前のpw.x計算と同じ並列のオプションを使わないとエラーとなる。
    > mpirun -np 6 pw.x < test.scf.in > test.scf.out &
    > mpirun -np 6 ph.x < test.ph.in > test.ph.out &

基準モードの計算(ラマン、赤外モード)

  • ラマン、IRの振動モードの計算が必要なだけなら、q=0だけの計算でいいので、epsil = .true.にする。フォノンの分散曲線を計算する場合は、qの範囲を指定する(ldisp = .true.でnq1,nq2,nq3で指定)。ph.xの出力には最初にq=0(gamma点)の結果がある(分散を計算させた場合、qの異なる結果が出力にはあるので注意)。
  • ph.xの計算では3つの並進運動の振動数(低周波数に出ている)がゼロになっていない。気にしないならそれでもよいが、きになるなら、dynmat.xでasr = 'simple'で計算すると並進運動がとゼロなるように再計算してくれる。TO-LO分裂も処理する。また、IRのcross sectionの計算も行われる。ラマン、IRの計算をする時はq=0のph.xの出力ファイル(.dyn1)を入力ファイルfildynとして指定する。
  • pw.x(ph.xも同様)では入力構造から点群を推定しているが、なかなか正しく認識させるのは難しい場合がある。なるべくQE標準のセルの取り方で、pw.xの最初の方の出力に認識された対称性が出てきて、点群も印刷されているので、出力を見て正しい点群になっているか確認する必要がある。なお、点群まで詳しく出力させるには、&controlのところで、verbosity = 'high'とする。
  • ph.xの出力では各振動モードについて点群の指標に基づいて、A1gなどと各ピークをassignしてくれるので便利。しかしdynmat.xでasrを補正すると、IRモードについては振動数が変わり、場合によってはph.xの出力リストの順番と入れ替わることがある。dynmat.xの出力ではモードがassignされていないので、対応は簡単でない。最終的には各モードの変位をチェックして、対応させるしかない? なお最適化を丁寧にやっておけば並進運動モードが最初から0に近いはずなので、ズレは大きくないので、補正しなくてもいいかも。ラマン活性モードについては振動数にあまり変化がないので(補正される並進運動はラマン活性モードと異なる指標を持つため)、ph.xの出力と直接対応できる。
  • ラマン散乱テンソルの計算はdynmat.xで行うが、ノルム保存(NC)擬ポテンシャルにしか対応していない。上記のように振動数だけが必要な場合は、この計算は必要ない。
  • GGAの場合、計算した振動数は実測よりも低波数側にずれる傾向。これは体積自体が少し過大評価されるため。私の経験では4-5%程度。必要なら圧力を上げて、実測体積に近いように調整するのも1つの手。
  • pbesol paw potentialで計算したstishovite(SiO2の高圧相でルチル構造)の振動モードの計算結果を下に示す。Ramanの実測値より数%低い周波数になっているが、よく再現されている。各振動モードの原子変位は*.0.dyn1の最後の方に出ているので、この図ではそれを使って、原子変位方向に矢印を付けている。この表示にはVestaのVector機能を使った。なお、stishoviteは圧力により変位型相転移をして、CaCl2構造になるが、その時のソフトモードはB1gが対応する。八面体がab面内で傾く振動であり、これに対応してセルがtetragonalからorthorhombicになる。この辺りは実験と計算でよく調べられている。下記は授業のデモンストレーション用に計算した。この場合、指標(B1gなど)の最後にgがつくのがRaman active, uがつくのがIR activeモードである。
    vibration_mode_stishovite.png

VESTAでの矢印表示用変換コード(ph2vesta.py)

  • 上記のstishoviteの場合はVESTA上で1個毎データを入力したが、振動モードが多いと大変なので、ph.xの出力から自動的に作成するプログラムをPythonで作成してみた。まだ、問題点が多いが、公開しておく(download)。Macのterminalでのみテストした。MacのOSデフォルトのPythonを使っている。
  • stishoviteの計算例を添付したので、とりあえずそのディレクトリーに移動して、
    >./ph2vesta.py stishovite
    とタイプしてみる。コードが実行されない時は実行権限の問題なので、chmod 755 ph2vesta.pyとタイプする。
  • 実行するとモードのリストが出るので、見たいモードを選んで、その番号を入力するとそのモードに対応するVESTAのファイルが作成され、VESTAに表示される。
  • arrowの色、長さ、太さは-1を入力すると変更できる(現在1種のみ)。なお長さのところをマイナスにすると矢印の向きを逆にできる。短時間で作ったのでエラー処理などはいい加減です。なお、stishoviteデータの10.2 cm-1以下の3つのモードは並進運動のモードとなります。
  • とりあえず、100以上モードがある計算データを1個づつVESTAに入力する必要はなくなった。

問題対応

  • FFT gridとの関係で対称性が正しく設定されない場合は(pw.xの出力でそう書かれている時)、ecutrhoの数値をいじるか、FFT gridを直接指定する(nr1,nr2,nr3)。
  • pw.xでは問題ないが、ph.xでset_irr_sym_newのエラーが出た時は、以前の関係する計算ファイルを削除するとエラーが出なくなった。理由は不明。
  • pw.xで入力構造の詳しい点群の対称性を出力させるにはverbosity = 'high' にして実行する。開始すぐに出力ファイルの先頭を眺めて、もし点群が間違えられているなら、すぐ計算を止めて、入力ファイルを直して、計算し直す。
  • どうしても本来の対称性よりも低くなってしまう場合は、とりあえず計算を行い、本来の点群と計算された点群で指標を対応させる。BilbaoサイトのRaman and Hyper-Raman scatteringの"Correlations Points"で対応関係を計算してくれる。
  • うまくいかないケース:monoclinicのbase center cellだと、そのままだと点群は正しいが、モード数が(原子数が2倍になるため)多くなる。かといってprimitive cellに直すと、triclinic cellとなり、今度は点群が正しくなくなる。
  • *.0.dyn1の出力は一部の情報が欠けているので、正しい格子定数と原子位置がちゃんと計算できない。点群の指標もこちらには出てこない。ph.xの標準の出力を同時に参照する必要がある。

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2017-06-26 (月) 17:33:11 (150d)