体積粘性係数の分子動力学計算 -

等方的な流体の場合、粘性係数には、shear viscosityとbulk viscosity(体積粘性係数)の2種類がある。水のような液体の場合、非圧縮性近似をするのが一般的なので、Navier-Stokes方程式で体積粘性係数がかかる項は0になる。気体の場合も、流速が音速よりずっと遅い場合には非圧縮近似はよい近似になっているとされている。固体の場合は、弾性が支配的なので、通常、粘性そのものが考慮されない。

Wikipedia情報によれば、体積粘性係数の値は多くの場合、あんまり知られてないらしい

Volume viscosity/Measurement
https://en.wikipedia.org/wiki/Volume_viscosity#Measurement

The volume viscosity of many fluids is inaccurately known, despite its fundamental role for fluid dynamics at high frequencies. The only values for the volume viscosity of simple Newtonian liquids known to us come from the old Litovitz and Davis review, see References. They report the volume viscosity of water at 15 °C is 3.09 centipoise.

そもそも、通常の流体運動の範疇では体積粘性係数項は消えてしまって観測量と関わってこない(というか、NS方程式から消えてるだけだけど)から、測定が難しいのは仕方ない。多分、体積粘性係数を考慮に入れるべき身近な現象は音の減衰だと思う。なので、体積粘性係数を測定する方法として、音を利用するのは自然に思える。

Bulk viscosity and compressibility measurement using acoustic spectroscopy
http://aip.scitation.org/doi/abs/10.1063/1.3095471?journalCode=jcp

という論文には"To date, measurement of ultrasound attenuation is the only known way of experimentally studying bulk viscosity."と書かれている(Section Bの冒頭)。この論文によれば、25度に於ける水の体積粘性係数は、2.4cPであるらしい(TABLE III)

医療分野では、流体ではない生体組織の弾性や粘性(の分布)を測定するエラストグラフィと総称される手法がいくつか存在し、知る限り、いずれの方法でも超音波を何らかの形で利用している。
Elastography
https://en.wikipedia.org/wiki/Elastography
elastographyという言葉は、名前の通り、弾性の計測を主要な目的としている。癌などの腫瘍組織は一般に周囲の組織より硬くなる(乳がんが、しこりとして検出できる理由でもある)ことから、病変検出に有用と考えられている。一方、音の伝播の方程式に粘性項を入れるのは、固体であっても流体と変わるところはないので、同じ方法で粘性係数を決めることも理屈上はできる。実際に、そのような研究は存在するものの、弾性に比べると粘性を決定する医学的意義はよく分からないので、広く行われているとは言いがたいっぽい。しかし、何はともあれ、固体・流体を問わず、実験的に体積粘性係数を知る方法は、超音波の伝播や減衰を調べる以外にないらしい

現在は21世紀なので、簡単な分子の場合は、実験によらず、分子動力学で計算するという道がありえる。原理的には、shear viscosityも体積粘性係数も、Green-Kubo公式で計算できるはずなので。shear viscosityの計算は色々やられていて、手法もいくつか存在する

[1]Determining the Bulk Viscosity of Rigid Water Models
http://pubs.acs.org/doi/pdf/10.1021/jp211952y

は、分子動力学で、shear viscosityとbulk viscosityを決める論文。水のモデルとしては、TIP4P/2005がよいらしい。これは、他の論文でも、精度高いと書いてある

というわけで、Gromacsを使って、上の論文と同じことができないか試した。Gromacsのバージョンは5.1.2(aptで入れたら、これになった)。マニュアルには、gmx viscosityコマンドで、上の論文の手法が実装されてるよって書いてあるけど、そんなコマンドは存在しない。というわけで、手動で頑張ってデータを解析するしかない

また、水のモデルは、デフォルトでは、TIP4P/2005が入ってない

GROMACS files for the TIP4P/2005 model
http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model

にあるのが、TIP4P/2005モデルらしいけど、面倒なので、普通のTIP4Pを使うことにした。

論文[1]によれば、

At the beginning, a MD simulation was performed at constant temperature and pressure conditions (NPT ensemble) for a period of 2 ns, in order to determine the density of liquid water as it is predicted by each water model.

The average energy ⟨E⟩ of the system at that density was computed from a second constant volume and temperature (NVT ensemble) simulation of 2 ns.

For the calculation of the bulk and shear viscosity at each temperature, we have performed a set of 100 independent constant volume and energy (NVE ensemble) simulations. For each trajectory, initial atomic momenta were sampled from a Maxwell distribution at temperature T.

らしいので、NPTアンサンブルで2ns→NVTアンサンブルで2ns→独立なNVE計算を100回(10psの平衡化の後、300psの計算をやると書いてある)という順番で計算するらしい。NVE計算では、長い時間やると、enerygy driftによって、エネルギー保存が成立しなくなるので短時間の計算を沢山やるということのよう。平衡化する場合、NVT→NPTの順が一般的な気がするけど、今回は逆

論文では、水256分子で計算しているので、適当に250分子を用意した。まあ、大体、以下のような感じのコマンドを叩く。エネルギー最小化は論文には書いてないので、特にやらなくてもいいかもしれない。NVE計算は、面倒なので100回やらずに30回しかやってない。

#-- We will get 250 molecules
gmx solvate -cs tip4p -o water.gro -box 2.1 2.0 2.0 -p conf/topol.top

#-- Energy minimization
gmx grompp -f conf/em1.mdp -o em1 -c water.gro -p conf/topol.top
gmx mdrun -s em1.tpr -deffnm em1 -v 
gmx grompp -f conf/em2.mdp -o em2 -pp em2 -po em2 -c em1 -t em1 -p conf/topol.top
gmx mdrun -s em2.tpr -deffnm em2

##gmx energy -f em1.edr -o min-energy.xvg
##gmx energy -f em2.edr -o min2-energy.xvg 

#-- NPT ensemble (to determine the density of liquid water)
gmx grompp -f conf/npt.mdp -o nptin -c em2.gro -p conf/topol.top -maxwarn 1
gmx mdrun -s nptin.tpr -deffnm nptout

#-- check density ,pressure, energy, temparature
echo -e "7\n8\n10\n14\n15\n" | gmx energy -f nptout.edr -o nptout.xvg

#-- NVT ensemble (to determine the average energy)
gmx grompp -f conf/nvt.mdp -o nvtin -c nptout -t nptout -p conf/topol.top
gmx mdrun -s nvtin.tpr -deffnm nvtout

#-- check energy
echo -e "7\n9\n11\n" | gmx energy -f nvtout.edr -o nvtout.xvg

#-- NVE ensemble
for i in `seq 1 30`;do
  gmx grompp -f conf/nve.mdp -o nvein -c nvtout -t nvtout -p conf/topol.top
  gmx mdrun -s nvein.tpr -deffnm nveout_$i
  echo -e "7\n8\n10\n" | gmx energy -f nveout_$i.edr -o nveout_$i.xvg    #--energy,temperature,pressure
done

使用したmdpファイルとtopol.topは、以下に置いた
https://gist.github.com/vertexoperator/aebb9676509e5eba608dfff6ff74bfc5

NPT計算の結果

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9772.19        4.7    154.428    -16.766  (kJ/mol)
Temperature                 298.147      0.062    11.6618   -0.23701  (K)
Pressure                   -4.63437          2    868.369   -8.01756  (bar)
Volume                      7.51647     0.0023   0.124593 -0.00247449  (nm^3)
Density                     995.296        0.3    16.3923   0.194343  (kg/m^3)

とかなった。密度が大体正しそうなので多分OK。温度は298.15K=25度を指定しているので、そのへんで動く。論文[1]のTable Iでは、25度でTIP4Pを使用した時、0.9953(g/cm^3)という値を得たと書いてあるので、妥当な数値っぽい

NVT計算の結果は

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -9755.22        3.2    134.505    -16.928  (kJ/mol)
Temperature                 298.143       0.11    11.1177   -0.10815  (K)
Pressure                   -386.763        9.7    777.601    18.4497  (bar)

まぁ、よさげ？

Gromacsには、tcafコマンドというのがあって、私が理解してない謎の手法で、shear viscosityを計算できるらしい。
gmx tcaf
http://manual.gromacs.org/programs/gmx-tcaf.html

それで、NVT計算の結果とかに対して、

gmx tcaf -f nvtout.trr -s nvtin.tpr

とかやると、visc_k.xvgに(一部ヘッダ省略)

@    title "Fits"
@    xaxis  label "k (nm\S-1\N)"
@    yaxis  label "\xh\f{} (10\S-3\N kg m\S-1\N s\S-1\N)"
@TYPE xy
@    s0 symbol 2
@    s0 symbol color 1
@    s0 linestyle 0
 3.113 0.51114
 3.269 0.493288
 3.269 0.48944
 4.514 0.386852
 4.514 0.387317
 4.514 0.398109
 4.514 0.389984
 4.623 0.392635
 4.623 0.385249
 5.573 0.329574
 5.573 0.330219
 5.573 0.326714
 5.573 0.325943
 6.226 0.296062
 6.538 0.28417
 6.538 0.284974

などと出力される(デフォルトでは、visc_k.xvgにkとetaのみの出力がされる)。マニュアルによると、これを適当にfittingすると、shear viscosityが計算できるらしい

以下のようなコードで、実際にfittingしてみた

import numpy as np
from scipy import stats

k,y=[],[]
for line in open("visc_k.xvg"):
    if line[0]=="@":continue
    if line[0]=="#":continue
    ls = line.strip().split()
    k.append(float(ls[0]))
    y.append(float(ls[1]))


x=np.array(k)**2
y=np.array(y)
r=stats.linregress(x,y)
print(r.intercept)

計算すると、例えば0.540536022685という値を得たので、0.54e-3(kg/(m s))という感じであるらしい。水の粘性係数は、25度で0.000894(Pa s)らしいので、2/3くらいの値になっているが、論文[1]の方法では、TIP4Pで計算した場合、25度で、(4.83±0.09)e-4 (Pa s)という値になったらしいし、まぁTIP4Pの粘性に関する精度はこんなものであるらしい。NPT,NVEの結果から、tcafで計算しても、同じような結果を得た。TIP4P/2005なら、もっとよい結果が得られるらしい

Green-Kubo公式で、体積粘性係数を計算するために
$B_{ACF}(t) = \langle \delta P(t) \delta P(0) \rangle$
を計算する。
After an equilibration period of 10 ps, a configuration with total energy equal to the ⟨E⟩ decided before was selected. For that configuration and for the next 300 ps, the instantaneous stress tensor elements were stored in the disk every 20 fs for further analysis.
とか書いてある(はNVT計算の結果を使うらしい)けど、に到達するのに、ものすごく時間がかかるっぽかったので無視した(なんか勘違いしてる?)

よくやるように、
$B_{ACF}{(t) = \langle \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} dt' \delta P(t+t') \delta P(t') \rangle$
で計算する。これが得られたら、論文の式(4)でfittingを行う。次のpythonコードが実装例

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


def func(t , C , tau_f , tau_s , omega , beta_f , beta_s):
   fast_term = (1-C)*np.exp(-np.power((t/tau_f) , beta_f))*np.cos(omega*t)
   slow_term = C*np.exp(-np.power(t/tau_s , beta_s))
   return (fast_term + slow_term)


BACF = []
for i in range(1,31):
   pressures = []
   for line in open("nveout_{}.xvg".format(i)):
        if line[0]=="@" or line[0]=="#":continue
        ls  = line.strip().split()
        assert(len(ls)==4)
        t,E,_,P = ls
        pressures.append( float(P) )
   Pbar = sum(pressures)/len(pressures)
   deltaP = np.array([P-Pbar for P in pressures])
   ACF = []  #--auto-correlation function
   for t in range(len(deltaP)):
      if t==0:
          sv = np.sum(deltaP*deltaP)/(len(deltaP))
      else:
          sv = np.sum(deltaP[t:]*deltaP[0:-t])/(len(deltaP) - t)
      ACF.append( sv )
   if len(BACF)==0:
      BACF = np.copy(ACF)
   else:
      BACF = BACF+ACF


dt = 0.02  #-- pico seconds
xdata = np.arange(len(BACF))*dt
ydata = BACF/30

popt , pcov = curve_fit(func , xdata , ydata/ydata[0])
plt.plot(xdata[:1000], func(xdata[:1000], *popt))
plt.show()

#--積分値(total relaxation time)
print( np.sum(func(xdata , *popt)) * dt)

あと、単位に注意する(粘性係数の単位をPa sとなるように合わせる)。上のdtはpsなので、10^{-12}(s)で、Gromacsが出力する圧力の単位は、bar(1bar = 10^5 Ps)。体積はNVT計算の平均値、温度も同様にして、計算する。上のコードで得られた諸々の数値をもとに、

(7.51647e-27/(1.38e-23 * 298.15))*(ydata[0]*1.0e10)*np.sum(func(xdata , *popt)) * dt*1.0e-12

で、求めたい体積粘性係数の値が出る。私の場合は、0.0018268112792585417という値を得た。論文のTABLE2によれば、298KでTIP4Pを使った場合、0.001173という値を得ている。オーダーは合っている。ずれは、アンサンブル数の少なさ、方法の若干の違いによるものと思われる。最初のほうで見た2.4cPという実測値と比較しても、悪くない。論文を信じれば、TIP4P/2005の方が精度はよいらしい。おしまい