3D-RISM Tutorial¶

A solute such as a molecule is treated by density functional theory (DFT) and the surrounding solution by 3D-RISM, and the two are solved self-consistently (3D-RISM-SCF). This page runs the "formaldehyde (HCHO) in water" and "H₂O in NaCl(aq)" calculations of the example PW/examples/RISM_example and looks at the three-dimensional solvent distribution and the solvation free energy. For the solution side (1D-RISM) see the 1D-RISM tutorial, and for the theory see 3D-RISM.

Setup¶

In addition to the 1D-RISM settings, specify the Lennard-Jones (LJ) parameters of each solute atom and how 3D-RISM is solved. The input of the HCHO example is:

 &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 selects the source of the LJ parameters for each solute species (here the OPLS-AA force field). As in the H₂O-in-NaCl(aq) example, the values can also be given directly with solute_epsilon (kcal/mol) and solute_sigma (Å).
3D-RISM is solved in the same periodic cell as the solute (here a 10 Å cube). Make the cell large enough that the solvent distribution is no longer affected by the solute at its boundary.
ecutsolv is the cutoff for the solvent correlation functions.

Run the example as follows. Each of the two calculations takes a little over ten seconds.

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

How 3D-RISM-SCF proceeds¶

At each SCF iteration, 3D-RISM is solved in the electrostatic potential of the current electron density, and the resulting solvent charge distribution is fed back into the next Kohn–Sham equation.

     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

In the HCHO example the SCF converges in 11 iterations, and 3D-RISM is solved at each of them. The 3D-RISM convergence threshold (convergence threshold) is tightened as the SCF approaches convergence, so that 3D-RISM is not solved more precisely than needed in the early SCF iterations.

Three-dimensional solvent distribution¶

Running pprism.x after the SCF writes the solvent distribution in CUBE format.

 &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

This produces the number density of each solvent site (HCHO.aq.3drism_rho#1_O.h2o.cube etc.), the solvent charge density and the potentials. Dividing the number density by its bulk value gives the three-dimensional distribution function \(g(\boldsymbol{r})\). The figure shows a cut through the molecular plane of HCHO.

Distribution of water O and H around HCHO

Water O (left): within about 3 Å of the molecule the solvent is excluded, and the first hydration shell surrounds this region.
Water H (right): water H atoms gather strongly on the side of the HCHO oxygen (red) (\(g>3\)), showing water hydrogen-bonded to the carbonyl oxygen. Almost no H gathers on the side of the HCHO hydrogens.

Solvation free energy¶

At the end of the output, the chemical potential of solvation (solvation free energy) is printed. The value is also included in the total-energy breakdown as 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

System	Solvent	Closure (kcal/mol)	GaussFluct (kcal/mol)
HCHO	water	+3.74	+0.80
H₂O	0.5 M NaCl(aq)	−4.24	−6.25

For H₂O in NaCl(aq) the breakdown by solvent species is also printed (Closure: H₂O −3.67, Na⁺ −1.35, Cl⁻ +0.78 kcal/mol). Closure and GaussFluct evaluate the solvation free energy from the same distribution with different approximate expressions, and their values generally differ. Choose between them according to the quantity and the system you want to compare.

Next¶

3D-RISM treats the solute and the solvent in the same three-dimensional periodic cell, so it cannot represent an interface where only one side faces the solution, such as an electrode surface. Electrode/electrolyte interfaces are handled in the ESM-RISM tutorial, which combines the ESM method with RISM in the Laue representation.