diff --git a/Readme.org b/Readme.org new file mode 100644 index 0000000..3216136 --- /dev/null +++ b/Readme.org @@ -0,0 +1,223 @@ +#+TITLE: dump2analysis +#+AUTHOR: 大窪 貴洋 (千葉大院工) ohkubo.takahiro@faculty.chiba-u.jp +#+LANGUAGE: ja +#+EMAIL: ohkubo.takahiro@faculty.chiba-u.jp +#+OPTIONS: toc:2 num:nil author:t creator:nil LaTeX:nil \n:nil ^:{} +#+SETUPFILE: https://amorphous.tf.chiba-u.jp/org-html-themes/org/theme-readtheorg.setup +#+MACRO: color @@html:$2@@ +#+MACRO: kakomi @@html:
$1
@@ + + +* はじめに +lammpsのdump commandで生成したtrajectryファイルを解析できます。ただし +以下のdumpコマンドでlammpstrjを出力する必要があります。 +#+BEGIN_EXAMPLE +dump 1 all custom 1 hoge.lammpstrj id type element x y z vx vy vz +dump_modify 1 element Ar sort id +#+END_EXAMPLE + +* 動径分布関数の計算 +配位数や原子間距離等の局所構造を表す動径分布関数の計算を計算します。動径分布関数$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からスタートすることに注意)) +#+BEGIN_EXAMPLE +dump2analysis -m gr -a 0,1,2,3 -b 8,9,10,11 -i hoge.lammpstrj -o hoge.gr +#+END_EXAMPLE + +PとLiのG(r)を計算する場合 +#+BEGIN_EXAMPLE +dump2analysis -m gr -x P -y Li -i hoge.lammpstrj -o hoge.gr +#+END_EXAMPLE + +「コンピュータシミュレーションによる物質化学」 p.47参照 + +* 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で計算する場合は以下のとおりです。 +#+BEGIN_EXAMPLE +dump2analysis -m msd -x Li --dt 2 -i hoge.lammpstrj -o hoge.msd +#+END_EXAMPLE + +* 角度分布の計算 + +S-P-S結合の角度分布は、以下のように計算します。--rcut_abと--rcut_bcで +結合の有無を判断する距離を指定する必要があります。結合原子を逐次、検索 +するので遅い。結合の組み換えがない場合は、自作した方が良いでしょう。 +#+BEGIN_EXAMPLE +dump2analysis -m angle -x S -y P -z S --rcut 3.0 -i hoge.lammpstrj -o huga.angle +#+END_EXAMPLE + +* 速度相関関数の計算 + 次の式にしたがって指定した原子の速度関数を計算します。 + $$ + 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ステップシフトさせながら計算する場 + 合。 + #+BEGIN_EXAMPLE + dump2analysis -m vcorr -i hoge.lammpstrj -o hoge.vcorr -x Li --dt 1 --tau 1000 --shift 100 + #+END_EXAMPLE + +** 回転相関関数の計算 + 次の式にしたがって指定した水分子の回転相関関数を計算します。 + $$ + 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ステップシフトさせながら計算する場 + 合。 + #+BEGIN_EXAMPLE + dump2analysis -m vcorr -i hoge.lammpstrj -o hoge.vcorr -x Li --dt 1 --tau 1000 --shift 100 + #+END_EXAMPLE + +* 空間分布の計算 + 空間をvoxcelで分割し、voxcel毎に指定した原子が存在する存在頻度、速 + 度、力のいずれかでcubeファイルを作成します。cubeファイルのパラメー + タは、--cube_typeで指定します。 + #+BEGIN_EXAMPLE + --cube_type 1 (指定原子の頻度) + --cube_type 2 (指定原子の速度ノルム) + --cube_type 3 (指定原子の力ノルム) + #+END_EXAMPLE + + 空間の分割数は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を調整します。 + + #+BEGIN_EXAMPLE + dump2analysis -m cube -x Li --cube_type 1 -i hoge.lammpstrj -o hoge.cube + #+END_EXAMPLE + + ※VMDのbugかと思いますがisosurfaceの原点が少しずれます。na, nb, nc + が大きいとほとんど問題になりませんが、気になる場合は手動で3行目の + offsetの値(1,2,3列のデータ)を調整してください。 + +* 有限サイズでの空間分布の計算 + 空間をvoxcelで分割し、指定した有限サイズの原子が占めるvoxelの密度分 + 布を求めます。半径は--radiusで与えます。radiusのデフォルト値は0.68 + \AA (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を調整します。 + + #+BEGIN_EXAMPLE + dump2analysis -m cube_radius --radius 0.68 -x Li -i hoge.lammpstrj -o hoge.cube + #+END_EXAMPLE + + ※VMDのbugかと思いますがisosurfaceの原点が少しずれます。na, nb, nc + が大きいとほとんど問題になりませんが、気になる場合は手動で3行目の + offsetの値(1,2,3列のデータ)を調整してください。 + +* jumpの計算 + 空間をvoxcelで分割し、voxcel毎に指定した原子が指定したステップの間 + に移動する距離でcubeファイルを作成します。移動する距離は--tauで指 + 定します。他は--nx, --ny, --nzはcubeファイルの計算と同じです。 + #+BEGIN_EXAMPLE + --tau 1 (jumpの時間ステップ) + #+END_EXAMPLE + + + #+BEGIN_EXAMPLE + dump2analysis -m cube_jump -x Li --tau 1 -i hoge.lammpstrj -o hoge.cube + #+END_EXAMPLE + + +* 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ステップのシフトで計算する場合。 + + #+BEGIN_EXAMPLE + dump2analysis -m vanHove -x Li --dt 5 --tau 100 --shift 10 -i hoge.lammpstrj -o hoge.cube + #+END_EXAMPLE + +** cube2slice.py +cubeファイルを任意の面でスライスしてコンタープロットします。 +#+BEGIN_EXAMPLE +-a スライスする面と垂直の軸を指定 a, b, cのいずれか +-s スライスする面の幅の始まりと終わりを指定する。0から1の範囲 +#+END_EXAMPLE + + diff --git a/angle.c b/angle.c index 65482fb..62d4d6e 100644 --- a/angle.c +++ b/angle.c @@ -78,7 +78,7 @@ void EstimateAngle(PARAM *param){ /* 距離の計算 */ for (s=0; ssteps; s++){ if (s%10 == 0) { - printf("Calculating Gr step: %d\r", s); + printf("Calculating Angle step: %d\r", s); fflush(stdout); } /* 原子の選択 */