Polarizable continuum model (PCM)

# Import necessary packages
from vampyr import vampyr3d as vp
import numpy as np
from scipy.special import erf

Polarizable continuum model (PCM)#

Polarizable continuum models have been used for almost half a century to describe solvation effects in quantum chemistry []. This method has proved itself an effective way to reduce the degrees of freedom necessary to compute solvent effects, as no solvent molecules are explicitly included in the calculation.

PCM treats the solvent as an isotropic continuum characterized by its permittivity, obtained from the solvent’s bulk properties. The solute is described at the quantum level, and their interaction is represented through an electrostatic potential supported on the so-called cavity of the solute. Below, we detail the mathematical formulation of this model and provide a practical example of how it can be implemented.

The cavity is the area where the solute is located and we define it using an interlocking spheres model where the spheres are centered in the atoms of the solute and the radii of the spheres is usually the vdW-radii.

In our case we define a cavity function as

(1)#\[\begin{align} C(\vec{r}) = \begin{cases}0 & \mathbf{r}\text{ inside cavity} \\ 1 & \mathbf{r} \text{ outside cavity}\end{cases} \end{align}\]

In our implementation we defined the boundary of the cavity to be smooth by use of the error function. Below we outline the mathemathical details of our implementation.

Each i-th sphere is defined by a center \(\mathbf{r}_i\) and a radius \(R_i\) and the cavity function is defined as

(2)#\[\begin{equation} C_i(\mathbf{r}) = \frac{1}{2} \left( 1 + \text{erf} \left( \frac{|\mathbf{r} - \mathbf{r}_i| - R_i}{ \sigma} \right) \right), \end{equation}\]

where erf is the error function and \(\sigma\) is a parameter that controls the width of the boundary.

The total cavity function is then defined as

(3)#\[\begin{equation} C(\mathbf{r}) = 1 - \prod_i \left(1-C_i(\mathbf{r})\right). \end{equation}\]

This has been implemented as a class below

class Cavity():
    
    def __init__(self, cav_coords, radii, width):
        self.centers = cav_coords # list of cavity centers. Normally, but not always, nucleus coordinates. 
        self.radii = radii  # list of cavity radii. Normally, but not always, nucleus radii.
        self.sigma = width  # width of the cavity boundary
   

    def __call__(self,  r):
        """
        r is a list of floats of length 3, can be a numpy array as well
        """
        r_vec = np.array(r)
        C = 1.0
        for i, r_i in enumerate(self.centers):
            r_vec_i = np.array(r_i)
            s_i = np.linalg.norm(r_vec_i - r_vec) - self.radii[i]
            O_i = (1.0/2.0)*(1 + erf(s_i/self.sigma))
            C_i = 1 - O_i
            C *= 1 - C_i
        C = 1.0 - C
        return C

The electrostatic potential \(V(\mathbf{r})\) satisfies the generalized Poisson equation (GPE)

(4)#\[\begin{equation} \nabla \cdot (\varepsilon(\mathbf{r}) \nabla V(\mathbf{r})) = -4\pi \rho(\mathbf{r}), \end{equation}\]

where \(\rho(\mathbf{r}\) is the total total density of the solute including both electronic, \(\rho_{el}\) and nuclear \(\rho_{nuc}\) densities,

(5)#\[\begin{equation} \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r}), \end{equation}\]

and \(\varepsilon(\mathbf{r})\) is the position dependent permittivity parametrized using the cavity of the solute where

(6)#\[\begin{align} \varepsilon(\mathbf{r}) = \begin{cases}\varepsilon_0 & \mathbf{r}~\text{inside}~\text{cavity} \\ \varepsilon_r & \mathbf{r}~\text{outside}~ \text{cavity}.\end{cases} \end{align}\]

where \(\varepsilon_0\) is the permittivity of free space and \(\varepsilon_r\) is the relative permittivity of the solvent.

We have defined the permittivity function as follows

(7)#\[\begin{equation} \varepsilon(\mathbf{r}) = \varepsilon_0 \exp\left(\log\left[\frac{\varepsilon_r}{ \varepsilon_0}\right]\left(1 - C(\mathbf{r})\right)\right). \end{equation}\]

This has been implemented as a class below

class Permittivity():
    def __init__(self, cavity, inside=1.0, outside=2.0):
        
        self.C = cavity # instance of cavity class
        self.inside = inside # permittivity of free space
        self.outside = outside # permittivity of solvent
    
    def __call__(self, r):
        C_eval = self.C(r)
        permittivity = self.inside*np.exp((np.log((self.outside/self.inside)))*(1.0 - C_eval))
    
        return permittivity

The electrostatic potential can be divided into two contributions, the solute-solvent interaction reaction potential \(V_R\) (also normally simply called reaction potential), and the vacuum potential \(V_{vac}\),

(8)#\[\begin{equation} V(\mathbf{r}) = V_R(\mathbf{r}) + V_{vac}(\mathbf{r}), \end{equation}\]

where

(9)#\[\begin{equation} \nabla^2 V_{vac}(\mathbf{r}) = -4\pi \rho(\mathbf{r}). \end{equation}\]

The energy contribution from the reaction potential, \(E_R\) is given as

(10)#\[\begin{equation} E_R =\frac{1}{2} \int d\mathbf{r}V_R(\mathbf{r})\rho(\mathbf{r}). \end{equation}\]

Since this energy is only dependent on the reaction potential, we want to directly solve for it.

Resolving the chain rule, substituting in equations 8 and 9, and collecting terms gives us the following expression for the reaction potential

(11)#\[\begin{equation} \nabla^2 V_R(\mathbf{r}) = -4\pi\left( \frac{\rho(\mathbf{r})}{\varepsilon(\mathbf{r})} -\rho(\mathbf{r})\right) - \frac{\nabla \varepsilon(\mathbf{r}) \cdot \nabla V(\mathbf{r})}{\varepsilon(\mathbf{r})} \end{equation}\]

We solve directly for the reaction potential by applying the Poisson operator \(P(\mathbf{r})\)

(12)#\[\begin{equation} V_R(\mathbf{r}) = P(\mathbf{r}) \star\left[ \rho_{eff}(\mathbf{r}) + \gamma_s(\mathbf{r}) \right] \end{equation}\]

where we did the following substitutions

(13)#\[\begin{align} \rho_{eff}(\mathbf{r}) = 4\pi \left(\frac{\rho(\mathbf{r})}{\varepsilon(\mathbf{r})} - \rho(\mathbf{r})\right) \\ ~\\ \gamma_s(\mathbf{r}) = \frac{\nabla \varepsilon(\mathbf{r}) \cdot \nabla V(\mathbf{r})}{\varepsilon(\mathbf{r})} \end{align}\]

which we call respectively the effective molecular density \(\rho_{eff}\) and the surface charge distribution \(\gamma_s\).

Note that the Reaction potential must be solved iteratively since it appears on both sides of equation 12 through the total electrostatic potential \(V\) which is part of \(\gamma_s\) as shown in equation 13.

Below is an implementation of the self consistent procedure for the reaction potential for a single positive charge inside a sphere cavity of unit radius. We start by defining several functions to compute the objects needed for the reaction potential

def constructChargeDensity(positions, charges, width_parameter):
    """Computes the charge density of a set of point charges using a Gaussian
    function. The Gaussian function is centered at the position of the charge
    and the width parameter is the standard deviation of the Gaussian.

    Args:
        positions (list of lists of floats): list of positions of the charges
        charges (list of floats): list of charges
        width_parameter (int, optional): Standard deviation of the Gaussian. Defaults to 1000.

    Returns:
        function: function that computes the charge density at a given position
    """
    charge_density = vp.GaussExp()
    for (pos, charge) in zip(positions, charges):
        beta = width_parameter
        alpha = (beta / np.pi)**(3.0/2.0)
        charge_density.append(vp.GaussFunc(beta=beta, alpha=alpha*charge, position=pos, poly_exponent=[0,0,0]))
    return charge_density


def computeGamma(Derivative_op, V, permittivity, epsilon):
    """Computes gamma_s for a given permittivity and potential.

    Args:
        Derivative_op (vp.ABGVDerivative): Derivative operator
        V (vp.FunctionTree): Potential used to compute gamma_s
        permittivity (vp.FunctionTree): permittivity function of the solvent
        epsilon (float): crop precision

    Returns:
        _type_: _description_
    """
    gamma =  vp.dot(vp.gradient(Derivative_op,V), vp.gradient(Derivative_op, permittivity)) * ( permittivity**(-1))
    gamma = gamma.crop(epsilon)
    return gamma.deepCopy()

Initialize the Multiresolution analysis

k = 5
L = [-10, 10]
MRA = vp.MultiResolutionAnalysis(order=k, box=L)
print(MRA)

================================================================
                    MultiResolution Analysis                    
----------------------------------------------------------------
 polynomial order      : 5
 polynomial type       : Interpolating
----------------------------------------------------------------
 total boxes           : 8
 boxes                 : [          2           2           2 ]
 unit lengths          : [  10.000000   10.000000   10.000000 ]
 scaling factor        : [   1.250000    1.250000    1.250000 ]
 lower bounds          : [ -10.000000  -10.000000  -10.000000 ]
 upper bounds          : [  10.000000   10.000000   10.000000 ]
 total length          : [  20.000000   20.000000   20.000000 ]
================================================================

Initialize all operators that will be used

epsilon = 1e-5
Projection_op = vp.ScalingProjector(mra=MRA, prec=epsilon)
Derivative_op = vp.ABGVDerivative(mra=MRA, a=0.0, b=0.0)
Poisson_op = vp.PoissonOperator(mra=MRA, prec=epsilon)

Compute the density of the solute

charge_coords = [[0.0000000000,    0.0000000000,    0.000000000]]
charges = [1.0]
charge_width = 1000.0
dens = Projection_op(constructChargeDensity(charge_coords, charges, width_parameter=charge_width))

Project the permittivity of the solvent

cav_coords = [[0.0000000000,    0.0000000000,    0.000000000]] # list of cavity centers. Normally, but not always, nucleus coordinates.
cav_radii = [1.0] # list of cavity radii. Normally, but not always, nucleus radii.
boundary_width = 0.2 # width of the cavity boundary

C = Cavity(cav_coords, cav_radii, boundary_width)
perm = Projection_op(Permittivity(C, inside=1.0, outside=2.0))

Construct the effective solute density as defined above and compute the vacuum potential \(V_{vac}\)

rho_eff = (4*np.pi)*((dens * (perm)**(-1))  - (dens))
V_vac =Poisson_op((4*np.pi)*(dens))

Compute the zero-th reaction potential before the iterative procedure

gamma_0 = computeGamma(Derivative_op, V_vac, perm, epsilon)
V_R = Poisson_op(rho_eff + gamma_0)

Iterate throught the application of the Poisson operator to compute the reaction potential

update = 1.0
E_r_old = 0.0
max_iter = 100
print(f"Iter.{' '*2}Norm{' '*12}Update{' '*10}Energy (a.u.){' '*3}Energy update (a.u.)")
print(f"{'-'*75}")

for i in range(max_iter):
    V_tot = V_vac + V_R
                
    # compute the surface charge distribution 
    gamma_s = computeGamma(Derivative_op, V_tot, perm, epsilon)

    # solve the generalized poisson equation for V_R 
    V_R_np1 = Poisson_op((rho_eff) + (gamma_s)) 
    dV_R = V_R_np1 - V_R
    
    update = dV_R.norm()
    
    V_R =  V_R + dV_R
    
    E_r = (1/2)*vp.dot(V_R, dens) # computing energy
    
    dE_r = E_r - E_r_old
    E_r_old = E_r
    
    print(f"{i}{' '*6}{V_R.norm():14.7e}  {update:14.7e}  {E_r:14.7e}  {dE_r:14.7e}") 
    
    if (update < epsilon):
        break
else:
    print("WARNING: SCRF did not converge")

Iter.  Norm            Update          Energy (a.u.)   Energy update (a.u.)
---------------------------------------------------------------------------
     5.4607414e+00   2.8853494e+00  -2.4138496e-01  -2.4138496e-01
     6.1262685e+00   6.6570352e-01  -2.6637517e-01  -2.4990213e-02
     6.0111386e+00   1.1516806e-01  -2.6215285e-01   4.2223193e-03
     6.0270817e+00   1.5949632e-02  -2.6272795e-01  -5.7509706e-04
     6.0252390e+00   1.8435902e-03  -2.6266231e-01   6.5638875e-05
     6.0254219e+00   1.8305719e-04  -2.6266876e-01  -6.4520209e-06
     6.0254060e+00   1.5950934e-05  -2.6266820e-01   5.5755347e-07
     6.0254072e+00   1.2399247e-06  -2.6266825e-01  -4.3037920e-08

The energy can then be computed as follows

(1/2)*vp.dot(V_R, dens)

-0.2626682475802428