チュートリアル06:ランダム力の時間相関関数

　本節ではFreeFlexを用いたランダム力の時間相関関数の計算方法を示す。

　今回は、水900 分子、Cl^-イオン1分子の系で、16個の初期構造を使用し、Cl^- の拡散係数を得ることを考える。初期構造の作成は、シェルのfor文を利用して、

> for (( i=0; i<=15; i++ )) do ii=`printf “%04d” $i`; ./add.exe add.nml $ii/system.ms; ./unlap.exe $ii/system.ms $ii/system2.ms; done

のようなコマンドを実行すればよい（add.nml や各分子のmsファイルはチュートリアルのディレクトリ内参照）。これで0000/ ~ 0015/ にunlap.exe による緩和後のmsファイルが得られる。

　初期構造が作成されたら実際のMDを走らせる。ここで注意したいのは、ランダム力の時間相関関数の計算では、注目分子の速度を固定する必要があることである。なぜなら、MDシミュレーションから得られる情報は注目分子に働く力であり、ランダム力とは一般化ランジュバン方程式を通じて

$F(t)=-{\frac {\partial U(x)}{\partial x}}-\int \limits _{0}^{\infty }\eta (\tau )v(t-\tau )d\tau +R(t)$

の関係があるためである。

$\langle R(t)R\rangle =\langle F(t)F\rangle -\langle F\rangle ^{2}$

が成り立つには、v=0 に固定し、v の時間依存性をなくさなければならない。このような計算をFreeFlex上で行うには、FreeFlex のネームリストに

=== input.nml =================================================================================
  .	
  .
  .
&addconstr
  condition = "FREEZE 901"
/
&monitor
  file = “monitor.dat”
  condition = “GROUP_F 901”
/
================================================================================================

のような指示を加える必要がある。 &addconstr ネームリストでは、分子を固定する命令、&monitor ネームリストでは分子に働く力を出力する命令をしている。詳しくは各ネームリストの説明を見ること。

　計算を実行すると、ディレクトリ内に monitor3.dat が作成される。ファイルの中身は、

=== monitor3.dat ================================================================================
&MONITOR
  .
  .
  .
/
  Fx(901), Fy(901), Fz(901)
     0.9573944741E-09    0.4891307747E-08    0.4926594522E-08
    -0.7770599862E-09    0.1329902147E-09    0.1808071790E-08
  .
  .
  .
================================================================================================

のようになっており、数字の羅列が x, y, z 方向に働く力を表わしている。このデータから tools/analysis.exe により時間相関関数を得ることができる。具体的は、

=== analysis.nml ================================================================================
&analysis
  skip = 5000
  method = "TCF"
/

&tcf
  Acolumn = 1
  Bcolumn = 1
  maxdata = 5000
  file = "tcf_xx.dat"
/
================================================================================================

のようなネームリストで解析手法を指定し、

./analysis.exe  analysis.nml  0*/monitor3.dat &

を実行することを実行することにより、cf_xx.dat に<F_x(t)F_X>のグラフを得ることができる。

拡散係数の計算について

拡散係数 $D$ はボルツマン定数 $k_{B}$ 、温度 $T$ を使って以下のように表すことができる。

$D=(k_{B}T)^{2}/\int _{0}^{T_{\rm {corr}}}\langle R(t)R\rangle dt$

$T_{\rm {corr}}$ は時間相関関数の周期である。ただし、 $t$ の大きいところでも相関が0に収束しないことがあり、それが誤差をもたらす要因にもなり得るので、実用的には以下のようにダンピングをかけて計算する方がよい。

$D=(k_{B}T)^{2}/\int _{0}^{T_{\rm {corr}}}\langle R(t)R\rangle \exp(-\alpha t^{2}/t_{0}^{2})dt$

Kikkawa(2016)では、 $\alpha =0.5$ 、 $t_{0}=T_{\rm {corr}}*0.5$ としている。以上の計算は例えば次に示すようなpythonスクリプトによって計算できる。ここでは $z$ 方向の力を配列に格納した後に、(1) $N(T_{\rm {corr}}=N\Delta t)$ で指定される周期分だけ時間相関関数を計算し、(2)時間相関関数にダンピングをかけたものを積分している。

=== tcf.py ================================================================================
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
file_force = "force.dat"
dt = 1.0e-15
N = 10000                 # T_corr = N * dt
n_skip = 0                # skip the n_skip lines
alpha = 0.5               # damping factor
KB = 1.3806485279e-23
T0 = 298.15
# Reading forces along z-axis
_, _, fz = np.loadtxt(fname=file_force, unpack=True, skiprows=14+n_skip)
# Calculating TCF
def calc_tcf(data, N):
    tcf = np.zeros(N)
    for j in range(N):
        s = 0
        for i in range(N):
            s += data[i]*data[i+j]
        tcf[j] = s
    return tcf
tcf = calc_tcf(fz, N)/N
time = np.linspace(0, dt*(N-1), N)
damp = np.exp(-time**2/(N*dt*0.5)**2*alpha)
zeta = integrate.simps(tcf*damp, time)/KB/T0
D = KB*T0/zeta
print("Diffusion coefficient: ", D)
# for plot
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(time, tcf*damp)
plt.show()
================================================================================================

チュートリアル06:ランダム力の時間相関関数

拡散係数の計算について

案内メニュー

検索