Solvation

Implicit solvation models are available for the calculation of the solvation free energies partitioned as

\[\Delta G_{\text{solv}} = \Delta G_{\text{polar}} + \Delta G_{\text{npol}} + \Delta G_{\text{shift}}\]

including the polar contribution \({\Delta G_{\text{polar}}}\) based on electrostatics, the non-polar contribution \({\Delta G_{\text{npol}}}\) based on cavity formation and dispersion, and a constant shift \({\Delta G_{\text{shift}}}\) depending on the thermodynamic state of initial gas and final liquid solution.

ALPB solvation model

The analytical linearized Poisson-Boltzmann (ALPB) model evaluates the polar contribution

\[\Delta G^{\text{ALPB}}_{\text{polar}} = - \frac{1}{2} \left(\frac{1}{\epsilon_{\text{in}}} - \frac{1}{\epsilon_{\text{out}}}\right) \frac{1}{1+\alpha\beta} \sum_{A,B} q_{A} q_{B} \left( \frac{1}{f(R_{AB, a_{A}, a_{B}})} + \frac{\alpha\beta}{\mathcal{A}_{\text{det}}} \right)\]

based on the ALPB constant \({\alpha}\) (set to 0.571214), the solute (\({\epsilon_{\text{in}}=1}\)) and solvent (\({\epsilon_{\text{out}}}\)) dielectric constants combined in \({\beta=\frac{\epsilon_{\text{in}}}{\epsilon_{\text{out}}}}\), atomic partial charges \({q_{A/B}}\), and the electrostatic size of the solute \({\mathcal{A}_{\text{det}}}\). [1] \({f(R_{AB}, a_{A}, a_{B})}\) is the interaction kernel with the Born radii \({a_{A/B}}\) and can take two forms, either

\[f(R_{AB}, a_{A}, a_{B}) = \left( R_{AB}^2 + a_{A} a_{B} \exp\left[-\frac{R_{AB}^2}{4 a_{A} a_{B}} \right] \right)^{\frac{1}{2}}\]

proposed by Still (default in GBSA), or the more recent P16 kernel (default for ALPB):

\[f(R_{AB}, a_{A}, a_{B}) = R_{AB} + \sqrt{a_{A} a_{B}} \left(1+\frac{1.028 R_{AB}}{16 \sqrt{a_{A} a_{B}}} \right)^{16}\]

For specific polar interactions, an atom-wise hydrogen bonding correction is introduced:

\[\Delta G^{\text{HB}}_{\text{polar}} = \sum \Delta G^{\text{HB}}_{\text{A}}\]

In addition to the polar contribution, the non-polar contribution is included with a cavity dispersion solvation term (CDS) based on the atomic surface tension \(\gamma_{A}\) and the solvent-accessible surface area (SASA) \(\sigma_{A}\):

\[\Delta G^{\text{CDS}}_{\text{npol}} = \sum_{A} \gamma_{A} \sigma_{A}\]

An additional empirical constant shift is applied to the solvation free energy. A solution state correction can be activated but is not included by default.

GBSA solvation model

The generalized Born solvation model (GBSA) is a simplified version of ALPB in the limit of an ideal conductor environment (\({\epsilon_{\text{out}}}\rightarrow \infty\) and \({\beta\rightarrow 0}\)). As for ALPB, CDS and a constant shift shift are applied, while a solution state correction can be activated (only if the solvent is specified by name).

Domain decomposition continuum solvation models

The domain decomposition solvation models implemented through the ddX library provide efficient formulations of the classical COSMO, CPCM, and PCM continuum models. The physical picture remains unchanged here, with the solute being embedded in a molecular cavity \(\Omega\) that is surrounded by a homogeneous polarizable continuum. The solvent response is represented by an apparent surface charge distribution \(\sigma\) on the cavity boundary \(\Gamma\), which generates the reaction field potential \(\phi^\sigma\).

In the COSMO model, the solvent is first approximated by an ideal conductor, \(\epsilon \rightarrow \infty\), such that the total electrostatic potential must vanish at the cavity boundary. Therefore, the reaction field potential exactly cancels the molecular electrostatic potential \(\phi^{\rho_0}\) on \(\Gamma\). Inside the cavity, the reaction field potential is harmonic, such that the following boundary value condition is fulfilled.

\[\begin{split}\begin{aligned} -\nabla^2 \phi^\sigma(\mathbf{r}) &= 0, \qquad \mathbf{r} \in \Omega, \\ \phi^\sigma(\mathbf{s}) &= -\phi^{\rho_0}(\mathbf{s}), \qquad \mathbf{s} \in \Gamma . \end{aligned}\end{split}\]

The central idea of the domain decomposition formulation is to rewrite this global boundary value problem as a set of coupled local problems. For this purpose, the molecular cavity is decomposed back into its initial union of atom-centered spheres,

\[\Omega = \bigcup_{A=1}^{N_\mathrm{at}} \Omega_A .\]

The surface of each atomic sphere is further separated into an exposed (e) and buried (interior, i) part,

\[\Gamma_A = \Gamma_A^\mathrm{e} \cup \Gamma_A^\mathrm{i} .\]

Thus, the boundary value problem can be rewritten as a set of \(N_\mathrm{at}\) coupled sphere-wise boundary value problems,

\[\begin{split}\begin{aligned} -\nabla^2 \phi_A(\mathbf{r}) &= 0, \qquad \mathbf{r} \in \Omega_A, \\ \phi_A(\mathbf{s}) &= -\phi^{\rho_0}(\mathbf{s}), \qquad \mathbf{s} \in \Gamma_A^\mathrm{e}, \\ \phi_A(\mathbf{s}) &= \frac{1}{|\mathcal{N}(A,\mathbf{s})|} \sum_{B \in \mathcal{N}(A,\mathbf{s})} \phi_B(\mathbf{s}), \qquad \mathbf{s} \in \Gamma_A^\mathrm{i}, \end{aligned}\end{split}\]

where the potential on the buried surface is obtained from the potentials of all intersecting spheres at the considered surface point. Here, \(\mathcal{N}(A,\mathbf{s})\) defines the list of neighbouring spheres of sphere \(A\) at surface point \(\mathbf{s}\). These local problems can be expressed in the known integral form by representing the reaction field potential on each sphere through a single-layer potential,

\[\phi_A(\mathbf{s}) = \int_{\Gamma_A} \frac{\sigma_A(\mathbf{s}')}{|\mathbf{s}-\mathbf{s}'|} \,d\mathbf{s}' .\]

The local surface charge distributions \(\sigma_A(\mathbf{s})\) are then expanded in real spherical harmonics,

\[\sigma_A(\mathbf{s}) = \sum_{\ell=0}^{\ell_\mathrm{max}} \sum_{m=-\ell}^{\ell} [X_A]_{\ell m} Y_{\ell m}(\mathbf{s}) ,\]

and integration is performed numerically on the unit sphere using Lebedev quadratures. This finally results in a linear system

\[\mathbf{L} \mathbf{X} = \mathbf{g} ,\]

where the vector \(\mathbf{X}\) contains the expansion coefficients of the local apparent surface charges, \(\mathbf{g}\) is the corresponding spherical harmonics expansion of the free-space potential generated by the solute’s charge distribution, and \(\mathbf{L}\) is a naturally block-sparse coupling matrix. This system is solved iteratively in a very efficient manner, for instance by Jacobi or DIIS procedures, which allows the solvation contribution to be computed with linear-scaling behavior.

In tblite, the solute charge density is represented by atom-centered point charges. The final solvation free energy is therefore obtained through a dot product of the vector containing the partial charges, \(\mathbf{\Psi}\), and the expansion coefficients,

\[\Delta G_\mathrm{polar} = \frac{1}{2} f(\epsilon,x) \langle \mathbf{X}, \mathbf{\Psi} \rangle ,\]

where the factor

\[f(\epsilon,x) = \frac{\epsilon - 1}{\epsilon + x}\]

scales the surface charges, or equivalently the energy, to account for the fact that the real dielectric permittivity of the solvent is finite. As established, it further distinguishes between COSMO, \(x = 0.5\), and CPCM, \(x = 0\).

In PCM, the continuum is not treated as an ideal conductor, leading to both single- and double-layer potential contributions. Addressing the general integral equation formalism (IEF) of PCM by the same sphere-wise framework used for ddCOSMO/ddCPCM, the ddPCM equation can be written as

\[\hat{\mathcal{R}}_\epsilon \hat{\mathcal{S}}\sigma = -\hat{\mathcal{R}}_\infty \phi^{\rho_0} ,\]

with the operators

\[\begin{split}\begin{aligned} \hat{\mathcal{R}}_\epsilon &= \frac{2\pi}{f(\epsilon)}\hat{1} - \hat{\mathcal{D}}, \\ \hat{\mathcal{R}}_\infty &= 2\pi\hat{1} - \hat{\mathcal{D}}, \end{aligned}\end{split}\]

and the single- and double-layer operators \(\hat{\mathcal{S}}\) and \(\hat{\mathcal{D}}\). The ddPCM model solves this problem by introducing an intermediate potential \(\phi_\epsilon\), which splits the problem into

\[\hat{\mathcal{R}}_\epsilon \phi_\epsilon = \hat{\mathcal{R}}_\infty \phi^{\rho_0}\]

and

\[\hat{\mathcal{S}}\sigma = -\phi_\epsilon .\]

While the second equation has the same structure as the single-layer problem solved in ddCOSMO/ddCPCM, the first equation stands out as the additional PCM problem.

For further reading, the reader is referred to comprehensive reviews [2][3] and the references cited therein.

Solution state correction

For solvation free energies, the state of the inital gas and final liquid solution can be changed with a solution state correction. By default no solution state correction is applied (gsolv, default), which is comparable with most other solvation models (SMD, COSMO-RS, …). For normal production runs, the option bar1mol should be used. For explicit comparisons with reference state corrected COSMO-RS, the reference option should be used (includes solvent-specific correction for infinite dilution). Solution state correction is available for the ALPB and GBSA solvation models.

Name	Definition
gsolv (default)	1 L of ideal gas and 1 L of liquid solution
bar1mol	1 bar of ideal gas and 1 mol/L liquid solution
reference	1 bar of ideal gas and 1 mol/L liquid solution at infinite dilution

Literature

[1]

Sebastian Ehlert, Marcel Stahn, Sebastian Spicher, and Stefan Grimme. Robust and efficient implicit solvation model for fast semiempirical methods. J. Chem. Theory Comput., 17(7):4250–4261, 2021. doi:10.1021/acs.jctc.1c00471.

[2]

Benjamin Stamm, Louis Lagardère, Giovanni Scalmani, Paolo Gatto, Eric Cancès, Jean-Philip Piquemal, Yvon Maday, Benedetta Mennucci, and Filippo Lipparini. How to make continuum solvation incredibly fast in a few simple steps: a practical guide to the domain decomposition paradigm for the conductor-like screening model. Int. J. Quantum Chem., 119(1):e25669, 2019. doi:10.1002/qua.25669.

[3]

Michele Nottoli, Michael F. Herbst, Aleksandr Mikhalev, Abhinav Jha, Filippo Lipparini, and Benjamin Stamm. Ddx: polarizable continuum solvation from small molecules to proteins. WIREs Comput. Mol. Sci., 14(4):e1726, 2024. doi:10.1002/wcms.1726.