In computational chemistry, a solvent model is a computational method that accounts for the behavior of solvated condensed phases.^[1][2][3] Solvent models enable simulations and thermodynamic calculations applicable to reactions and processes which take place in solution. These include biological, chemical and environmental processes.^[1] Such calculations can lead to new predictions about the physical processes occurring by improved understanding.

Solvent models have been extensively tested and reviewed in the scientific literature. The various models can generally be divided into two classes, explicit and implicit models, all of which have their own advantages and disadvantages. Implicit models are generally computationally efficient and can provide a reasonable description of the solvent behavior, but fail to account for the local fluctuations in solvent density around a solute molecule. The density fluctuation behavior is due to solvent ordering around a solute and is particularly prevalent when one is considering water as the solvent. Explicit models are often less computationally economical, but can provide a physical spatially resolved description of the solvent. However, many of these explicit models are computationally demanding and can fail to reproduce some experimental results, often due to certain fitting methods and parametrization. Hybrid methodologies are another option. These methods incorporate aspects of implicit and explicit aiming to minimize computational cost while retaining at least some spatial resolution of the solvent. These methods can require more experience to use them correctly and often contain post-calculation correction terms.^[4]

Implicit models

Implicit solvents or continuum solvents, are models in which one accepts the assumption that implicit solvent molecules can be replaced by a homogeneously polarizable medium as long as this medium, to a good approximation, gives equivalent properties.^[1] No explicit solvent molecules are present and so explicit solvent coordinates are not given. Continuum models consider thermally averaged and usually isotropic solvents,^[3] which is why only a small number of parameters can be used to represent the solvent with reasonable accuracy in many situations. The main parameter is the dielectric constant (ε), this is often supplemented with further parameters, for example solvent surface tension. The dielectric constant is the value responsible for defining the degree of polarizability of the solvent. Generally speaking, for implicit solvents, a calculation proceeds by encapsulating a solute in a tiled cavity (See the figure below). The cavity containing the solute is embedded in homogeneously polarizable continuum describing the solvent. The solute's charge distribution meets the continuous dielectric field at the surface of the cavity and polarizes the surrounding medium, which causes a change in the polarization on the solute. This defines the reaction potential, a response to the change in polarization. This recursive reaction potential is then iterated to self-consistency. Continuum models have widespread use, including use in force field methods and quantum chemical situations. In quantum chemistry, where charge distributions come from ab initio methods (Hartree-Fock (HF), Post-HF and density functional theory (DFT)) the implicit solvent models represent the solvent as a perturbation to the solute Hamiltonian. In general, mathematically, these approaches can be thought of in the following way:^[3][5][6][7]

${\displaystyle {\hat {H))^{\mathrm {total} }(r_{\mathrm {m} })={\hat {H))^{\mathrm {molecule} }(r_{\mathrm {m} })+{\hat {V))^{\text{molecule + solvent))(r_{\mathrm {m} })}$

Note here that the implicit nature of the solvent is shown mathematically in the equation above, as the equation is only dependent on solute molecule coordinates ${\displaystyle (r_{\mathrm {m} })}$. The second right hand term ${\displaystyle {\hat {V))^{\text{molecules + solvent))}$ is composed of interaction operators. These interaction operators calculate the systems responses as a result of going from a gaseous infinitely separated system to one in a continuum solution. If one is therefore modelling a reaction this process is akin to modelling the reaction in the gas phase and providing a perturbation to the Hamiltonian in this reaction.^[4]

${\displaystyle Q(m)=Q_{\mathrm {cavity} }+Q_{\mathrm {electrostatic} }+Q_{\mathrm {dispersion} }+Q_{\mathrm {repulsion} ))$

${\displaystyle G=G_{\mathrm {cavity} }+G_{\mathrm {electrostatic} }+G_{\mathrm {dispersion} }+G_{\mathrm {repulsion} }+G_{\text{thermal motion))}$

Top: Four interaction operators generally considered in the continuum solvation models. Bottom: Five contributing Gibbs energy terms from continuum solvation models.^[5]

The interaction operators have a clear meaning and are physically well defined. 1st - cavity creation; a term accounting for the energy spent to build a cavity in the solvent of suitable size and shape as to house the solute. Physically, this is energy cost of compressing the solvents structure when creating a void in the solvent. 2nd term - electrostatic energy; This term deals with the polarization of the solute and solvent. 3rd term - quantum mechanical dispersion energy; can be approximated using an averaging procedure for the solvent charge distribution.^[5] 4th term - an approximation for the quantum mechanical exchange repulsion; given the implicit solvent this term can only be approximated against high level theoretical calculations.

These models can make useful contributions when the solvent being modelled can be modelled by a single function i.e. it is not varying significantly from the bulk. They can also be a useful way to include approximate solvent effects where the solvent is not an active constituent in the reaction or process. Additionally, if computer resources are limited, considerable computational resources can be saved by evoking the implicit solvent approximation instead of explicit solvent molecules. Implicit solvent models have been applied to model the solvent in computational investigations of reactions and to predict hydration Gibbs energy (Δ_hydG).^[8] Several standard models exist and have all been used successfully in a number of situations. The Polarizable continuum model (PCM) is a commonly used implicit model and has seeded the birth of several variants.^[5] The model is based on the Poisson-Boltzmann equation, which is an expansion of the original Poisson's equation. Solvation Models (SMx) and the Solvation Model based on Density (SMD) have also seen wide spread use. SMx models (where x is an alphanumeric label to show the version) are based on the generalized Born equation. This is an approximation of Poisson's equation suitable for arbitrary cavity shapes. The SMD model solves the Poisson-Boltzmann equation analogously to PCM, but does so using a set of specifically parametrised radii which construct the cavity.^[9] The COSMO solvation model is another popular implicit solvation model.^[10] This model uses the scaled conductor boundary condition, which is a fast and robust approximation to the exact dielectric equations and reduces the outlying charge errors as compared to PCM.^[11] The approximations lead to a root mean square deviation in the order of 0.07 kcal/mol to the exact solutions.^[12]

Explicit models

Explicit solvent models treat explicitly (i.e. the coordinates and usually at least some of the molecular degrees of freedom are included) the solvent molecules. This is a more intuitively realistic picture in which there are direct, specific solvent interactions with a solute, in contrast to continuum models. These models generally occur in the application of molecular mechanics (MM) and dynamics (MD) or Monte Carlo (MC) simulations, although some quantum chemical calculations do use solvent clusters. Molecular dynamics simulations allow scientists to study the time evolution of a chemical system in discrete time intervals. These simulations often utilize molecular mechanics force fields which are generally empirical, parametrized functions which can efficiently calculate the properties and motions of large systems.^[6][7] Parametrization is often to a higher level theory or experimental data. MC simulations allow one to explore the potential energy surface of a system by perturbing the system and calculating the energy after the perturbation. Prior criteria are defined to aid the algorithm in deciding whether to accept the newly perturbed system or not.

In general, force field methods are based on similar energy evaluation functionals which usually contain terms representing the bond stretching, angle bending, torsions and terms for repulsion and dispersion, such as the Buckingham potential or Lennard-Jones potential. Commonly used solvents, such as water, often have idealized models generated. These idealized models allow one to reduce the degrees of freedom which are to be evaluated in the energy calculation without a significant loss in the overall accuracy; although this can lead certain models becoming useful only in specific circumstances. Models such as TIPXP (where X is an integer suggesting the number of sites used for energy evaluation)^[13] and the simple point charge model (SPC) of water have been used extensively. A typical model of this kind uses a fixed number of sites (often three for water), on each site is placed a parametrized point charge and repulsion and dispersion parameter. These models are commonly geometrically constrained with aspects of the geometry fixed such as the bond length or angles.^[14]

Advancements around 2010 onwards in explicit solvent modelling see the use of a new generation of polarizable force fields, which are currently being created. These force fields are able to account for changes in the molecular charge distribution. A number of these force fields are being developed to utilise multipole moments, as opposed to point charges, given that multipole moments can reflect the charge anisotropy of the molecules. One such method is the Atomic Multipole Optimised Energetics for Biomolecular Applications (AMOEBA) force field.^[15] This method has been used to study the solvation dynamics of ions.^[1] Other emerging polarizable forcefields which have been applied to condensed phase systems are; the Sum of Interactions between Fragments ab initio computed (SIBFA)^[16] and the Quantum Chemical Topology Force Field (QCTFF).^[17] Polarizable water models are also being produced. The so-called charge on spring (COS) model gives water models with the ability to polarize due to one of the interaction sites being flexible (on spring).^[18]

Hybrid models

Hybrid models, as then name suggests, are in the middle between explicit and implicit models. Hybrid models can usually be considered closer to one or other implicit or explicit. Mixed quantum mechanics and molecular mechanics models,(QM/MM) schemes, can be thought of in this context. QM/MM methods here are closer to explicit models. One can imagine having a QM core treatment containing the solute and may be a small number of explicit solvent molecules. The second layer could then comprise MM water molecules, with a final third layer of implicit solvent representing the bulk. The Reference Interaction Site Model (RISM) can be thought of being closer to implicit solvent representations. RISM allows the solvent density to fluctuate in a local environment, achieving a description of the solvent shell behaviour.^[1][2][5]

QM/MM methods enable a section of the system to be calculated using quantum mechanics, for example the active site in a biological molecule, whilst the rest of the system is modeled using MM force fields. By continuing to a third layer with an implicit solvent the bulk water effect can be modeled more cheaply than using all explicit solvent molecules. There are many different combinations that can be used with the QM/MM technique. Alternatively, a few explicit solvent molecules can be added to a QM region and the rest of the solvent treated implicitly. Previous work has shown mixed results upon the addition of explicit solvent molecules to an implicit solvent. One example added up to three explicit water molecules to a QM calculation with an implicit COSMO water model. The results suggest that using either implicit or explicit solvent alone provide a good approximation to experiment, however, the mixed models had mixed results and possibly some dependence on the number of added explicit solvent molecules.^[19]

RISM, a classical statistical mechanics methodology, has it roots in the integral equation theory of liquids (IET). By statistically modeling of the solvent, an appreciation of the dynamics of the system can be acquired. This is more useful than a static model as the dynamics of the solvent can be important in some processes. The statistical modeling is done using radial distribution function (RDF). RDF are probabilistic functions which can represent the probability of locating solvent atoms/molecules in a specific area or at a specific distance from the reference point; generally taken as the solute molecule. As the probability of locating solvent atoms and molecules from the reference point can be determined in RISM theory, solvent shell structure can be directly derived.^[20]

The molecular Ornstein-Zernike equation (MOZ) is the starting point for RISM calculations.^[5] Within the MOZ equations a solvated system can be defined in 3D space by three spatial coordinates (r) and three angles (Θ). Using relative RDF's the MOZ equations for the solvated system can define the total correlation function h(r - r';ʘ - ʘ'). The equations have a high dimensionality (6D).

${\displaystyle h(r-r';\Theta -\Theta ')=g(r-r';\Theta -\Theta ')-1}$

${\displaystyle h(r;\Theta )}$ is the total correlation function, ${\displaystyle g(r;\Theta )}$ is the radial distribution function accounting for the direct effects of one molecule on another separated by r.^[5]

It is a common approximation to assume spherical symmetry, allowing one to remove the orientational (angular) degrees of freedom. The MOZ equation splits the total correlation function in two. First the direct correlation function c(r), concerned with the effect of one particle on one other over the distance r. The second, the indirect correlation function, accounts for the effects of a third particle in a system. The indirect correlation function is given as the direct correlation function between the first and the third particles ${\displaystyle c(r_{1,3})}$ in addition to the total correlation function between the second and third particles ${\displaystyle h(r_{2,3})}$.^[21]

${\displaystyle h(r)=c(r_{1,2})+\int \mathrm {d} r_{3}\,c(r_{1,3})\rho (r_{3})h(r_{2,3})}$

Ornstein-Zernike equation with the assumption of spherical symmetry. ρ is the liquid density, r is the separating distance, h(r) is the total correlation function, c(r) is the direct correlation function.

h(r) and c(r) are the solutions to the MOZ equations. In order to solve for h(r) and c(r), another equation must be introduced. This new equation is called a closure relation. The exact closure relation is unknown, due to the so-called bridge functions exact form being unclear, we, therefore, must introduce approximations. There are several valid approximations, the first was the HyperNetted Chain (HNC), which sets the unknown terms in the closure relation to zero. Although appearing crude the HNC has been generally quite successfully applied, although it shows slow convergence and divergent behaviour in some cases.^[22] A modern alternative closure relation has been suggested the Partially Linearised HyperNetted Chain (PLHNC) or Kovalenko Hirata closure.^[23] The PLHNC partially linearises the exponential function if it exceeds its cutoff value. This causes a much more reliable convergence of the equations.^[4]

${\displaystyle h_{\alpha }(r)={\begin{cases}\mathrm {e} ^{-\beta U(r)+T(r)}-1&({\text{when))-\beta \upsilon _{a}(r)+h_{a}(r)-c_{a}(r)\leq 0)\\-\beta U(r)+T(r)&({\text{when))-\beta \upsilon _{a}(r)+h_{a}(r)-c_{a}(r)>0)\end{cases))}$

The PLHNC closure, where ${\displaystyle \beta ={\frac {1}{k_{B}T))}$ and ${\displaystyle U(r)}$ is the interaction potential, a typical interaction potential is shown below. T(r) is the indirect correlation function, as it is the difference of the total and the direct correlation functions.

${\displaystyle U(r)=4\epsilon \left[\left({\frac {\sigma _{1)){r_{12))}\right)^{12}-\left({\frac {\sigma _{2)){r_{12))}\right)^{6}\right]+{\frac {Q_{1}Q_{2)){r_{12))))$

There are various approximations of the RISM equations. Two popular approximations are 3D RISM and 1D RISM.^[1] There are known deficiencies in these approximate RISM models. 3D RISM makes a poor estimation of the cavity creation term. 1D RISM has been found to not be properly accounting for the spatial correlations of the solvent density around the solute. However, both methods are quick to calculate, 1D RISM can be calculated in a matter of seconds on a modern computer, making it an attractive model for high through put computation. Both 3D RISM and 1D RISM have had correction schemes proposed which make the predictions reach a comparable level of accuracy to traditional implicit and explicit models.^[22][24][25]

The COSMO-RS model is another hybrid model using the surface polarization charge density derived from continuum COSMO calculations to estimate the interaction energies with neighbored molecules. COSMO-RS is able to account for a major part of reorientation and strong directional interactions like hydrogen bonding within the first solvation shell. It yields thermodynamically consistent mixture thermodynamics and is often used in addition to UNIFAC in chemical engineering applications.

Applications to QSAR and QSPR

Quantitative Structure–Activity Relationships (QSAR)/Quantitative Structure–Property Relationships (QSPR), whilst unable to directly model the physical process occurring in a condensed solvent phase, can provide useful predictions of solvent and solvation properties and activities; such as the solubility of a solute.^{[26][27][28][4]} These methods come in a varied way from simple regression models to sophisticated machine learning methods. Generally, QSAR/QSPR methods require descriptors; these come in many different forms and are used to represent physical features and properties of a system of interest. Descriptors are generally single numerical values which hold some information about a physical property.^[29] A regression model or statistical learning model is then applied to find a correlation between the descriptor(s) and the property of interest. Once trained on some known data these model can be applied to similar unknown data to make predictions. Typically the known data comes from experimental measurement, although there is no reason why similar methods can not be used to correlate descriptor(s) with theoretical or predicted values. It is currently debated whether if more accurate experimental data was used to train these models whether the prediction from such models would be more accurate.^[30]

More recently the rise of deep learning has provided many methods to generate embedded representations of molecules.^[31][27] Some of these methods have also been applied to solvation properties such as solubility prediction ^[32]

^[27]