3D-RISM チュートリアル¶

分子などの溶質を密度汎関数理論（DFT）で、まわりの溶液を3D-RISMで扱い、両者を自己無撞着に解きます（3D-RISM-SCF）。例題 PW/examples/RISM_example の「水中のホルムアルデヒド（HCHO）」と「NaCl水溶液中のH₂O」を計算し、溶媒の3次元分布と溶媒和自由エネルギーを見ていきます。溶液側の準備（1D-RISM）は1D-RISMチュートリアル、理論は3D-RISMを参照してください。

準備¶

1D-RISMの設定に加えて、溶質の各原子のLennard-Jones（LJ）パラメータと、3D-RISMの解き方を指定します。HCHOの例題の入力は次のとおりです。

 &control
    calculation = 'scf'
    trism       = .true.
 /
 &rism
    nsolv    = 1
    closure  = 'kh'
    tempv    = 300.0  ! Kelvin
    ecutsolv = 144.0  ! Rydberg
    ! Lennard-Jones
    solute_lj(1) = 'opls-aa'
    solute_lj(2) = 'opls-aa'
    solute_lj(3) = 'opls-aa'
 /
ATOMIC_SPECIES
 C  -1.0  C.pbe-van_ak.UPF
 O  -1.0  O.pbe-van_ak.UPF
 H  -1.0  H.pbe-van_ak.UPF
CELL_PARAMETERS angstrom
   10.00000  0.00000  0.00000
    0.0000  10.00000  0.00000
    0.0000   0.00000 10.00000
ATOMIC_POSITIONS {angstrom}
 C  5.000  5.609  5.000
 O  5.000  4.425  5.000
 H  4.075  6.191  5.000
 H  5.925  6.191  5.000
K_POINTS {gamma}
SOLVENTS {g/cm^3}
 H2O  1.0  H2O.spc.MOL

solute_lj は溶質の原子種ごとのLJパラメータの出典を選びます（ここではOPLS-AA力場）。NaCl水溶液中のH₂Oの例題のように、solute_epsilon（kcal/mol）と solute_sigma（Å）で値を直接与えることもできます。
3D-RISMは溶質と同じ周期的なセル（ここでは10 Åの立方体）の中で解きます。セルは、溶媒の分布が溶質の影響を受けなくなる程度に大きくとります。
ecutsolv は溶媒の相関関数を表すカットオフです。

例題は次のように実行できます。2つの計算はそれぞれ十数秒で終わります。

cd q-e/PW/examples/RISM_example
./run_example_3D-RISM

3D-RISM-SCFの流れ¶

SCFの各反復で、その時点の電子密度がつくる静電ポテンシャルのもとで3D-RISMを解き、得られた溶媒の電荷分布を次のKohn–Sham方程式に戻します。

     3D-RISM Calculation

     convergence threshold    = 5.012E-03
     iter. #     1  RMS(g-h-1)= 8.113E-01  nbox=  0
     iter. #    11  RMS(g-h-1)= 2.247E-02  nbox= 10
     iter. #    21  RMS(g-h-1)= 9.501E-03  nbox= 10
     iter. #    25  RMS(g-h-1)= 4.317E-03  nbox= 10

     End of 3D-RISM calculation

HCHOの例ではSCFが11回の反復で収束し、そのたびに3D-RISMが解かれています。3D-RISMの収束判定（convergence threshold）は、SCFが収束に近づくにつれて厳しくなります。SCFの初期に3D-RISMを必要以上に精密に解かないための工夫です。

溶媒の3次元分布¶

SCFの後に pprism.x を実行すると、溶媒の分布をCUBE形式で書き出せます。

 &inputpp
    prefix  = 'HCHO.aq'
    outdir  = './tmp/'
    lpunch  = .true.
 /
 &plot
    iflag         = 3  ! 3D plot
    output_format = 6  ! Gaussian's CUBE
 /

pprism.x < HCHO.aq.post.in > HCHO.aq.post.out

溶媒の各サイトの数密度（HCHO.aq.3drism_rho#1_O.h2o.cube など）、溶媒の電荷密度、ポテンシャルのファイルができます。数密度をバルクの値で割ると、3次元の分布関数 \(g(\boldsymbol{r})\) になります。HCHOの分子面で切った断面を示します。

HCHOのまわりの水のOとHの分布

水のO（左）：分子のまわり約3 Åは溶媒が入れない領域で、その外側を第1水和殻が取り囲んでいます。
水のH（右）：HCHOの酸素（赤）の側に、水のHが強く集まっています（\(g>3\)）。水がカルボニル酸素に水素結合している様子です。HCHOの水素の側では、Hはほとんど集まりません。

溶媒和自由エネルギー¶

出力の最後に、溶媒和の化学ポテンシャル（溶媒和自由エネルギー）が書かれます。この値は、全エネルギーの内訳の solvation energy (RISM) にも含まれています。

     Chemical potential of solvation
     H2O        Closure      0.374101E+01 kcal/mol
                GaussFluct   0.799804E+00 kcal/mol
     Total      Closure      0.374101E+01 kcal/mol
                GaussFluct   0.799804E+00 kcal/mol
     ...
     solvation energy (RISM)   =       0.01193133 Ry

系	溶媒	Closure（kcal/mol）	GaussFluct（kcal/mol）
HCHO	水	+3.74	+0.80
H₂O	0.5 M NaCl水溶液	−4.24	−6.25

NaCl水溶液中のH₂Oでは、溶媒種ごとの内訳も出力されます（Closure：H₂O −3.67、Na⁺ −1.35、Cl⁻ +0.78 kcal/mol）。Closure と GaussFluct は、同じ分布から異なる近似式で評価した溶媒和自由エネルギーで、両者の値は一般に異なります。どちらを使うかは、比べたい量と系に応じて選びます。

次へ¶

3D-RISMは溶質と溶媒を同じ3次元周期セルで扱うため、電極表面のように片側だけが溶液に接する界面は表せません。電極/電解液界面は、ESM法とLaue表示のRISMを組み合わせたESM-RISMチュートリアルで扱います。