dump2sq
Table of Contents
はじめに
複数の時間ステップからなる時間ステップの構造データから、動径分布関数 \(g(r)\) 、積算配位数 \(CN(r)\) 、構造因子 \(S(Q)\) を計算するプログラムです。 Triclinicなセルでも計算できます。原子ペア毎の\(S(Q)\) を計算した後、原子 散乱因子または中性子散乱長を使って、X-rayと中性子回折実験的で観測され る \(S(Q)\) を計算します。\(g(r)\) は、セルサイズ \(L/2\) を超えても計算でき ます。\(L/2\sqrt{3}\) まで確保しても良さそうです。産総研の土田さんにアド バイスいただきました。
https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
1ステップ毎に\(g(r)\) と\(S(Q)\) を計算して最後に時間平均しています。よっ て、NPT計算でセルサイズが変わったとしても問題なく計算できます。計算コ ストのほとんどは、原子間距離の計算です。
入力ファイル
拡張子で入力ファイルのformatを判断しています。
LAMMPS
dumpコマンドで生成したファイルで、拡張子は.trajectryです。1ス テップ目で原子のならび決めて、2ステップ目以降は、同じ原子の並びとし て計算しています。時間ステップで原子の並びが変わると、正しい結果に ならないので注意してください。データの並びを決め打ちしています。そ のため、以下のようにdumpとdump_modifyで元素記号の出力し、idでsortし ておく必要があります。
dump 1 all custom 10 hoge.lammpstrj id type element x y z dump_modify 1 element Si O sort id
XYZ
xyz形式で、拡張子は.xyzまたは.10です。す。1行目に原子数、2行目に Step数と、格子定数を記述しておきます。alpha, beta, gammaがない場合 は、90, 90, 90とします。3行目以降にデータがならびます。2ステップ目 以降も同じ形式でデータを並べます。
420 Step 225837 18.4116484 18.4116484 18.4116484 90 90 90 Si 10.560073 17.438902 6.614775 O 17.812448 0.110128 0.3462461 . . .
出力ファイル
4つのファイルに拡張子を追加して出力します。basenameは、指定しなけれ ば、インプットファイルと同じになります。
積算配位数
拡張子はcnです。原子ペア毎の \(CN(r)\) が保存されています。
動径分布関数
拡張子はgrです。原子ペア毎の \(g(r)\) と全原子 \(G(r)\) が保存されています。
構造因子
拡張子はsqです。原子ペア毎の \(S(Q)\) とX線と中性子で観測される \(S(Q)\) が保存されています。
インストール
使い方
並列計算のために、環境変数をセットしてスレッド数を指定しておく。
export OMP_NUM_THREADS=4
example以下のglass.lammpstrjを引数で渡して実行してみます。
dump2sq example/glass.lammpstrj
デフォルトで、\(g(r)\) を計算する\(r\) の最大値は、もっとも小さいセル長 さの半分としています。プログラムが実行されると、次のように出力されて、 4つのファイルが出力されます。
basename: example/glass Maximum cell length = 17.284100 Rmax = 8.642050 Start Gr and Sq calculation... Calculation 0 steps... Calculation 10 steps... Calculation 20 steps... Calculation 30 steps... Calculation 40 steps... Calculation 50 steps... Calculation 60 steps... Calculation 70 steps... Calculation 80 steps... Calculation 90 steps... Total: 99 steps example/glass.coeff was generated. example/glass.cn was generated. example/glass.gr was generated. example/glass.sq was generated.
\(g(r)\) の\(r\) ステップや最大値、\(Q\) の範囲をコマンドラインで指定できます。詳細 は、ヘルプを参照してください。ヘルプで表示される[ ]の値はデフォルト 値です。GRmaxは、セルサイズを超えても計算できますが、物理的な意味は ないでしょう。
./dump2sq -h
usage: dump2sq [options] inputfile
inputfile format
LAMMPS: *.lammpstrj
FEMTECK: *.10
options:
-h --help show this help message
--SQmin=%f Q minimum [0.2]
--SQmax=%f Q maximum [30]
--SQdel=%d Q delta [0.02]
--GRmin=%f r minimum [0.025]
-r --GRmax=%f r maximum [L(minium cell length)*0.5]
-s --GRdel=%d r delta [0.05]
-o --output=%s Output basename [basename]
tools
coeffploter.py
元素を指定して、原子散乱因子と中性子散乱長をプロットする。引数に元 素名を複数していできる。
gr2sq.py
Qの刻みを変えて再計算するために、dump2sqで生成したgrファイルから S(Q)を計算しなおす。結晶のパターンを計算するために、\(2\theta\) で\(Q\) を指定して計算することもできる。\(\theta\) と\(Q\) の関係は以下のとおり。 defaultで波長は、CuのK alpha線の1.54183 Åとしている。 \[ Q = \frac{4\pi}{\lambda}\sin(\theta) \] 計算した$S(Q)$はpngとxlsxを生成する
cn2Iq.c
dump2sqで生成された*.cnと*.coeffから$I(q)$を計算するCプログラム。 openmpでスレッド並列する。f*fは足していない。sincのwindow関数をかけ る。コンパイルは、
gcc -fopenmp -O3 -Wall cn2Iq.c -o cn2Iq
実行は、*.cnと*.coeffを渡す。I(q)は標準出力に出力されうｒ．
./cn2Iq fort.cn fort.cn.coeff
複数のcnファイルやxslxに保存したい場合は、cn2Iq.pyを使う。
cn2Iq.py
subprocessでCのcn2Iqを呼び出して、I(Q)のpngとxlsxを生成する。Qの範囲 またはthetaの範囲を指定できる。