#+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:<font color="$1">$2</font>@@
#+MACRO: kakomi @@html:<div style="margin: 1em auto; padding: 1em; width: 100%; border: 2px dotted #ff0000;background-color: #ff;border-radius: 1px">$1</div>@@


* はじめに
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