In response to the Comment by Scalmani and Frisch, we clarify certain claims made in the context of our “switching/Gaussian” discretization procedure. Furthermore, an explanation is proposed to explain observed similarities between this technique and the “continuous surface charge” method introduced by Scalmani and Frisch.

Recently, we described a technique1,2 for discretizing the integral equations that appear in certain “apparent surface charge” dielectric continuum solvation models, commonly known as polarizable continuum models (PCMs).3 Our switching/Gaussian (SWIG) approach1,2 generalizes a technique that was introduced by York and Karplus4 (YK) in the context of the conductor-like screening model (COSMO), in order to ensure that finite-element discretization of the solute/continuum interface does not introduce discontinuities in the solute's potential energy surface. The SWIG procedure extends the YK ideas to a more general class of PCMs that are known in the literature as either the “integral equation formalism” (IEF-PCM),3,5–7 or else the “surface and simulation of volume polarization for electrostatics” [SS(V)PE] method.8–10 

Subsequent to our initial report of the SWIG methodology,1 Scalmani and Frisch (SF) reported their own “continuous surface charge” (CSC) extension of the YK technique.11 In a Comment12 directed at our follow-up paper detailing the SWIG approach,2 SF point out that Ref. 2 also reports a “subSWIG” model that is not a full implementation of the CSC method, and that conclusions based on the subSWIG model are therefore not transferable to the CSC method. Below, we clarify what was claimed in our paper,2 and also comment on similarities and differences between the SWIG and CSC approaches.

To understand this discussion, a bit of context is required. The original YK procedure has two essential features: a switching function that smoothly attenuates the contribution of each surface discretization point,

$\vec{s}_i$
si⁠, as that point passes through a buffer region that surrounds the solute cavity; and Gaussian blurring of the surface charges that are located at the points
$\vec{s}_i$
si
. We have shown that Gaussian blurring is essential when a switching function is employed, because the switching function allows the surface points to approach one another more closely than they otherwise would. This can lead to undesirable (albeit continuously differentiable) oscillations in the energy gradient, which are largely eliminated by Gaussian blurring.1 

To generalize the YK technique to IEF-PCM/SS(V)PE and related methods, one must decide how to discretize the integral operators

$\hat{D}$
D̂ and
$\hat{D}^\dagger$
D̂
, which do not appear in COSMO. (The operator
$\hat{D}^\dagger$
D̂
generates the normal negative of the outward-directed electric field at the cavity surface.2,3,6,8) The off-diagonal matrix elements Dij can be obtained in a straightforward way from the off-diagonal matrix elements of the surface Coulomb operator,
$\hat{S}$
Ŝ
, via the relation3 

\begin{equation}D_{ij} = \vec{n}_j {\cdot }\frac{\partial S_{ij}}{\partial \vec{r}_j} .\end{equation}
Dij=nj·Sijrj.
(1)

Diagonal matrix elements are more problematic. The matrix elements Sii (which represent the Coulomb self-energies of the surface elements) have traditionally been obtained by numerical integration of the Coulomb potential,13–15 although they can be evaluated analytically within the SWIG approach.1,2 In any case, Eq. (1) is not helpful in defining Dii. Often Dii is simply set to zero,15 but the accuracy improves15,16 if Dii is defined by means of an exact sum rule,16 

\begin{equation}D_{ii} = -\frac{1}{a_i}\biggl ( 2\pi + \sum _{j \ne i} D_{ij} \, a_j \biggr ) .\end{equation}
Dii=1ai2π+jiDijaj.
(2)

Here, ai represents the area associated with the ith surface element.

In a preliminary report of our extension of the YK scheme, we set Dii = 0.1 This is justified in the limit ai → 0, but this choice was also pragmatic, as we encountered numerical problems when the definition in Eq. (2) was used instead. Subsequently, we presented a mathematical proof that Eq. (2) does not hold when a switching function is used.2 A variety of numerical problems, including nonvariational solvation energies, are encountered when this sum rule is used to define Dii within SWIG-type approaches.2 Numerical problems are avoided by setting Dii = 0, but this choice degrades the accuracy of the predicted solvation energies. Instead, we chose the definition

\begin{equation}D_{ii} = -\frac{S_{ii}}{2R_I} ,\end{equation}
Dii=Sii2RI,
(3)

which is exact for a spherical cavity of radius RI.5 (We take RI to be the radius of the atomic sphere on which the point

$\vec{s}_i$
si resides.2) Using this definition, the SWIG implementation of IEF-PCM/SS(V)PE appears to be free of numerical artifacts. Numerical tests in Ref. 2 show that solvation energies obtained using SWIG discretization are in good agreement with those obtained using straightforward pointwise discretization of the solute cavity, for which Eq. (2) is valid.

Returning now to the contention at hand, SF claim12 that we have reported an incorrect implementation of their CSC method. What is in fact reported in our paper,2 for comparison to our preferred SWIG approach, is a well-defined, alternative method in which the sum rule in Eq. (2) is substituted in place of the definition in Eq. (3). We therefore chose to refer to this alternative method as the “subSWIG” discretization procedure, and this terminology is used consistently in labeling the data (see Sec. V of Ref. 2). Unfortunately, in a few instances we did make statements that may imply an equivalence between CSC and subSWIG, where in fact none exists. In particular, we concede that the expression given for the CSC value of Dii (Table III in Ref. 2) as well as that for Sii [Eq. (4.5) in Ref. 2] are incorrect and misleading, as the summations over surface grid points in these expressions should have been restricted to a single atomic sphere.

At the same time, this in no way affects the veracity or importance of the subSWIG comparisons summarized above. Although the subSWIG approach was not intended for realistic chemical applications, results obtained using this model are nevertheless important in view of the widespread use of Eq. (2) in the PCM literature.3,10,15,17–19 SubSWIG results confirm the formal proof that the PCM solvation energy is no longer variational when Eq. (2) is used in conjunction with a switching function. Other numerical artifacts are also encountered using subSWIG, including anomalous vibrational frequencies and catastrophic failures to conserve energy in molecular dynamics simulations.2 These artifacts are avoided in our preferred SWIG approach,2 and also, it appears, in the CSC approach.12 

In the CSC discretization method,11 the matrix elements Dii are defined by considering the atomic spheres individually, independent of how they intersect to form the cavity surface. The sum rule in Eq. (2) is applied to each individual atomic sphere (in the context of the Born ion model), using exact expressions for Dij and Sij that are the same as those that would be used in our SWIG procedure, if the cavity were spherical. Thus, the CSC method first disassembles the cavity surface into atomic spheres in order to define matrix elements, then reassembles it to perform the PCM calculation. (This is the “separation between model and cavity” emphasized by SF.11,12) Numerical results,12 using several of the same numerical tests that we used in Ref. 2, appear to be in reasonable agreement with SWIG results. This makes some sense in view of the fact that the definition in Eq. (3) is obtained from a model of a strictly spherical cavity surface.5 

In summary, the subSWIG model studied in Ref. 2 is not equivalent to the CSC discretization procedure, yet the former plays an important role in the development of smooth PCMs, as it demonstrates the incompatibility of switching functions with the widely used sum rule in Eq. (2). Our SWIG discretization method2 avoids the numerical problems associated with this sum rule. Calculations performed thus far suggest that the SWIG approach, while formally distinct from the CSC method,11 is numerically quite similar, and an explanation for this observation is suggested herein.

1.
A. W.
Lange
and
J. M.
Herbert
,
J. Phys. Chem. Lett.
1
,
556
(
2010
).
2.
A. W.
Lange
and
J. M.
Herbert
,
J. Chem. Phys.
133
,
244111
(
2010
).
3.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
,
Chem. Rev.
105
,
2999
(
2005
).
4.
D. M.
York
and
M.
Karplus
,
J. Phys. Chem. A
103
,
11060
(
1999
).
5.
B.
Mennucci
,
E.
Cancès
, and
J.
Tomasi
,
J. Phys. Chem. B
101
,
10506
(
1997
).
6.
J.
Tomasi
,
B.
Mennucci
, and
E.
Cancès
,
J. Mol. Struct.: THEOCHEM
464
,
211
(
1999
).
7.
E.
Cancès
and
B.
Mennucci
,
J. Chem. Phys.
114
,
4744
(
2001
).
8.
D. M.
Chipman
,
J. Chem. Phys.
112
,
5558
(
2000
).
9.
D. M.
Chipman
,
Theor. Chem. Acc.
107
,
80
(
2002
).
10.
D. M.
Chipman
and
M.
Dupuis
,
Theor. Chem. Acc.
107
,
90
(
2002
).
11.
G.
Scalmani
and
M. J.
Frisch
,
J. Chem. Phys.
132
,
114110
(
2010
).
12.
G.
Scalmani
and
M. J.
Frisch
,
J. Chem. Phys.
134
,
117101
(
2011
).
13.
A.
Klamt
and
G.
Schüürmann
,
J. Chem. Soc., Perkin Trans.
2
,
799
(
1993
).
14.
D. M.
Chipman
,
J. Chem. Phys.
110
,
8012
(
1999
).
15.
M.
Cossi
,
G.
Scalmani
,
N.
Rega
, and
V.
Barone
,
J. Chem. Phys.
117
,
43
(
2002
).
16.
E. O.
Purisima
and
S. H.
Nilar
,
J. Comput. Chem.
16
,
681
(
1995
).
17.
M.
Cossi
,
N.
Rega
,
G.
Scalmani
, and
V.
Barone
,
J. Chem. Phys.
114
,
5691
(
2001
).
18.
G.
Scalmani
,
V.
Barone
,
K. N.
Kudin
,
C. S.
Pomelli
,
G. E.
Scusceria
, and
M. J.
Frisch
,
Theor. Chem. Acc.
111
,
90
(
2004
).
19.
C. S.
Pomelli
, in
Continuum Solvation Models in Chemical Physics
, edited by
B.
Mennucci
and
R.
Cammi
(
Wiley
,
Chichester, UK
,
2007
), pp.
49
63
.