Quantum Espresso(QE)のプリ/ポスト処理プログラム

目次

1 はじめに

Quantum espresso はフリーの平面波擬ポテンシャル法による第一原理計算プ ログラムです。(http://www.quantum-espresso.org/) このプログラムを使っ て結晶やガラス等の固体物質の第一原理計算を行うことができます。第一原理 計算の大まかな手順は以下のとおりです。

(1) 結晶データ等の初期構造からインプットファイルの雛形をつくる
(2) インプットファイルを編集する
(3) プログラムの実行と収束の確認
(4) 計算から求めた波動関数や原子座標を評価

このうち(1)-(3)はプリ処理用のプログラムを使います。(4)はポスト処理プロ グラムを使って解析することになります。

構造や電子分布の確認、可視化はVESTA, VMD, xcrysdenを使います。VESTAは GUIで使いやすいのですが、静的な構造だけでMD計算等から得られる複数の構 造を扱えない欠点があります。またあまりに原子数が多すぎると重くて使いづ らいです。xcrysdenは古いプログラムで開発も止まっているようでが、QEの inputとoutputから直接、構造を可視化できます。VMDは、正直、インターフェー スがショボくて使い辛いですが、MDや原子数の制限は実質ありません。また自 由度がもっとも高くレイトレーシングな綺麗な表現が可能です。またすべての ファイル操作や画像表現をスクリプトで記述できるので、原子を動かしながら 各ステップの画像を連番のファイル名で保存しffmpegを使って画像を連結する ことで動画を作成することができます。

可視化用プログラムのリンクは以下のとおりです。全部フリー。

  1. xcrysden http://www.xcrysden.org/ (siliconからXを飛ばして使ってください)
  2. vesta http://jp-minerals.org/vesta/jp/ (各自のPCにインストールして使える)
  3. vmd http://www.ks.uiuc.edu/Research/vmd/ (各自のPCにインストールして使える)

2 プログラムのインストール

プログラムはpythonとC言語で作成しています。pythonが適切にインストール されていて*.pyはパスの通ったところに配置すれば動きます。*.cはコンパイ ルせねばなりません。コンパイルの方法は適当なディレクトリに*.cと Makefileを置いて

make

すれば、HOGE.cに対応した実行ファイルが生成されます。iccが使える場合は、 icc用のMakefileを-fで指定してください。

make -f Makefile.icc

プリ/ポスト処理プログラム一式をtgzしたもの qe-pre-post.tgz
(研究室内からのみアクセス可です。プログラムに興味があれば大窪に連絡してください)

3 プリ処理

3.1 cif2qein.py

cifファイルからQEのインプットファイルを生成します。またnx, ny, nzで  super cellを作成できます。y方向(b軸)に2つのsuper cellを作成したいなら ば、つぎのようにします。

# cif2qein.py -ny=2  hoge.cif

nxとnzも同じように指定できます。プログラムが正しく実行するとhoge.inが 生成されます。QEINPUT.inは、qe用のインプットファイルの雛形になります。 擬ポテンシャルファイルの指定等、インプットファイルを正しく編集せねば 計算できません。少なくとも雛形に書いているパラメータを意味をすべて理 解してから計算を行ってください。また擬ポテンシャル用のファイル等は、 マニュアルでダウンロードする必要があります。

3.2 make_qein.py

cifファイルからQEのインプットファイル(構造最適化とNMR計算)job投入ス クリプトを生成します。入力するcifファイルは、vestaで対称性を消して exportしたものです。Vestaでcifファイルを開いたらEdit->Edit Data->Unit cellでRemove symmetryした後にcifファイルをexportして保存 してください。

make_qein.py  hoge.cif 

これで、hoge.cifからhoge.qeinhoge.shが生成され、擬ポテンシャルファ イルの候補が次のように標準出力に表示されます。

----------------------------------------------------------------------------
posible pseudo potential files in ~/projects/pseudo/kjpaw/
----------------------------------------------------------------------------
-Li(6.941) wfc=103.0 rho=410.0
*Li.pbe-s-kjpaw_psl.1.0.0.UPF
 Li.pbe-sl-kjpaw_psl.1.0.0.UPF
-P(30.973762) wfc=34.0 rho=171.0
*P.pbe-n-kjpaw_psl.1.0.0.UPF
 P.pbe-nl-kjpaw_psl.1.0.0.UPF
-S(32.066) wfc=39.0 rho=181.0
*S.pbe-n-kjpaw_psl.1.0.0.UPF
 S.pbe-nl-kjpaw_psl.1.0.0.UPF

defaultで、擬ポテンシャルと周期表から元素データを以下のファイルから読 み込みます。make_qein.pyにベタ打ちしています。

ppdir = "~/projects/pseudo/kjpaw/"
eperiodic = '~/.emacs.d/site-lisp/misc/eperiodic.el'

3.3 updateqein.py

QEのlogファイルとインプットファイルからk点やecutwfc等の収束条件、また は、calculationを変えたインプットファイルを生成します。また、nmr計算 用のインプットファイルを生成します。

4 ポスト処理

4.1 qeout2trans

QEのMD計算から得られたlogファイルから必要なデータを抜き出します。可視 化用cif(最終ステップのみ)、解析用の*.lammpstrj、統計情報の*.statistics を生成します。データが大きすぎるときは、skipを指定してデータを読み飛ば すようすると良いです。-oで出力ファイルのbasenameを指定します。vc-relax でCELL_PARAMETERSが三角行列になってなくてもプログラムは正しく動きます。

qeout2trans -i hoge.log -o huga

引数が不足しているとプログラムは実行できません。プログラムを実行して表 示されるヘルプを参照すればオプションの渡し方がわかると思います。

なお以下のソースコードをqeout2trans用に修正しています。修正したコード で生成したpw.xを使わないとqeout2transが正しく動作しません。

  • PW/src/dynamics_module.f90
  • PW/src/output_tau.f90

4.2 dmol2trans

dmolのlogファイルをlammpstrjに変換するプログラム。qeout2transと使い方は同じ。

dmol2trans -i test.outmol -o test

4.3 qerelax2convert.py

inputやlogファイルをcifファイルに変換するプログラム。格子定数も出力します。

hoge.logを変換してhoge.lammpstrj(lammpstrj)とhoge.dat(格子定数の変化)それからhoge.cif(最終ステップの構造)を生成する場合。

qerelax2convert.py hoge.log -o hoge

4.4 dump2analysis

qeout2transで生成した*.lammpstrjの解析(動径分布G(r)や平均二乗変位msd等の計算)を行います。 もし、lammpsで生成したdumpの結果も解析できます。ただし以下のdumpコマンドでlammpstrjを出力する必要があります。

dump        1 all custom 1 hoge.lammpstrj id type element x y z vx vy vz
dump_modify 1 element Ar sort id

4.4.1 動径分布関数の計算

配位数や原子間距離等の局所構造を表す動径分布関数の計算を計算します。動径分布関数\(g(r)\) は平均数密度\(\rho_0\) からのずれとして定義されます。 \[ g(r) = \frac{\left< n(r+\Delta r) \right>}{4\pi r^2 \Delta r}\cdot \frac{1}{\rho_0} \]

「コンピュータシミュレーションによる物質化学」 p.50参照

プログラムの実行は、計算modeをgr(-m gr)として行います。計算対象の原子 ペアは、元素記号や原子タイプ、複数のidを指定して行うことができます。id はlammpstrjの数字です。VMDで原子をピックアップすると0から始まるidが表 示されるので間違えないようしてください。また–rcutで指定した距離内の平均結合距離を計算します。

ID=0,1,2,3とID=8,9,10,11の(r)を計算する場合(IDは0からスタートすることに注意))

dump2analysis -m gr -a 0,1,2,3 -b 8,9,10,11 -i hoge.lammpstrj -o hoge.gr

PとLiのG(r)を計算する場合

dump2analysis -m gr -x P -y Li -i hoge.lammpstrj -o hoge.gr

「コンピュータシミュレーションによる物質化学」 p.47参照

4.4.2 MSDの計算

原子の平均二乗変位\(MSD(t)\) の計算をします。 \[ MSD(t) = \left< \left[\vec{r_i(0)} - \vec{r_i(t)}\right]^2 \right> \] \(MSD(t)\) は自己拡散係数と以下の関係があります。 \[ D = \frac{1}{6t} MSD(t) \] よって\(MSD(t)\) が直線になったところの傾きから\(D\) を求めることができます。

「コンピュータシミュレーションによる物質化学」 p.50参照

計算対象する原子は、元素名、原子タイプ、IDで選択します。–dtで1ステッ プの時間(fs)を指定する必要があります。またMSDのスタートステップをシフ トしながら計算させる場合は、–tauと–tau_shiftを指定します。–tau_shiftが 短いと相関が残って結果がおかしくなるので注意が必要です。

LiのMSDを1ステップ2fsで計算する場合は以下のとおりです。

dump2analysis -m msd -x Li --dt 2 -i hoge.lammpstrj -o hoge.msd

4.4.3 角度分布の計算

S-P-S結合の角度分布は、以下のように計算します。–rcut_abと–rcut_bcで 結合の有無を判断する距離を指定する必要があります。

dump2analysis -m angle -x S -y P -z S --rcut 3.0 -i hoge.lammpstrj  -o huga.angle

4.4.4 速度相関関数の計算

次の式にしたがって指定した原子の速度関数を計算します。 \[ C(t) = \vec{v(0)}\cdot \vec{v(t')} \] また拡散係数と速度相関関数は以下の関係があります。 \[ D = \frac{1}{3}\int_0^\infty C(t) dt \] 「コンピュータシミュレーションによる物質化学」 p.50参照

さらに速度相関関数のフーリエ変換よりIRスペクトルを計算できます。 \[ I(\nu) = \int_0^\infty \exp(-2\pi it ) C(t) dt \]

–dtで1ステップの時間(fs)を指定する必要があります。MSDの計算と同様、 に–tauと–shiftで相関関数を計算するステップ数とシフトを指定して統 計を稼ぐことができます。十分に相関関数が減衰するよう–tauの指定する 必要があります。

1000ステップ分の速度相関関数を100ステップシフトさせながら計算する場 合。

dump2analysis -m vcorr -i hoge.lammpstrj -o hoge.vcorr -x Li --dt 1 --tau 1000 --shift 100

4.4.5 回転相関関数の計算

次の式にしたがって指定した水分子の回転相関関数を計算します。 \[ C(t) = \left< \frac{3\cos^2\theta -1 }{2}\right> \] \(\cos \theta\) は水の永久双極子モーメントの向きです。 また回転相関関数を時間積分すると^2H NMRの緩和時間と関係する回転相関時間が求まります。 \[ \tau_R = \frac{1}{3}\int_0^\infty C(t) dt \] 日本化学会編「溶液の分子論的描像」 p.99参照 「分子シミュレーション入門」 p.101参照

–dtで1ステップの時間(fs)を指定する必要があります。MSDの計算と同様、 –tauと–shiftで相関関数を計算するステップ数とシフトを指定して統計 –を稼ぐことができます。十分に相関関数が減衰するよう–tauの指定する –必要があります。

1000ステップ分の速度相関関数を100ステップシフトさせながら計算する場 合。

dump2analysis -m vcorr -i hoge.lammpstrj -o hoge.vcorr -x Li --dt 1 --tau 1000 --shift 100

4.4.6 空間分布の計算

空間をvoxcelで分割し、voxcel毎に指定した原子が存在する存在頻度、速 度、力のいずれかでcubeファイルを作成します。cubeファイルのパラメー タは、–cube_typeで指定します。

--cube_type 1  (指定原子の頻度)
--cube_type 2  (指定原子の速度ノルム)
--cube_type 3  (指定原子の力ノルム)

空間の分割数はa,b,c軸方向の分割数を–na, –nb, –ncで指定します。微 小体積の1辺の長さは、a/na(a軸方向)になるので、以下の式を満たせば voxcelは、1辺の長さが同じ平行6面体になります。 \[ \frac{a}{na} = \frac{b}{nb} = \frac{c}{nc} \]

なおcubeファイルの原子位置は全ステップの平均座標になります。VMDで isosurfaceの表示は、Drawing methodでisosurfaceを選択した後、 isovalueとRangeを調整します。色を変えたい場合は、Coloring methodで Volumeを選択しTrajectory->Color Scale Data Rangeを調整します。

dump2analysis -m cube -x Li --cube_type 1 -i hoge.lammpstrj -o hoge.cube

※VMDのbugかと思いますがisosurfaceの原点が少しずれます。na, nb, nc が大きいとほとんど問題になりませんが、気になる場合は手動で3行目の offsetの値(1,2,3列のデータ)を調整してください。

4.4.7 有限サイズでの空間分布の計算

空間をvoxcelで分割し、指定した有限サイズの原子が占めるvoxelの密度分 布を求めます。半径は–radiusで与えます。radiusのデフォルト値は0.68 Å (Li+のイオン半径)です。各原子と全voxelの中心距離を求めて、指定 した半径より小さい場合、原子が占めるvoxelとしてカウントします。(計 算にかなり時間がかかります)。

空間の分割数はa,b,c軸方向の分割数を–na, –nb, –ncで指定します。微 小体積の1辺の長さは、a/na(a軸方向)になるので、以下の式を満たせば voxcelは、1辺の長さが同じ平行6面体になります。 \[ \frac{a}{na} = \frac{b}{nb} = \frac{c}{nc} \]

なおcubeファイルの原子位置は全ステップの平均座標になります。VMDで isosurfaceの表示は、Drawing methodでisosurfaceを選択した後、 isovalueとRangeを調整します。色を変えたい場合は、Coloring methodで Volumeを選択しTrajectory->Color Scale Data Rangeを調整します。

dump2analysis -m cube_radius --radius 0.68 -x Li -i hoge.lammpstrj -o hoge.cube

※VMDのbugかと思いますがisosurfaceの原点が少しずれます。na, nb, nc が大きいとほとんど問題になりませんが、気になる場合は手動で3行目の offsetの値(1,2,3列のデータ)を調整してください。

4.4.8 jumpの計算

空間をvoxcelで分割し、voxcel毎に指定した原子が指定したステップの間 に移動する距離でcubeファイルを作成します。移動する距離は–tauで指 定します。他は–nx, –ny, –nzはcubeファイルの計算と同じです。

--tau 1  (jumpの時間ステップ)
dump2analysis -m cube_jump -x Li --tau 1 -i hoge.lammpstrj -o hoge.cube

4.4.9 van Hove関数の計算

\(N\) 粒子の jump時間\(t\) と移動距離の分布であるvan Hove関数を計算します。\(t\) は–tauで指定します。–shiftで統計を稼ぎくことができます。 \[ G(r,t) = \frac{1}{N}\sum_{i=1}^{N} \left< \delta\left(r - |r_i(t) - r_i(0)| \right) \right> \]

5fsでのvan Hove関数 \(G(r, 500fs)\) を10ステップのシフトで計算する場合。

dump2analysis -m vanHove -x Li --dt 5 --tau 100 --shift 10 -i hoge.lammpstrj -o hoge.cube

4.5 cube2slice.py

cubeファイルを任意の面でスライスしてコンタープロットします。

-a スライスする面と垂直の軸を指定 a, b, cのいずれか
-s スライスする面の幅の始まりと終わりを指定する。0から1の範囲

5 その他

5.1 submit_ito.py

ITO上でmdを連続して行うために、目的とするjobがrunまたはqueしてい るか確認する。runをqueもしていなかったら、logファイルを日付にmvして jobをsubmitする。外部からcron+sshで自動実行することを想定しています。 中身を読んで自分用に改造して使うこと。

5.2 qelog2in.py

relax作業中に計算がストップしてしまった時にrestartさせるためのpython プログラム。logファイルから最終ステップの格子と座標データを読み込ん で、inputファイルを更新します。元のinputファイルとlogファイルは、バッ クアップとしてファイルの後ろに3つの数字をつけてコピーします。

6 TIPS

6.1 .vmdrc

vmdの初期設定での表示は、しょぼいので qe-pre-post.tgz に含まれている.vmdrcを$HOMEに置いて初期の描画 方法を指定すると良いと思います。VESTAのマネをして原子サイズや色を設定 しています。

6.2 vmdで動画

MDの結果からvmdで連番の画像ファイルを作成して動画を作成できます。vmd はtcl/tk言語でスクリプト処理できます。

連番の画像ファイルを生成できたら以下のffmpegで動画を作成します。pptにも貼り付けれます。

ffmpeg -i %04d.png -qscale 0 -f mp4 -vcodec libx264 -pix_fmt yuv420p out.m4v

fps(frames/sec)は defaultで30くらいです。なので30秒の動画をつくりたかっ たら画像ファイルを900毎準備すればよいことになります。

著者: 大窪 貴洋

Created: 2019-12-03 火 19:34

Validate