3D-RISM: 3D reference Interaction Site Model¶

Introduction

The three-dimensional reference interaction site model (3D-RISM) provides the solvent structure in the form of a 3D site distribution function, \(g_{\gamma}^{UV}(r)\), for each solvent site, \(\gamma\).

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 based on a molecular-mechanics type description of the solvent. The method yields all of the features available by using other solvation approaches.

The 3D-RISM implementation has been revised and supplemented by a new option that uses the fitted electron density in the interaction between solvent and solute directly. Without the use of point charges, numerical stability has been improved. Using this option, 3D-RISM can now be used for structure optimization, (numerical) frequencies, vertical excitation energies and magnetic properties. The code has also been parallelized. Details can be found in Ref. [1]. Starting from AMS2022, new point charge based options are available to be used in combination with structure optimization, numerical frequencies and vertical excitation energies.

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

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, partial molar volume and compressibility. Starting from AMS2023, enthalpic and entropic decomposition of the solvation free energy can be obtained approximately using analytic expression (which are exact for CHRGLVL=MM or pure Lennard-Jones fluids), also including the macroscopic dependence of the dielectric constant on temperature and solvent density. [9] 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 solvation and in gas phase.

Input

When performing 3D-RISM simulations, each atom in the ATOMS block in the AMS part of the input must have two Lennard-Jones parameters specified, adf.SigU and adf.EpsU, for example:

System
   ATOMS
     C     0.00    0.00    0.00     adf.SigU=3.50    adf.EpsU=0.066    adf.ChgU=-0.4000
     ...
   END
End

They can be obtained from a Lennard-Jones force-field parameter set. Atom centered point charges adf.ChgU can optionally be set as well but will only be used in combination with CHRGLVL=MM (see below).

All other 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=weight nAtoms=NSiteTypes
      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 preceding 1D-RISM calculation of the solvent(s) to create the bulk susceptibility function. Even though all technical RISM1D sub-keys have reasonable defaults, the FLUIDPARAM sub-key will most likely require attention because its default values for density and dielectric constant are only applicable if the solvent is water. Thus, you may need to change these values when modeling a different solvent. 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 for each component 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. Usually, only the Udens part should be changed. Then follow the actual solvent parameters. In principle, each solvent consists of multiple atoms and functional groups (so called sites). For example, in 3D-RISM terms, methanol consists of 3 sites: CH₃ , O, and H. Each such site has a set of three parameters, shared between all sites 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 site types that follow. The first line for each site type contains, in this order: number of equivalent sites of this type, \(z_\alpha\) (the point charge assigned) , \(\sigma_\alpha\) , \(\epsilon_\alpha\) (the Lennard-Jones parameters assigned), three coordinates for the first site of this type. If there is more than one site of this type then the coordinates for the 2nd and other sites 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. Using Udens=MOL, the specific density needs to be given in molar fractions, and should therefore be one 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.064605
    2     0.4238   1.000    0.0560      -0.816495  0.00000  0.512747
                                         0.816495  0.00000  0.512747
  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 at least be twice as large as the molecule and the BOXGRID values must be even numbers. Due to the 3D-FFT routines, the largest prime factor of the BOXGRID values should be 7 for optimal performance. The size/np ratio defines the grid spacing in each direction and this should be not larger than 0.5 angstrom for energy calculations and not larger than 0.35 angstrom for structure optimizations. When the BOXSIZE and BOXGRID keywords are not specified the grid is calculated based on the GRIDSTEP and BOXEXTENT values, see below. The number of grid points along the Cartesian coordinate i will then be determined as \(N_i = (range(i) + 2*BOXEXTENT)/GRIDSTEP + 1\), where range(i) is the span of Cartesian coordinate i of the solute molecule. The number of grid points may be increased to be a product of prime numbers not larger than 7.

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

! optional RISM1D subkeys with their default values:
RISM1D Theory=RISM Closure=KH NOSPEEDUP
  FLUIDPARAM temper=298. DielConst=78.497 UTotDens=1/A3 0.03333
  OUTPUT PrintLev=5 File=solvent
  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
  DODEPS alpha=-1.0 kappa=-1.0 dedtp=-1.0 dedpt=-1.0
  Restart1D=restartfile
SUBEND

RISM1D

Theory - Version of 1D-RISM used: RISM - extended RISM, DRISM - dielectrically consistent RISM (should be used for solvents with a static dipole moment);

Closure - Closure for the 1D problem: KH - Kovalenko-Hirata, HNC - hypernetted chain, PSEn - partial series expansion of HNC up to order n (n=1 is identical to KH), PY - Perkus-Yevik, KGK - Kobryn-Gusarov-Kovalenko;

NOSPEEDUP - Switch off the automatic adjustment of the MDIIS step width.

FLUIDPARAM

Temper - temperature;

DielConst - dielectric constant (only used when theory=DRISM);

UTotDens - units for total density: G/CM3, KG/M3, 1/A3 are valid units followed by the density value.

OUTPUT

PrintLev - print level;

File - common base name for output files;

OutList - which information will be written to additional output files (no additional information will be written by default): G - distribution function, C - direct correlation function, H - total correlation function, U - interaction function, X - “.xvv” file (Only required up to ADF2019), T - thermodynamics, D - density and temperature derivatives of all correlation functions (automatically set, when DERIV option is chosen in SOLUTE block).

GRID

Nr - the size of the 1D-RISM grid, must be a power of 2 (max. value: 16384);

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;

Abdcor - switching parameter for dielectric correction (only used if theory=DRISM).

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.

DODEPS

If this block is set, the temperature and density dependence of the dielectric constant is included in the calculation of temperature and density derivatives (OutList=D);

Alpha - isobaric thermal expansion coefficient of the solvent (in 1/K, not set by default!);

Kappa - isothermal compressibility of the solvent (in 1/Pa, not set by default!);

dedtp - temperature derivative of the dielectric constant at constant pressure (in 1/K, not set by default!);

dedpt - pressure derivative of the dielectric constant at constant temperature (in 1/Pa, not set by default!).

Restart1D - use the given restart file (of the RISMDATA type) to restart the 1D-RISM calculation.

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=ZLM
  GRIDSTEP=0.33
  BOXEXTENT=12.0
SUBEND

Outlist - output requested: M - detailed output for the excess chemical potential , F - some corrections (GF, PC+, UC, …) to the excess chemical potential , G - distribution function, C - direct correlation function, H - total correlation function;

Closure - closure for the 3D problem: KH – Kovalenko-Hirata, HNC - hypernetted chain approximations, PSEn - partial series expansion of HNC up to order n (n=1,2,…; n=1 is identical to KH), KGK - Kobryn-Gusarov-Kovalenko;

RBC - Use the current closure in combination with the repulsive bridge correction (incompatible with analytical gradients);

UCA - The value of the a-parameter in the universal correction (a*V + b, mandatory to compute the UC correction);

UCB - The value of the b-parameter in the universal correction (a*V + b, mandatory to compute the UC correction);

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

Outfile - name of the output text files;

Output - print level;

CHRGLVL - which charges computed by ADF to use. This can be ZLM - use the electron density (default), CHELPG - use point charges fitted to the electron density in analogy to the CHELPG procedure, RESP - same as CHELPG but with a harmonic restraint to keep the charges small, MM - use constant point charges given in the System block, MDCq, MDCd, MDCm, or EXMDC (labeled EXACT in older versions). The MDCq, MDCd, MDCm or EXMDC levels can not be combined with analytical gradients;

OLDZLM - reproduce the exact ZLM results obtained with ADF2019 and ADF2020;

RESPWEIGHT - pre-factor of the harmonic restraint in the RESP-like procedure, default value: 0.05;

VDWTYPE - functional form of the van der Waals interaction requested: 12-6 (default), 9-6;

SPHINT - calculate the spherical averaged distribution function around the given atom: 1,2,… - atom in input order , 0 - Center of mass;

SOLVPOS - find the most likely positions of a solvent site: 1,2,… - site in input order , 0 - all sites;

OLDCODE - run the old version of the code (ADF2019 and older), should only be used by experts;

DERIV - compute the temperature and density derivative of the excess chemical potential analytically (only exact for CHRGLVL=MM or pure Lennard-Jones fluids). When DODEPS is set in the RISM1D block, enthalpy and entropy of solvation are calculated as well;

NOSPEEDUP - switch off the automatic adjustment of the MDIIS step width;

GRIDSTEP - Grid spacing, in Angstrom;

BOXEXTENT - Extent of the grid, in Angstrom, beyond the space occupied by atoms in each direction.

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

The following, compiled set of solvent parameters is taken from the literature (see Ref. [1] and references therein for these and additional solvents).

OPLS-AA Parameters for common atoms and atom group

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

References