First-principles modeling of chemistry in mixed solvents: Where to go from here?

Solvents give rise to drastic effects on chemical kinetics, thermodynamics, and reaction mechanisms.¹ Mixed solvents are appealing because their polarity and Brønsted acidity and basicity can be continuously modulated to optimize reaction outcomes.^2–7 For a recent example, tuning the concentration of hydroxyl groups in aqueous organic solvent mixtures tremendously increases rates of various acid-catalyzed reactions up to 1000 fold.⁸ The synergistic effects of complicated solvent environments are clearly worth studying, but it remains unclear how to best translate practical calculation schemes used for single-component (pure) solvents (e.g., for predicting pK_as,⁹ standard redox potentials,¹⁰ and hydricities¹¹) to make insightful predictions in complicated mixed solvent environments.

In this Perspective, solvation thermodynamics, which are more thoroughly discussed in other references,^12–14 will be approximated as arising mostly from electrostatic, dispersion, and configurational entropy contributions. Figure 1 graphically represents key differences between pure and mixed solvents acting on a system having a net dipole. For charged or polar solutes in polar or ionic solvents, the largest contributions to solvation free energies are electrostatic, and these terms have been well treated historically by existing theories. A significant amount of electrostatic solvation free energy is gained through the reorientation of solvent molecules to cancel the electric field of the solute. Under solvation in a pure solvent, the solvent density remains homogeneously distributed apart from heterogeneities due to formation of direct contacts with the solute [Fig. 1(a)]. However, in a mixed solvent, the relative densities of the solvent components can change, and this allows for charge migration as well as segregation of molecules having different dipole moments to result in complex double-layer structures [Fig. 1(b)].

Beyond electrostatic terms, electronic fluctuations of polarizable molecules lead to dispersion forces between the solvent and the solute. Traditionally, most solvation models have treated dispersion forces as short-range interactions, but recent work has shown that electronic fluctuations can also drive long-range and highly collective interactions that are not adequately modeled as short-range two-body forces.^15,16

Finally, flexible solute molecules may possess significant configurational entropy arising from their nuclear motion, and this can either increase or decrease (in much the same way as dispersion energies) depending on how collective degrees of freedom of the solute are coupled to solvent fluctuations. These three categories of solvation free energies are furthermore all coupled, crossing over in particular when mesoscale electrostatic effects alter the solute electronic structure.

Once the physical components of solvation free energies are accounted for, computational procedures can be used to predict reaction energies, acidity constants, standard redox potentials, and other quantities.^17,18 Using a generalized thermodynamic cycle, as shown in Scheme 1, one can compute the free energy change for a chemical process in any solvating liquid phase as the lower horizontal leg, $ΔGliq*$⁠, as a function of an (electro)chemical environment (e.g., pH and an applied bias, ϕ).¹⁸ The superscript ^* denotes the liquid standard state of 1M. The upper horizontal leg of the cycle is the gas phase free energy difference for the same reaction (using a standard state of 1 bar, denoted by superscript °). An assortment of quantum chemical methods from Kohn–Sham density functional theory (DFT) to correlated wave function theories with corresponding ideal gas, rigid rotor, and harmonic oscillator approximations may be convenient in this setting.

The solvation free energy contributions shown by the vertical legs of this thermodynamic cycle (Scheme 1) are normally assumed to be challenging to determine precisely, as well as accurately. For instance, solvation free energies can be determined using implicit, explicit, or hybrid solvation models, but those calculated energies can be sensitive to solute geometries and simulation parameters.

Scheme 1 also requires knowledge about the proton solvation free energy in the solvent, $ΔGS*(H+)$⁠, as well as a reference (i.e., absolute) free energy for an electron in a reference electrode interacting with that solvent, G(e⁻). The latter can be derived using the former, and both are known in water and some other pure solvents through semiempirical schemes.^19,20 Note that $ΔGS*(H+)$ and G(e⁻) are generally not known in mixed solvents. In addition, there are discrepancies in the literature due to whether the reported proton solvation free energy and its corresponding absolute electrode potential are correctly accounting for the phase potential for the solvent, as discussed later.^17,21 In the traditional derivation of Trasatti,²² the absolute proton solvation free energy used (−1088 kJ/mol) led to an absolute reduction potential of −4.44 V for the standard hydrogen electrode. Another widely used value for the proton solvation free energy in water from Tissandier et al.,²³ −1105 kJ/mol, is a value that we believe represents a real solvation energy that includes a negative valued phase potential contribution. Kelly et al.²⁰ calculated a similar value to derive an absolute reduction potential of −4.28 V, a value also obtained by Isse and Gennaro¹⁹ with $ΔGS*(H+)=−1101$ kJ/mol and Fermi–Dirac instead of Boltzmann statistics.

The generalized thermodynamic cycle in Scheme 1 allows for the practical determination of changes in Gibbs free energies for liquid phase reactions. Since gas-phase free energy calculations can be treated at reasonably high accuracy with appropriate quantum chemistry methods, the remaining challenges in using this cycle are the determination of solvation free energies, as well as the important terms for $ΔGS*(H+)$ and G(e⁻). We now discuss issues that obfuscate the determination of these values.

In a practical sense, the phase potential of a chemically balanced reaction in a single phase will cancel so long as real or absolute solvation free energies are used consistently.³² Being unsure whether a phase potential is included in a solvation free energy can lead to errors in properties (e.g., pK_as and redox potentials) calculated from the thermodynamic cycle in Scheme 1. Furthermore, quantifying the phase potential is difficult and lacks clear consensus.^33–35 There have been numerous classical and quantum mechanical estimates for water ranging from −1.10 V to +3.63 V. Investigation of other pure and mixed solvent phase potentials would require carefully crafted classical force fields or significant computational cost for molecular dynamics with quantum chemistry. There is currently little empirical or computational benchmark information for the phase potentials of mixed solvents with air, indicating a need for the further work in this field. In the absence of accurate phase potential calculations or measurements, the best (but unsatisfactory) approximation is to refer to any known solvents hypothesized to have a comparable surface structure to the system at hand.

Preferential solvation is the local concentration variation in a mixed solvent that occurs when one chemical species has a different affinity for the solute than others. In the example of Fig. 1(b), this situation leads to formation of alternating layers. For complex species or for strong concentration gradients, this phenomenon can become highly intractable, even leading to local demixing and the formation of droplets or other structures. In the special case of monatomic ions around charged or polar solutes, Poisson–Boltzmann–Debye–Hückel theory is available to rewrite the Coulomb electrostatic potential with an additional ion screening term, solved in spherical symmetry to give ϕ ∝ e^κ/r/ϵr (where ϵ is the dielectric constant and 1/κ is the Debye length over which screening by mobile charges becomes significant), rather than ϕ ∝ 1/ϵr as in the absence of salt.³⁶ Poisson–Boltzmann–Debye–Hückel theory is powerful, but it neglects shape effects and nonelectrostatic forces, thus being mainly useful for monatomic cations, and mainly far from the solute surface. Additions to the theory to treat the modification of ion density in a short range by finite radius of ions and the presence of nonelectrostatic interactions exist and can be highly effective.³⁷ However, the theory remains more useful for long-range electrostatics but less so for steric and dispersion effects, which are usually in the short range.

Other theoretical frameworks provide a preferential solvation parameter on a general basis for a solute X in a mixture of i and j, δ_Xi. This parameter is defined as the difference between the local and bulk mole fractions (x) of a solvent component i around the solute X, $δXi=xXiL−xi$⁠. For example, Kirkwood–Buff (KB) theory relates preferential solvation to integrals of the radial distribution function in the form of a KB integral (KBI). Such integrals usefully distill the affinity between X and i to a scalar greater than one for two species that associate, or less than one for two species that repel each other.^38,39 Inferring δ_Xi can be achieved by using the concentration dependence of experimental data such as isothermal compressibility, partial molar volumes, or excess molar Gibbs free energy of mixing. Additionally, the standard molar Gibbs free energy of transfer (ΔG_tr) of the solute from a reference solvent i to a mixture of i and j is required to obtain the KBI via the inversion of KB theory.⁴⁰ These properties and their derivatives must be determined with high accuracy, such that experiments may be difficult or expensive,⁴¹ so the alternative of molecular dynamics is sometimes used to determine the KBI. Convergence of integrals over correlation functions using molecular dynamics is, however, difficult due the large weight given to long-range parts of the correlation function (which are also the slowest to sample),^42,43 ensemble effects, and systematic force field errors that may be present in medium and long ranges. The preferential solvation parameter is directly dependent on the derivatives of these properties, thus motivating the need for highly accurate parameterizations.

The quasi-lattice quasi-chemical (QLQC) method is another way to predict δ_Xi.^44–46 The QLQC method also requires information typically acquired through experiments such as ΔG_tr and excess molar Gibbs free energy of mixing of components i and j at equimolar composition. Additionally, the QLQC method requires the selection of a lattice parameter, Z, which defines the number of neighbors around the solute and every solvent molecule. Z is conceptually related to the coordination number of the solute (and solvent molecules), but it mostly serves as a fitting parameter for the excess Gibbs energy of mixing curve. This adds emphasis to quality experimental data and can make a priori prediction of δ_Xi difficult. The QLQC method typically uses the same Z for all compositions of the same mixed solvent; this would be valid if the coordination number is independent of composition, which is not always the case. There are ways to have concentration dependence in Z, but this is generally avoided because it greatly complicates calculations and reduces intuitive understanding. In addition, the assumption of ideal mixing of the neighboring molecules in the quasi-lattice sites is inherently incompatible with preferential solvation. With all of these characteristics, the QLQC method is typically restricted to small solutes with sufficient experimental data.⁴⁷ 

In summary, the above models and others^48–50 are useful and powerful for predicting preferential solvation effects for various solutes in mixed solvents. However, each of them has limitations due to computational cost, narrow range of well-treated cases, or requirement for specific empirical parameters. One might be tempted to average properties of the pure solvent components linearly to obtain a reasonable approximation to a mixed solvent system, as discussed in some textbooks on the subject.^46,51 This approach is enticing because it avoids the effort of predicting preferential solvation and seems plausible, but in general, mixed solvent properties exhibit nonlinearities and solute-solvent properties are dominated by the local solvent structure. Figure 2 shows the Gibbs free energy required to transfer Na⁺ or CH₃COO^– from pure water to a mixture of water/acetonitrile at varying mole fractions.^52,53 Both ions deviate from the linear interpolation line between pure water and pure acetonitrile. Na⁺ even has favorable transfer Gibbs free energy at acetonitrile mole fractions ranging from 10% to 50% compared to transfer from pure water to pure acetonitrile. Thus, simple averaging schemes not only are sometimes unphysical approximations but may also produce significant errors that arise from preferential solvation effects.

Explicit solvation methods involve placing the solute(s) in a large box surrounded by solvent molecules, where forces on each atom are calculated to determine the trajectories for the entire system. This straightforward technique over time has become a ubiquitous treatment for mixed solvent systems. Evaluation of the forces involves classical force fields, quantum chemical methods, and combinations of both. However, each method is dependent either on accurate parameterizations or on computationally expensive molecular dynamics with quantum chemistry. Useful information about mixed solvent systems is available from explicit modeling, but other less costly methods should be considered.

An implicit solvation model is an alternative to explicit methods. In this case, the solvent is treated as a dielectric continuum, possibly with electrostatic screening by mobile ions. Simplifying the solvent to a homogeneous medium requires extensive experimental data for significant parameterization. Furthermore, implicit solvation models calculate only a single parameter per lattice point (the electrostatic potential or its gradient), meaning that complex ordering of co-solvents, often important near the solute, cannot be captured. Forfeiting the capability to account for preferential solvation makes mixed solvents challenging for implicit solvation models, especially for polar or ionic chemical species. However, some implicit solvation models [e.g., the Solvation Model based on Density (SMD) and the Integral Equation Formalism (IEF)] have seen success with ionic liquids, but so far applications have been limited to small, neutral molecules.^54,121

Improved theoretical foundations are needed to physically and accurately allow for predictions of solutes in mixed solvent environments without the high computational cost of explicit quantum solvation. Two such models, COnductor-like Screening MOdel for Real Solvents (COSMO-RS) and reference interaction site model (RISM), are briefly summarized as they have been demonstrated to model mixed solvents with reasonable accuracy and computational cost.

The COnductor-like Screening MOdel for Real Solvents (COSMO-RS) is classified as a hybrid solvation model.^55–59 Every COSMO-RS calculation starts from the COSMO implicit solvation model cavity with an infinite dielectric (perfect conductor) reference state (ϵ = ∞). The local interaction energy, $EintA,B$⁠, of two molecules contacting each other in a liquid is quantified in COSMO-RS from the COSMO polarization charge densities, σ and σ′, from their respective COSMO surfaces. COSMO-RS sums nonideal screening from unbalanced σ and σ′ contributions (so-called electrostatic misfits) and hydrogen bonding by using surface segments with a_cont contact surface area to account for $EintA,B$⁠.

By summing all energy corrections over all molecular surface contacts, COSMO-RS obtains a total energy correction for that specific molecular configuration. This energy correction then needs to be expanded into an ensemble to represent a liquid. Instead of employing molecular dynamics or Monte Carlo procedures, a simpler statistical thermodynamics model of an ensemble of independently pairwise interacting surface segments is used.

Overall, COSMO-RS is a valuable tool for chemical engineering thermodynamics.⁵⁸ The COSMO-RS scheme significantly simplifies the problem of lengthy explicit simulations and produces solvation free energies in pure and mixed solvents. A self-consistent generalization of COSMO-RS, Direct COSMO-RS (DCOSMO-RS), has been implemented in several quantum chemistry packages.⁶⁰ 

While COSMO-RS has several significant improvements over COSMO and other implicit solvation models, it still experiences some shortcomings. COSMO-RS also requires parameterizations for nonelectrostatic contributions from cavitation, dispersion, and repulsion,⁶¹ but robust parameterizations in turn require extensive databases. Thus, approximate parameters are generally used in practice even though they might not be suitable for highly accurate computational investigations. COSMO-RS was originally developed and parameterized for the neutral solute and solvents and incorrectly assigned misfit charges while modeling charged solutes.⁶² However, recent developments have shown promising success with ionic liquids⁶³ and electrolytes.⁶⁴ In addition, an extension called COSMO-RSC was developed to relieve some of the error associated with highly localized charges; however, these corrections are not consistently implemented in quantum chemistry packages and thus should be used with caution.^21,65

The reference interaction site model (RISM) is another hybrid method, based on the integral equation theory of liquids. It was first developed in 1972 by Chandler and Anderson,⁶⁶ but it has gone through several iterations to make it an intriguing solvation model for the study of future problems. Extensive literature on the RISM already exists,^67–71 so only a brief summary is given here.

In the 1D-RISM, a single radial distribution for each solvent site around the solute is determined. In the so-called 3D-RISM, a three-dimensional density distribution $gγr$ for each interaction site (γ) of each solvent species present is determined. A value of $gγr$ greater or less than one shows solvent site enrichment or depletion, respectively. Thermodynamic quantities, such as entropy, enthalpy, and pressure, are then determined as integrals (single, double, or higher) over the 3D density fields of the various interacting sites.

An advantage of working with 3D density maps is that distributions of specific atoms or functional groups of the solvent(s) in relation to the solute can be plotted directly. Unlike implicit solvation models, the 3D-RISM directly accounts for nontrivial ordering in the solvent, such as electric double-layers or even more complex density variations, as illustrated with a peptide fibril structure⁷² in Fig. 3.

The 3D-RISM connects the 3D total correlation h_αγ(r), an emergent many-body property, and the direct two-body correlation functions c_αγ(r). The two-body correlation is directly determined by the susceptibility and can be calculated to a good approximation for small fluctuations directly from the pairwise interaction potential −βu_γα(r). Such two-body correlations are also observable with scattering experiments made at low concentrations such that many-body effects are limited. Conversely, the theory can also find intermolecular potentials based on h_αγ(r) from scattering experiments at high concentrations.

By this point, it should be clear that the 3D-RISM can be highly effective for analyzing the role of mixed solvents.^76–78 The computational expense of the 3D-RISM depends on which closure relation is chosen, as well as the parameters (e.g., Lennard-Jones) of the given system. As the solvent increases in complexity (i.e., the number of co-solvents or number of sites treated per solvent), the number of pair correlations and therefore cost of a self-consistency iteration rise as the square of the number of sites. Furthermore, the number of density field grids needed increases per site and thus requires significant computer memory resources. The 3D-RISM is therefore quite elegant but also very challenging for complex and flexible solvents. These challenges can be reduced if complicated co-solvents are treated as being present only in the low concentration limit so that solvent-solvent site interactions can be ignored within the problematic species.⁷⁷ The required fineness of the grid increases with the sharpness of concentration gradients; therefore, small or strongly interacting solvent/co-solvent sites can dramatically increase the computational cost. The generalization of the method to adaptive grids would appear to be a potential avenue for progress. Strong interactions and high concentrations also increase the importance of three-body and higher terms in the mapping from pair to total correlations, thus increasing the likelihood of slow or failed convergence depending on the closure.

Because the 3D-RISM can use pairwise interaction potentials as an input, the user is at liberty to supply the input data at arbitrary levels of quantum mechanical accuracy, although of course many-body electronic interactions cannot enter via the two-body parameterization process. With all of these technical aspects aside, the 3D-RISM is convenient to use with conventional DFT calculations, and it also has overlapping characteristics as an algorithm, so it has been integrated into some quantum chemistry packages such as Amsterdam Density Functional (ADF).⁷⁹ 

Another technique commonly used to model solvents is cluster-continuum or mixed implicit/explicit solvation modeling. Explicit solvent molecules are added around the solute and then this cluster is placed in an implicit solvation model, such that the blue region described as a “large volume with weak order” in Fig. 1(a) is treated as a continuum dielectric, while the inner region is modeled classically or via quantum chemistry. This procedure is advantageous for several reasons, including balanced computational costs and retention of local solvent environment characteristics.

Several challenges still remain, and there is a lack of widely accepted best practices for cluster-continuum modeling. For example, there is no clear consensus on choosing the number and location of the explicit solvent molecules. Often, chemical intuition has guided the placement of solvent molecules, but this is limited by a priori understanding of the chemical process.⁸⁰ In this section, we present our interpretation for modeling chemistry in a complicated, mixed solvent environment.

Partitioning the contributions of the inner and outer solvation shells around a solute is fundamental in these cluster-continuum solvation schemes. Quasi-chemical theory (QCT) represents the most robust method to partition the system and calculate the solvation free energy of any solute. Only a brief introduction will be given here as there is already substantial literature on the subject.^27,81–84

Here, the excess chemical potential is broken into three terms. The first being concerned with the ligand-association chemical reaction in an ideal gas phase is denoted by the (0) superscript. An advantage of the QCT formulation is that calculation of $Kn(0)$ may be accomplished with gas-phase electronic structure calculations for the free energies of the solute, individual ligand, and solute-ligand cluster. The overall availability of ligands for the reaction is represented by the bulk density of the ligand, ρ_L. This association reaction is designated to occur within a volume chosen to represent the inner shell.

The second term accounts for the thermal fluctuations of having n ligands occupying the inner shell, defined by considering the probability of observing that particular case, $pXn$⁠. Evaluating this probability involves molecular dynamics, usually with quantum chemistry, with adequate sampling of the solute surrounded by the solvent. In some circumstances, classical force fields can be used, but then results become dependent on the quality of parameterization.⁸⁶ 

The outer shell (bulk) contribution to the solvation free energy of X is captured in the third term. Essentially, it is the free energy contribution of solvating the solute-ligand complex and removing the ligands from the outer solvation environment. A standard implicit solvation model, while an approximation, has effectively evaluated $μXLn(ex)$ and $nμL(ex)$ in previous applications.^27,87–98

In QCT, selection of the inner solvation shell is typically done manually. While any shape can be used, a sphere with the radius λ centered around X is commonly used. The selection of λ is generally chosen to capture the innermost-solvation shell defined by a neighborship analysis of solvent molecules that contribute to the first radial distribution function peak;^27,99,100 however, any choice of λ is valid in QCT. Ultimately, the choice of λ dictates the probabilities of observing n ligands in the inner shell, and any case of n ligands for the inner solvation shell can be evaluated in QCT.

Some care toward selecting λ and n can deliver tremendous computational feasibility and minimize errors for several reasons. Implicit solvation models depend on the definition of the cavity to calculate solvation free energies, which are sensitive to the radii being used. When a complete solvation shell is defined, the free energy of solvation by the implicit solvation model is less sensitive to specific radius definitions. This trend suggests choosing λ that captures complete inner-shell occupancies for convenient error cancellation. While any number of complete solvation shells could be treated, the simplest inner shell volume contains only the nearest neighboring ligands that are continuously interacting with the solute. This “no split occupancies rule” means λ should be selected by a neighborship decomposition analysis,⁹⁹ where distant neighbors that encompass many shells (multimodal) are excluded.^27,92,96,100

As previously mentioned, $pXn$ must be determined by the observation frequency. Choosing n that is rarely observed would require long simulation times.⁹³ When molecular dynamics with quantum chemistry is used, the higher computational cost motivates choosing the most probable n within λ for the solute, $n¯$⁠. Applying QCT to the case of $n¯$ would have $pXn$ near one, so the term can be dropped altogether. The case of $n¯$ with QCT can reduce the computational cost, but it is not always known. While $n¯$ can be determined from explicit simulations and neighborship analyses, it can become computationally impractical for large solutes without vetted, highly accurate force fields. This makes direct implementation into computational chemistry challenging.

QCT has produced predictions of solvation free energies in excellent agreement with experimental values for hydration of cations and anions,^{27,87–93,95–98} as well as other molecules,^101,102 and of free energies of ion binding to proteins and solvation by nonaqueous solvents.^{27,83,86,94,95,100,103–106} In addition, QCT predictions of pK_as and redox potentials agree well with experimental values.^107,108 While QCT is completely generalized for any solute and solvent, its application to mixed solvents has not yet been published. Nevertheless, these characteristics of QCT demonstrate its robustness and theoretical rigor for calculating absolute solvation free energies.

The advantage of QCT is the separation of solvation free energies into components from inner-shell and outer-shell solvent molecules, permitting exploitation of the known chemical physics analyses of clusters. Issues that can be treated in cluster analyses include proper overlap repulsions, polarizability, charge transfer, London dispersion interactions, n-body ligand interactions generally, distortion of flexible ligands, and zero-point motion of the cluster. Large displacements in the vibrational motions of the clusters can be treated with molecular dynamics, as done recently for anions using quantum chemical simulations, and earlier for hydrogen using classical force field simulations.^98,102 Current implementations of QCT require significant effort from manual evaluation of terms, with molecular dynamics simulations and neighborship analyses used to guide selection of inner and outer solvation shells and determine $n¯$ for computational feasibility and to minimize errors.

When modeling chemical processes (e.g., reaction mechanisms with several intermediates), QCT may become impractical due to the high computational cost of molecular dynamics or requirement of accurate force field parameters for all structures (i.e., reactants, products, intermediates, and transition states) to determine $pXn$ or $n¯$⁠. However, theoretical and conceptual foundations of QCT are valuable for developing future solvation schemes. Here, we discuss a novel, automatable procedure for studying solvent effects and computing solvation free energies.¹⁰⁹ The key feature is the integration of the QCT framework twice (ligand free energy and ion solvation free energy) into a standard thermodynamic cycle in combination with a machine learning analysis of local solvent environments. This dual-QCT approach achieves additional error cancellation by balancing QCT terms in such a way that evaluating $pXn$ or seeking $n¯$ and using harmonic approximations are unnecessary. Additional advantages are automation and eliminating potentially cost-prohibitive molecular dynamics with quantum chemistry and neighborship analyses.

The central idea is to generate many candidate clusters of different structures and a number of solvent molecules (sizes). Then, clusters that satisfy specific criteria are used in the cluster cycle to calculate the solute solvation free energy. In past applications of the cluster cycle, solvent molecules are manually placed or carved out of an explicit simulation. This route proves to be tedious or computationally inefficient. Alternatively, an algorithm can be used to automatically generate the prospective clusters. One such example is ABCluster, which employs the artificial bee colony algorithm to find global and local minima of atomic or molecular clusters.^111,112 For each size, thousands of clusters can be rapidly generated and initially ranked according to their energy with classical force fields. Then, the lowest (e.g., five) can be further optimized at 0 K with quantum chemical methods such as DFT. To achieve an adequate representation of the local solvent environment with this multistep filtering procedure, the Boltzmann-weighted average of same sized clusters is taken. Note that the cluster cycle requires solvation free energies for solvent and solute-solvent clusters, so the above procedure is performed both in the absence and in the presence of the solute.

This still leaves the main challenge of determining the appropriate cluster size unsolved. One could compare to experimental data to solve this problem; however, this would not be different from current practices of ad hoc selection of the number of solvent molecules. Furthermore, having experimental solvation energies for any and all chemical species and environments is an unfeasible endeavor; there need to be separate selection criteria independent of experimental data.

As mentioned before, the local solvent environment is the largest contributor to solvation energetics, so this is a natural choice for selection criteria. To avoid requiring a priori information about the local solvent environment, a separate technique can be used to compare clusters of different sizes. Comparing the structures of molecules, or clusters, is not novel and is commonly done through connectivity and distances as in root mean square displacement (RMSD). Calculating RMSDs, however, requires manual alignment of the systems being compared that becomes impossible when comparing clusters containing different numbers of atoms. Over the last few years, several fingerprint functions have been used to characterize the environment of atoms, but one function that fits all does not apply here. Each is designed to differentiate between specific system characteristics. For example, some fingerprint functions characterize a system by elemental characteristics or bond angles which would not be able to differentiate local solvent environments. Hence, selecting an appropriate fingerprint function is paramount for the critical evaluation of these microsolvated clusters.

One promising algorithm is the smooth overlap of atomic positions (SOAP) kernel that enables comparison between clusters that is invariant to translations, rotations, and permutations.^113,114 The SOAP kernel evaluates the local environment by placing a Gaussian function centered on each atom. Then, the SOAP gathers information from where there is overlap of Gaussian functions (local atomic neighbor densities) and reports high-dimensional pair-similarity data. Essentially, this quantifies the similarities and differences of the local solvent environment of different clusters.

After employing the SOAP kernel, one has detailed information about the local solvent environments but needs an intuitive way to compare them. A data dimensionality reduction and visualization scheme called the sketch-map can be used to compare and visualize the high-dimensional pairwise data from the SOAP.^115,116 The sketch-map algorithm optimizes a nonlinear objective function that finds a metric, typically Euclidean distance, which best describes the high-dimensional SOAP metrics. Basically, it takes the SOAP data and converts them to a 2D space for human readability. There is an abundance of data that can be gleaned from a sketch-map analysis. One potentially confusing aspect is the axes that have no physical meaning; they could even be flipped without changing the information from the SOAP. Conclusions about the similarity of the local solvent environment can only be made qualitatively from the distances between points. In other words, points that are near each other (clustered) represent similar local solvent environments.

An example of a sketch-map analysis on the solvation free energy of Na⁺ in water is shown in Fig. 4(a).¹⁰⁹ The dots are color coded based on the number of solvent molecules (in this case water) included in the cluster. To analyze the sketch-map, the proximity of dots is examined. For example, the dark slate blue dots (all having four water molecules) in the upper middle of the sketch-map have some space between them and are not very close to each other. This means in the five local minima, clusters generated with four water molecules have different local solvent environments. This communicates that the Gaussian functions placed on the atoms in each cluster overlapped in different amounts and locations representing dissimilarity in the SOAP. To the right in the sketch-map, the dark green dots cluster more for the eight water molecule case signifying the structures and local solvent environment begin to look similar.

The sketch-map places the 12 water clusters very close to each other and even overlapping with the 16-water molecule case. This has important significance for determining the correct cluster size for the cluster cycle. When additional solvent molecules are added, they start populating outer solvation shells and are not captured within the solute’s Gaussian function. Namely, the local solvent environment starts to look the same. Once this occurs, the appropriate cluster size (12 water molecules for Na⁺) is used with the cluster cycle to calculate the solvation free energy to avoid errors from the implicit solvation model’s sensitivity to radius definitions due to incomplete solvation shells. The solvation free energies calculated from the cluster cycle using n water molecules for Na⁺ are shown in Fig. 4(b). Starting from zero water molecules (around −30 kcal/mol error), the cluster cycle solvation energies begin to approach the experimental value of −101.3 kcal/mol.²³ Using the sketch-map selection criteria of 12 water molecules provided the most accurate solvation free energy.

The solvation technique discussed here attempts to capture the following three criteria in an automatable fashion. First, information of the local solvent environment needs to be retained. Preferential solvation must be accounted for to ensure accurate solvation energetics in mixed solvents. Second, explicit simulations should not be required. While substantial information can be gained, the computational cost of explicit simulations and/or their reliance on parameterization prohibits their routine use. Third, no a priori information about the system should be required. For example, knowing the coordination number of the solute drastically reduces the number of applicable systems or demands explicit simulations.

This approach has already been successfully tested in multiple systems. Basdogan and Keith were able to correctly identify an elusive reaction mechanism for the acid-catalyzed Morita–Baylis–Hillman reaction in methanol using this technique.¹¹⁷ Additionally, a complete analysis was also performed for 11 different ionic species’ solvation free energies ranging from 2+ to 2− charges that agreed with experimental values.¹⁰⁹ 

In summary, we reaffirm that the intent of this Perspective is to stimulate development, benchmarking, and application of models used for studying mixed solvents. Over the years, there have been several challenges of predicting solvation free energies for drug and pesticide molecules in water (e.g., SAMPL), and in those, the COSMO-RS model has performed exceedingly well.^118,119 Some computationally intensive RISM implementations have also performed comparably well,¹²⁰ while cluster-continuum modeling has largely not been used except for single ion and small molecule solvation free energies. However, as explained above, modeling intermediates in chemical reactions (especially charged intermediates in electrochemical reactions) requires computational treatments with one of the above procedures. While COSMO-RS and RISM have shown strong performances in pure solvents, it is our perspective that cluster-continuum modeling with QCT and facilitated with the SOAP/sketch-map would be a useful treatment for studies in mixed solvents if and when COSMO-RS and RISM fail. The dual-QCT approach discussed above should make predicting solvation free energies of challenging systems tractable in mixed solvents without explicit simulations.

In this Perspective, we demonstrated that mixed solvents exhibit several theoretical and computational complexities that require attention. We discussed two known, but computationally challenging, characteristics of mixed solvents: the phase potential and preferential solvation. Phase potentials could impact computational studies with certain solvation treatments (e.g., cluster-continuum modeling) or inconsistent use of absolute or real solvation free energies. However, we primarily focused on current solvation models’ potential to capture preferential solvation. Explicit solvation methods should produce the correct solvation free energies in mixed solvents provided sufficient sampling; however, this confidence comes at a significant computational cost of molecular dynamics with quantum chemistry or high quality force fields. Computationally cheaper hybrid methods (e.g., COSMO-RS and RISM) can more efficiently treat mixed solvent environments, but these methods have some known pitfalls, particularly with charged intermediates, that may limit their widespread use.

It is our perspective that additional solvation techniques should require limited a priori information without the need for ab initio molecular dynamics or well-parameterized force fields. Cluster-continuum modeling, when rooted in QCT, has been useful for modeling charged intermediates, but the conventional QCT approach also requires explicit solvation to properly evaluate solvation free energies. We show here that the dual-QCT approach (leveraging global optimization techniques and SOAP/sketch-map analyses) should provide a suitable framework to account for preferential solvation in mixed solvents.

We acknowledge support from the R. K. Mellon Foundation, the U.S. National Science Foundation (Grant Nos. CBET-1653392 and CBET-1705592), the Luxembourg National Research Fund, Sandia’s Laboratory-Directed Research and Development (LDRD) program, and the Center for Integrated NanoTechnology (CINT). Sandia National Laboratories is a multi-mission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE’s NNSA under Contract No. DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in this paper do not necessarily represent the views of the U.S. DOE or the U.S. government.