Recall that PCM electrostatics calculations require the solution of the set of linear equations given in Eq. (11.2), to determine the vector $\mathbf{q}$ of apparent surface charges. The precise forms of the matrices $\mathbf{K}$ and $\mathbf{R}$ depend upon the particular PCM (Table 11.3), but in any case they have dimension ${N}_{\mathrm{grid}}\times {N}_{\mathrm{grid}}$, where ${N}_{\mathrm{grid}}$ is the number of Lebedev grid points used to discretize the cavity surface. Construction of the matrix ${\mathbf{K}}^{-1}\mathbf{R}$ affords a numerically exact solution to Eq. (11.2), whose cost scales as $\mathcal{O}({N}_{\mathrm{grid}}^{3})$ in CPU time and $\mathcal{O}({N}_{\mathrm{grid}}^{2})$ in memory. This cost is exacerbated by smooth PCMs, which discard fewer interior grid points so that ${N}_{\mathrm{grid}}$ tends to be larger, for a given solute, as compared to traditional discretization schemes. For QM solutes, the cost of inverting $\mathbf{K}$ is usually negligible relative to the cost of the electronic structure calculation, but for the large values of ${N}_{\mathrm{grid}}$ that are encountered in MM/PCM or QM/MM/PCM jobs, the $\mathcal{O}({N}_{\mathrm{grid}}^{3})$ cost of inverting $\mathbf{K}$ is often prohibitively expensive.

To avoid this bottleneck, Lange and Herbert^{Herbert:2016a} have developed
an iterative conjugate gradient (CG) solver for Eq. (11.2) whose cost
scales as $\mathcal{O}({N}_{\mathrm{grid}}^{2})$ in CPU time and $\mathcal{O}({N}_{\mathrm{grid}})$ in memory. A
number of other cost-saving options are available, including efficient
pre-conditioners and matrix factorizations that speed up convergence of the CG
iterations, and a fast multipole algorithm for computing the electrostatic
interactions.^{Li:2009} Together, these features lend themselves to a
solution of Eq. (11.2) whose cost scales as $\mathcal{O}({N}_{\mathrm{grid}})$ in
both memory and CPU time, for sufficiently large systems.^{Herbert:2016a}
Currently, these options are available only for C-PCM, not for SS(V)PE/IEF-PCM.

Listed below are job control variables for the CG solver, which should be
specified within the *$pcm* input section. Researchers who use this feature
are asked to cite the original SWIG PCM
references^{Lange:2010b, Lange:2011b} as well as the reference for the CG
solver.^{Herbert:2016a}

Solver

Specifies the algorithm used to solve the PCM equations.

INPUT SECTION: *$pcm*

TYPE:

STRING

DEFAULT:

INVERSION

OPTIONS:

INVERSION
Direct matrix inversion
CG
Iterative conjugate gradient

RECOMMENDATION:

Matrix inversion is faster for small solutes because it needs to be performed
only once in a single-point calculation. However, the CG solver (which must be
applied at each SCF iteration) is recommended for large MM/PCM or
QM/MM/PCM calculations.

CGThresh

The threshold for convergence of the conjugate gradient solver.

INPUT SECTION: *$pcm*

TYPE:

INTEGER

DEFAULT:

6

OPTIONS:

$n$
Conjugate gradient converges when the maximum residual is less than ${10}^{-n}$.

RECOMMENDATION:

The default typically affords PCM energies on par with the precision of matrix
inversion for small systems. For systems that have difficulty with SCF
convergence, one should increase $n$ or try the matrix inversion solver. For
well-behaved or very large systems, a smaller $n$ might be permissible.

DComp

Controls decomposition of matrices to reduce the matrix norm for the CG Solver.

INPUT SECTION: *$pcm*

TYPE:

INTEGER

DEFAULT:

1

OPTIONS:

0
Turns off matrix decomposition
1
Turns on matrix decomposition
3
Option 1 plus only stores upper half of matrix and enhances gradient evaluation

RECOMMENDATION:

None

PreCond

Controls the use of the pre-conditioner for the CG solver.

INPUT SECTION: *$pcm*

TYPE:

None

DEFAULT:

Off

OPTIONS:

No options. Specify the keyword to enable pre-conditioning.

RECOMMENDATION:

A Jacobi block-diagonal pre-conditioner is applied during the conjugate
gradient algorithm to improve the rate of convergence. This reduces the number
of CG iterations, at the expense of some overhead. Pre-conditioning is
generally recommended for large systems.

NoMatrix

Specifies whether PCM matrices should be explicitly constructed and stored.

INPUT SECTION: *$pcm*

TYPE:

None

DEFAULT:

Off

OPTIONS:

No options. Specify the keyword to avoid explicit construction of PCM matrices.

RECOMMENDATION:

Storing the PCM matrices requires $\mathcal{O}({N}_{\mathrm{grid}}^{2})$ memory. If this is prohibitive,
the NoMatrix option forgoes explicit construction of the PCM matrices, and
instead constructs the matrix elements as needed, reducing the memory requirement to
$\mathcal{O}({N}_{\mathrm{grid}})$ at the expense of additional computation.

UseMultipole

Controls the use of the adaptive fast multipole method in the CG solver.

INPUT SECTION: *$pcm*

TYPE:

None

DEFAULT:

Off

OPTIONS:

No options. Specify the keyword in order to enable the fast multipole method.

RECOMMENDATION:

The fast multipole approach formally reduces the CPU time to $\mathcal{O}({N}_{\mathrm{grid}})$, but
is only beneficial for spatially extended systems with several
thousand cavity grid points. Requires the use of NoMatrix.

MultipoleOrder

Specifies the highest multipole order to use in the FMM.

INPUT SECTION: *$pcm*

TYPE:

INTEGER

DEFAULT:

4

OPTIONS:

$n$
The highest order multipole in the multipole expansion.

RECOMMENDATION:

Increasing the multipole order improves accuracy but also adds more
computational expense. The default yields satisfactory performance in common
QM/MM/PCM applications.

Theta

The multipole acceptance criterion.

INPUT SECTION: *$pcm*

TYPE:

FLOAT

DEFAULT:

0.6

OPTIONS:

$n$
A number between zero and one.

RECOMMENDATION:

The default is recommended for general usage. This variable determines when
the use of a multipole expansion is valid. For a given grid point and box
center in the FMM, a multipole expansion is accepted when $r/d\le $
Theta, where $d$ is the distance from the grid point to the box center
and $r$ is the radius of the box. Setting Theta to one will accept
all multipole expansions, whereas setting it to zero will accept none. If not
accepted, the grid point’s interaction with each point inside the box is
computed explicitly. A low Theta is more accurate but also more
expensive than a higher Theta.

NBox

The FMM boxing threshold.

INPUT SECTION: *$pcm*

TYPE:

INTEGER

DEFAULT:

100

OPTIONS:

$n$
The maximum number of grid points for a leaf box.

RECOMMENDATION:

The default is recommended. This option is for advanced users only.
The adaptive FMM boxing algorithm divides space into smaller and smaller boxes
until each box has less than or equal to NBox grid points. Modification
of the threshold can lead to speedup or slowdown depending on the molecular system and
other FMM variables.

A sample input file for the linear-scaling QM/MM/PCM methodology can be found in the $QC/samples directory, under the name QMMMPCM_crambin.in. This sample involves a QM/MM description of a protein (crambin) in which a single tyrosine side chain is taken to be the QM region. The entire protein is immersed in a dielectric using C-PCM with SWIG discretization.