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.

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. [300], with applications in Ref. [301]. The theory of 3D-RISM-KH in combination with DFT can be found in Refs. [302-304, 349]. A combination of 3D-RISM-KH with FDE (frozen-density embedding) can be found in Ref. [305].

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”: CH₃ , 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

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 [433, 434, 435]. These parameters are kindly provided by Leonardo Costa.