Theory

Introduction

Cavity and solute density

ddX considers vdW- and SAS-cavities that can be defined by a set of \(M\) balls \(\Omega_j\) each of which is located at \(\textbf{x}_j\in\mathbb R^3\) with radius \(r_j\in\mathbb R\). The region occupied by the solute molecule is given by

\[ \Omega = \cup_{j=1}^M \Omega_j. \]

The charge density of the solute molecule is denoted by \(\rho\) and it can consist of point-charges (partial charges) for classical models or nuclear point-charges with electronic densities for quantum-mechanical models. Within ddX, the solute molecule is completely described by the cavity \(\Omega\) and the charge distribution \(\rho\). The methods used in ddX assume that the charge \(\rho\) is supported in \(\Omega\), i.e. \(\rho(\textbf{x})=0\) for \(\textbf{x}\in\mathbb R^3\backslash\Omega\).

Fundamental solutions and free-space potentials

The free-space potential, i.e. the potential generated by \(\rho\) in vacuum, is then defined by

\[ \Phi(\textbf{x}) = \int_{\Omega} \rho(\textbf{y}) \, G_0(|\textbf{x}-\textbf{y}|) d\textbf{y} \qquad\textrm{ and satisfies }\qquad -\Delta \Phi = 4\pi \rho, \textrm{ in } \mathbb R^3, \]

where \(G_0(r) = \frac{1}{4\pi r} \) denotes the fundamental solution of the Laplace-operator \(\mathcal L_\Delta=-\Delta\). We also introduce the Yukawa-potential \(G_\kappa(r) = \frac{e^{-\kappa r}}{4\pi r} \), which is the fundamental solution of the operator \(\mathcal L_\kappa = -\Delta+\kappa^2\).

Models

The three methods COSMO, PCM and LPB differ by the model of the physical between the solute and the implicit/bulk solvent.

Linearized Poisson-Boltzmann model (LPB)

The LPB-model is the most general one of the three models and the (implicit) solvent is assumed to be homogeneous continuum with relative dielectric permittivity \(\varepsilon_s\) and Debye Hückel constant \(\kappa_s\). Then the LPB-equations are given by: find the potential \(V\) solution to

\[ -\nabla\cdot (\varepsilon\nabla V) + \kappa^2 V = 4\pi \rho, \qquad \textrm{in } \mathbb R^3, \]

where

\[ \varepsilon(\textbf{x}) = \begin{cases} 1 & \textrm{if } \textbf{x}\in\Omega, \\ \varepsilon_s & \textrm{otherwise}, \end{cases} \qquad\textrm{ and }\qquad \kappa(\textbf{x}) = \begin{cases} 0 & \textrm{if } \textbf{x}\in\Omega, \\ \sqrt{\varepsilon_s}\kappa_s & \textrm{otherwise}. \end{cases} \]

ddLPB is based on the following coupled equations: find the reaction-potential \(W=V-\Phi\) and the auxilliary potential \(V_e\) that satisfy

\begin{align*} -\Delta W &= 0, && \textrm{in }\Omega, \\ -\Delta V_e + \kappa_s^2 V_e^2 &= 0, &&\textrm{in }\Omega, \end{align*}

with (nonlocal) boundary conditions

\begin{align*} W + \Phi &= V_e, && \textrm{on }\partial\Omega, \\ V_e + S_\kappa\big(\tfrac{1}{\varepsilon_s} \partial_n W - \partial_n V_e\big) &= -S_\kappa\big(\tfrac{1}{\varepsilon_s} \partial_n \Phi\big), && \textrm{on }\partial\Omega, \end{align*}

where \(S_\kappa\) denotes the single-layer potential with respect to the operator \(\mathcal L_\kappa\) defined by

\begin{align*} (S_\kappa\sigma)(\textbf{s}) &= \int_{\partial\Omega} \sigma(\textbf{s'}) \, G_\kappa(|\textbf{s}-\textbf{s'}|) \, d\textbf{s'}. \end{align*}

Note that \(\partial_n=\frac{\partial}{\partial n}\) stands for the normal derivative.

Polarizable Continuum Model (PCM)

The PCM is the particular case of LPB with \(\kappa_s=0\). ddPCM relies on the following equivalent Integral-Equation formulation: find the charge density \(\sigma\) and the intermediate potential \(\Phi_\varepsilon\) that satisfy

\begin{align*} R_\varepsilon \Phi_\varepsilon &= R_\infty \Phi, &&\textrm{on }\partial\Omega, \\ S \sigma &= -\Phi_\varepsilon, &&\textrm{on }\partial\Omega, \end{align*}

where

\[ R_\varepsilon = 2\pi \frac{\varepsilon_s+1}{\varepsilon_s-1} I - D, \]

and \(S\) and \(D\) denote the single and double-layer boundary operators respectively with respect to the Laplace-operator \(\mathcal L_\Delta\) defined by

\begin{align*} \forall \textbf{s}\in\partial\Omega: \qquad (S\sigma)(\textbf{s}) &= \int_{\partial\Omega} \sigma(\textbf{s'}) \, G_0(|\textbf{s}-\textbf{s'}|) \, d\textbf{s'}, \\ \forall \textbf{s}\in\partial\Omega: \qquad (D\sigma)(\textbf{s}) &= \int_{\partial\Omega} \sigma(\textbf{s'}) \, \frac{\partial}{\partial n(\textbf{s'})} \, G_0(|\textbf{s}-\textbf{s'}|) \, d\textbf{s'}. \end{align*}

Note that \(n(\textbf{s'})\) denotes the outward-pointing unit normal at \(\textbf{s'}\in\partial\Omega\).

COnductor-like Screening MOdel (COSMO)

The COSMO is the particular case of PCM with \(\varepsilon_s=+\infty\), i.e. the solvent is assumed to be a perfect conductor. In consequence, \(V=0\) in \(\mathbb R^3\backslash \Omega\) and the problem reduces to

\[ -\Delta V = 4\pi \rho, \qquad \textrm{in } \Omega, \]

which,for \(W=V-\Phi\), can be reformulated as

\begin{align*} -\Delta W &= 0, && \textrm{in } \Omega, \\ W &= - \Phi, && \textrm{on } \partial\Omega. \end{align*}

The corresponding Integral-Equation formulation is: find the charge density \(\sigma:\partial\Omega\to\mathbb R\) such that

\[ S\sigma = - \Phi, \qquad \textrm{on } \partial\Omega, \]

and \(\sigma\) relates to \(W\) through

\[ \forall \textbf{x}\in\mathbb R^3: \qquad W(\textbf{x}) = \int_{\partial\Omega} \sigma(\textbf{y}) \, G_0(|\textbf{x}-\textbf{y}|) \, d\textbf{y}. \]

General remarks

1. As \(\varepsilon_s\to\infty\), the solution to the PCM becomes the COSMO solution. The first PCM equation yields \(\Phi_\varepsilon \to \Phi\) and the second PCM equation corresponds to the Integral-Equation formulation of COSMO.
2. lpb different formulation

In all three models, the solvation energy \(E_s\) is defined as

\[ E_s=\frac12 \int_{\Omega} \rho(\textbf{x}) \, W(\textbf{x}) \, d\textbf{x}, \]

where \(W=V-\Phi\) denotes the reaction-potential and the force-vector \(F_i\) acting on the \(i\)th atom corresponding to the solvation energy is given by

\[ F_i = -\nabla_{\textbf{x}_i} E_s. \]

Discrete equations in a general framework - overview

The resulting linear system for all three methods can be written in the general form

\[ \tilde{\operatorname{L}} \tilde{\operatorname{X}} = \tilde{\operatorname{F}}, \]

and the (approximate) energy can be computed as

\[ E_s = \frac12 \langle \tilde \Psi, \tilde{\operatorname{X}} \rangle, \]

where \(\langle \cdot, \cdot \rangle\) denotes the scalar product of the two enclosed vectors.

ddCOSMO:

\[ \tilde{\operatorname{L}} = \begin{pmatrix} \operatorname{L} & 0 \\ 0 & 0 \end{pmatrix} \qquad \tilde{\operatorname{X}} = \begin{pmatrix} \operatorname{X} \\ 0 \end{pmatrix} \qquad \tilde{\operatorname{F}} = \begin{pmatrix} -\Phi \\ 0 \end{pmatrix} \qquad \tilde{\Psi} = f(\varepsilon)\begin{pmatrix} \Psi \\ 0 \end{pmatrix} \]

Note: The notation used in the ddCOSMO-literature, see e.g. [1,2,7], uses \(\operatorname{g} = -\Phi\). The references [1–9] also give the definition of the matrix \(\operatorname{L}\) and the vectors \(\operatorname{g} = -\Phi\), \(\Psi\) as well as the scalar function \(f(\varepsilon)\).

The subroutine for matrix-vector multiplications of the \({\operatorname{L}}\) with a vector is

> lx(...)     in ddx_operators.f90

ddPCM:

\[ \tilde{\operatorname{L}} = \begin{pmatrix} \operatorname{L} & \operatorname{I}_d \\ 0 & \operatorname{R}_\varepsilon \end{pmatrix} \qquad \tilde{\operatorname{S}} = \begin{pmatrix} \operatorname{X} \\ \Phi_\varepsilon \end{pmatrix} \qquad \tilde{\operatorname{F}} = \begin{pmatrix} \operatorname{0} \\ \operatorname{R}_\infty\Phi \end{pmatrix} \qquad \tilde{\Psi} = \begin{pmatrix} \Psi \\ 0 \end{pmatrix} \]

Note: The notation used in the references [10,11], uses \(\sigma = \operatorname{X}\). The references [10–13] also give the definition of the matrices \(\operatorname{R}_\varepsilon\), \(\operatorname{R}_\infty\) and the vector \(\operatorname{g} = -\Phi\). \(\operatorname{I}_d\) denotes the identity matrix.

The subroutine for matrix-vector multiplications of the \(\operatorname{R}_\varepsilon\) and \(\operatorname{R}_\infty\) with a vector is

> repsx(...)     in ddx_operators.f90
> rinfx(...)     in ddx_operators.f90

ddLPB:

\[ \tilde{\operatorname{L}} = \tilde{\operatorname{T}} = \begin{pmatrix} \operatorname{A} & 0 \\ 0 & \operatorname{B} \end{pmatrix} + \tilde{\operatorname{C}} \qquad \tilde{\operatorname{C}} = \begin{pmatrix} \operatorname{C}_1 & \operatorname{C}_2 \\ \operatorname{C}_1 & \operatorname{C}_2 \end{pmatrix} \qquad \tilde{\operatorname{X}} = \begin{pmatrix} \operatorname{X} \\ \operatorname{X}_e \end{pmatrix} \qquad \tilde{\operatorname{F}} = \begin{pmatrix} \operatorname{F}_0-\Phi \\ \operatorname{F}_0 \end{pmatrix} \qquad \tilde{\Psi} = \begin{pmatrix} \Psi \\ 0 \end{pmatrix} \]

Note: The notation used in the ddLPB-literature, see e.g. [14,15], uses \(\operatorname{X}_r = \operatorname{X}\) and \(\operatorname{G}_0 = -\Phi\). These references also provide the definition of the matrices \(\operatorname{A}\) (which is a scaled verion of \(\operatorname{L}\)), \(\operatorname{B}\), \(\operatorname{C}_1\), \(\operatorname{C}_2\) as well as the vectors \(\operatorname{F}_0\), \(\operatorname{G}_0 = -\Phi\).

The subroutine for matrix-vector multiplications of the \(\tilde{\operatorname{T}}\), \(\operatorname{A}\), \(\operatorname{B}\) and \(\tilde{\operatorname{C}}\) with a vector is

> tx(...)     in ddx_operators.f90
> ax(...)     in ddx_operators.f90
> bx(...)     in ddx_operators.f90
> cx(...)     in ddx_operators.f90

Gradient computations

The derivative of the energy with respect to a parameter \(\lambda\) is computed using the adjoint method

\begin{align*} E^\lambda_s &= \frac12 \langle \tilde \Psi^\lambda, \tilde{\operatorname{X}} \rangle + \frac12 \langle \tilde \Psi, \tilde{\operatorname{X}}^\lambda \rangle = \frac12 \langle \tilde \Psi^\lambda, \tilde{\operatorname{X}} \rangle + \frac12 \langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle \end{align*}

where \(\tilde{\operatorname{S}}\) is solution to the adjoint linear system

\[ \tilde{\operatorname{L}}^\top \tilde{\operatorname{S}} = \tilde{\Psi}, \]

and

\[ \tilde{\operatorname{H}} = \tilde{\operatorname{F}}^\lambda - \tilde{\operatorname{L}}^\lambda \tilde{\operatorname{X}}. \]

ddCOSMO:

\[ \tilde{\operatorname{L}}^\top = \begin{pmatrix} \operatorname{L}^\top & 0 \\ 0 & 0 \end{pmatrix} \qquad \tilde{\operatorname{S}} = \begin{pmatrix} \operatorname{S} \\ 0 \end{pmatrix} \qquad \tilde{\operatorname{H}} = \begin{pmatrix} -\Phi^\lambda - \operatorname{L}^\lambda \operatorname{X} \\ 0 \end{pmatrix} \]

The subroutine for matrix-vector multiplications of the \(\tilde{\operatorname{L}}^\top\) with a vector is

> lstarx(...)     in ddx_operators.f90

The term \(\langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle\) then writes

\[ \langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle = - \langle \operatorname{S}, \Phi^\lambda\rangle - \langle \operatorname{S}, \operatorname{L}^\lambda\operatorname{X} \rangle \]

The subroutine for the contraction of the differentiated matrix \(\operatorname{L}^\lambda\) with two vectors is

> contract_grad_l(...)     in ddx_gradients.f90

ddPCM:

\[ \tilde{\operatorname{L}}^\top = \begin{pmatrix} \operatorname{L}^\top & 0 \\ \operatorname{I}_d & \operatorname{R}_\varepsilon^\top \end{pmatrix} \qquad \tilde{\operatorname{S}} = \begin{pmatrix} \operatorname{S} \\ -\operatorname{Y} \end{pmatrix} \qquad \tilde{\operatorname{H}} = \begin{pmatrix} -\operatorname{L}^\lambda\operatorname{X} \\ \operatorname{R}^\lambda(\Phi-\Phi_\varepsilon)+\operatorname{R}_\infty\Phi^\lambda \end{pmatrix} \]

Note: The minus sign in front of \(\operatorname{Y}\) appears for consistency with the notation introduced in [12].

The subroutine for matrix-vector multiplications of the \(\operatorname{R}_\varepsilon^\top\) and \(\operatorname{R}_\infty^\top\) with a vector is

> repsstarx(...)     in ddx_operators.f90
> rinfstarx(...)     in ddx_operators.f90

Remarks:

1. There holds \(\operatorname{S}=\operatorname{R}_\varepsilon^\top\operatorname{Y}\).
2. Since \(\operatorname{R}_\infty-\operatorname{R}_\varepsilon=-\frac{4\pi}{\varepsilon_s-1}\) is constant there holds \(\operatorname{R}^\lambda:=\operatorname{R}^\lambda_\varepsilon=\operatorname{R}^\lambda_\infty\).
3. As consequence of 1. and 2., it is advantageous to write

   \[ \langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle = - \langle \operatorname{S}, \operatorname{L}^\lambda\operatorname{X} \rangle + \langle \operatorname{Y}, \operatorname{R}^\lambda(\Phi_\varepsilon-\Phi) \rangle - \langle \operatorname{Q}, \Phi^\lambda \rangle \qquad \operatorname{Q} := \operatorname{S} - \frac{4\pi}{\varepsilon_s-1} \operatorname{Y}. \]

The subroutine for the contraction of the differentiated matrix \(\operatorname{R}^\lambda\) with two vectors is

> contract_gradr(...)     in ddx_gradients.f90

ddLPB:

\[ \tilde{\operatorname{L}}^\top = \tilde{\operatorname{T}}^\top = \begin{pmatrix} \operatorname{A}^\top & 0 \\ 0 & \operatorname{B}^\top \end{pmatrix} + \tilde{\operatorname{C}}^\top \qquad \tilde{\operatorname{C}}^\top = \begin{pmatrix} \operatorname{C}_1^\top & \operatorname{C}_1^\top \\ \operatorname{C}_2^\top & \operatorname{C}_2^\top \end{pmatrix} \qquad \tilde{\operatorname{X}} = \begin{pmatrix} \operatorname{S} \\ \operatorname{S}_e \end{pmatrix} \qquad \tilde{\operatorname{H}} = \begin{pmatrix} \operatorname{F}_0^\lambda-\Phi^\lambda - \operatorname{A}^\lambda - \operatorname{C}_1^\lambda \operatorname{X} - \operatorname{C}_2^\lambda \operatorname{X}_e \\ \operatorname{F}_0^\lambda - \operatorname{B}^\lambda - \operatorname{C}_1^\lambda \operatorname{X} - \operatorname{C}_2^\lambda \operatorname{X}_e \end{pmatrix} \]

The subroutine for matrix-vector multiplications of the \(\tilde{\operatorname{T}}^\top\), \(\operatorname{A}^\top\), \(\operatorname{B}^\top\) and \(\tilde{\operatorname{C}}^\top\) with a vector is

> tstarx(...)     in ddx_operators.f90
> astarx(...)     in ddx_operators.f90
> bstarx(...)     in ddx_operators.f90
> cstarx(...)     in ddx_operators.f90

The term \(\langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle\) then writes

\[ \langle \tilde{\operatorname{S}}, \tilde{\operatorname{H}} \rangle = \langle \operatorname{S}, \operatorname{F}_0^\lambda-\Phi\lambda^\rangle + \langle \operatorname{S}_e, \operatorname{F}_0^\lambda\rangle - \langle \operatorname{S}, \operatorname{A}^\lambda\operatorname{X} \rangle - \langle \operatorname{S}_e, \operatorname{B}^\lambda\operatorname{X}_e \rangle - \langle \tilde{\operatorname{S}}, \tilde{\operatorname{C}}_1^\lambda\tilde{\operatorname{X}} \rangle. \]

The subroutine for the contraction of the differentiated matrices \(\operatorname{A}^\lambda\), \(\operatorname{B}^\lambda\) and \(\tilde{\operatorname{C}}^\lambda\) with two vectors is

> contract_grad_a(...)     in ddx_gradients.f90
> contract_grad_b(...)     in ddx_gradients.f90
> contract_grad_c(...)     in ddx_gradients.f90