In a recent paper,1 Lange and Herbert discuss some formal aspects of the theory of the polarizable continuum models (PCM) of solvation2 and they describe the surface integral discretization scheme (SWIG/ISWIG) employed in their implementation of the model. They compare such discretization scheme with a few others previously published in the literature.3–5 In particular, they implemented elements of the scheme described by us within the context of the continuous surface charge (CSC) formalism5,6 for PCM in order to assess its quality and performance against the SWIG scheme. They refer to this modified CSC scheme as “subSWIG” and, at times, as “CSC/subSWIG.” However, the implementation of the element of our approach, designated subSWIG in the paper by Lange and Herbert, as described in Ref. 1, is not an actual implementation of our CSC approach, and the conclusions based on analysis of this alternative model are not transferable to the CSC approach. Specifically, the discretization scheme used within our CSC formalism smoothly dissociates bonds, does not produce spurious features in simulated vibrational spectra, and conserves energy during molecular dynamic (MD) simulations.

The PCM solvation models require the evaluation of surface integrals6,7 over the solute–solvent interface Γ, usually assumed to be the boundary of a molecular cavity dug within a polarizable dielectric or conductor medium. Such integrals are evaluated using numerical techniques involving the discretization of the surface in many small elements, each representing a portion of the surface, where the integral kernels are evaluated. The technical details of such discretization schemes are typically of no interest for the average user of the solvation models, but they are of critical importance to ensure a smooth dependence of the free energy functional upon the solute geometry and other parameters in the model. Lange and Herbert's discretization scheme was first introduced in a previous letter8 by the same authors and, just like the scheme we adopted in the CSC formalism,5 is based on the work of York and Karplus,9 published over ten years ago, which represents the true breakthrough toward the definition and actual implementation of a smooth free energy functional.

In Ref. 1 Lange and Herbert did not implement correctly the content of Sec. III D of our CSC paper,5 where we describe the optimal choice of the so-called self-field factorsgi. This is one of the six quantities associated with each surface element arising from the discretization of the solute–solvent interface Γ. The other five quantities associated with each surface element are the position

$\mathbf {s}_{i}$
si⁠, the weight ai, the outward normal direction to the surface
$\hat{\mathbf {n}}_{i}$
n̂i
, the self-potential factorfi, and the exponent ζi of the Gaussian basis function used in the expansion of the apparent surface charge (ASC) density. Both self-factors are defined on the unit sphere as soon as the appropriate Lebedev grid has been selected to fulfill user-defined criteria (typically a given number of points per atomic sphere or per Å2 of surface area) and before the spheres are intersected to form the surface of the molecular cavity. In addition to the dependence on the particular Lebedev grid being used, the self-factors depend only on the radius R of the sphere according to the simple scaling relations fi(R) = Rfi(1) and gi(R) = gi(1). The self-potential factor are defined by evaluating the exact Born energy using the discretized surface of the unit sphere, while the self-field factors are defined using a sum rule, discussed by Purisima et al.,10 which descends directly from Gauss's law and involves the integral of the normal electric field on a spherical surface. More than once in Sec. III D of Ref. 5 we state that the results being developed pertain to the specific case of a single sphere. In the last two paragraph of Sec. III D we explain how we assume that the values of the self-factors are transferable from the single sphere to the final cavity. Therefore, for each sphere the values of fi and gi are set once the integration grid has been chosen, are scaled to the actual value of the radius of the atomic sphere, and are used for the final cavity discretization without recomputing them. In the spirit of the “separation between model and cavity,” which is a core concept in the CSC formalism, the self-factors are intrinsic properties of the surface elements that are defined within some discretization scheme, and are used to describe the singular portion of the discretized
${\hat{\mathcal S}}$
Ŝ
,
${\hat{\mathcal D}}$
D̂
, and
$\hat{\mathcal {D}}^{\star }$
D̂
integral operators, i.e., the
$\mathbf {S}$
S
,
$\mathbf {D}$
D
, and
$\mathbf {D}^{\star }$
D
matrices. In Ref. 1, it appears that the sums in Eqs. (4.2) and (4.5) are extended to all surface elements of the final molecular surface, after the intersection of the atomic spheres has been carried out, and that is not consistent with what we described in Ref. 5. Therefore, also the “CSC” column of Table III is incorrect. The definition of Dii in Eq. (4.2) has been previously used by us in the context of an older discretization scheme,11 which was not based on the work of York and Karplus,9 i.e., employed point-charges to represent the ASC density. That scheme has been superseded as described in Refs. 5 and 6 exactly because of its many limitations, some of them correctly pointed out by Lange and Herbert in Sec. IV B of Ref. 1.

In conclusion we would like to show the actual results obtained with the CSC approach, which is implemented in the GAUSSIAN 09 program.12 Most likely very similar results would be obtained by a modified subSWIG approach using a definition of Dii consistent with the CSC one. Figure 1 shows the behavior of the relevant quantities as the bond of Na–Cl in water is broken. To obtain these data we used exactly the same setup used by Lange and Herbert and therefore Fig. 1 can be directly compared with Figs. 1 and 2 of Ref. 1. No discontinuities are observed in the energy gradient, and all the eigenvalues of the

$\mathbf {Q}$
Q matrix have the correct sign. In Fig. S1 of the supplementary material13 the simulated IR spectrum of liquid glycerol is plotted, as obtained by us using analytical energy second derivatives, and the reader should compare it with Fig. 7 of Ref. 1 and note the absence of any spurious features in the high wavenumber region. Finally, Figs. 2 and S2 show results from a MD simulation of glycine in water that we carried out under very similar conditions to the ones reported in Sec. VI B 1 of Ref. 1. In particular, we used the same level of theory, solvation model, temperature, time step, and overall duration. Despite the limited number of time steps, our simulation does not present any apparent energy conservation issue, as shown in the top panel of Fig. 2, which should be compared with Fig. 9 of Ref. 1. Moreover, the bottom panel of Fig. 2 shows the polarization energy profile for the overall simulation and an expanded view of the region where the transition between neutral and zwitterionic form occurs, and again there is no apparent issue with our method.

FIG. 1.

Behavior of the CSC discretization scheme on the dissociation of NaCl in water. The explanation of the different plots can be found in the captions to Figs. 1 and 2 of Ref. 1. The SWIG discretization scheme behaves almost exactly such as the CSC scheme of Ref. 5.

FIG. 1.

Behavior of the CSC discretization scheme on the dissociation of NaCl in water. The explanation of the different plots can be found in the captions to Figs. 1 and 2 of Ref. 1. The SWIG discretization scheme behaves almost exactly such as the CSC scheme of Ref. 5.

Close modal
FIG. 2.

Top panel: Energy conservation profile for ab initio MD simulation of glycine in water under the same conditions used in Ref. 1. Bottom panel: Partial profile of the energy of polarization during the same MD simulation showing the region where the transition between the neutral and the zwitterionic form occurs.

FIG. 2.

Top panel: Energy conservation profile for ab initio MD simulation of glycine in water under the same conditions used in Ref. 1. Bottom panel: Partial profile of the energy of polarization during the same MD simulation showing the region where the transition between the neutral and the zwitterionic form occurs.

Close modal

Therefore (i) Lange and Herbert's implementation of an element of our CSC approach in the subSWIG model is not complete and, as demonstrated here, produces results inconsistent with our CSC approach, (ii) all the implied criticism of our approach made on the basis of the performance of their “subSWIG model” is without merit and, (iii) their SWIG/ISWIG does not appear to provide results which are significantly different from the ones obtained with our CSC approach.

1.
A. W.
Lange
and
J. M.
Herbert
,
J. Chem. Phys.
133
,
244111
(
2010
).
2.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
,
Chem. Rev.
105
,
2999
(
2005
).
3.
H.
Li
and
J.
Jensen
,
J. Comput. Chem.
25
,
1449
(
2004
).
4.
P.
Su
and
H.
Li
,
J. Chem. Phys.
130
,
074109
(
2009
).
5.
G.
Scalmani
and
M. J.
Frisch
,
J. Chem. Phys.
132
,
114110
(
2010
).
6.
F.
Lipparini
,
G.
Scalmani
,
B.
Mennucci
,
E.
Cances
,
M.
Caricato
, and
M. J.
Frisch
,
J. Chem. Phys.
133
014106
(
2010
).
7.
E.
Cances
, “
Integral equation approaches for continuum models
,” in
Continuum Solvation Models in Chemical Physics
, edited by
B.
Mennucci
and
R.
Cammi
(
Wiley
,
New York
,
2007
), Chap. 1.2, pp.
29
48
.
8.
A. W.
Lange
and
J. M.
Herbert
,
J. Phys. Chem. Lett.
1
,
556
(
2010
).
9.
D.
York
and
M.
Karplus
,
J. Phys. Chem. A
103
,
11060
(
1999
).
10.
E.
Purisima
and
S.
Nilar
,
J. Comput. Chem.
16
,
681
(
1995
).
11.
M.
Cossi
,
G.
Scalmani
,
N.
Rega
, and
V.
Barone
,
J. Chem. Phys.
117
,
43
(
2002
).
12.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
 et al, GAUSSIAN 09, Revision B.01, Gaussian, Inc., Wallingford, CT,
2009
.
13.
See supplementary material at http://dx.doi.org/10.1063/1.3567489 for the additional figures discussed in the text.

Supplementary Material