3DRISM: 3D reference Interaction Site Model¶
Introduction
The threedimensional reference interaction site model (3DRISM) 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 3DRISM, one can study chemical reactions, including reaction coordinates and transition state search, with the molecular solvation described from the first principles based on a molecularmechanics type description of the solvent. The method yields all of the features available by using other solvation approaches.
The 3DRISM 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, 3DRISM 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 3DRISMKH in ADF can be found in Ref. 2, with applications in Ref. 3. The theory of 3DRISMKH in combination with DFT can be found in Refs. 4 5 6 7 . A combination of 3DRISMKH with FDE (frozendensity embedding) can be found in Ref. 8.
Similar to explicit solvent simulations, 3DRISM 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 LennardJones 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 socalled 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) ]_3DRISM [ (Total Bonding Energy) + (Internal Energy)  T * (Entropy) ]_GasPhase
 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) ]_3DRISM [ (Total Bonding Energy)]_GasPhase
However, a formally accurate calculation should include the difference between the thermal corrections from frequency calculations produced by ADF in the SCF calculation with 3DRISM solvation and in gas phase.
Input
When performing 3DRISM simulations, each atom in the ATOMS block in the AMS part of the input must have two LennardJones 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 LennardJones forcefield 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 3DRISMrelated 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 subblock contains general parameters for the preceding 1DRISM calculation of the solvent(s) to create the bulk susceptibility function. Even though all technical RISM1D subkeys have reasonable defaults, the FLUIDPARAM subkey 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 subblock the subblock itself must be specified, even if empty. See below for complete description of the RISM1D subblock.
The SOLVENT subblock can be repeated for each component if the solvent is a mixture. Each SOLVENT subblock 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 subblock 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 3DRISM terms, methanol consists of 3 sites: CH_{3} , 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 LennardJones 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 subblock 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 monocomponent solvents. Using Udens=MOL, the specific density needs to be given in molar fractions, and should therefore be one for monocomponent 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 subblock specifies 3DRISM parameters for your molecule. The BOXSIZE and BOXGRID subkeys 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 3DFFT 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 3DRISM keys for the RISM1D and SOLUTE subblocks 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.e12
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 1DRISM 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  KovalenkoHirata, HNC  hypernetted chain, PSEn  partial series expansion of HNC up to order n (n=1 is identical to KH), PY  PerkusYevik, KGK  KobrynGusarovKovalenko;
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 1DRISM 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 1DRISM calculation.
SOLUTE ArbitrarySoluteName
outlist=M closure=KH xvvfile=solvent.xvv outfile=rism3d
Nis=10 DELOZ=0.5 TOLOZ=2.0e6
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 – KovalenkoHirata, HNC  hypernetted chain approximations, PSEn  partial series expansion of HNC up to order n (n=1,2,…; n=1 is identical to KH), KGK  KobrynGusarovKovalenko;
RBC  Use the current closure in combination with the repulsive bridge correction (incompatible with analytical gradients);
UCA  The value of the aparameter in the universal correction (a*V + b, mandatory to compute the UC correction);
UCB  The value of the bparameter in the universal correction (a*V + b, mandatory to compute the UC correction);
Xvvfile  name of the file with the results of the 1DRISM 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  prefactor of the harmonic restraint in the RESPlike procedure, default value: 0.05;
VDWTYPE  functional form of the van der Waals interaction requested: 126 (default), 96;
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 LennardJones 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 3DRISM 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).
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
O 
0.8476 
3.166 
0.1554 
0.000000 
0.000000 
0.064605 
H 
0.4238 
1.000 
0.0560 
0.000000 
0.816495 
0.512747 
0.000000 
0.816495 
0.512747 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
H 
0.435 
1.000 
0.0560 
0.00000 
0.00000 
0.00000 
O 
0.700 
3.070 
0.1700 
0.00000 
0.00000 
0.94500 
CH3 
0.265 
3.775 
0.2070 
1.35136 
0.00000 
1.39716 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
H 
0.435 
1.000 
0.0560 
0.000000 
0.000000 
0.000000 
O 
0.700 
3.070 
0.1700 
0.000000 
0.000000 
0.945000 
CH2 
0.265 
3.775 
0.2070 
1.356103 
0.000000 
1.398746 
CH3 
0.000 
3.905 
0.1750 
1.342751 
0.000000 
2.928687 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
H 
0.435 
1.000 
0.0560 
0.000000 
0.000000 
0.000000 
O 
0.700 
3.070 
0.1700 
0.000000 
0.000000 
0.945000 
CH 
0.265 
3.850 
0.0800 
1.356103 
0.000000 
1.398746 
CH3 
0.000 
3.910 
0.1600 
1.578230 
1.268465 
2.224915 
1.578230 
1.268465 
2.224915 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
CH3 
0.150 
3.775 
0.2070 
1.458000 
0.000000 
0.000000 
C 
0.280 
3.650 
0.1500 
0.000000 
0.000000 
0.000000 
N 
0.430 
3.200 
0.1700 
1.157000 
0.000000 
0.000000 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
O 
0.459 
2.930 
0.2800 
0.000000 
0.000000 
1.530000 
S 
0.139 
3.560 
0.3950 
0.000000 
0.000000 
0.000000 
CH3 
0.160 
3.810 
0.1600 
0.000000 
1.716700 
0.541300 
1.665300 
0.416800 
0.541300 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
O 
0.424 
2.960 
0.2100 
0.000000 
0.000000 
1.208737 
C 
0.300 
3.750 
0.1050 
0.000000 
0.000000 
0.013263 
CH3 
0.062 
3.910 
0.1600 
0.000000 
1.286301 
0.798425 
0.000000 
1.286301 
0.798425 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
C 
0.364 
3.400 
0.1094 
0.000000 
0.000000 
0.814053 
H 
0.218 
2.293 
0.0157 
0.926635 
0.000000 
1.406798 
0.926635 
0.000000 
1.406798 

Cl 
0.036 
3.564 
0.2550 
0.000000 
1.451416 
0.177892 
0.000000 
1.451416 
0.177892 
Atom 
\(z_\alpha\) / a.u. 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
X/Å 
Y/Å 
Z/Å 
C 
0.386 
3.400 
0.1094 
0.000000 
0.000000 
0.463135 
H 
0.266 
2.115 
0.0157 
0.000000 
0.000000 
1.563010 
Cl 
0.040 
3.564 
0.2550 
1.607217 
0.475149 
0.067114 
0.392118 
1.629465 
0.067114 

1.215099 
1.154317 
0.067114 
OPLSAA 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.
Atom 
Example 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
Reference 
C (SP3) 
methane 
3.73 
0.294 

C (SP3) 
neopentane 
3.80 
0.050 

C (SP2) 
isobutene 
3.75 
0.105 

C (SP) 
acetonitrile 
3.40 
0.099 

C (arom) 
benzene 
3.55 
0.07 

C 
chloroform 
4.10 
0.05 

H 
O—H 
0.7 
0.046 

H 
hydrocarbons 
2.42 
0.03 

O 
water 
3.166 
0.1554 

O 
alcohol 
3.07 
0.17 

O 
sulfoxide 
2.93 
0.28 

N (SP) 
acetonitrile 
3.3 
0.099 

S 
sulfoxide 
3.56 
0.395 

Cl 
chloroform 
3.40 
0.300 
Atom 
Example 
\(\sigma_\alpha\) / Å 
\(\epsilon_\alpha\) / kcal/mol 
Reference 
CH (SP3) 
isobutane 
3.850 
0.080 

CH (SP2) 
2butene 
3.800 
0.115 

CH (arom) 
benzene 
3.750 
0.110 

CH2 (SP3) 
nbutane 
3.905 
0.118 

CH2 (SP2) 
1butene 
3.850 
0.140 

CH3 
hydrocarbon 
3.775 
0.207 

CH3 
acetonitrile 
3.6 
0.38 

CH3 
sulfoxide 
3.81 
0.16 

NH2 
amine 
3.3 
0.17 
References
 1(1,2)
M. Reimann and M. Kaupp, Evaluation of an Efficient 3DRISMSCF Implementation as a Tool for Computational Spectroscopy in Solution, Journal of Physical Chemistry A 124, 7439 (2020)
 2
S. Gusarov, T. Ziegler and A. Kovalenko, SelfConsistent Combination of the ThreeDimensional RISM Theory of Molecular Solvation with Analytical Gradients and the Amsterdam Density Functional Package, Journal of Physical Chemistry A 110, 6083 (2006)
 3
D. Casanova, S. Gusarov, A. Kovalenko and T. Ziegler, Evaluation of the SCF Combination of KSDFT and 3DRISMKH; Solvation Effect on Conformational Equilibria, Tautomerization Energies, and Activation Barriers, Journal of Chemical Theory and Computation 3, 458 (2007)
 4
A. Kovalenko and F. Hirata, Potentials of mean force of simple ions in ambient aqueous solution. II. Solvation structure from the threedimensional referenceinteraction site model approach, and comparison with simulations, Journal of Chemical Physics 112, 10403 (2000)
 5
A. Kovalenko and F. Hirata, Selfconsistent description of a metalwater interface by the KohnSham density functional theory and the threedimensional referenceinteraction 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. Threedimensional referenceinteraction site model approach, Journal of Chemical Physics 112, 10391 (2000)
 7
A. Kovalenko, Threedimensional RISM theory for molecular liquids and solidliquid 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 169275.
 8
J.W. Kaminski, S. Gusarov, A. Kovalenko and T.A. Wesolowski, Modeling solvatochromic shifts using the orbitalfree embedding potential at statistically mechanically averaged solvent density, Journal of Physical Chemistry A 114, 6082 (2010)
 9
M. Reimann and M. Kaupp Reaction Entropies in Solution from Analytical ThreeDimensional Reference Interaction Site Model Derivatives with Application to Redox and SpinCrossover Processes, Journal of Physical Chemistry A 126, 3708 (2022)
 10(1,2,3,4,5,6,7,8,9,10)
W.L. Jorgensen, J.D. Madura and C.J. Swenson, Optimized intermolecular potential functions for liquid hydrocarbons, Journal of the American Chemical Society 106, 6638 (1984)
 11(1,2,3,4,5,6,7,8,9,10,11,12,13)
A.E. Kobryn and A. Kovalenko, Molecular theory of hydrodynamic boundary conditions in nanofluidics, Journal of Chemical Physics 129, 134701 (2008)
 12(1,2,3)
O. Acevedo and 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)