GIPAW出力の処理 の変更点


#author("2021-03-30T08:37:09+09:00;2021-03-26T10:06:00+09:00","default:masami","masami")
#author("2021-04-17T16:58:57+09:00;2021-03-26T10:06:00+09:00","default:masami","masami")
#image(http://www.misasa.okayama-u.ac.jp/~masami/images/OUOperatingGrant.png,right,20%)
*QE-GIPAW出力(NMR化学シフト)の処理 [#g0eb40d6]
-Quantum EspressoのGIPAWコードを使ってNMR化学シフトやEFG(電場勾配)が計算できる(ESRのパラメータも)。普段それを利用して、化学シフトの予測計算を行っているが、席が数個だと出力からコピー&ペーストすれば済むが、最近は独立な席が数十個あるケースを計算しているので、化学シフト値を取り出すだけでもかなり大変。そこで、Pythonで簡単に化学シフト値のリストを取り出せるようにした。これは簡単な処理。
-実測スペクトルと比較する場合には、計算された化学シフト値にGaussian分布をかけて、スペクトルとして比較した方が分かりやすい。この部分の計算もPythonで書いてみた。結果はcsvファイルに出力する。計算時にFWHMを引数で指定できるようにした。これを使って、tridymiteの3つの多形のSUP{SIZE(9){29}}Si MAS NMRスペクトルを予測した結果をプロットしてみた。プロット自体はPlot2 Proを使っている。横軸は出力そのままで、標準物質による補正をまだしていない。青い線で示した多形(PO10)は独立なSi席が実に80もある。これの計算は18 core & 64 GB RAMのIntel i9 Linux PC(私物)で行ったが、pw+gipawで2週間かかった。
#image(http://www.misasa.okayama-u.ac.jp/~masami/images/tridymite-29Si-NMR.png,right,40%)
-赤い線で示した多形(MC)には実測の報告があって、それと大体一致する。一方、緑で示した多形(MX-1)も実測がいくつか報告されている。今回の計算はそれらとの一致がそれほど良くないが、MX-1はインコメンシュレート相とされており、そこに原因がありそうである。今回は報告されているsuperstructureを使って計算している。青で示した多形(PO10)ついては実測自体がなされていないようだ。また、以下に書いている構造との相関からは、どうもこのPO10の計算に少し問題があるかもしれないことが分かってきた。
*化学シフトと局所構造の相関をみる
-SUP{SIZE(9){29}}Si MAS NMRの場合は、化学シフトとSi原子周りの局所構造との相関がよく調べられている。そして経験的な関係式も多く発表されている。Quantum Espressoの構造最適化結果(一度VESTAで読めるようにして、xtlファイルで出力)を処理して、局所構造パラメータを計算するコードをPythonで作成した。それと上記の化学シフトを取り出したものを使って相関をみた。
#image(http://www.misasa.okayama-u.ac.jp/~masami/images/shift_vs_local_structure.png,right,40%)
-この図はMC多形の場合の相関で、横軸にSherriff and Grundy (Nature, 332, 819, 1988)の式4を使ったものである。これはSi-O-Si角、Si-Oの中点へのもう1つのSiからの距離などの局所構造パラメータから計算される。この式自体は一般のケイ酸塩鉱物に使えるのであるが、今回は対象がSiO2組成だけなので、式4のbond valenceの寄与は省略している(第2近接が全てSiなので)。非常に良い相関が得られた。もっともより簡単なcos(Si-O-Si角)を横軸に使ってもかなり良い相関が得られるが。この場合は実験との相関もよい。
-これにMX-1の計算した相関を重ねると、ほぼ一致した。一方、PO10の結果を重ねたところ、半分くらいは同じ直線に乗るが、残りは最大3 ppmくらい下側にずれる。どうも計算結果に問題があるようである。PO10は原子数が多くて計算時間が非常に長いため、cutoff, K mesh等を振って計算結果が正しいかの確認できてなかったが、構造の相関の方からダメ出しが出た。一部合うところ、バラつきが不思議ではあるが、現在K meshを増やして再計算中(cutoffsは3相で全く同じなので多分そのせいではなさそう)。しかし時間がかかり、GIPAW計算は1.5ヶ月かかる予定で、途中で定期点検の停電が来てしまう。GIPAWは(pw.xも)restartが可能なので、停電があってもいいのだが。
-これにMX-1の計算した相関を重ねると、ほぼ一致した。一方、PO10の結果を重ねたところ、半分くらいは同じ直線に乗るが、残りは最大3 ppmくらい下側にずれる。どうも計算結果に問題があるようである。PO10は原子数が多くて計算時間が非常に長いため、cutoff, K mesh等を振って計算結果が正しいかの確認できてなかったが、構造との経験的な相関の方からダメ出しが出た。一部は合っているところが不思議ではあるが、現在K meshを増やして再計算中(cutoffsは3相で全く同じなので多分そのせいではなさそう)。しかし時間がかかり、GIPAW計算は1.5ヶ月かかる予定である。定期点検の停電が途中で来たので、一度止めて、restartした。GIPAWは(pw.xも)restartが可能なのである。
-GIPAWも途中で停止する時はprefix.EXITという空のファイル名を同じディレクトリに作成する。この作成にはtouchコマンドを使う。停止して、prefix.EXITが消去される。GIPAWの入力ファイルでは、restartモードがデフォルトなので、(restart_modeをfrom_scratchとしていない限り)同じ入力ファイルを使って再実行すると、K pointsの途中から計算を始めてくれる。