3D-RISM: 3D reference Interaction Site Model

Introduction

The three-dimensional reference interaction site model with the closure relation by Kovalenko and Hirata (3D-RISM-KH) provides the solvent structure in the form of a 3D site distribution function, \(g_{\gamma}^{UV}(r)\), for each solvent site, \(\gamma\). Note that the use of 3D-RISM as implemented in ADF is an expert option.

Warning

We have no workflow that consistently produces accurate results.

It enables, at modest computational cost, the calculations of thermodynamics, electronic properties and molecular solvation structure of a solute molecule in a given molecular liquid or mixture. Using 3D-RISM, one can study chemical reactions, including reaction coordinates and transition state search, with the molecular solvation described from the first principles. The method yields all of the features available by using other solvation approaches. The 3D-RISM part of ADF has not been parallelized.

Details on the implementation of 3D-RISM-KH in ADF can be found in Ref. [1], with applications in Ref. [2]. The theory of 3D-RISM-KH in combination with DFT can be found in Refs. [5] [6] [7] [4]. A combination of 3D-RISM-KH with FDE (frozen-density embedding) can be found in Ref. [3].

Similar to explicit solvent simulations, 3D-RISM properly accounts for chemical peculiarities of both solute and solvent molecules, such as hydrogen bonding and hydrophobic forces, by yielding the 3D site density distributions of the solvent. Moreover, it readily provides, via analytical expressions, all of the solvation thermodynamics, including the solvation free energy potential, its energetic and entropic decomposition, and partial molar volume and compressibility. The expression for the solvation free energy (and its derivatives) in terms of integrals of the correlation functions follows from a particular approximation for the so-called closure relation used to complete the integral equation for the direct and total correlation functions.

Solvation free energy

The solvation free energy can be calculated as:
Solvation Free Energy =
[ (Total Bonding Energy) + (Internal Energy)
- T * (Entropy) + (Excess Chemical Potential) ]_3D-RISM
- [ (Total Bonding Energy) + (Internal Energy) - T * (Entropy) ]_Gas-Phase
If one assumes that the internal energy and the vibrational and rotational entropy of the solute molecule is the same in solution and in gas phase, then this simplifies to:
Solvation Free Energy =
[ (Total Bonding Energy) + (Excess Chemical Potential) ]_3D-RISM
- [ (Total Bonding Energy)]_Gas-Phase

However, a formally accurate calculation should include the difference between the thermal corrections from frequency calculations produced by ADF in the SCF calculation with 3D-RISM-KH solvation and in gas phase.

Input

When performing 3D-RISM simulations, each atom in the ATOMS block must have two parameters specified, SigU and EpsU, for example:

ATOMS
  C     0.00    0.00    0.00     SigU=3.50    EpsU=0.066
  ...
END

The SigU and EpsU parameters have the same meaning as Sigma_alpha and Eps_alpha for atoms of the solvent in the SOLVENT sub-block below. They can be obtained from a Lennard-Jones force-field parameter sets.

All 3D-RISM-related input keys are contained in a RISM input block. Below, only the mandatory keywords are shown. Optional keywords are described in the next section.

RISM title
  RISM1D
    FLUIDPARAM temper=298. DielConst=78.497 UTotDens=1/A3 0.03333
  SUBEND
  SOLVENT ArbitrarySolventName
    UNITS uWeight=g/mol  ULJsize=A  ULJenergy=kcal/mol Ucoord=A Udens=1/A3
    PARAMETERS Weight NAtomTypes
      N1    Z_alpha1   Sigma_alpha1   Eps_alpha1   X1_1 Y1_1 Z1_1
                                                   X1_2 Y1_2 Z1_2
                                                   ...  ...  ...
      N2    Z_alpha2   Sigma_alpha2   Eps_alpha2   X2_1 Y2_1 Z2_1
                                                   X2_2 Y2_2 Z2_2
                                                   ...  ...  ...
      ...
      DENSPE=density
  SUBEND
  SOLUTE ArbitrarySoluteName
    BOXSIZE   sizeX  sizeY  sizeZ
    BOXGRID    npX    npY    npZ
  SUBEND
END

The RISM1D sub-block contains general parameters for the 1D-RISM calculation of the solvent(s). Even though all RISM1D sub-keys have reasonable defaults, the FLUIDPARAM sub-key deserves a special attention because its default values are only applicable if the solvent is water. Thus, you may need to change some of these values when modeling a different solvent, at least the dielectric constant and the density. Note that even when using all default values from the RISM1D sub-block the sub-block itself must be specified, even if empty. See below for complete description of the RISM1D sub-block.

The SOLVENT sub-block can be repeated if the solvent is a mixture. Each SOLVENT sub-block contains parameters for one solvent. First, each solvent has a name, which is specified on the SOLVENT keyword’s line and is arbitrary. The first line in the SOLVENT sub-block must contain the UNITS key. You should leave it at the default values. Then follow the actual solvent parameters. In principle, each solvent consists of multiple atoms and functional groups. For simplicity, we will call each of them an atom. For example, in 3D-RISM therms, methanol consists of 3 “atoms”: CH3 , O, and H. Each such atom has a set of three parameters, shared between all atoms of the same type, and the coordinates. These parameters follow the PARAMETERS keyword. The line with the PARAMETERS keyword itself must specify the molecular weight of the solvent and the number of atom types that follow. The first line for each atom type contains, in this order: number of atoms of this type, \(z_\alpha\) , \(\sigma_\alpha\) , \(\epsilon_\alpha\) , three coordinates for the first atom of this type. If there is more than one atom of this type then the coordinates for the 2nd and other atoms follow on subsequent lines. The SOLVENT sub-block is concluded by the specific density of this solvent, by default, in molecules per cubic angstrom. This number should be equal to the total density for mono-component solvents.

For example, the SOLVENT block for water would typically look as:

SOLVENT water
  UNITS      uWeight=g/mol  ULJsize=A  ULJenergy=kcal/mol Ucoord=A Udens=1/A3
  PARAMETERS Weight=18.015   nAtoms=2
    1    -0.8476   3.166    0.1554       0.000000  0.00000  0.000000
    2     0.4238   1.000    0.0460      -0.816490  0.00000  0.577359
                                         0.816490  0.00000  0.577359
  DENSPE=0.03333
SUBEND

The SOLUTE sub-block specifies 3D-RISM parameters for your molecule. The BOXSIZE and BOXGRID sub-keys specify dimensions of the simulation box, in Angstrom, and the number of points of grid in each direction. The box should be twice as large as the molecule and the BOXGRID values must be a power of 2. The size/np ratio defines the grid spacing in each direction and this should be not larger than 0.5 angstrom.

The optional 3D-RISM keys for the RISM1D and SOLUTE sub-blocks are listed below together with their defaults.

RISM1D Theory=RISM Closure=KH
! optional RISM1D subkeys with their default values:
  FLUIDPARAM temper=298. DielConst=78.497 UTotDens=1/A3 0.03333
  OUTPUT PrintLev=5 File=solvent OutList=GXT
  GRID   Nr=8192  dR=0.025 MaxRout=100.0 MaxKout=0.0
  MDIIS  N=20 Step=0.5 Tolerance=1.e-12
  ELSTAT LRsmear=1.0 Adbcor=0.5
  ITER  Ksave=-1 Kshow=1 Max=10000
SUBEND
OUTPUT

PrintLev - print level;

File - common base name for output files;

OutList - which information will be written to output files: G - distribution function, C - direct correlation function, H - total correlation function, T - thermodynamics.

GRID

Nr - the size of the 1D-RISM grid, must be a power of 2;

dR - mesh size in Angstrom;

MaxRout - plot range in direct space;

MaxKOut - plot range in reciprocal space.

MDIIS

N - number of vectors in the DIIS space;

Step - step size;

Tolerance - convergence criterion.

ELSTAT
LRsmear - smearing parameters for coulomb potential;
ITER

Ksave - save the current solution every Ksave steps (0 - do not save);

Kshow - print convergence progress every Kshow steps;

Max - maximum number of iterations.

FLUIDPARAM

Temper - temperature;

DielConst - dielectric constant;

UTotDens - units for total density: g/cm3, kg/m3, 1/A3 are valid units followed by the density value.

SOLUTE ArbitrarySoluteName
  outlist=M  closure=KH  xvvfile=solvent.xvv outfile=rism3d
  Nis=10  DELOZ=0.5 TOLOZ=2.0e-6
  Ksave=-1  Kshow=1  Maxste=10000
  Output=4
  CHRGLVL=MDCq
SUBEND

Outlist - output requested: G - distribution function, C - direct correlation function, H - total correlation function

Closure - closure for the 3D problem: KH – Kovalenko-Hirata, HNC - hypernetted chain approximations, PY - Perkus-Yevik

Xvvfile - name of the file with the results of the 1D-RISM calculation specified in the RISM1D keyword above, with .xvv appended to it;

Outfile - name of the output text file;

Output - print level;

CHRGLVL - which charges computed by ADF to use. This can be MDCq (default), MDCd, MDCm, or EXACT. The Nis, DELOZ, and TOLOZ have the same meaning for 3D-RISM as parameters of the MDIIS keyword of the RISM1D block. Likewise, Ksave, Kshow, and Maxste are analogous to the parameters of the ITER key in RISM1D.

Parameters for some solvents

Table 5 Water Weight=18.015 nAtoms=3
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
O -0.8476 3.166 0.1554 -0.0646 0.0 0.0
H 0.4238 0.7 0.046 0.5127 0.8165 0.0
H 0.4238 0.7 0.046 0.5127 -0.8165 0.0
Table 6 Methanol Weight=32.042 nAtoms=3
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
CH3 0.265 3.775 0.207 -1.43 0.0 0.0
O -0.7 3.07 0.17 0.0 0.0 0.0
H 0.435 0.7 0.046 0.2998 0.8961 0.0
Table 7 Ethanol
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
CH3 0.0 3.775 0.207 -1.9028 -1.4551 0.0
CH2 0.265 3.905 0.118 -1.43 0.0 0.0
O -0.7 3.07 0.17 0.0 0.0 0.0
H 0.435 0.7 0.046 0.2998 0.8961 0.0
Table 8 Acetonitrile Weight=41.0520 nAtoms=3
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
CH3 0.269 3.6 0.38 1.46 0.0 0.0
C 0.129 3.4 0.099 0.0 0.0 0.0
N -0.398 3.3 0.099 -1.17 0.0 0.0
Table 9 Dimethylsulfoxide Weight=78.135 nAtoms=4
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
CH3 0.160 3.81 0.16 0.0 1.7167 -0.5413
CH3 0.160 3.81 0.16 -1.6653 -0.4168 -0.5413
S 0.139 3.56 0.395 0.0 0.0 0.0
O -0.459 2.93 0.28 0.0 0.0 1.53
Table 10 Isopropanol Weight=60.09502 nAtoms=5
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
O -0.700 3.0700 0.1700 1.32393 0.14660 -0.01452
H 0.435 0.7000 0.0460 1.84501 0.01837 -0.79238
CH 0.265 3.8500 0.0800 0.84520 1.49406 -0.00612
CH3 0.000 3.7750 0.2070 1.38816 2.18026 1.24895
CH3 0.000 3.7750 0.2070 -0.68415 1.45191 -0.02015
Table 11 Benzene
Atom \(z_\alpha\) / a.u. \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol X/Å Y/Å Z/Å
C -0.115 3.55 0.07 1.4 0.0 0.0
C -0.115 3.55 0.07 0.7 -1.2124 0.0
C -0.115 3.55 0.07 -0.7 -1.2124 0.0
C -0.115 3.55 0.07 -1.4 0.0 0.0
C -0.115 3.55 0.07 -0.7 1.2124 0.0
C -0.115 3.55 0.07 0.7 1.2124 0.0
H 0.115 2.42 0.03 2.48 0.0 0.0
H 0.115 2.42 0.03 1.24 -2.1477 0.0
H 0.115 2.42 0.03 -1.24 -2.1477 0.0
H 0.115 2.42 0.03 -2.48 0.0 0.0
H 0.115 2.42 0.03 -1.24 2.1477 0.0
H 0.115 2.42 0.03 1.24 2.1477 0.0

OPLS-AA Parameters for common atoms and atom group

The table below contains sigma and epsilon parameters for some of the most common solvent groups collected from ref. [8] [9] [10]. These parameters are kindly provided by Leonardo Costa.

Table 12 All atoms model
Atom Example \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol Reference
C (SP3) methane 3.73 0.294 433
C (SP3) neopentane 3.80 0.050 433
C (SP2) isobutene 3.75 0.105 433
C (SP) acetonitrile 3.40 0.099 434
C (arom) benzene 3.55 0.07 434
C chloroform 4.10 0.05 435
H O—H 0.7 0.046 434
H hydrocarbons 2.42 0.03 434
O water 3.166 0.1554 434
O alcohol 3.07 0.17 434
O sulfoxide 2.93 0.28 434
N (SP) acetonitrile 3.3 0.099 434
S sulfoxide 3.56 0.395 434
Cl chloroform 3.40 0.300 435
Table 13 United atom model
Atom Example \(\sigma_\alpha\) / Å \(\epsilon_\alpha\) / kcal/mol Reference
CH (SP3) isobutane 3.850 0.080 433
CH (SP2) 2-butene 3.800 0.115 433
CH (arom) benzene 3.750 0.110 433
CH2 (SP3) n-butane 3.905 0.118 433
CH2 (SP2) 1-butene 3.850 0.140 433
CH3 hydrocarbon 3.775 0.207 433
CH3 acetonitrile 3.6 0.38 434
CH3 sulfoxide 3.81 0.16 434
NH2 amine 3.3 0.17 434

References

[1]S. Gusarov, T. Ziegler, and A. Kovalenko, Self-Consistent Combination of the Three-Dimensional RISM Theory of Molecular Solvation with Analytical Gradients and the Amsterdam Density Functional Package, Journal of Physical Chemistry A 110, 6083 (2006)
[2]D. Casanova, S. Gusarov, A. Kovalenko, and T. Ziegler, Evaluation of the SCF Combination of KS-DFT and 3D-RISM-KH; Solvation Effect on Conformational Equilibria, Tautomerization Energies, and Activation Barriers, Journal of Chemical Theory and Computation 3, 458 (2007)
[3]J.W. Kaminski, S. Gusarov, A. Kovalenko, T.A. Wesolowski, Modeling solvatochromic shifts using the orbital-free embedding potential at statistically mechanically averaged solvent density, Journal of Physical Chemistry A 114, 6082 (2010)
[4]A. Kovalenko and F. Hirata, Potentials of mean force of simple ions in ambient aqueous solution. II. Solvation structure from the three-dimensional reference-interaction site model approach, and comparison with simulations, Journal of Chemical Physics 112, 10403 (2000)
[5]A. Kovalenko and F. Hirata, Self-consistent description of a metal-water interface by the Kohn-Sham density functional theory and the three-dimensional reference-interaction site model, Journal of Chemical Physics 110, 10095 (1999)
[6]A. Kovalenko and F. Hirata, Potentials of mean force of simple ions in ambient aqueous solution. I. Three-dimensional reference-interaction site model approach, Journal of Chemical Physics 112, 10391 (2000)
[7]A. Kovalenko, Three-dimensional RISM theory for molecular liquids and solid-liquid interfaces, In Molecular Theory of Solvation; Hirata, Fumio, Ed.; Understanding Chemical Reactivity (series); Mezey, Paul G., Series Ed.; Kluwer Acadamic Publishers: Dordrecht, The Netherlands, 2003; Vol. 24, pp 169-275.
[8]W.L. Jorgensen, J.D. Madura, C.J. Swenson, Optimized intermolecular potential functions for liquid hydrocarbons, Journal of the American Chemical Society 106, 6638 (1984)
[9]A.E. Kobryn, A. Kovalenko, Molecular theory of hydrodynamic boundary conditions in nanofluidics, Journal of Chemical Physics 129, 134701 (2008)
[10]O. Acevedo, W.L. Jorgensen, Influence of Inter- and Intramolecular Hydrogen Bonding on Kemp Decarboxylations from QM/MM Simulations, Journal of the American Chemical Society 127, 8829 (2005)