QEでの振動モード計算メモ の履歴(No.1)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
mkanzaki.sakura.ne.jp/codes
(2023/12/17) ラマン強度計算にforsteriteの例を追加。これはノルム保存のpbesol擬ポテンシャルで計算したもの。
(2023/12/8) 敢えて（プログラムをちょっとイジって）pbe(GGA)を使う場合、ノルム保存でpbe, pbesolのポテンシャルがpslibraryにはありません。Pseudo Dojo(http://www.pseudo-dojo.org/)にはあるので、計算したい方はそこからDLしたらいいようです。ノルム保存(ONCVPSP v0.4.1)＋pbesol+standardで少し計算してみましたが、pw-mt_fhiポテンシャル(LDA)と大きな違いはないようでした。注意して使う必要がありますが。
(2023/11/24) ph.x結果からgif動画を作るph2gif.pyを公開。それに応じて内容更新。
(2023/11/16) ph2vesta.pyを更新。オプションを追加。「a」を入力すると全モードのVESTAファイルを一度に作成するようにした。この場合、VESTAは開かない。ファイル等更新。
(2023/11/14) ph2vesta.pyを更新。今回から予め結晶構造のVESTAファイルを用意しておいて（pwout2xtl.pyでpw.xの出力から作れます）、それに変位の矢印を追加するように変更した。なので新しく３番目の引数としてVESTAファイル名が必要です。予めVESTAファイルを最終的に見たいようにイジってから処理するので、こちらの方が便利だと思います。それに応じてph2vesta.zipも内容を更新。
(2023/04/02) 「計算と実測スペクトルが合わない時」を追加。
(2023/02/10) enstatite多形の低周波数振動モード動画を最後に追加。
(2022/12/12) またph2vesta.pyにバグ発見。E, Tの縮重したモードがある時にそれらをまとめて示す場合と個別に示す場合があって（対称心があるかどうかで違うようだ）、それにうまく対応してなかったので場合によってエラーが出る。改訂版をアップロードした。
基準振動モードの振動数だけが知りたい場合にはph.xで(lraman = .false.)とすると、使用する擬ポテンシャルに制約はありません。IR強度も計算できます。しかし、ラマン強度（ラマンテンソル、偏光解消度など）を計算したい場合は使用する擬ポテンシャルに制約があります。最後の方をご覧ください。振動数計算の手順は、まず普通のscfまたは構造最適化計算をpw.xを使って実行する。原子位置等が安定位置からずれていると振動数がおかしくなることがあるので、構造最適化計算(vc-relax)を行った方がよいが、色々な事情で観測されている構造をそのまま使ってscf計算する場合もある。できるだけ収束条件などは精度が高くなるようにする。続いてph.xで格子振動モード等の計算を行う。ph.xの入力ファイルでは、prefixで指定する名前、ディレクトリ等はpw.xの入力ファイルと同じにする。この時にqをどの範囲で計算するかを指定する（これは分散曲線書かせたい場合で、Raman, IRではq=0の計算をする）。さらに必要に応じて、dynmat.x, matdyn.xなどを実行する。振動数だけでいいならdynmat.xを必ずしも使う必要はないが、ASR(acoustic sum rule)を使う場合やIRやRamanの散乱断面積計算には実行する必要がある。
なお、Quantum EspressoのPHonon/examplesディレクトリにph.xの計算例があるので、ph.xの入力ファイルの書き方などはそれらを参考にする。PHonon/DOCディレクトリにはマニュアル等が置いてあるので参考にする。
pw.x、ph.xの計算を並列で行う時は、同じ並列(コア/CPU)数のオプションを使わないとエラーとなる。dynmat.xも同様。dynmat.xはすぐ終わる計算なので並列の必要は全くないが、そうしないとエラーが出る。
> mpirun -np 8 pw.x < test.scf.in > test.scf.out &
なお、xtlファイルを作るためにpw.xの後に以下を実行する。
> ./pwout2xtl.py test.scf.out
これでtest.scf.out.xtlができるので、それをVESTAで読み込み、VESTA形式(test.scf.out.vesta)でセーブする。
> mpirun -np 8 ph.x < test.ph.in > test.ph.out & > mpirun -np 8 dynmat.x < test.dm.in > test.dm.out &
振動モードの振動数が必要なだけなら、q=0だけの計算でいいので、ldisp = .false.にする。デフォルトではfalseなので指定しなくてもよい。フォノン分散曲線を計算する場合は、qの範囲を指定する(ldisp = .true.でnq1,nq2,nq3で指定)。分散曲線計算させた場合にはph.xの出力の最初にq=0(gamma点)の結果がある（qの異なる結果が出力には沢山出ているのでそれをラマン、IRの振動数と勘違いしないように）。
ph.xの計算出力では３つの並進運動の振動数(低周波数に出ている３つ)が厳密にゼロにはなっていない。気にする必要もないが、気になるなら、dynmat.xでasr = 'crystal'で計算すると並進運動がゼロなるように再計算してくれる。TO-LO分裂も処理する（その場合はepsil=.true.にする）。また、IRの散乱断面積の計算も行われる（これはlraman = .false.でも）。
pw.x(ph.xも同様)では入力構造から点群を推定しているが、なかなか正しく認識させるのは難しい場合がある。なるべくQuantum Espresso標準のセルの取り方を使って、pw.xの出力の最初の方に認識された対称性が出てきて、推定された点群も出力されているので、出力を見て正しい点群になっているか確認する必要がある。なお、pw.xで点群まで詳しく出力させるには、&controlのところで、verbosity = 'high'とする必要がある。この問題については、空間群で構造を指定するとほぼ解決できそうである。空間群で構造を指定する方法はこちらに書いた--> QuantumESPRESSO_空間群入力。私はCrystal Open Database (COD)から取ってきたcifファイルをVESTAで表示させて、必要なら空間群を標準な設定に変更して、構造をVESTA形式で出力。それをpw.xの入力ファイルに変換するPythonコードを使って変換。大体はそれでpw.x, ph.xがうまく計算ができている。しかし、空間群を使った場合でおかしくなることもたまにあるので、その場合はibravでセルを指定するしかない。
ph.xの出力では各振動モードについて点群の指標に基づいて、A1gなどと各ピークをassignしてくれる。しかしdynmat.xでasrを補正すると、IRモードについては振動数がかなり変わり、場合によってはph.xの出力リストの順番と入れ替わることがある。dynmat.xの出力ではモードがassignされていないので、対応は簡単でない。最終的には各モードの変位をチェックして、対応させるしかない？ なお最適化を丁寧にやっておけば並進運動モードが最初から0に近いはずなので、ズレは大きくないので、補正しなくてもいいかも。ラマン活性モードについては振動数にあまり変化がないので（補正される並進運動はラマン活性モードと異なる指標を持つため）、ph.xの出力と直接対応できる。なお、asrは結晶の場合は並進だけだが、分子だと回転にも対応する必要があって、それに応じてasr = ''のところを変える必要がある。結晶の場合は、'simple'か'crystal'を使う。
なおassignしてくれる指標であるが、ほとんどの場合正しいが、まれにちょっとおかしい場合がある。BilbaoのSAMでも指標を求めて同じになるかチェックした方がいいかもしれない。
GGAの場合(pbe-*.UPF)、計算した振動数は実測よりも低波数側にずれる傾向。これは体積自体が少し過大評価されるため。私の経験では4-5%程度。pbeより、pbesolの方が差が小さくなる。必要なら圧力を上げて、実測体積に近いように調整するのも1つの手。
pbesol-PAW (*.pbesol-n-kjpaw_psl.UPF)で計算したstishovite(SiO2の高圧相でルチル構造)の振動モードの計算結果（表）を下に示す。ラマン実測値より数％低い周波数になっているが、まずまず再現されている。実測値は私の実測データであるが、文献値と1 cm-1以内で合っている。各振動モードの原子変位は*.dyn（このファイル名はfildynで指定できる）の最後の方に出ているので、それを使って原子変位方向に矢印を付けた図を下に示した。この表示（矢印）にはVESTAのVector機能を使った（下記参照）。
| モード | 実測シフト | 計算シフト |
| B1g | 231 cm-1 | 217 cm-1 |
| Eg | 589 cm-1 | 571 cm-1 |
| A1g | 754 cm-1 | 733 cm-1 |
| B2g | 966 cm-1 | 927 cm-1 |
なお、stishoviteは圧力によって変位型の相転移をして、CaCl2構造になるが、その時のソフトモードにはB1gが対応する。下の図でも分かるように、B1gは八面体がab面内で回転する振動であり（繋がっている八面体は反対方向に回転）、これに対応してセルがtetragonalからorthorhombicになる。つまりこのモードの振動方向が転移する時に原子が移動する方向と一致している。そのためこのモードの復元力は転移点に近づくにつれ弱くなり、ソフト化する（ソフトモード）。この辺りは実験と計算でよく調べられている。下記は授業のデモンストレーション用に計算した。この場合、A1g, B1g, B2g, EgがRaman active、A2u, EuがIR active、A2gとB1uはどちらもinactive。圧力を変えて構造最適化計算するとこの転移を計算で再現することができる。原子数も少なくて計算が早いので都合よく、私はよくデモや演習で使う。この時に重要なことは、stishovite構造のまま圧力を上げても、酸素が鞍点(saddle point)にあるので、自発的にCaCl2構造に転移しないことである。なぜなら鞍点では原子に働く力がゼロとなるためで、鞍点から少しずれると力がゼロでなくなる。なので酸素位置(x or/and y)を最初に少しずらせておく必要がある。
振動モード計算をした場合、各モードで原子がどう変位するのかを知りたいところである。そのためにVESTAのVector表示機能を使って、変位方向に矢印を表示させてみた。上記のstishoviteの場合はVESTA上で1個毎データ（変位成分）を手入力したが、振動モードが多いと大変なので、ph.xの出力から自動的に作成するプログラムをPythonで作成してみた。なお、振動に伴う変位はph.xが出力するファイルに書かれているので、それを読んで使っている。このpythonコード（例もまとめたzipファイル）はdownloadできます。Macのterminalで普段テストしています。現在はanacondaでインストールしたpythonを使ってますが、特に特別なモジュールはこのコードでは使ってません。うまく動かない場合は最初の行で実際のPythonのパスを設定してやる必要があります。zipファイルには計算例が２つ付属していて、stishoviteとquartzである。なお、１つのモードの出力を例としてつけている。さらにdynmat.xの出力も含まれている。
例えばstishoviteの計算例でテストする場合は、（３つのファイルを指定にするように2023/11/14現在なっています）
>./ph2vesta.py stishovite.ph.out stishovite.dyn stishovite.scf.out.vesta
とタイプしてみる。コードが実行されない時は実行権限の問題なので、chmod 755 ph2vesta.pyとタイプする。またはパスが通っていない。Macだと.bashrcの$PATHにこのコードのあるパスを設定しておくとよい。
実行するとモードのリストが出るので、見たいモードを選んで、その番号を入力するとそのモードに対応するVESTAのファイルが作成され、VESTAに表示される（Macの場合）。Mac以外の場合は最後のVESTAを呼び出すところ2行をコメントアウトしておくか、それらで使えるように変更する（例を示している）。
arrowの色、長さ、太さは0を入力すると変更できる（現在arrowは１種のみ）。なお長さのところをマイナスにすると矢印の向きを逆にできる。短時間で作ったのでエラー処理などはいい加減です。なお、stishoviteデータの10.2 cm-1以下の３つのモードは並進運動モードである（振動モードではない）。
とりあえず、沢山モードがある計算データを１個づつVESTAのVectorとして入力する苦行はなくなった。下図はその使用例。これはAlPO4 moganite相の振動モードの１つである。高温相の構造で振動計算を行なったら、周波数が負のモードが１つだけ得られた。これがソフトモードである。もちろん周波数が負（正確には虚数）になることはないので、これはこの変位に対して復元力が働かないことを意味している（構造が鞍点にあるため）。これらの変位方向に動かしていくと低温相の構造が得られる。この変位と低温相構造の対応は綺麗で、最初にこの図が出てきた時は感動した。この図は上記のph2vesta.pyを使ってVESTA用ファイルを作成して、さらにVESTA上で多面体表示などの設定を行なったものである（現在のph2vesta.pyだとそういう設定は実行前にやっておけばよい）。
実のところQuantum Espressoはあまりラマン強度計算のための環境が整っていない。どうもVASPを使った方がいいそうだが、これを見られている方はQuantum Espressoユーザーだと思うし、私はVASP使ってないので…
まず知っておくべきことは、現在のところ(〜ver. 7.2)ノルム保存(Norm Conservative) pseudopotentialにしか対応していないことである。メタル・半導体もダメ（計算させることはできる）。また、LDAのみ対応で、PAW, GGA, ultrasoftポテンシャルを使うとエラーで止まる。ただし、GGAについては何とか対応させることが可能である（以下のexample15のworkaround参照）。ノルム保存のポテンシャルは最近のQuantum Espressoのpseudopotentialリストでデフォルトでは出てこないが、こちらの左側のswitch tableのHartwigensen-Goedecker-Hutter PP tableかFHI PP tableをクリックすると出てくるpseudopotentialがノルム保存のもの。私のSiO2多形の計算では、Hartwigensen-Goedecker-Hutter PPは全然よくなかった。FHI PPはまともな結果が出る。なお、LDAを使う場合はpz(pw)とあるものを、GGAの場合はpbeとあるものをダウンロードする。なお適当なものがない場合は、Quantum Espressoのld1.xを使うと必要なpseudopotentialを作ることもできるが…私はずっと以前pawがまだ少なかった頃、自分でpawを作っていた。GIPAW計算ではpawしか使えないため。自分で作った場合は、それが正しいかチェックするための計算が結構大変である。
それらの対応をして、ph.xでlraman = .true.にして、振動数の時と同じように計算すればph.xの出力にラマンテンソルが、そしてdynmat.xの出力の最後にラマン強度（正確にはactivity）と偏光解消度が出ているはずである。なお、dynmat.xの計算は粉末測定に対応したもの（非偏光）になっている。
少しSiO2鉱物でやってみたが、ノルム保存のポテンシャルを使うと構造最適化で密度が実際よりもかなり高くなることが見られた。構造最適化はしないで、実測構造やpbesol-kjpawで最適化した構造を使ってscf計算するだけの方がいいかもしれない。比較的マシだったのが、FHIの*.pw-mt_fhi.UPFだった。
stishoviteの場合の入力・出力ファイル等を添付している（ページ一番下のところのqe-phonon.zip）。これはSi.pw-mt_fhi.UPF、O.pw-mt_fhi.UPFを使ってscf計算して、ph.xでlraman = .true.の計算、そしてdynmat.xでラマン強度を得るための入力ファイルと出力ファイルを圧縮したもの（Ver. 7.0で計算）。正確に比較してないが、強度は測定値と大体対応しているようである。振動周波数は上記の表よりさらに低い側にずれている。最も強度の高いA1gモード（全対称モード）では、偏光解消度がゼロに近いことなどもdynmat.xの出力結果から分かる。使っているpseudopotentialは、こちらからダウンロードできる。
石英で計算した例（青）を実測ラマンスペクトル（赤）と一緒に示す。波数が実測からちょっとずれるので、1.05を計算振動数にかけている。ガウス曲線を使ってピークらしくしている。最近石英のラマンスペクトルを測定してないので（低周波数領域のスペクトルはある）、少し前のデータを掘り出してきた。大体合っている。この場合TO-LO分裂を計算していないが（計算は可能）、石英は対称心がないので、TO-LO分裂が生じる。いくつか計算で再現されていないピークはTO-LO分裂したEモードと思われる。例えば500 cm-1くらいのピークがそう。
もう１つの計算例でforsterite(Mg2SiO4)を示します。これはGGAでも計算できるようにして（下のexample15のworkaroundを参照）、pbesolを使って計算したもの。ただし、PSlibraryのpbesolはノルム保存型ではないので、Pseudo Dojoからダウンロードしたノルム保存型でpbesolな擬ポテンシャルを使ってます。赤は合成のforsteriteペレットで測定したラマンスペクトル。周波数の補正をしてないので、計算結果はピーク位置が全体的に低周波数側にずれてますが、ピークを大体対応させられることが分かります。
先の石英のように多くの場合は計算と実測は結構合う。しかし合わない場合がたまにある。合わない理由はいくつか考えられるが、最近経験したケースでは、構造解析で得られた構造が実は平均構造で、実際にはより対称性の低い構造の間を移動していて、時間平均した結晶構造を計算で使ったために不一致が生じた。このようなケースは高温相で生じやすい。X線はどうしても時間・空間的に平均構造を与えるが、ラマンでは我々の時間スケールに比べて瞬間的な構造に対応したスペクトルを示すので、一致しないことが起こる。このような場合、計算した振動モードの低周波数モードをチェックすることが有用である。ソフトモードらしきものが見つかる場合がある。もし有ったら、その振動モードの原子変位方向に原子を移動させておいて、構造最適化を行う。変位させた構造が安定化したら、そのモードはソフトモードであり、低温側で相転移が期待される。
Quantum EspressoのPHonon/examplesディレクトリのexample15にラマン強度計算の例があり、CO2分子とZnO結晶の計算を行うことができる。しかし実はこの例ではlraman = .true.がコメントアウトされているので、このまま計算しても赤外の方は問題ないが、ラマンは計算されない（Ver. 7.2でもこのままだった）。また、lramanのコメントアウトを外しても、ZnOの場合にはメタル対応の計算となっており、メタルでラマン強度計算するとエラーが出る。さらに使われているpsuedopotentialがどちらもGGAで、これもエラーを出す（公式にはLDAのみ対応）。なので、この例を実行したい場合はZnOの場合はメタル対応の計算をやめる（それでいいのか疑問ではあるが）。LDAでノルム保存のpseudopotentialを使うことが必要。
なお、GGAの場合でも、このサイトのようにPHonon/PHディレクトリのphq_setup.f90の以下の行を直してやって（1を0にするだけ）、
変更前
IF (xclib_dft_is('gradient').and.(lraman.or.elop)) call errore('phq_setup', &
'third order derivatives not implemented with GGA', 1)
変更後
IF (xclib_dft_is('gradient').and.(lraman.or.elop)) call errore('phq_setup', &
'third order derivatives not implemented with GGA', 0)
その後にmake phでphだけコンパイルし直すとGGAでもエラーで止まらなくなる。third order derivativesを使わなくても結果はあまり変わらないとの話であるが、実測と比べるとか、LDAの結果と比べて特におかしくないかのチェックが必要である。最近いくつか計算してみているが、今のところ問題はないようである。
なので、lramanのコメントアウトを外して、pseudopotentialをLDAに変えるか、GGA対応にする上記の変更を行う。例で使われているポテンシャルはGGAであるがノルム保存でもあるので、GGA対応策を取ればそのまま使えるはず。ZnOの場合はメタル扱いを止める。これらを実施すればZnOのラマン強度が一応は計算されるようになる。ノルム保存でGGA(pbe, pbesol)な擬ポテンシャルはPseudo Dojoで入手できる。
ちなみにここにこのexample15を実行した例が出ているが、これは例のまま計算しているだけなので、ラマン強度は当然計算されておらず、示されたラマンスペクトルにはTHz単位の周波数がラマン強度のところに入っているようだ（なので直線的に強度が変わっている）。当然実測スペクトルとは強度が全く合いません。
上記のstishoviteの例も参考にしてください。
振動モードを動画的に表示するコードph2gif.pyを作りました。VESTAファイルの原子座標を振動モードの変位で少しづつずらせたものを沢山作成して、それをVESTAから画像pngファイルとして出力させます（これにはVESTAを外部から制御する機能を利用）。できたpngファイルをpillowというライブラリーを使ってgif動画を作成します。pillowのインストールが必要です（pipかcondaで可能）。download なお、これは例も含んだzipファイルです。使い方はQE用pythonコードの方に書いてます。
下はbrucite Mg(OH)2のOHの伸縮モードで、上がラマンアクティブなA1gモード、下が赤外アクティブなA2uモード。対向するOHの位相が２つのモードで逆になります。上側はOHの向きが逆で同位相なので、双極子モーメントが相殺されますが、下では相殺されないので赤外アクティブになります。分極率の方は上では振動により単調な増減が起こり、ラマンアクティブになります。下側でも分極率の増減自体はありますが、ちょうど変位ゼロを中心に対称的に増加が起こるために、変位ゼロでの分極率の微分がゼロになって、ラマンは不活性となります。
これらのモードは基本的にはOH伸縮振動だけですが、それらは同じ波数にはならずに60 cm-1ほど離れています。これはOH振動が相互作用をするためで、振動数が非常に近くなると分裂が起こります（correlation field splitting）。bruciteの場合にもこれが生じていて、振動数が分裂していることになります。ただ、それぞれラマンのみ赤外のみで観察されるので、２本のピークが同時に見られることはありませんが、そういう効果が存在することを意識しないと、ピークの圧力依存性などの解釈を誤ることになります。
これは現在研究している異極鉱hemimorphite（Zn4(Si2O7)(OH)2•H2O）の高圧相の振動モードです。約2.5 GPaで現れる高圧相について計算された低周波数に出てくるA1モードで、ソフトモードと思われます。実際、ラマン分光法でもソフトモードが観察されました（論文作成中）。この場合はSi2O7と４つのZnO4が組み合わさったユニットがそれぞれ逆方向に回転する運動になります。中央のユニットの反時計方向の回転をさらに進めると、低圧相の構造が得られます。
エンスタタイト(enstatite)は常圧常温で３つの多形(orthoenstatite, low-clinoenstatite, protoenstatite)が観察できます。天然岩石で出現するのはほとんどがorthoenstatite（実際には多少鉄を含む）で、low-clinoenstatiteは火山岩が急冷したような場合などちょっと特殊な状況で出現します。protoenstatiteは実験的には一部急冷ができるのですが、天然では報告例がないと思います。これら３相はラマンではprotoenstatiteとそれ以外は比較的簡単に区別できます。それは670 cm-1くらいのピークがprotoenstatiteでは１つなのですが、orthoenstatiteとclinoenstatiteでは２つになることで簡単に分かります。一方、orthoenstatiteとclinoenstatiteはかなりスペクトルが似ているので、ラマンスペクトルでの判別は結構難しいところがあります。ただ、低周波数側(80-120 cm-1)が測れている場合は、それぞれに特徴的なピークがあって、それで簡単に判別できます。orthoenstatiteでは82 cm-1に結構強いピークが、low-clinoenstatiteでは118 cm-1にピークが、protoenstatiteでは105 cm-1にピークが出ます。これはBruno Reynardさんらが提案している方法です(B. Reynard et al., Journal of European Ceramic Society, 28, 2459-2462, 2008)。
そういう低周波数のモードがどういう振動なのか気になったので、計算してgif動画を作ってみました。計算された振動数はほぼ2 cm-1以内で実測のラマンピーク位置と一致しました（みんな全対称のAgモードです）。
