1D-RISM チュートリアル¶
Quantum ESPRESSOのRISM計算では、最初に1D-RISMで溶媒どうしの相関（溶液そのものの構造）を求め、その結果を使って3D-RISMやESM-RISM（Laue-RISM）で溶質のまわりの溶媒分布を求めます。このページでは、例題 PW/examples/RISM_example の「NaCl水溶液中のH2O」の計算のうち、最初に行われる1D-RISMの部分を取り出して見ていきます。理論は1D-RISMを参照してください。
準備¶
RISMは &control の trism = .true. で有効になり、&rism 名前リストと SOLVENTS カードで溶媒を指定します。1D-RISMに関係する部分は次のとおりです。
&control
calculation = 'scf'
trism = .true.
/
&rism
nsolv = 3 ! number of solvent species
closure = 'kh' ! closure equation (Kovalenko-Hirata)
tempv = 300.0 ! temperature (K)
ecutsolv = 144.0 ! cutoff for the solvent (Ry)
! 1D-RISM's setting
starting1d = 'zero'
rism1d_conv_thr = 1.0e-8
rism1d_maxstep = 10000
mdiis1d_size = 20
mdiis1d_step = 0.5
/
SOLVENTS {mol/L}
H2O -1.0 H2O.spc.MOL
Na+ 0.5 Na+.aq.MOL
Cl- 0.5 Cl-.aq.MOL
SOLVENTSの各行は「溶媒の名前・濃度・MOLファイル」です。濃度に負の値を与えると、MOLファイルに書かれた密度が使われます（水は1 g/cm3）。この例はNa+とCl−が0.5 mol/L、つまり0.5 MのNaCl水溶液です。- MOLファイルは溶媒分子の構造と古典力場（各原子の電荷とLennard-Jonesパラメータ）を与えます。水はSPCモデル（水素のLJパラメータを修正したもの）です。例題の
MOL/ディレクトリにあり、pseudo_dirに置いて読み込ませます。
例題は次のように実行できます。この計算は数十秒で終わります。
cd q-e/PW/examples/RISM_example
./run_example_3D-RISM
ここで示す結果は、例題に付属する参照出力と全エネルギーが \(10^{-6}\) Ry以内で一致することを確認しています。
1D-RISMの計算の流れ¶
出力の冒頭に、1D-RISMの設定と収束の様子が書かれます。
1D-RISM info
------------
closure equation = KH
temperature = 300.0000 kelvin
coulomb smearing radius = 2.0000 bohr
number of solvent sites = 5
number of grids = 5008
...
iter. # 1 RMS(g-h-1)= 1.104E-03 nbox= 0 T= 2400.00K
iter. # 101 RMS(g-h-1)= 5.513E-05 nbox= 5 T= 711.11K
iter. # 201 RMS(g-h-1)= 6.915E-05 nbox= 4 T= 316.05K
iter. # 273 RMS(g-h-1)= 9.291E-09 nbox= 20 T= 300.00K
convergence has been achieved in 273 iterations
溶媒のサイトは水のO・H、Na+、Cl−の5種類です。反復の途中で温度 T が2400 Kから目標の300 Kまで下がっていきます。相互作用が相対的に弱い高温から解き始めて徐々に冷やすことで、イオンを含む溶液でも安定に収束させています。
出力ファイル *.1drism¶
1D-RISMで得られたサイト間の対分布関数 \(g(r)\) が prefix.1drism に書き出されます。
# Pair distribution function
# Gvv(r)
#
H2O:H2O H2O:Na+ ...
r/Angstrom O/h2o:O/h2o O/h2o:H/h2o H/h2o:H/h2o O/h2o:Na ...
1列目が距離 \(r\)（Å）、2列目以降がサイトの組ごとの \(g(r)\) です。\(g(r)=1\) はバルクの平均密度、1より大きいところはその距離にサイトが集まりやすいことを表します。
水の構造¶
- O–O：2.96 Åに高いピーク（\(g\approx2.4\)）があり、これが水分子の第1水和殻です。
- O–H（異なる分子どうし）：1.8 Å付近の小さなピークが水素結合に対応します。
- 遠方ではどの曲線も1に近づき、一様な液体になります。
イオンの水和¶
- Na+：Na+–Oのピーク（2.32 Å、\(g\approx4.0\)）がNa+–H（3.06 Å）より内側にあります。陽イオンには水の酸素が向いています。
- Cl−：Cl−–Hのピーク（2.22 Å）がCl−–O（3.49 Å）より内側にあります。陰イオンには水の水素が向いています。
イオンの電荷の符号によって水分子の向きが逆になる様子が、分布関数から読み取れます。
溶媒和の化学ポテンシャル¶
1D-RISMの後に、溶液中での各溶媒種の溶媒和の化学ポテンシャル（過剰化学ポテンシャル）が出力されます。
Chemical potential of solvation
H2O H2O Closure -0.153796E+01 kcal/mol
...
Total Closure -0.229340E+01 kcal/mol
GaussFluct -0.583926E+01 kcal/mol
Na+ H2O Closure -0.446607E+02 kcal/mol
Na+ Closure 0.175759E+06 kcal/mol
Cl- Closure -0.175796E+06 kcal/mol
Total Closure -0.820464E+02 kcal/mol
| 溶媒種 | Closure（kcal/mol） | GaussFluct（kcal/mol） |
|---|---|---|
| H2O | −2.29 | −5.84 |
| Na+ | −82.0 | −84.3 |
| Cl− | −76.4 | −79.2 |
Closureは閉包式（ここではKH）に対応する式、GaussFluctはガウスゆらぎ（GF）近似の式で評価した値です。- イオンどうしの成分（Na+–Na+、Na+–Cl− など）は個別には大きな値になりますが、互いに打ち消し合います。意味を持つのは合計（
Total）です。