Advantage, limitations, and performance insight of commonly used computation chemistry methods for solvation free energy estimation

Advantage, limitations, and performance insight of commonly used computation chemistry methods for solvation free energy estimation Chengze Yang Department of Molecular Engineering, University of Chicago, Chicago, IL, USA Abstract The accurate calculation of free energy is a cornerstone of modern computational chemistry, essential for predicting the spontaneity of chemical reactions and the stability of molecular systems. Solvation free energy ΔG𝑠𝑜𝑙𝑣 is a fundamental thermodynamic quantity that describes the free energy change associated with transferring a solute from the gas phase into a solvent. It is particularly useful for estimating reaction rates and pathways of organic reactions, most of which occur in the liquid phase, and has important applications in organic synthesis, biochemistry, and drug discovery.1 The field has evolved significantly, moving from foundational statistical mechanics methods like Free Energy Perturbation (FEP) and Thermodynamic Integration (TI) to more sophisticated multiscale and hybrid approaches.2,3 A central, persistent challenge is the delicate balance between achieving high accuracy from quantum mechanical descriptions of electronic structure and performing the extensive configurational sampling required to capture entropic effects.4 This report provides an overall review of the theoretical frameworks, implementations, and predictive accuracy of commonly used methods for calculating solvation free energy ΔG𝑠𝑜𝑙𝑣 , assessing alchemical methods, quantum chemistry methods, hybrid QM/MM approaches, while also addressing practical considerations. Finally, we examining the transformative role of emerging methods like machine learning, which are poised to accelerate these calculations, making rigorous, large-scale free energy simulations an option and better balancing the cost-accuracy trade-off that has long constrained the field.5,6 Introduction: The Gibbs free energy, denoted as ΔG, is the fundamental thermodynamic quantity that determines the spontaneity and equilibrium of a chemical process at certain temperature and pressure.7 It is defined by the equation ΔG=ΔH−TΔS, which highlights the critical interplay between the system's change in enthalpy (ΔH) and its change in entropy (ΔS). Solvation free energy ΔG𝑠𝑜𝑙𝑣 is a fundamental thermodynamic property defined as the change in Gibbs free energy when a molecule is transferred from the ideal gas phase to a solvent at a defined temperature and pressure.8 In the context of computational chemistry, the primary challenge is complicated by the need to reconcile two disparate domains: the high-level, first principles quantum mechanical (QM) or various approximation description of electronic structure required for phenomena like bond breaking and formation, and the statistical mechanical requirement for comprehensive configurational sampling to account for entropic contributions arising from molecular motion.9 The inherent trade-off between the accuracy of the electronic structure method and the computational cost of achieving adequate sampling remains the central theme and principal bottleneck in the field of computational free energy.10,11 In practice, the landscape of solvation modeling is broadly divided into two primary theoretical classes: implicit and explicit solvent models. Implicit models, such as the Polarizable Continuum Model (PCM) and Generalized Born (GB)/Poisson-Boltzmann (PB) methods, approximate the solvent as a continuous, polarizable medium.12,13 While computationally efficient, these methods are limited by core approximations, including linear and local dielectric response, which can lead to inaccuracies.12 Conversely, explicit solvent models treat each solvent molecule separately, offering a more physically realistic depiction of solute-solvent interactions at the cost of significant computational expense. Alchemical free energy methods, including Free Energy Perturbation (FEP) and Thermodynamic Integration (TI), have become the standard for rigorous calculations in explicit solvent.13 These methods are widely implemented in molecular dynamics (MD) software such as AMBER, CHARMM, and GROMACS.12,13,14 To address this trade-off, modern approaches often combine techniques. For example, advances of hybrid method between QM and implicit solvation model have led to specialized tools such as AMSOL and COSMO, which provide a well-balanced alternative.15,16 COSMO quantum chemistry–based solvation model achieves high predictive accuracy with minimal empirical parameterization, offering a computationally efficient option for certain applications.4 This review article provides an overview of the most commonly used computation chemistry methods for calculating solvation free energy, from the foundational theories, statistical mechanical principles to the emerging approaches. Beginning with a discussion of implicit and explicit models, grounded in foundational theories of force fields, alchemical and thermodynamic integration methods, and statistical mechanics approaches such as Monte Carlo (MC) and Molecular Dynamics (MD) for configuration generation.17,9,18 It then transitions to quantum chemistry approaches, including Density Functional Theory (DFT), and various post HF methods with an emphasize of their performance in free energy prediction.19 The review subsequently moves to multiscale QM/MM methods, which aim to balance computational accuracy with system size.20 Finally, the report concludes by examining the emerging frontier of machine learning and its transformative implications for the future of free energy calculations.5,6 Implicit Solvent Models These models treat the solvent as a continuous, homogeneously polarizable medium characterized by a dielectric constant.21 This approximation significantly reduces the number of particles and degrees of freedom in a simulation, as solvent molecules can constitute over 90% of the atoms in a system.22 The Polarizable Continuum Model (PCM) and the Conductor-like Screening Model (COSMO) are widely used, with COSMO offering a fast and robust approximation that reduces outlying charge errors compared to PCM.21 Other popular models include the Generalized Born (GB) and Poisson-Boltzmann (PB) models, which are particularly prevalent in biomedical research.23 The general mathematical framework was described as the following. Potential energy is expressed as the sum of the gas-phase energy of the solute and a solvation free energy term, V𝑡𝑜𝑡𝑎𝑙 = V𝑠𝑜𝑙𝑢𝑡𝑒 + ΔG𝑠𝑜𝑙𝑣𝑒.24 The term ΔG𝑠𝑜𝑙𝑣𝑒 is usually composed of Polar Solvation Energy (ΔG𝑝𝑜𝑙𝑎𝑟), which can be estimated by the Poisson-Boltzmann (PB) equation or the Generalized Born (GB) model, and non-polar solvation energy (ΔG𝑛𝑜𝑛𝑝𝑜𝑙𝑎𝑟) can be approximated by Solvent Accessible Surface Area (SASA).25,26,27 The derivation of passion function can be expressed as the following, ∇E = −∇2φ𝑖(𝑟) = ρ(r)𝑓 ε𝑖 , which lead to the general formula for Poisson-Boltzmann (PB) models −∇2ε𝑖 (𝑟)φ𝑖(𝑟) + k2(𝑟)φ𝑖(𝑟) = 4π ρ𝑠𝑜𝑙𝑢𝑡𝑒(𝑟). Generalized Born (GB) Model is a faster, analytical approximation to the PB equation. ΔG𝑝𝑜𝑙𝑎𝑟 = 1 2 ∑q𝑖q𝑗 f𝐺𝐵(r𝑖𝑗, α𝑖 , α𝑗 ) 𝑖,𝑗 , where q𝑖 and q𝑗 are the partial charges of atoms i and j, r𝑖𝑗,is the distance, and f𝐺𝐵 is a mathematical function that depends on effective Born radii, α𝑖 𝑎𝑛𝑑 α𝑗 . ΔG𝑛𝑜𝑛𝑝𝑜𝑙𝑎𝑟 = γ × SASA.25,26,27 However, this efficiency comes at a cost to accuracy. Implicit models fail to account for specific, nonelectrostatic interactions, such as hydrogen bonds or specific ion-solvent interactions.23 The assumption of a linear dielectric response can break down in regions with strong electric fields, such as around nucleic acids or highly charged proteins, where it can induce unrealistically strong polarization.23 Furthermore, the assumption of a local dielectric response, where polarization at a point depends only on the field at that same point, neglects the complex, non-local interactions of real molecular solvents, and also lack the physical representation of entropic effects.23 Finally, the definition of the dielectric interface between the solute and solvent is often ambiguous, leading to a high sensitivity to parameters.23 Explicit Solvent Models In contrast to implicit models, explicit solvent models treat each solvent molecule as a discrete particle with its own coordinates and degrees of freedom, providing a physically resolved description of the solvent environment.28 This approach can capture direct, specific interactions between the solute and solvent, such as hydrogen bonds and complex entropic effects.9 The general mathematical framework for explicit model are the following: V𝑡𝑜𝑡𝑎𝑙 = V𝑏𝑜𝑛𝑑𝑒𝑑 + V𝑛𝑜𝑛−𝑏𝑜𝑛𝑑𝑒𝑑. For bonded interaction, Vbonded = ∑kb( r − r0)2 bonds + ∑kθ( θ − θ0)2 angle + ∑kϕ[1 + cos(nϕ − δ)] dihedral , where r, θ, ϕ are the bond length, angle, and dihedral angle, respectively, kb, kθ, kϕ are the force constants, and r0,θ0 𝑎𝑚𝑑 δ are the equilibrium values.30 V𝑛𝑜𝑛−𝑏𝑜𝑛𝑑𝑒𝑑 = ∑qiqj 4πε0 εr rij, i<j + ∑( Aij rij 12i<j - Bij rij 6 ), where the first term is the Coulomb potential and second term is the Lennard-Jones potential.31 Explicit models are generally more accurate for providing a physically resolved description of the solvent and its role in solvation.28 The primary challenge with explicit solvent models is their high computational cost. The large number of solvent molecules required to model a bulk solution leads to a significant computational burden. 23 This makes explicit solvent simulations prohibitively expensive for high-throughput applications.29 Furthermore, while explicit models are theoretically more rigorous, the accuracy of these models is fundamentally constrained by the accuracy of the energy functions they employ and tied to the quality of the underlying empirical force fields, which are susceptible to polarization.32,33 Hybrid and Discrete-Continuum Models A practical strategy to balance accuracy and computational cost is to combine implicit and explicit solvation approaches. Common implementation of this hybrid strategy is the discrete–continuum model, in which the first solvation shell of solvent molecules is treated explicitly at the atomic level, while the remaining bulk solvent is represented as a continuum.34 The total solvation free energy, ΔG𝑡𝑜𝑡𝑎𝑙, is typically expressed as: ΔG𝑡𝑜𝑡𝑎𝑙 = ΔGsolv explicit + ΔGsolv 𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑢𝑚 - ΔGsolv overlap, where ΔGsolv explicit is free energy contribution from the explicitly treated term, ΔGsolv 𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑢𝑚 is free energy contribution from the implicit bulk, and ΔGsolv overlap is a correction to avoid overlap in regions between explicit and implicit treatments. Force Field: Force fields represent a pragmatic compromise between computational efficiency and physical accuracy, as they sacrifice the explicit treatment of electronic structure that quantum mechanical methods provide. Classical force fields decompose the total potential energy into bonded and non-bonded contributions. The bonded terms typically include harmonic potentials for bond stretching and angle bending, along with periodic torsional potentials for dihedral rotations. The non-bonded interactions consist of electrostatic and van der Waals contributions. Force field parameters are usually derived through a combination of quantum mechanical calculations on small molecular fragments and empirical fitting to reproduce experimental observables such as densities, heats of vaporization, and solvation free energies.35 This parameterization strategy seeks to balance transferability across chemical space with accuracy for specific systems. The AMBER family of force fields employs the restrained electrostatic potential (RESP) charge model and has undergone continuous refinement.36 CHARMM force fields optimize interaction parameters and, in newer versions, explicitly consider polarization effects through the use of lone pairs and virtual sites.37 The OPLS family focuses on reproducing liquid-phase properties, making it particularly relevant for solvation free energy calculations.38 For small organic molecules, the Generalized AMBER Force Field (GAFF) and the CHARMM General Force Field (CGenFF) provide broad coverage of chemical space and transferable parameters. Since water molecules constitute the majority of atoms in explicit solvent simulations, the choice of water model is especially important. Widely used models such as TIP3P, SPC/E, OPC, and TIP4P adopt different strategies for representing water’s charge distribution and consequently exhibit distinct physical properties, including density, dielectric constant, diffusion coefficient, and hydrogen-bonding behavior.39 While modern force fields can achieve good accuracy for neutral organic molecules, larger errors persist for charged species and molecules with complex electronic structure.40 Fixed-charge models cannot capture environmental response for ions and highly polar molecules due to polarization effect. Various approaches have been implemented, including induced dipole models (AMOEBA), fluctuating charge schemes (CHARMM-FQ), and Drude oscillator methods (CHARMM-Drude).37,41 These models allow the charge distribution to respond to the local electric field. Furthermore, force field errors propagate through alchemical transformation pathways, such as FEP and TI calculations amplifying underlying inaccuracies.42 Thus, the parameterization plays an important role in calculation process and the refinement of parameters might be used to address those issues. Monte Carlo and Molecular Dynamics Simulations The accurate calculation of solvation free energy requires not only an appropriate potential energy function but also adequate sampling of the configurational space to capture the statistical mechanical ensemble properties. Two primary simulation paradigms dominate this task: Monte Carlo (MC) and Molecular Dynamics (MD), and both methods serve as the foundation for the free energy perturbation and thermodynamic integration techniques discussed earlier. Monte Carlo Methods Monte Carlo methods generate configurations of a molecular system by making random moves that are accepted or rejected based on statistical criteria. The most widely used algorithm is the Metropolis-Hastings algorithm, which ensures that the generated configurations follow the correct statistical distribution.43 The core principle of MC sampling is based on importance sampling, where configurations are generated with a probability proportional to their Boltzmann weight: P(r𝑁) ∝ exp(−β U(r𝑁)), where β = 1 K𝐵 𝑇, U(r𝑁) is the potential energy of configuration r𝑁, and K𝐵 is the Boltzmann constant. The Metropolis acceptance criterion determines whether a proposed move from state i to state j is accepted: P𝑎𝑐𝑐 (i → j) = min (1, exp [−β (U (rj N) – U (ri N))]). In the context of solvation free energy calculations, MC methods are especially valuable for efficiently sampling discrete changes in chemical identity (e.g., alchemical transformations) and overcoming large energy barriers through specialized move sets, such as particle insertion/deletion or conformational swaps.44 Molecular Dynamics Methods Molecular Dynamics simulations solve Newton's equations of motion for all atoms in the system, generating trajectories that sample phase space according to the microcanonical (NVE), canonical (NVT), or isothermalisobaric (NPT) ensemble when appropriate thermostats and barostats are applied.45 The fundamental equation governing MD simulations is: mid2ri dt2 = -∇iU(r𝑁) = Fi, where Fi is the force on atom i, mi is its mass, d2ri dt2 is its acceleration, and U(r𝑁) is the total potential energy. This produces time-resolved trajectories that naturally capture dynamical processes such as solvation shell reorganization, diffusion, and conformational transitions. MD simulations are the backbone of explicit solvent free energy calculations, providing both structural and dynamical insight into solvation phenomena. Foundational Principles: Theoretical statistical and thermodynamics principles. The computational landscape for modeling solvation is grounded in statistical mechanics and thermodynamics, which provide the link between molecular interactions and macroscopic free energies. In free energy simulations, these principles are applied directly in methods such as Free Energy Perturbation (FEP) and Thermodynamic Integration (TI). Both rely on ensemble averages derived from statistical mechanics to evaluate free energy differences between two states. FEP uses exponential averaging of energy differences, while TI integrates ensemble-averaged derivatives of the Hamiltonian with respect to a coupling parameter. In both cases, accurate sampling of configurations and a reliable potential energy function are essential to prevent error propagation. Free Energy Perturbation (FEP) Free Energy Perturbation (FEP) is a well-established method in computational chemistry based on the principles of statistical mechanics. It is designed to compute the free energy difference, ΔF, between two states, A and B, from molecular dynamics (MD) or Metropolis Monte Carlo (MC) simulations. The central pillar of FEP is the Zwanzig equation, presented as: ΔF(A→B) = FB − FA = −KBT ln〈exp(−EB − EA KBT)〉 , where T is the temperature, KB is the Boltzmann constant, and the angular brackets denote an ensemble average over a simulation run for state A.46 This equation provides a link between the microscopic energy fluctuations of a system and its macroscopic thermodynamic properties. A primary application of FEP is in alchemical transformations, where a molecule is computationally mutated into another.47 The core characteristic is their reliance on non-physical intermediate states that bridge the gap between two physical end states, allowing the simulation to bypass difficult-to-sample regions. FEP has a significant limitation: the calculations only converge properly when the free energy difference between the two states is small.44 This practical constraint means that large transformations, such as the alchemical change of one molecule into a very different one, must be broken down into a series of smaller, independent windows. Thermodynamic Integration (TI) Thermodynamic Integration (TI) is an alternative to FEP that is often more robust for complex transformations.48 It calculates the free energy difference by integrating the derivative of the Hamiltonian with respect to a coupling parameter, λ, as the system is smoothly transitioned from an initial state (A, where λ=0) to a final state (B, where λ=1). With new potential energy function defined as: U(λ) = U(λ) + λ (UB − UA). The free energy of this system is defined as: F (N, V, T, λ) = −KBT ln〈exp(−U(λ) KBT)〉. This approach provides a rigorous thermodynamic path, and the free energy difference is obtained by integrating the ensemble-averaged derivative of the potential energy along this path.49 ΔF(A→B) = ∫∂F(λ)) 𝜕λ dλ = ∫〈∂U(λ)) 𝜕λ 〉 dλ = ∫〈UB(λ) − UA(λ) 〉 dλ Hybrid and Enhanced Sampling Methods The limitations of FEP and TI, particularly their focus on a single perturbation between two end states, have driven the development of more advanced methods. Umbrella Sampling (US) is a well-known enhanced sampling technique that computes a potential of mean force (PMF), which is a free energy map along a specific reaction coordinate.50 A more recent development is Replica-Exchange Enveloping Distribution Sampling (RE-EDS), which belongs to the broader EDS family of methods. A major advantage of RE-EDS is its ability to calculate multiple free energy differences in a single simulation by defining a reference state51, which eliminates the need to explicitly define intermediate states between end points. The development of these methods demonstrates a progression driven by the need for sampling efficiency. The initial challenges of FEP, which requires sufficiently small perturbations and closely spaced λ-windows to converge when end states are very different, led to refinements such as multi-window FEP and the parallel development of TI as an alternative integration scheme. These approaches were subsequently enhanced by techniques like RE-EDS, which can handle multiple simultaneous end states and improve sampling efficiency through its enveloping reference state approach. Quantum Chemistry Approaches: First Principles for Reaction Energetics While alchemical methods are powerful for calculating free energy differences, other approaches leverage the full power of quantum mechanics to compute free energy from first principles. Density Functional Theory (DFT) is one of the most widely used quantum mechanical methods; when combined with statistical thermodynamics, it can provide thermochemical quantities such as enthalpy, entropy, and Gibbs free energy.52 Beyond DFT, higher-level post-Hartree-Fock (post-HF) methods, including Møller–Plesset perturbation theory (MP2), Coupled Cluster (CCSD, CCSD(T)), and Configuration Interaction (CI), allow for a more accurate treatment of electron correlation effects, which are often critical for quantitatively predicting reaction energetics.53 The QM/MM Approach For large and complex systems, neither a full QM treatment nor a purely classical molecular mechanics (MM) simulation is adequate. The QM/MM method offers a powerful solution by combining the strengths of both.54,55 The core concept of a QM/MM simulation is to divide the system into two regions. A small, chemically critical region, where bonds are being broken or formed, is treated with a high-level QM method. The larger surrounding environment—the bulk of the solvent—is modeled using a computationally efficient MM force field. 54,55 The total energy of the system is described by a combined Hamiltonian that accounts for the energy of each region and their interactions. The total energy of the system (E𝐐𝐌/𝐌𝐌 ) is expressed as an additive sum of three terms: E𝐐𝐌/𝐌𝐌 = E𝐐𝐌 (R𝐐𝐌) + E𝐌𝐌 (R𝐌𝐌) + E𝐐𝐌/𝐌𝐌 (R𝐐𝐌, R𝐌𝐌), where E𝐐𝐌(R𝐐𝐌) is the energy of the QM subsystem calculated by solving the electronic Schrödinger equation, and the Hamiltonian is not isolated influenced by the electrostatic field generated by the MM atoms. E𝐌𝐌 (R𝐌𝐌) is the energy of the MM subsystem, calculated with a classical force field. E𝐐𝐌/𝐌𝐌 (R𝐐𝐌, R𝐌𝐌) is the coupling term that describes the interaction between the QM and MM regions, and E𝐐𝐌/𝐌𝐌 (R𝐐𝐌, R𝐌𝐌) = E𝐐𝐌/𝐌𝐌 non−bonded + E𝐐𝐌/𝐌𝐌 bonded, where E𝐐𝐌/𝐌𝐌 non−bonded part handles the long-range interactions, primarily electrostatics and van der Waals forces between the QM and MM atoms, and E𝐐𝐌/𝐌𝐌 bonded is needed when a covalent bond is modified. The QM region's electronic structure is sensitive to the surrounding MM environment, and ensuring a consistent and accurate description of this interaction is complex. The computational cost of QM/MM simulations is considerably higher than pure MM simulations, and this increased cost directly limits the amount of conformational sampling that can be performed. QM/MM simulations have become an invaluable tool for analyzing enzyme catalysis, drug design, and protein engineering.55 Moreover, instead of being a part of the simulation, QM can also be used as a source of high-quality data to define the parameters for the classical force field equations. Emerging Methods: Machine Learning in Free Energy Calculations The recent rise of machine learning (ML) presents a new data driven approach to solve the perennial problem by dramatically accelerating the core calculations. Machine-learned potentials (MLPs) are a promising new class of methods that are trained on a set of reference QM calculations to reproduce a system's potential energy and forces.5 Once trained, a well-parameterized MLP can predict these QM properties with "near-QM level of accuracy" but with much lower computational costs. In hybrid QM/MM simulations, the expensive QM calculation can be replaced by an efficient MLP energy prediction, allowing for more extensive sampling. This direct substitution circumvents the primary computational bottleneck of the QM/MM approach without large scarification on accuracy.56 By providing a fast and accurate surrogate for the DFT potential energy surface, MLPs make the extensive configurational sampling required for rigorous free energy calculations feasible. Moreover, ML can be trained on different data set for direction free energy prediction or used in a hybrid model in a creative way for higher efficiency.6 However, the unique architectures of ML models might limit its transferability and accuracy, its full promise hinges on continued methodological development to ensure the robustness and reliability of these new, ML-accelerated workflows. General Pipeline for Solvation Free Energy Calculation The computational pipeline for solvation free energy calculation diverges significantly based on whether an explicit or implicit solvent model is employed, representing a fundamental trade-off between atomic-level accuracy and computational efficiency. The explicit solvent approach, which begins by placing the solute molecule within a box of explicit solvent molecules, assigning a force field. This assembly then undergoes extensive equilibration through energy minimization, heating, and pressurization to achieve a stable, realistic thermodynamic state. The solute-solvent interactions are gradually turned on or off using a coupling parameter (λ) across many independent simulation windows. Methods like Thermodynamic Integration (TI) or Free Energy Perturbation (FEP) are used within Molecular Dynamics (MD) or Monte Carlo (MC) simulations to sample the system's configurations along this alchemical pathway. The final analysis stage involves integrating the data to compute the free energy change and testing its statistical error.57 In stark contrast, the implicit solvent model offers a streamlined and computationally inexpensive workflow by foregoing the simulation of individual solvent molecules. Here, the solvent is treated as a continuous, polarizable medium characterized by a dielectric constant. It requires only the solute's structure and a few parameters like the solvent's dielectric constant and a probe radius. There is no need for system equilibration or lengthy sampling. Instead, the solvation free energy is calculated directly. The total energy is typically decomposed into a polar component, calculated by solving an electrostatic equation like Poisson-Boltzmann or Generalized Born, and a non-polar component, estimated from the solute's solvent-accessible surface area (SASA). Analysis therefore focuses on the sensitivity of the result to the chosen parameters.58 Hybrid or mixed solvation models combine the strengths of both explicit and implicit approaches. In these models, the solute is surrounded by a limited number of explicit solvent molecules to capture specific, shortrange interactions such as hydrogen bonding or coordination, while the surrounding bulk solvent is treated as a continuous, polarizable medium characterized by a dielectric constant. The explicit region might be equilibrated and sampled using MD or MC simulations. Meanwhile, the implicit region contributes to the solvation free energy through electrostatic and nonpolar terms treated as implicit models.59 average unsigned error (1.10 kcal/mol) using OPLS_2005, outperforming AM1-BCC/GAFF (R² = 0.87, AUE = 1.17 kcal/mol) and CHelpG/CHARMm-MSI (R² = 0.72, AUE = 1.88 kcal/mol).68 Polar compounds, particularly those capable of strong hydrogen bonding, generally showed larger prediction errors. Table 9. Performance of MD/FEP (OPLS_2005) methods on solvation free energies simulation of 239 neutral small molecules System / Dataset Method AUE (kcal/mol) Correlation Coefficient R to experimental data Notes / Reference Solvation free energies for 239 neutral small molecules OPLS_2005/OPLS_2005 1.10 0.94 Shivakumar et al. AM1-BCC/GAFF 1.17 0.87 CHelpG/CHARMm-MSI 1.88 0.72 Comparing the implicit and explicit model for accuracy The accuracy of predicting solvation free energy has been compared between the implicit and explicit model by few researchers. Computational expensive explicit model should offer more accurate and precise simulations for a wider range of objects as we all believe. Indeed, explicit model do offer more precise simulation for most of the case. However, it might not represent all the cases, particularly for solvation energy prediction tasks with heavy parameterization, which in turn heavily influence the performance of the prediction. Some researcher found that the performance of implicit model actually performs better in solvation energy prediction in some cases. While some argues that the improvement of accuracy would be hard for explicit model once approach the threshold of 1 kcal/mol error limit, due to the inherent mechanism required for interaction parameterization of the explicit model. Thus, when balancing the accuracy and computational cost, implicit indeed is of great choice for certain cases. Pande et al. conducted a blind test to assess the predictive accuracy of two distinct computational methods for calculating the aqueous solvation free energies of 17 challenging, polyfunctional small molecules.69 The study compared Poisson-Boltzmann (PB) implicit solvent models against alchemical free energy calculations in explicit solvent (using GAFF force field, TIP3P water, and multiple charge models). Results demonstrated that explicit solvent alchemical calculations provided the slightly highest accuracy, with the best performance coming from the AM1-BCC charge model.69 The implicit solvent PB method performs calculation significantly faster. Both methods faced challenges with specific chemical classes, such as esters and benzamides, revealing limitations in the underlying force fields and parametrizations, also high-lighted the choice of charges models.69 Table 10. Performance of implicit and explicit methods on solvation free energies prediction of 17 mall molecules System / Dataset Method RMS Error (kcal/mol) Notes / Reference 17 Challenging Small Molecules Explicit Solvent: AM1-BCC v2.6A24/GAFF 1.33 ± 0.05 Pande et al. Explicit Solvent: AM1-BCC v1/GAFF 1.53 ± 0.05 Explicit Solvent: AM1-BCC v1/GAFF (Merck-Frosst) 1.71 ± 0.05 Explicit Solvent: RESP HF/6-31G*/GAFF 2.05 ± 0.05 Implicit Solvent: PB (ZAP-9 Radii, AM1-BCC v1) 1.87 ± 0.03 Implicit Solvent: PB (Bondi Radii, AM1-BCC v1) 2.57 ± 0.03 Chen et al. investigated the accuracy of implicit and explicit solvent models for predicting the change in solvation free energy along the reaction coordinate of the Menschutkin reaction in aqueous solution.70 The study compared popular quantum mechanical (QM) implicit solvent models (SMD, SM12, COSMO-RS) against a molecular mechanical (MM) explicit solvent model based on Free Energy Perturbation (FEP) with the CHARMM General Force Field (CGenFF) and TIP3P water model. Contrary to the common assumption that explicit solvent models are inherently more accurate, the standard FEP(MM) approach performed the worst, significantly over-stabilizing the charged transition state in water due to the use of fixed, nonpolarizable Lennard-Jones parameters.70 In contrast, the QM-based implicit solvent models provided more balanced and reasonable estimates.70 The accuracy of the explicit solvent model was dramatically improved by applying end-state corrections, which demonstrating that performance of explicit solvent models requires efforts for accurate parameterization. Table 11. Performance of implicit and explicit methods on prediction of solvation free energy change of Menschutkin reaction in aqueous solution. System / Dataset Method Resulting ΔG (kcal/mol) Notes / Reference Menschutkin Reaction SMD (M06-2X) 26.5 Chen et al. COSMO-RS (BP/TZP) 27.3 SM12 (M06-2X) 29.1 FEP(MM) w/ MM→QM Correction 26.2 FEP(MM) / CGenFF 17.3 Experiment (NH₃ + CH₃I) >23.5 Hybrid Models Hybrid solvation models aim to balance the efficiency of implicit methods with the accuracy of explicit representations. A common implementation is the discrete–continuum model, in which the first solvation shell is treated explicitly to capture short-range, directional interactions such as hydrogen bonding, while the surrounding bulk solvent is modeled as a polarizable continuum to account for long-range electrostatics. QM/MM-based hybrid solvation models extend this concept by treating the solute or reactive site quantum mechanically (QM) to capture electronic effects and chemical reactivity, while the solvent is handled at the molecular mechanics (MM) level with explicit molecules and/or implicit continuum for the bulk. This approach enables accurate treatment of electronic polarization, charge transfer, and local solvation effects without the full computational cost of a fully quantum simulation. König et al. evaluated a novel hybrid quantum mechanics/molecular mechanics (QM/MM) approach for predicting hydration free energies in the SAMPL4 blind challenge.71 The method uses molecular mechanics (MM) to efficiently sample conformational space while employing quantum mechanics (QM) to evaluate potential energies, with free energies determined through Non-Boltzmann Bennett (NBB) reweighting. This approach, termed QM-NBB, addresses a key limitation of pure QM methods by incorporating solute entropy through MD sampling. MM with implicit solvent (MM-GB), MM with explicit TIP3P water (MM-TIP3P), thermodynamic perturbation from MM to QM/MM states (QM/MM-TP), and the full QM/MM-NBB method were tested for predictive performance. Results showed progressive improvement in accuracy when incorporating QM corrections.71 MM-TIP3P achieved better performance than MM-GB, and QM/MM-NBB outperformed both MM methods and the simpler QM/MM-TP approach. Additional analyses using the SMD implicit solvent model demonstrated that reweighting MM trajectories with QM energies can substantially improve predictions, particularly for molecules with multiple protonation states or conformers.71 However, the computational cost of QM/MM-NBB is approximately 50 times higher than pure MM simulations. Table 12. Performance of hybrid QM/MM methods on prediction of hydration free energy for SAMPL4 21 molecules System / Dataset Method RMSD (kcal/mol) Reference SAMPL4 blind subset (21 molecules) MM-GB (implicit) 2.8 König et al. MM-TIP3P (explicit) 2.3 QM/MM-TP (B3LYP/6-31G(d)) 2.0 QM/MM-NBB (B3LYP/6-31G(d)) 1.6 M06-2X/SMD single conformer 1.42 M06-2X/cc-pVTZ/SMD single 1.01 Selected molecules (1, 22, 23, 24) M06-2X/SMD-NBB 1.2 QM/MM-NBB 4.3 SAMPL4 overall Best submission (#561) 1.0 Median submission 1.9 Average (outliers removed) 2.2 Machine Learning Models and Data Driven Methods Data driven methods like AI and machine learning have offered another alternative for solvation free energy modeling. Due to the inherent nature of this data driven method, large range of application could be made. For example, depending on the specific target of prediction, machine learning and neural network could be trained to predict target value of solvation free energy either based on experimental tested value or high-level quantum chemical calculation. Moreover, instead of directly predicting the solvation free energy value, ML and ANN (Artificial neural network) could be train to predicted force field parameter values, deviation of computational chemistry value to the experimental values, or various parameter values in the computational approximation process, which could then be combined with traditional computation chemistry method to offer a fast, low-cost prediction with relative acceptable accuracy. Borhani et al. developed hybrid structure-property relationship (QSPR) models—using Partial Least Squares (PLS) and Multivariate Linear Regression (MLR)—to predict the Gibbs free energy of solvation across a broad range of solute/solvent pairs.72 These hybrid models combine experimental descriptors for solvents with quantum mechanical descriptors for solutes, enabling predictions for 295 solutes in 210 solvents. The MLR model uses only five descriptors: three for solutes (polarizability, LUMO energy, electrostatic acidity) and two for solvents (heat of vaporization, octanol-water partition coefficient).72 The PLS model incorporates all 21 descriptors. Performance metrics demonstrated comparable accuracy to computationally expensive quantum mechanical methods, offering fast predictions without requiring new experimental data or intensive calculations for each compound.72 Hybrid QSPR models combining experimental solvent properties with theoretical solute descriptors provide an efficient alternative to purely computational methods, achieving accuracy comparable to continuum solvation models. Table 13. Performance of QSPR models and various implicit based methods on prediction of free energy of solvation System / Dataset Method RMSE (kcal/mol) MUE (kcal/mol) Notes / Reference MLR (5 descriptors) 0.59 (train), 0.55 (test) 0.44 Borhani et al. 295 solutes, 210 solvents (1777 pairs) PLS (6 latent variables) 0.52 (train), 0.55 (test) 0.43 318 solutes, 91 solvents (2346 pairs) SMD (various levels) — 0.6–1.0 Marenich et al. COSMO-RS — 0.48 Klamt et al. 51 solutes in 3 solvents (77 pairs) SMD/X3LYP 1.11 0.83 Zanith & Pliego SM8/B3LYP 1.08 0.79 MLR (this work) 0.59 0.46 Tested on same dataset PLS (this work) 0.71 0.59 Chung et al. developed three complementary methods to predict Abraham solute parameters and solvation properties at 298 K: a group contribution method (Solute-GC) based model, a machine learning model for solute parameters (Solute-ML), and a direct machine learning model (Direct-ML).73 The Direct-ML is a blackbox machine learning model predictor that predict solvation free energy directly, Solute-ML provides solute parameters which later used to estimate solvation free energy, and Solute-GC which relates parameters to chemical substructures for free energy prediction. The models were trained on extensive compiled databases containing 8,366 solute parameters, 20,253 solvation free energies, and 6,322 solvation enthalpies. Solute-GC uses atom-centered functional groups with ring strain corrections.73 Both ML models employ directed message passing neural networks (D-MPNN) from the Chemprop architecture, with Solute-ML predicting five Abraham parameters that are then used with linear solvation energy relationships (LSERs), while Direct-ML directly predicts solvation free energy and enthalpy from solvent-solute SMILES pairs.73 Performance was evaluated on carefully constructed test sets using random and substructure-based splits to assess generalization to unseen solutes.73 Direct-ML consistently outperformed the other methods, achieving accuracy comparable to advanced quantum chemistry approaches. However, the three methods showed different error patterns across chemical substructures and solvents. Averaging predictions from multiple models reduced errors by outliers, demonstrating the value of data driven ML approach and potential of model combination.73 Table 13. Performance of ML models and GC methods on prediction of free energy of solvation System / Dataset Method Split Type MAE (kcal/mol) RMSE (kcal/mol) Notes / Reference ΔG solvation free energy Direct-ML Random 0.40 0.73 Chung et al. Substructure 0.89 1.32 (20,253 pairs, 5,991 solutes) Solute-ML Random 0.48 0.95 Via LSER Substructure 1.01 1.45 Solute-GC Random 0.63 1.06 Substructure 1.18 1.66 2-model average Random 0.40 0.76 SoluteML + DirectML 3-model average Random 0.42 0.77 All three models Li et al. demonstrated that machine learning can successfully parametrize polarizable force fields using exclusively quantum mechanical calculation data, without any experimental input during training.74 The study focused on methanol using the AMOEBA force field, addressing whether QM-trained models could accurately predict condensed phase properties (density and heat of vaporization) across temperatures. The methodology involved two stages: first, optimizing 44 electrostatic parameters from 4,943 methanol dimer configurations; second, optimizing 10 van der Waals parameters from 1,250 molecular clusters extracted from MD simulations at various thermodynamic states.74 Three QM levels were tested, MP2/6-31G(d,p), DFMP2(fc)/jul-cc-pVDZ, and DFMP2(fc)/jul-cc-pVTZ. The author mention that training on molecular clusters rather than dimers was essential—dimer-only training failed to capture many-body dispersion effects, yielding inaccurate bulk properties.74 The cluster-based approach effectively incorporated many-body effects into parameters through the training process and the best model (DFMP2(fc)/jul-cc-pVTZ with offset) matched or exceeded AMOEBA's performance despite being trained without experimental data.74 It accurately predicted density (0.781 vs 0.786 g/mL) and heat of vaporization (8.98 vs 8.95 kcal/mol) at 298 K, and maintained accuracy across -5°C to 60°C.74 The study establishes that ML strategies can generate physically accurate force fields from pure QM data, provided another option for computational chemistry calculation. Table 14. Performance of ML models predicted force field for the calculation of heat of vaporization Property Method QM Level Density (g/mL) ΔH vap (kcal/mol) Reference/Note Methanol 298 K Experiment — 0.786 8.95 Li et al. AMOEBA — 0.774 9.17 Empirically fitted ML/GA optimized MP2/6-31G(d,p) 0.401 6.54 No offset ML/GA optimized MP2/6-31G(d,p) 0.759 9.44 With offset ML/GA optimized DFMP2(fc)/jul-ccpVDZ 0.568 7.23 No offset ML/GA optimized DFMP2(fc)/jul-ccpVDZ 0.763 9.43 With offset ML/GA optimized DFMP2(fc)/jul-ccpVTZ 0.686 8.58 No offset ML/GA optimized DFMP2(fc)/jul-ccpVTZ 0.781 8.98 Best model Conclusion The accurate calculation of solvation free energy remains a central challenge in computational chemistry, fundamentally defined by the trade-off between the quantum mechanical electronic structure accuracy of the solute description and the computational cost of extensive configurational sampling for entropic contributions. This review has provided the theoretical foundations, computational implementations, and predictive performance of the major approaches currently employed in the field—from foundational alchemical methods like FEP and TI, through implicit and explicit solvent models, to hybrid QM/MM strategies and emerging machine learning techniques. The landscape of available methods presents researchers with a complex set of trade-offs. Implicit solvent models, such as SMD, COSMO-RS, and PB/GB, offer high computational efficiency and are often sufficiently accurate for neutral organic molecules, with errors frequently around 1-2 kcal/mol.1,62 However, they systematically struggle with charged or ionic species and specific, directional interactions like hydrogen bonding due to their inherent approximations of a continuous, lacking solute-solvent interactions, linearly responding dielectric medium. Explicit solvent models, simulated via Molecular Dynamics or Monte Carlo with force fields like GAFF, OPLS-AA, and CHARMM, provide a more physically detailed description and generally achieve higher accuracy, with errors for neutral molecules often approaching a seemingly practical limit of ~1 kcal/mol.40 Their performance, however, is critically dependent on the quality of the force field parameterization, particularly the charge model, and they incur a significantly higher computational cost. However, the "1 kcal/mol error barrier" appears difficult to solve systematically with fixed-charge force fields, suggesting fundamental limitations in the classical approximation for highly polar and charged systems. Hybrid models, including discrete-continuum and QM/MM methods, strategically bridge this gap by combining the accuracy of explicit treatment for critical regions with the efficiency of implicit models for the bulk solvent. While promising, they introduce their own complexities in parametrization and coupling, and their cost can be substantial, especially for QM/MM. After all, computational expense does not automatically guarantee superior accuracy. Several studies demonstrate that well-parameterized implicit models can outperform poorly parameterized explicit approaches, particularly when fixed Lennard-Jones parameters fail to capture polarization effects in charged transition states.58 The quality of underlying physical approximations and parameterization often matters more than the level of theory alone. ML offers a paradigm shift, not only through direct, fast property prediction but also by creating highly accurate, QM-informed force fields and surrogates that accelerate traditional sampling methods. ML-based approaches now achieve accuracy around MAE of 0.4-0.7 kcal/mol for free energy prediction, while dramatically reducing computational demands.74 The demonstrated success of training polarizable force fields exclusively on QM data,75 and the effectiveness of ensemble ML models that combine multiple prediction strategies, suggest that data-driven methods will play an increasingly central role in solvation modeling. However, challenges in transferability and robustness remain. Continued progress is crucial—particularly in improving the treatment of charged species, refining polarization models for greater force field accuracy, and developing robust benchmarks to evaluate emerging machine learning architectures. At present, no single approach universally suits solvation free energy calculations; the optimal choice depends on the system under study, the required accuracy, and the computational resources available. Looking ahead, advances in polarizable force fields, sophisticated hybrid schemes, and machine learning methodologies hold the promise of enabling solvation free energies to be computed with unprecedented speed and reliability, greater generalizability, even for charged and highly complex chemical systems. Reference: 1. Mobley, D. L.; Guthrie, J. P. FreeSolv: A Database of Experimental and Calculated Hydration Free Energies, with Input Files. J. Comput.-Aided Mol. Des. 2014, *28* (7), 711–720. 2. Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, *22* (8), 1420–1426. 3. Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, *3* (5), 300–313. 4. Klamt, A. The COSMO and COSMO-RS Solvation Models. WIREs Comput. Mol. Sci. 2011, *1* (5), 699–709. 5. Noé, F.; Tkatchenko, A.; Müller, K.-R.; Clementi, C. Machine Learning for Molecular Simulation. Annu. Rev. Phys. Chem. 2020, *71*, 361–390. 6. Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine Learning Force Fields. Chem. Rev. 2021, *121* (16), 10142–10186. 7. Atkins, P.; de Paula, J. Physical Chemistry, 8th ed.; W.H. Freeman and Company: New York, 2006; pp 103–109. 8. Ben-Naim, A. Solvation Thermodynamics; Plenum Press: New York, 1987. 9. Chipot, C.; Pohorille, A. (Eds.) Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer Series in Chemical Physics, Vol. 86; Springer: Berlin, Heidelberg, 2007. 10. Christ, C. D.; Fox, T. Accuracy Assessment and Automation of Free Energy Calculations for Drug Design. J. Chem. Inf. Model. 2014, *54* (1), 108–120. 11. Cournia, Z.; Allen, B.; Sherman, W. Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations. J. Chem. Inf. Model. 2017, *57* (12), 2911–2937. 12. Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, *105* (8), 2999–3094. 13. Bashford, D.; Case, D. A. Generalized Born Models of Macromolecular Solvation Effects. Annu. Rev. Phys. Chem. 2000, *51*, 129–152. 14. Van Der Spoel, D. et al. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, *26* (16), 1701– 1718. 15. Cramer, C. J.; Truhlar, D. G. AMSOL: A General-Purpose Quantum Mechanical Solvation Model. In Quantum Mechanical Simulation Methods for Studying Biological Systems; Springer Berlin Heidelberg: Berlin, Heidelberg, 1996; pp 51–61. 16. Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, No. 5, 799– 805. 17. Leach, A. R. Molecular Modelling: Principles and Applications; Prentice Hall: Harlow, England, 2001. 18. Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. 19. Kohn, W. Nobel Lecture: Electronic Structure of Matter—Wave Functions and Density Functionals. Rev. Mod. Phys. 1999, *71* (5), 1253–1266. 20. Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, *48* (7), 1198–1229. 21. Fowles, D. J.; McHardy, R. G.; Ahmad, A.; Palmer, D. S. Accurately Predicting Solvation Free Energy in Aqueous and Organic Solvents beyond 298 K by Combining Deep Learning and the 1D Reference Interaction Site Model. Digit. Discov. 2023, 2 (1), 177–188. 22. Zhang, J.; Zhang, H.; Wu, T.; Wang, Q.; van der Spoel, D. Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents. J. Chem. Theory Comput. 2017, 13 (3), 1034–1043. 23. “Solvation Model Background — APBS 3.1.3 Documentation.” APBS (version 3.1.3). Accessed September 13, 2025. 24. Roux, B.; Simonson, T. Implicit Solvent Models. Biophys. Chem. 1999, *78* (1–2), 1–20. 25. Davis, M. E.; McCammon, J. A. Electrostatics in Biomolecular Structure and Dynamics. Chem. Rev. 1990, *90* (3), 509–521. 26. Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, *112* (16), 6127–6129. 27. Eisenberg, D.; McLachlan, A. D. Solvation Energy in Protein Folding and Binding. Nature 1986, *319* (6050), 199–203 28. “Solvent model.” Wikipedia. Accessed September 13, 2025. https://en.wikipedia.org/wiki/Solvent_model 29. Zhang, J.; Zhang, H.; Wu, T.; Wang, Q.; van der Spoel, D. Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents. J. Chem. Theory Comput. 2017, 13 (3), 1034–1043. 30. Leach, A. R. Molecular Modelling: Principles and Applications; Prentice Hall: Harlow, England, 2001. 31. Jones, J. E. On the Determination of Molecular Fields. —II. From the Equation of State of a Gas. Proc. R. Soc. London, Ser. A 1924, *106* (738), 463–477. 32. Ponder, J. W.; Case, D. A. Force Fields for Protein Simulations. Adv. Protein Chem. 2003, *66*, 27–85. 33. Baker, C. M. Polarizable Force Fields for Molecular Dynamics Simulations of Biomolecules. WIREs Comput. Mol. Sci. 2015, *5* (2), 241–254. 34. Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, *105* (8), 2999–3094. 35. Vanommeslaeghe, K.; Guvench, O.; MacKerell, A. D., Jr. Molecular Mechanics. Curr. Pharm. Des. 2014, *20* (20), 3281–3292. 36. Cornell, W. D. et al. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, *117* (19), 5179–5197. 37. Savelyev, A.; MacKerell, A. D., Jr. All-Atom Polarizable Force Field for DNA Based on the Classical Drude Oscillator Model. J. Comput. Chem. 2014, *35* (16), 1219–1239. 38. Jorgensen, W. L.; Tirado-Rives, J. The OPLS [Optimized Potentials for Liquid Simulations] Potential