Hirano: ページの作成:「　本節ではFreeFlexを用いたランダム力の時間相関関数の計算方法を示す。　今回は、水900 分子、Cl^-イオン1分子の系で…」

2026-05-26T02:35:59Z

ページの作成:「　本節ではFreeFlexを用いたランダム力の時間相関関数の計算方法を示す。　今回は、水900 分子、Cl-イオン1分子の系で…」

新規ページ

本節では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シミュレーションから得られる情報は注目分子に働く力であり、ランダム力とは一般化ランジュバン方程式を通じて

<math>F(t)=-\frac{\partial U(x)}{\partial x}-\int \limits_{0}^{\infty}\eta(\tau)v(t-\tau)d\tau +R(t)</math>

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

<math> \langle R(t)R \rangle = \langle F(t)F \rangle - \langle F \rangle^2 </math>

が成り立つには、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 に<Fx(t)FX>のグラフを得ることができる。

== 拡散係数の計算について ==
拡散係数<math>D</math>はボルツマン定数<math>k_B</math>、温度<math>T</math>を使って以下のように表すことができる。

<math> D=(k_BT)^2/ \int^{T_{\rm corr}}_{0}\langle R(t)R \rangle dt </math>

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

<math> D=(k_BT)^2/ \int^{T_{\rm corr}}_{0}\langle R(t)R \rangle \exp(-\alpha t^2/t_0^2) dt </math>

Kikkawa(2016)では、<math>\alpha=0.5</math>、<math>t_0 = T_{\rm corr}*0.5</math>としている。
以上の計算は例えば次に示すようなpythonスクリプトによって計算できる。
ここでは<math>z</math>方向の力を配列に格納した後に、(1)<math>N(T_{\rm corr}=N\Delta t)</math>で指定される周期分だけ時間相関関数を計算し、(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:ランダム力の時間相関関数 - 版の履歴

Hirano: ページの作成:「 本節ではFreeFlexを用いたランダム力の時間相関関数の計算方法を示す。 今回は、水900 分子、Cl-イオン1分子の系で…」

Hirano: ページの作成:「　本節ではFreeFlexを用いたランダム力の時間相関関数の計算方法を示す。　今回は、水900 分子、Cl^-イオン1分子の系で…」