12.2.10 Poisson Boundary Conditions

The most general formulation of Poisson’s equation, for a dielectric function rather than a dielectric constant, is

$$\widehat{\mathbf{\nabla }}\mathbf{\cdot }\left[\epsilon (𝐫)\widehat{\mathbf{\nabla }}{\phi }_{\mathrm{tot}}(𝐫)\right]=-4\pi {\rho }_{\mathrm{sol}}(𝐫).$$

(12.13)

The solute’s charge density will be separated into nuclear and electronic components,

$${\rho }_{\mathrm{sol}}(𝐫)={\rho }_{\mathrm{nuc}}(𝐫)+{\rho }_{\mathrm{elec}}(𝐫).$$

(12.14)

The total electrostatic potential ${\phi }_{\mathrm{tot}}$ is comprised of the solute’s electrostatic potential ${\phi }_{\mathrm{elec}}+{\phi }_{\mathrm{nuc}}$, which is obtained from the density in Eq. (12.14), plus the polarization potential induced in the medium:

$${\phi }_{\mathrm{tot}}(𝐫)={\phi }_{\mathrm{elec}}(𝐫)+{\phi }_{\mathrm{nuc}}(𝐫)+{\phi }_{\mathrm{pol}}(𝐫).$$

(12.15)

To compute ${\rho }_{\mathrm{elec}}$, the electronic contribution to the electrostatic potential in vacuum is first evaluated on a Cartesian grid,

$${\phi }_{\mathrm{elec}}({𝐫}_i)=\sum _{\mu \nu }^{N_{\mathrm{basis}}}{P}_{\mu \nu }⟨g_{\mu }(𝐫)\left|\frac{1}{|𝐫-{𝐫}_i|}\right|g_{\nu }(𝐫)⟩.$$

(12.16)

Here P is the one-electron density matrix, $g_{\mu }$ and $g_{\nu }$ are atom-centered Gaussian basis functions, and the integration is over $𝐫$. The electronic part of the solute’s charge density is then obtained from

$${\rho }_{\mathrm{elec}}(𝐫)=-\frac{1}{4\pi }{\widehat{\nabla }}^2{\phi }_{\mathrm{elec}}(𝐫),$$

(12.17)

where in practice the Laplacian of ${\phi }_{\mathrm{elec}}$ is evaluated using an eighth-order finite difference scheme.¹⁸⁸ The nuclear charge density ${\rho }_{\mathrm{nuc}}$ is described classically,

$${\rho }_{\mathrm{nuc}}(𝐫)=-\sum _A^{N_{\mathrm{atoms}}}{Z}_A\delta \left(𝐫-{𝐑}_A\right),$$

(12.18)

but to avoid numerical problems the nuclear charges are smeared out over one Cartesian grid voxel in practice. Operationally, this is achieved by adding $Z_A/dV$ to the grid point nearest ${𝐑}_A$, where $dV$ is the volume of a Cartesian voxel. The exact nuclear electrostatic potential in vacuum is computed according to

$${\phi }_{\mathrm{nuc}}(𝐫)=-\frac{1}{4\pi }\sum _A^{N_{\mathrm{atoms}}}\frac{Z_A}{|𝐫-{𝐑}_A|}.$$

(12.19)

The Dirac delta function in Eq. (12.18) has dimensions of (volume)${}^{-1}$ and therefore incorporates the aforementioned factor of $d{V}^{-1}$

Having computed the electronic and nuclear charge densities and the electrostatic potentials in vacuum, one has in hand an exact solution of Eq. (12.13) for the case where $\epsilon (𝐫)=1$. The induced polarization charge density and corresponding electrostatic potential are then obtained iteratively, following closely the procedure outlined in Refs. 263 and 37. Equation (12.13) is rewritten to appear as Poisson’s equation in vacuum with an additional source charge density:

(12.20)

where the polarization charge density is

$${\rho }_{\mathrm{pol}}(𝐫)={\rho }_{\mathrm{iter}}(𝐫)+\left(\frac{1-\epsilon (𝐫)}{\epsilon (𝐫)}\right){\rho }_{\mathrm{sol}}(𝐫).$$

(12.21)

The iterative charge density ${\rho }_{\mathrm{iter}}(𝐫)$ takes the form

$${\rho }_{\mathrm{iter}}^{(k+1)}(𝐫)=\frac{\widehat{\mathbf{\nabla }}\ln \epsilon (𝐫)\mathbf{\cdot }\widehat{\mathbf{\nabla }}{\phi }_{\mathrm{tot}}^{(k)}(𝐫)}{4\pi }$$

(12.22)

at iteration $k+1$. This function is nonzero only in transition regions where the dielectric is being interpolated from vacuum to solvent. Note also that the second term in Eq. (12.21) vanishes wherever $\epsilon (𝐫)=1$. The quantity ${\rho }_{\mathrm{iter}}(𝐫)$ is iterated to convergence by performing updates ${\rho }_{\mathrm{iter}}^{(k)}(𝐫)\to {\rho }_{\mathrm{iter}}^{(k+1)}(𝐫)$ using Eq. (12.22), taking ${\phi }_{\mathrm{tot}}^{(0)}(𝐫)$ to be the solute’s electrostatic potential in vacuum. Each time ${\rho }_{\mathrm{iter}}(𝐫)$ is updated, a new polarization charge density is generated and Eq. (12.20) is solved to obtain a new total electrostatic potential. To stabilize the iterative updates of ${\rho }_{\mathrm{iter}}(𝐫)$, we use a damping procedure

(12.23)

rather than using Eq. (12.22) as written. Following Refs. 263 and 37, we use $\eta =0.6$.

Once the total charge density and electrostatic potential are iterated to self-consistency, the free energy of solvation,

$$G_{\mathrm{pol}}=\frac{1}{2}\int 𝑑\text{𝐫}{\phi }_{\mathrm{pol}}(\text{𝐫}){\rho }_{\mathrm{sol}}(𝐫),$$

(12.24)

is computed, where

$${\phi }_{\mathrm{pol}}(\text{𝐫})={\phi }_{\mathrm{tot}}(\text{𝐫})-{\phi }_{\mathrm{sol}}(𝐫).$$

(12.25)

Since ${\phi }_{\mathrm{tot}}$ and ${\rho }_{\mathrm{tot}}$ are computed self-consistently at each SCF iteration, the polarization effects must be incorporated into the Fock matrix. This is accomplished by adding the correction term $\mathrm{\Delta }\widehat{F}=\delta G_{\mathrm{pol}}/\delta {\rho }_{\mathrm{elec}}$. The matrix form of this correction in the atomic orbital basis is

$$\mathrm{\Delta }{F}_{\mu \nu }=\int 𝑑\text{𝐫}{\phi }_{\mathrm{pol}}(\text{𝐫})g_{\mu }(\text{𝐫})g_{\nu }(\text{𝐫}).$$

(12.26)

For additional details and a full description of the algorithm, see Ref. 188.

Nonequilibrium solvation within the framework of PCMs is discussed in Section 12.2.2.3, and that methodology^{635, 1039} has been adapted for vertical ionization processes in conjunction with Poisson boundary conditions.^{189, 188} A state-specific approach has been implemented that requires solution of a nonlinear Schrödinger equation, Eq. (12.6). For the reference state, corresponding to $i=0$ in Eq. (12.6), this is equivalent to solving Eq. (12.13) for the induced polarization potential and computing the solvation free energy and Fock matrix correction according to Eqs. (12.24) and (12.26), respectively. For the ionized state, denoted by a subscript $i$, the solvent polarization is partitioned into fast and slow components. Within the Marcus partitioning scheme,¹⁰³⁹ Eq. (12.20) becomes¹⁸⁸

$$-\frac{1}{4\pi }{\widehat{\nabla }}^2{\phi }_{\mathrm{tot},i}(𝐫)={\rho }_{\mathrm{sol},i}(𝐫)+{\rho }_{\mathrm{pol},0}^{\mathrm{slow}}(𝐫)+{\rho }_{\mathrm{pol},i}^{\mathrm{fast}}(𝐫),$$

(12.27)

where the fast polarization charge density of the ionized state, ${\rho }_{\mathrm{pol},i}^{\mathrm{fast}}$, is given by

$${\rho }_{\mathrm{pol},i}^{\mathrm{fast}}(𝐫)={\rho }_{\mathrm{iter},i}(𝐫)+\left(\frac{1-{\epsilon }_{\mathrm{opt}}(𝐫)}{{\epsilon }_{\mathrm{opt}}(𝐫)}\right)\left[{\rho }_{\mathrm{sol},i}(𝐫)+{\rho }_{\mathrm{pol},0}^{\mathrm{slow}}(𝐫)\right]$$

(12.28)

and the fast polarization potential, ${\phi }_{\mathrm{pol},i}^{\mathrm{fast}}$, is

$${\phi }_{\mathrm{pol},i}^{\mathrm{fast}}(\text{𝐫})={\phi }_{\mathrm{tot},i}(\text{𝐫})-{\phi }_{\mathrm{sol},i}(𝐫)-{\phi }_{\mathrm{pol},0}^{\mathrm{slow}}(𝐫).$$

(12.29)

Subscripts $i$ and 0 refer to the ionized state and the reference state, respectively. The quantity ${\epsilon }_{\mathrm{opt}}$ in Eq. (12.28) is the solvent’s optical dielectric constant, and ${\rho }_{\mathrm{iter},i}$ is computed using Eq. (12.23) but with a dielectric function constructed using ${\epsilon }_{\mathrm{opt}}$ and with ${\phi }_{\mathrm{tot},i}^{\mathrm{fast}}$ substituted in place of ${\phi }_{\mathrm{tot}}$. For the Marcus partitioning scheme, the slow component of the induced solvent polarization charge density is computed as the three-dimensional (volume charge rather than surface charge) analogue of Eq. (12.8).

A calculation using Poisson boundary conditions is requested by setting SOLVE_PEQ = TRUE in the $rem section. For computational efficiency, solvent effects from the Poisson equation solver are not incorporated until the error in the vacuum SCF calculation falls below a specified threshold, taken to be 10${}^{-3}$ a.u. by default and controllable using the $rem variable PEQ_SWITCH.

Note:  When SOLVE_PEQ = TRUE, the keywords SYM_IGNORE and NO_REORIENT are also set to TRUE. Poisson’s equation is discretized on a three-dimensional grid and the energy need not be rotationally invariant unless the grid spacing is sufficiently small.

Further job control for the Poisson equation solver (PEqS) routines is set up in two sections: $peqs_grid and $peqs. The former is used to define the Cartesian discretization grid, as follows:

$peqs_grid
   DimX <Nx> <Xmin> <Xmax>
   DimY <Ny> <Ymin> <Ymax>
   DimZ <Nz> <Zmin> <Zmax>
$end

The Cartesian grid must be orthorhombic and Nx, Ny, and Nz refer to the number of grid points for the $x$, $y$, and $z$ axes.

Note:  A fourth-order multigrid algorithm is used to accelerate convergence of the conjugate gradient PEqS routine that solves Eq. (12.20).¹⁸⁸ This requires that the number of grid points be odd, with the additional constraint that the quantities $N_x\mathrm{-}\mathrm{1}$, $N_y\mathrm{-}\mathrm{1}$, and $N_z\mathrm{-}\mathrm{1}$ must all be divisible by 8.

Values in the last two columns of the $peqs_grid section specify the minimum and maximum values of $x$, $y$, and $z$, in units of Ångstroms. The length of the grid in the $x$ direction is $L_x=x_{\max }-x_{\min }$ and the grid spacing is $\mathrm{\Delta }x=L_x/(N_x-1)$. The volume element discussed in the context of Eq. (12.18) is $dV=\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z$.

The $peqs section controls other aspects of a PEqS calculation, including construction of the solute cavity and other required parameters for interfacial or nonequilibrium solvation. The format is:

$peqs
   <Keyword>  <parameter/option>
$end

Available keywords are described below.

In conventional PCMs the solute cavity is a rigid two-dimensional surface constructed from a union of atomic spheres, the radii of which are generally take to be equal to the van der Waals radius ($r_{\mathrm{vdW}}$) scaled by a factor of 1.2; see Section 12.2.3.³⁷³ This cavity is rigid in the sense that once the atomic coordinates are specified, it remains unchanged during the SCF cycles. (The isodensity cavity discussed in Section 12.2.5 is an exception, but is not available for the PEqS method.) Furthermore, there is an abrupt and discontinuous change in the dielectric constant at the solute/continuum boundary. A three-dimensional analogue of this cavity, using continuous and differentiable spherically-symmetric error functions, has been implemented for PEqS calculations, following the procedure in Ref. 263. The dielectric function, which depends parametrically on the nuclear positions ${𝐑}_A$, is

$$\epsilon (𝐫;\{{𝐑}_A\})=\left({\epsilon }_{\mathrm{solvent}}-1\right)\left\{\prod _A^{N_{\mathrm{atoms}}}h(d_A,\mathrm{\Delta };|𝐫-{𝐑}_A|)\right\},$$

(12.30)

where

$$h(d_A,\mathrm{\Delta };|𝐫-{𝐑}_A|)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{|𝐫-{𝐑}_A|-d_A}{\mathrm{\Delta }}\right)\right].$$

(12.31)

Equation (12.30) smoothly interpolates between $\epsilon =1$ and $\epsilon ={\epsilon }_{\mathrm{solvent}}$, over a length scale of $\approx 4\mathrm{\Delta }$. The value of $\mathrm{\Delta }$ is specified using the RigidScale keyword. The quantity $d_A$ in Eqs. (12.30) and (12.31) sets radius of the atomic sphere for atom $A$. By default, $d_A=1.2r_{\mathrm{vdW}}$ but this can be controlled as described below.

Setting the SoluteCavity keyword to Spherical requests the construction of a spherical cavity around the solute, and the dielectric is smoothly interpolated from vacuum to solvent using a hyperbolic tangent function:¹⁸⁹

$$\epsilon (r)=\frac{1}{2}\left\{({\epsilon }_{\mathrm{solvent}}+1)+({\epsilon }_{\mathrm{solvent}}-1)\text{tanh}\left[\alpha (r-r_{\mathrm{mid}})\right]\right\}.$$

(12.32)

Here, $r_{\mathrm{mid}}$ is the distance where the dielectric assumes the value $\epsilon (r_{\mathrm{mid}})=({\epsilon }_{\mathrm{solvent}}+1)/2$. The value of $r_{\mathrm{mid}}$ is taken to be the sum of the sphere radius, $R$, and half of the interpolation length, $L$: $r_{\mathrm{mid}}=R+L/2$. The sphere radius and interpolation length are controlled with the SphereRadius and InterpolLength keywords described below. The parameter $\alpha$ controls the sharpness of the switching process, with $\alpha =4/L$ by default.

By setting Interface = True, the dielectric function on the Cartesian grid is further modified to mimic a solvent/vacuum interface. First, the solute cavity is constructed as specified by the SoluteCavity keyword as discussed above. Then, the dielectric function is smoothly interpolated in the $z$ direction across the Gibbs dividing surface ($z\equiv z_{\mathrm{GDS}}$) using the following hyperbolic tangent switching function:^{189, 188}

$$\epsilon (z)=\frac{1}{2}\left\{({\epsilon }_{\mathrm{solvent}}+1)+(1-{\epsilon }_{\mathrm{solvent}})\text{tanh}\left[\beta (z-z_{\mathrm{GDS}})\right]\right\}.$$

(12.33)

This interpolates the dielectric function over the interface length $L_{\mathrm{interface}}$, centered at the Gibbs dividing surface, $z\equiv z_{\mathrm{GDS}}$. Both of these values are controlled by the keywords InterfaceLength and GibbsDS, respectively. Similar to the parameter $\alpha$ used in the spherical cavity construction, the parameter $\beta =4/L_{\mathrm{interface}}$ controls the sharpness of the interpolation across the Gibbs dividing surface.

The following example computes the solvation free energy of a water molecule immersed in water. The Cartesian grid is cubic with a side length of 15.0 Å and 73 grid points in each direction. A rigid vdW cavity is used based on scaled vdW radii for the atomic spheres. The dielectric is not set, so takes the default value of 78.39.

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          wB97X-V
   BASIS             6-31+G*
   SOLVE_PEQ         TRUE
   PEQ_SWITCH        0
$end

$peqs
   SOLUTECAVITY RIGIDVDW
$end

$peqs_grid
   DimX 73 -7.50 7.50
   DimY 73 -7.50 7.50
   DimZ 73 -7.50 7.50
$end

The next example illustrates calculation of the solvation free energy of a water molecule at a water/vacuum interface, with the Gibbs dividing surface placed at $z=0.50$ Å. The length of the interface is set to $L_{\mathrm{interface}}=2.75$ Å  and the dielectric is interpolated from bulk solvent to vacuum in the positive $z$ direction across the Gibbs dividing surface. The Cartesian grid, solute cavity, and solvent dielectric are the same as in the previous example.

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          wB97X-V
   BASIS             6-31+G*
   SOLVE_PEQ         TRUE
   PEQ_SWITCH        0
$end

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$peqs
   SOLUTECAVITY RIGIDVDW
   INTERFACE TRUE
   INTERFACELENGTH 2.75
   GIBBSDS 0.50
   INTERFACEDIRECTION POSITIVE
$end

$peqs_grid
   DimX 73 -7.50 7.50
   DimY 73 -7.50 7.50
   DimZ 73 -7.50 7.50
$end

The final example illustrates a nonequilibrium (NonEquilJob = True) solvation calculation for the vertical ionization of ${\mathrm{H}}_2{\mathrm{O}}^{-}$ in bulk water. This is a compound job that first calculates the equilibrium solvation free energy of the anionic state (NonEquilState = Reference), then computes the nonequilibrium energy correction for the ionized state (NonEquilState = Ionized). The two jobs are separated by “@@@”. The Cartesian grid and solvent dielectric are the same as the previous examples, but the solute cavity is chosen to be spherical with a radius of 2.0 Å. The Marcus scheme (NonEquilPartition = Marcus) is used to partition the solvent response into fast and slow components. Since this is the default method, the NonEquilPartition keyword is omitted.

$molecule
   -1 2
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          HF
   BASIS             6-31++G*
   SOLVE_PEQ         TRUE
   PEQ_SWITCH        0
$end

$peqs
   SOLUTECAVITY   SPHERICAL
   SPHERERADIUS   2.00
   NONEQUILJOB    TRUE
   NONEQUILSTATE  REFERENCE
$end

$peqs_grid
   DimX 73 -7.50 7.50
   DimY 73 -7.50 7.50
   DimZ 73 -7.50 7.50
$end

@@@

$molecule
   0 1
   O     0.053004    -0.020947    -0.034784
   H     0.003424     0.185855     0.910594
   H    -0.844842     0.146674    -0.358413
$end

$rem
   SCF_CONVERGENCE   5
   THRESH            14
   EXCHANGE          HF
   BASIS             6-31++G*
   SOLVE_PEQ         TRUE
   PEQ_SWITCH        0
$end

$peqs
   SOLUTECAVITY   SPHERICAL
   SPHERERADIUS   2.00
   NONEQUILJOB    TRUE
   NONEQUILSTATE  IONIZED
$end

$peqs_grid
   DimX 73 -7.50 7.50
   DimY 73 -7.50 7.50
   DimZ 73 -7.50 7.50
$end