Each of the implicit solvation models described above is designed for
isotropic, bulk solvation—conditions that can be qualitatively described
using a scalar dielectric constant, $\epsilon $. For an anisotropic
environment, such as a liquid/vapor interface, a more general approach is
to solve Poisson’s equation with a spatially-varying dielectric
*function*, $\epsilon (\mathbf{r})$. Such an approach can be used, for
example, to model the air/water interface by describing certain regions
of space using $\epsilon =1$ (air) and other regions using $\epsilon =78$
(water).^{Coons:2016, Coons:2018} The atomistic region is described at the
SCF level. Construction of the dielectric function $\epsilon (\mathbf{r})$
is based on a solute cavity, somewhat analogous to the solute cavity in a PCM
calculation but where the dielectric changes smoothly (rather than abruptly)
from $\epsilon =1$ inside the cavity (in the atomistic QM region) to
$\epsilon ={\epsilon}_{\mathrm{solvent}}$ outside, in the continuum solvent.
Options for modeling a liquid/vapor interface are also available; see
Section 11.2.10.3. In contrast to PCMs, which must
approximate the “volume polarization” due to penetration of the tails of the
QM charge density beyond the solute cavity, this effect is described exactly
under Poisson boundary conditions. The price to be paid for lifting this
approximation is that charge densities must be discretized and integrated over
three-dimensional space rather than simply a two-dimensional solute cavity
surface.

Q-Chem’s implementation of Poisson boundary conditions was introduced in Ref. Coons:2016 and refined in Ref. Coons:2018. Users of this methodology are asked to cite at least Ref. Coons:2018.

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 (\mathbf{r})\widehat{\mathbf{\nabla}}{\phi}_{\mathrm{tot}}(\mathbf{r})\right]=-4\pi {\rho}_{\mathrm{sol}}(\mathbf{r}).$$ | (11.13) |

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

$${\rho}_{\mathrm{sol}}(\mathbf{r})={\rho}_{\mathrm{nuc}}(\mathbf{r})+{\rho}_{\mathrm{elec}}(\mathbf{r}).$$ | (11.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. (11.14), plus the polarization potential induced in the medium:

$${\phi}_{\mathrm{tot}}(\mathbf{r})={\phi}_{\mathrm{elec}}(\mathbf{r})+{\phi}_{\mathrm{nuc}}(\mathbf{r})+{\phi}_{\mathrm{pol}}(\mathbf{r}).$$ | (11.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}}({\mathbf{r}}_{i})=\sum _{\mu \nu}^{{N}_{\mathrm{basis}}}{P}_{\mu \nu}\u27e8{g}_{\mu}(\mathbf{r})\left|\frac{1}{|\mathbf{r}-{\mathbf{r}}_{i}|}\right|{g}_{\nu}(\mathbf{r})\u27e9.$$ | (11.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 $\mathbf{r}$. The electronic part of the solute’s charge density is then obtained from

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

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

$${\rho}_{\mathrm{nuc}}(\mathbf{r})=-\sum _{A}^{{N}_{\mathrm{atoms}}}{Z}_{A}\delta \left(\mathbf{r}-{\mathbf{R}}_{A}\right),$$ | (11.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 ${\mathbf{R}}_{A}$, where $dV$ is the volume of a Cartesian voxel. The exact nuclear electrostatic potential in vacuum is computed according to

$${\phi}_{\mathrm{nuc}}(\mathbf{r})=-\frac{1}{4\pi}\sum _{A}^{{N}_{\mathrm{atoms}}}\frac{{Z}_{A}}{|\mathbf{r}-{\mathbf{R}}_{A}|}.$$ | (11.19) |

The Dirac delta function in Eq. (11.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. (11.13) for the case where $\epsilon (\mathbf{r})=1$. The induced polarization charge density and corresponding electrostatic potential are then obtained iteratively, following closely the procedure outlined in Refs. Fisicaro:2016 and Andreussi:2012. Equation (11.13) is rewritten to appear as Poisson’s equation in vacuum with an additional source charge density:

$\begin{array}{cc}\hfill {\widehat{\nabla}}^{2}{\phi}_{\mathrm{tot}}(\mathbf{r})& =-4\pi \left[{\displaystyle \frac{{\rho}_{\mathrm{sol}}(\mathbf{r})}{\epsilon (\mathbf{r})}}+{\displaystyle \frac{\widehat{\mathbf{\nabla}}\mathrm{ln}\epsilon (\mathbf{r})\mathbf{\cdot}\widehat{\mathbf{\nabla}}{\phi}_{\mathrm{tot}}(\mathbf{r})}{4\pi}}\right]\hfill \\ & =-4\pi \left[{\rho}_{\mathrm{sol}}(\mathbf{r})+{\rho}_{\mathrm{pol}}(\mathbf{r})\right],\hfill \end{array}$ | (11.20) |

where the polarization charge density is

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

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

$${\rho}_{\mathrm{iter}}^{(k+1)}(\mathbf{r})=\frac{\widehat{\mathbf{\nabla}}\mathrm{ln}\epsilon (\mathbf{r})\mathbf{\cdot}\widehat{\mathbf{\nabla}}{\phi}_{\mathrm{tot}}^{(k)}(\mathbf{r})}{4\pi}$$ | (11.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. (11.21) vanishes wherever $\epsilon (\mathbf{r})=1$. The quantity ${\rho}_{\mathrm{iter}}(\mathbf{r})$ is iterated to convergence by performing updates ${\rho}_{\mathrm{iter}}^{(k)}(\mathbf{r})\to {\rho}_{\mathrm{iter}}^{(k+1)}(\mathbf{r})$ using Eq. (11.22), taking ${\phi}_{\mathrm{tot}}^{(0)}(\mathbf{r})$ to be the solute’s electrostatic potential in vacuum. Each time ${\rho}_{\mathrm{iter}}(\mathbf{r})$ is updated, a new polarization charge density is generated and Eq. (11.20) is solved to obtain a new total electrostatic potential. To stabilize the iterative updates of ${\rho}_{\mathrm{iter}}(\mathbf{r})$, we use a damping procedure

$\begin{array}{cc}\hfill {\rho}_{\mathrm{iter}}^{(k+1)}(\mathbf{r})& ={\displaystyle \frac{\eta}{4\pi}}\left[\widehat{\mathbf{\nabla}}\mathrm{ln}\epsilon (\mathbf{r})\mathbf{\cdot}\widehat{\mathbf{\nabla}}{\phi}_{\mathrm{tot}}^{(k)}(\mathbf{r})\right]+(1-\eta ){\rho}_{\mathrm{iter}}^{(k)}(\mathbf{r})\hfill \\ & =\eta {\rho}_{\mathrm{iter}}^{(k+1)}+(1-\eta ){\rho}_{\mathrm{iter}}^{(k)}(\mathbf{r})\hfill \end{array}$ | (11.23) |

rather than using Eq. (11.22) as written. Following Refs. Fisicaro:2016 and Andreussi:2012, 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 \mathit{d}\text{\mathbf{r}}{\phi}_{\mathrm{pol}}(\text{\mathbf{r}}){\rho}_{\mathrm{sol}}(\mathbf{r}),$$ | (11.24) |

is computed, where

$${\phi}_{\mathrm{pol}}(\text{\mathbf{r}})={\phi}_{\mathrm{tot}}(\text{\mathbf{r}})-{\phi}_{\mathrm{sol}}(\mathbf{r}).$$ | (11.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 \mathit{d}\text{\mathbf{r}}{\phi}_{\mathrm{pol}}(\text{\mathbf{r}}){g}_{\mu}(\text{\mathbf{r}}){g}_{\nu}(\text{\mathbf{r}}).$$ | (11.26) |

For additional details and a full description of the algorithm, see Ref. Coons:2018.

Nonequilibrium solvation within the framework of PCMs is discussed in
Section 11.2.2.3, and that
methodology^{Mewes:2015a, You:2015} has been adapted for vertical ionization
processes in conjunction with Poisson boundary
conditions.^{Coons:2016, Coons:2018} A state-specific approach has been
implemented that requires solution of a nonlinear Schrödinger equation,
Eq. (11.6). For the reference state, corresponding to $i=0$ in
Eq. (11.6), this is equivalent to solving Eq. (11.13)
for the induced polarization potential and computing the solvation free energy
and Fock matrix correction according to
Eqs. (11.24) and (11.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,^{You:2015}
Eq. (11.20) becomes^{Coons:2018}

$$-\frac{1}{4\pi}{\widehat{\nabla}}^{2}{\phi}_{\mathrm{tot},i}(\mathbf{r})={\rho}_{\mathrm{sol},i}(\mathbf{r})+{\rho}_{\mathrm{pol},0}^{\mathrm{slow}}(\mathbf{r})+{\rho}_{\mathrm{pol},i}^{\mathrm{fast}}(\mathbf{r}),$$ | (11.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}}(\mathbf{r})={\rho}_{\mathrm{iter},i}(\mathbf{r})+\left(\frac{1-{\epsilon}_{\mathrm{opt}}(\mathbf{r})}{{\epsilon}_{\mathrm{opt}}(\mathbf{r})}\right)\left[{\rho}_{\mathrm{sol},i}(\mathbf{r})+{\rho}_{\mathrm{pol},0}^{\mathrm{slow}}(\mathbf{r})\right]$$ | (11.28) |

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

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

Subscripts $i$ and 0 refer to the ionized state and the reference state, respectively. The quantity ${\epsilon}_{\mathrm{opt}}$ in Eq. (11.28) is the solvent’s optical dielectric constant, and ${\rho}_{\mathrm{iter},i}$ is computed using Eq. (11.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. (11.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.

SOLVE_PEQ

Perform a solvation free energy calculation on a Cartesian grid using Poisson
equation boundary conditions.

TYPE:

STRING

DEFAULT:

False

OPTIONS:

TRUE/FALSE

RECOMMENDATION:

None.

PEQ_SWITCH

Inclusion of solvent effects begins when the SCF error falls below 10${}^{-\mathrm{PEQ}\mathrm{\_}\mathrm{SWITCH}}$.

TYPE:

INTEGER

DEFAULT:

3

OPTIONS:

$n$
Corresponding to 10${}^{-n}$

RECOMMENDATION:

Use the default unless solvent effects need to be incorporated earlier in the SCF procedure.

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. (11.20).^{Coons:2018} 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}_{\mathrm{max}}-{x}_{\mathrm{min}}$
and the grid spacing is $\mathrm{\Delta}x={L}_{x}/({N}_{x}-1)$. The volume element discussed
in the context of Eq. (11.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.

PEQMaxIter

Sets the maximum number of iterations used in the conjugate gradient solver.

INPUT SECTION: *$peqs*

TYPE:

INTEGER

DEFAULT:

500

OPTIONS:

User-defined.

RECOMMENDATION:

Use the default unless the calculation fails to converge.

PEQSolverThresh

The electrostatic potential is considered converged when the error falls below
10${}^{-\mathrm{PEQSolverThresh}}$.

INPUT SECTION: *$peqs*

TYPE:

INTEGER

DEFAULT:

5

OPTIONS:

$n$
Corresponding to 10${}^{-n}$

RECOMMENDATION:

Use the default unless a tighter convergence criterion is desired, at greater
computational cost.

PolarIterScale

Specifies the mixing parameter $\eta $ that is used in
Eq. (11.23) to stabilize iterative
solution for the polarization charge density.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

0.6

OPTIONS:

$\eta $
Desired value of the mixing parameter (unit-less).

RECOMMENDATION:

Use the default, which was tested in Refs. Fisicaro:2016
and Andreussi:2012,
unless the calculation proves difficult to converge.

BatchSize

Evaluate electrostatic potential integrals in batches over small parts of the
Cartesian grid, rather than computing all integrals in a single batch.

INPUT SECTION: *$peqs*

TYPE:

INTEGER

DEFAULT:

5000

OPTIONS:

$n$
Corresponding to a batch size of $n$ grid points.

RECOMMENDATION:

For large grids ($\ge {10}^{6}$ grid points), a batch size $n\approx {N}_{x}{N}_{y}{N}_{z}/5$ is recommended.

SolventDielectric

Sets the dielectric value for the solvent.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

78.39

OPTIONS:

$\epsilon $
Desired value of $\epsilon $ (dimensionless).

RECOMMENDATION:

The default value corresponds to water at 25${}^{\circ}$C.

OpticalDielectric

Sets the optical dielectric value of the solvent.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

1.7778

OPTIONS:

${\epsilon}_{\mathrm{opt}}$
Desired value of ${\epsilon}_{\mathrm{opt}}$ (dimensionless).

RECOMMENDATION:

The optical dielectric is equal to the square of the index of refraction. The
default value corresponds to water at 25${}^{\circ}$C.

SoluteCavity

Specifies the type of solute cavity, which determines the form of the
dielectric function $\epsilon (\mathbf{r})$.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

Vacuum

OPTIONS:

Vacuum
Perform calculation in vacuum, $\epsilon (\mathbf{r})\equiv 1$.
RigidVDW
Create a cavity using spherically symmetric error functions.
Spherical
Create a single spherical cavity around the solute.

RECOMMENDATION:

None.

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 11.2.3.^{Herbert:2016a}
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 11.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. Fisicaro:2016. The
dielectric function, which depends parametrically on the nuclear positions
${\mathbf{R}}_{A}$, is

$$\epsilon (\mathbf{r};\{{\mathbf{R}}_{A}\})=\left({\epsilon}_{\mathrm{solvent}}-1\right)\left\{\prod _{A}^{{N}_{\mathrm{atoms}}}h({d}_{A},\mathrm{\Delta};|\mathbf{r}-{\mathbf{R}}_{A}|)\right\},$$ | (11.30) |

where

$$h({d}_{A},\mathrm{\Delta};|\mathbf{r}-{\mathbf{R}}_{A}|)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{|\mathbf{r}-{\mathbf{R}}_{A}|-{d}_{A}}{\mathrm{\Delta}}\right)\right].$$ | (11.31) |

Equation (11.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. (11.30) and (11.31) sets radius of the atomic sphere for atom $A$. By default, ${d}_{A}=1.2{r}_{\mathrm{vdW}}$ but this can be controlled as described below.

RigidScale

Sets the length scale on which the error function employed in the rigid vdW
cavity construction interpolates the dielectric value from vacuum to solvent.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

0.265 Å

OPTIONS:

$\mathrm{\Delta}$
Specifies the desired value (in Å); see Eq. (11.31).

RECOMMENDATION:

Use the default value, which was tuned in Ref. Fisicaro:2016 so
that errors in small-molecule
solvation energies were $\approx 1$ kcal/mol.

VDWType

Specifies details for the rigid vdW cavity construction.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

Scaled

OPTIONS:

Unscaled
${d}_{A}={r}_{\mathrm{vdW}}$
Scaled
${d}_{A}={r}_{\mathrm{vdW}}\times scale$
Shifted
${d}_{A}=({r}_{\mathrm{vdW}}+shift)\times scale$

RECOMMENDATION:

None. The values of $scale$ and $shift$ are set with the VDWScale and
VDWShift keywords.

VDWScale

Sets the empirical scale factor applied to the atomic van der Waals radius.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

1.2

OPTIONS:

$scale$
Specifies the desired dimensionless scaling factor for the atomic vdW radii.

RECOMMENDATION:

Use the default value.

VDWShift

Adjusts the center of the spherically-symmetric error functions when
constructing the rigid vdW cavity.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

0.0 Ångstroms

OPTIONS:

$shift$
Specifies the desired shift (in Å).

RECOMMENDATION:

None. If VDWType is set to Shifted, the vdW scale factor is
set to 1.0 by default. This can be adjusted using
the VDWScale keyword.

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:^{Coons:2016}

$$\epsilon (r)=\frac{1}{2}\left\{({\epsilon}_{\mathrm{solvent}}+1)+({\epsilon}_{\mathrm{solvent}}-1)\text{tanh}\left[\alpha (r-{r}_{\mathrm{mid}})\right]\right\}.$$ | (11.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.

SphereRadius

Sets the radius of the spherical solute cavity.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

No default.

OPTIONS:

$R$
Desired spherical cavity radius (in Å).

RECOMMENDATION:

None. See the Supporting Information of Ref. Coons:2016 for more information.

InterpolLength

Sets the length scale on which the dielectric is smoothly interpolated from
vacuum to solvent.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

No default.

OPTIONS:

$L$
Desired interpolation length (in Å).

RECOMMENDATION:

None.

InterpolScale

For a given interpolation length ($L$), InterpolScale ($\alpha $) sets
the sharpness of the dielectric transition from vacuum to solvent.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

$\alpha $ = $4/L$

OPTIONS:

$\alpha $
Desired interpolation scale factor (in Å${}^{-1}$).

RECOMMENDATION:

Use the default unless a broader (smaller $\alpha $) or narrower (larger
$\alpha $) transition region is desired.

Interface

Perform a solvation calculation at a solvent/vacuum interface.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

False

OPTIONS:

True
Modify the dielectric function to simulate an interface between solvent and vacuum.
False
Perform a solvation calculation in bulk solvent.

RECOMMENDATION:

The user will also need to specify the length scale on which the dielectric is
smoothly interpolated from the bulk solvent value to vacuum, and the location
of the Gibbs dividing surface (see below).

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:^{Coons:2016, Coons:2018}

$$\epsilon (z)=\frac{1}{2}\left\{({\epsilon}_{\mathrm{solvent}}+1)+(1-{\epsilon}_{\mathrm{solvent}})\text{tanh}\left[\beta (z-{z}_{\mathrm{GDS}})\right]\right\}.$$ | (11.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.

InterfaceLength

Sets the length scale over which the dielectric function is smoothly
transitioned from bulk solvent to vacuum in the $z$ direction

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

None

OPTIONS:

${L}_{\mathrm{interface}}$
Desired interface length (in Å).

RECOMMENDATION:

This sets the value $\beta =4/{L}_{\mathrm{interface}}$ in Eq. (11.33).
See the Supporting Information of Ref. Coons:2016 for a full
description of the solvent/vacuum interface construction.

GibbsDS

Sets the location of the Gibbs dividing surface, in the $z$ direction.

INPUT SECTION: *$peqs*

TYPE:

FLOAT

DEFAULT:

None

OPTIONS:

${z}_{\mathrm{GDS}}$
Desired location of the Gibbs dividing surface (in Å).

RECOMMENDATION:

Consult the literature. One such way to determine this value is to compute a
density profile of the solvent as a function of $z$ and set this location to be
where the solvent density has decreased to 50% of the bulk value. Usually
${z}_{\mathrm{GDS}}\approx {L}_{\mathrm{interface}}/2$.

NonequilJob

Obtain the nonequilibrium free energy of solvation for a vertical ionization process.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

False

OPTIONS:

True
Compute the nonequilibrium free energy for a vertical ionization process.
False
Compute the equilibrium solvation free energy.

RECOMMENDATION:

None.

NonequilPartition

Specifies the manner in which the solvent response is partitioned into fast
and slow components.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

Marcus

OPTIONS:

Marcus
Employ the Marcus partitioning scheme.
Pekar
Employ the Pekar partitioning scheme.

RECOMMENDATION:

Use the default. Although the fast and slow solvation responses are different
between the two approaches, the total solvation free energy is the
same,You:2015 but the Pekar scheme is computationally more
expensive than the Marcus scheme.

NonequilState

Specifies the state of interest for a nonequilibrium vertical ionization.

INPUT SECTION: *$peqs*

TYPE:

STRING

DEFAULT:

Reference

OPTIONS:

Reference
The reference (initial) state, from which an electron will be removed.
Ionized
The final (ionized) state.

RECOMMENDATION:

None. Both values will be needed in a compound input job, in order to compute
the non-equilibrium response to vertical ionization;
see Example 11.2.10.4.

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