Cavity-Free Continuum Solvation: Implementation and Parametrization in a Multiwavelet Framework

We present a multiwavelet-based implementation of a quantum/classical polarizable continuum model. The solvent model uses a diffuse solute–solvent boundary and a position-dependent permittivity, lifting the sharp-boundary assumption underlying many existing continuum solvation models. We are able to include both surface and volume polarization effects in the quantum/classical coupling, with guaranteed precision, due to the adaptive refinement strategies of our multiwavelet implementation. The model can account for complex solvent environments and does not need a posteriori corrections for volume polarization effects. We validate our results against a sharp-boundary continuum model and find a very good correlation of the polarization energies computed for the Minnesota solvation database.


INTRODUCTION
Continuum solvation models have been used in quantum chemistry for half a century. 1−4 Their use is motivated by the need to simulate the effect of a large solvent environment on a molecular solute, keeping at the same time the computational cost to a minimum.
Several models and flavors have throughout the years been developed. Common to essentially all such models are two basic assumptions: 1) the solvent degrees of freedom can be conveniently described in terms of a continuum, parametrized using macroscopic properties of the solvent; 2) the quantum system is confined inside a cavity and the solute−solvent interaction is described in terms of functions (charge density/ potential) supported on the cavity surface. Whereas the former assumption is a physical one, giving a prescription for the underlying physical laws, 5 the latter is a convenient mathematical formulation, which reduces the computational cost transforming a three-dimensional problem in the whole space to a two-dimensional one on the boundary of the molecular cavity. Despite the convenience, a sharp boundary between neighboring molecules assumes that no electronic density is present beyond the cavity surface. This is not physically sound, because electronic densities of solute and solvent in reality overlap. Initially, this issue has been dealt with by simple renormalization procedures: 3 more elaborate corrections have later been proposed, 6−8 and for the Integral Equation Formalism (IEF) formulation of the polarizable continuum model (PCM) it can be shown that a first-order correction is already included in the model. 9 A full account of this issue is, however, not practical in terms of a surface model, and the ever increasing basis sets employed in routine calculations, including very diffuse functions, aggravate the problem further by allowing more and more of the electron density to "escape" the cavity.
Neglecting electronic charge overlap between solute and solvent does not only impact the electrostatic energy: excitation energies depend on the charge distribution in the excited states, which is invariably more diffuse than in the ground state, and other interaction terms, such as the repulsion energy, depend explicitly on the overlap between solute and solvent densities. 10 The parametric description of the cavity surface also presents challenges, not only from a formal point of view to define the correct cavity boundary, 2,3 but also from a technical standpoint, especially for larger molecules. The development of stable cavity generators is still an active area of research. 11 −20 In recent years, several real-space methods for quantum chemistry have been developed, 21−26 and with these, the treatment of solvation as a three-dimensional problem has become a feasible alternative. The advantage is a seamless integration with the quantum mechanical implementation: the electrostatic potential is no longer computed in a vacuum but in the generalized dielectric medium with a position-dependent permittivity. Several real-space codes have so far adopted this strategy. 27−32 Another advantage of this approach is an increased flexibility: no constraints are placed on the form of the permittivity function, and complex environments consisting of surfaces, droplets, membranes, can be treated without the need of ad-hoc implementations, which are often limited to a handful of special cases. 33−35 In this contribution, we will present our implementation, which makes use of a multiwavelet (MW) framework 36−39 to solve both the Kohn−Sham (KS) equations of density functional theory (DFT) 40,41 and the Generalized Poisson Equation (GPE) 27 for the solvent reaction potential. We will also show a set of benchmark calculations to showcase the implementation's theoretical correctness, parametrization, and flexibility. MWs constitute a basis that can give accurate results up to a user-defined precision, thanks to an automatic adaptive refinement. 39 Our implementation is included in the opensource MW computational chemistry software package MRChem. 26 The combination of MW-based KS-DFT and GPE solver provides a methodology for the assessment of solvent effects with controlled precision.

THEORY
In the theoretical framework adopted in this work, molecules are described through quantum mechanics, whereas the solvent is modeled as a classical entity, described by macroscopic properties. The two subsystems are connected by the solute− solvent interaction, which describes the mutual polarization of the two subsystems. 2,3 Such an interaction is described by classical electrostatics. In almost all implementations, the quantum and the classical problem are solved with very different methods: the most widely used approach makes use of Boundary Element Method (BEM) 42 techniques to solve the electrostatic problem (environment) and Gaussian Type Orbital (GTO) bases 43,44 to describe the quantum problem. The use of Multiwavelets offers a unique opportunity to treat both problems with the same tools and methods. We will here recap the basic concepts of multiresolution analysis (MRA) and how it is employed to solve the electrostatic and the quantum problem.

Multiresolution Analysis and Multiwavelets.
MRA is a mathematical framework that considers a space spanned by a basis of functions with self-similarity and regularity properties. 45 In practice, all basis functions are constructed by simple translation and dilation of a small set of starting functions ϕ(x): The core idea of MRA is that the space spanned by the basis functions at a given scale n is a subspace of those at scale n + 1. Such a ladder of spaces can be extended indefinitely, and its limit is by construction dense in L 2 . Successive refinements thus provide a systematic strategy to reach completeness, with a handful of predefined functions. This is in stark contrast with traditional GTO methods, where extending a basis requires a complete reparametrization of the basis set, atom by atom. The wavelet functions are obtained by taking the difference between two consecutive scaling spaces, and they convey information about the error incurred at each scale n due to neglecting the refinement at scale n + 1, see Figure 1 for a 1-dimensional illustration.
As long as the fundamental properties of self-similarity and completeness are preserved, the choice of a specific basis set can be guided by numerical considerations to obtain compact representation of functions and efficient application of operators.
Alpert's Multiwavelets 36 constitute a practical realization of MRA by considering a set of polynomial functions (e.g., Legendre or Interpolating polynomials) defined on an interval. The main advantages of Multiwavelets are the simplicity of the original basis (a polynomial set) and the disjoint support (basis functions are zero outside their support node). 37 The latter enables adaptive refinement of functions to minimize the storage needs and the computational overhead. The extension to three-dimensional functions is obtained by tensor-product methods, and operators are efficiently applied in a separated form. 38 Multiwavelets are an ideal framework to deal with integral operators, and this allows both the KS equations for the quantum system 40 and the Poisson equation for the solvent polarization 27 to be solved within the same formalism, once the equations are converted from the conventional differential form to the appropriate integral form. Functions are projected/ computed on an adaptive grid to guarantee the requested precision. All operations (operator applications, algebraic manipulations) are defined within the requested precision, in such a way that the developer can easily implement new algorithms with little effort 46 and the end-user only needs to specify the requested precision. 26,47−49 For details about how to solve the KS equations within a MW framework, we refer to the literature. 39  GPE, we will expose the derivation and the implementation details in the remainder of this section.

Electrostatics of Continuous Media.
Any material is a bound aggregate of nuclei and electrons: at microscopic level these charged particles obey the microscopic Maxwell equations. We are, however, interested in the macroscopic behavior of the material in the presence of external sources of charge ρ(r) and current j(r). Following Jackson, 5 we can perform a spatial average to arrive at the macroscopic Maxwell equations: These equations are expressed in terms of the usual electric and magnetic fields, E and B, and additionally the displacement D and magnetization H fields appear as a result of the spatial averaging. In the quasistatic limit, the electric field has zero curl and can thus be written in terms of a scalar potential function: where V is the electrostatic potential. To relate the external sources to the potential it is first necessary to relate the fields E and D with a constitutive relation, 5,50 which is, in general, a nonlinear and space-time nonlocal relationship between the fields. For linear and local continuous media the constitutive relation is where the permittivity ε(r) is a position-dependent, rank-3 symmetric tensor. Upon inserting the constitutive relation into the first of Maxwell's equations, we obtain the GPE: In the following, we will further specialize to the isotropic case ε(r) = ε(r)I, with I the rank-3 identity: We remark that the permittivity is still position-dependent, in contrast to the usual PCM treatment. The solution to eq 5 can be partitioned as where V ρ is the electrostatic potential in a vacuum and V R is the reaction potential. The polarization energy is then defined as We write the reaction potential as a functional of the charge density: the functional dependence is linear. 9 2.3. The Quantum-Classical Coupling. Our quantum mechanical treatment of the system will be based on KS-DFT. For an N-electron system coupled with a classical polarizable continuum environment, the KS-DFT f ree energy 2 functional 51,52 reads:   (8) The molecular charge density is separated into electronic and nuclear components: E xc [ρ e , ∇ρ e ] is a GGA exchange-correlation functional, and the nuclear-electron potential is defined as ζ is a scalar factor influencing the portion of exact exchange included in the energy. The 1-body reduced density matrix (RDM) and electronic density function appear in the energy expression: The minimum is found by constrained optimization, to enforce idempotency and normalization of the RDM:  (12) and leads to the variational condition: 51,53 where the effective one-electron Fock operator appears:

Solving the Generalized Poisson Equation.
The solution to the GPE is a function supported on the entire space Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article 3 . Apparent surface charge formulations of continuum solvation models do not solve eq 5 directly, but rather reformulate it as a boundary integral equation and solve it by boundary-element discretization. The apparent surface charge, supported on the closed solute−solvent boundary, is the sought-after quantity to compute the polarization energy. 9 Such a procedure is generally based on two underlying assumptions: (1) the charge density is entirely contained inside the cavity boundary, and (2) the permittivity is unitary inside the cavity and constant outside the cavity, with a jump condition that defines the electrostatic potential and field across the cavity boundary. With a real-space approach both assumptions can be relaxed and the equation can be solved directly. We recap here the procedure outlined by Fosso-Tande and Harrison. 27 We rewrite eq 5 in terms of the Laplacian of the potential V: The second term on the right-hand side contains both the gradient of the permittivity and the gradient of the potential. When the permittivity is not constant, the equation cannot be solved in one step by inversion of the Laplacian, i.e., by convolution of the right-hand side with the Laplacian's Green's function. An iterative strategy must be employed instead. Let us then define the effective charge: = eff (16) and the polarization function: such that eq 15 becomes We can now formally solve eq 15 in terms of the Laplacian's Green's function: However, both the polarization energy in eq 7 and the solute− solvent interaction term in the Fock operator are expressed in terms of the reaction potential, rather than the total electrostatic potential. By making use of the partition of V in eq 6 and recalling that which can be formally inverted using the Poisson kernel: We stress that γ is a function of V = V ρ + V R and eq 21 must therefore be solved iteratively.

IMPLEMENTATION
In this section we present details about our specific choice of parametrization for the permittivity and how we compute the electrostatic potential between solute and solvent. We also show how we couple this to a standard self-consistent field (SCF) optimization procedure.

The Permittivity Function Parametrization.
We partition space into two regions: a cavity containing the solute, and the remainder. The cavity surface is defined as the union set of a collection of interlocking spheres centered on the nuclei. Their radii are parametrized by using the corresponding van der Waals radii times a factor. This factor is often set to either 1.1 or 1.2, 2 but it might vary, e.g., depending on the charge of the solute. For standard continuum models the cavity boundary is the support of the electrostatic problem for the solute−solvent interaction. In the current model it serves as a support to define the parametrization of the positiondependent ε(r). In Section 4 the appropriate parametrization of the cavity for the present model will be discussed.
Following Fosso-Tande and Harrison, we write the permittivity as a function of the molecular cavity function: 27 The exponential parametrization proves convenient in light of the definition of γ in eq 17, which lets us define its gradient using the cavity function, C(r), only. The molecular cavity function is constructed as follows. For each sphere α centered at r α with radius R α , we can measure the signed normal distance of any point in space as Given s α (r), we define a smoothed boundary of the sphere as where σ is a user-defined smoothing parameter: C α approaches the Heaviside step function as σ → 0. The molecular cavity function is then a product of all N spheres: see Figure 2 for an example. The log-derivative of the permittivity in eq 17 is then: requiring evaluation of the gradient of the cavity function. For interlocking-spheres cavities, a closed-form analytical expression is available, see Appendix A, and is implemented in our code. Note, however, that, in a real-space, multiwavelet framework, we can compute this gradient by direct application of the derivative operator, 54 which allows one to use more complex or even numerical definitions of the boundary, e.g., as isodensity surfaces.
3.2. The Self-Consistent Reaction Field. The selfconsistent reaction field (SCRF) is the iterative procedure to solve the GPE for any given molecular density. At convergence, the iterations produce the reaction potential V R , which can be directly employed in the solution of the KS-DFT equations.
Algorithm 1 shows the iterative procedure implemented to solve the GPE within the SCF iterations. The input parameters at iteration n are the charge density ρ [n] , the permittivity ε(r), a A straightforward implementation of the microiterations suffers from slow convergence of the reaction potential, thus adding a significant prefactor to each SCF iteration. We use the Krylov-accelerated inexact Newton (KAIN) method, 41 which is a convergence acceleration technique, similar to Pulay's DIIS 55 and Anderson's mixing. 56 At each microiteration i, the updated reaction potential [ + ] V i R 1 is constructed as a linear combination, with constraints, of N previous iterates. The KAIN history length N impacts both convergence and memory: N = 5 is generally a good compromise between fast convergence (fewer iterations) and acceptable memory footprint. The KAIN acceleration is combined with an adaptive threshold to improve the convergence rate of the microiterations: instead of converging the reaction potential to the same predefined threshold ϵ used for the orbitals, we make use of a threshold, δ, chosen to be the norm of the orbital update in the parent SCF macroiteration. δ is thus updated during the SCF procedure. There are two parameters that affect the convergence pattern of the reaction potential, V R : 1. The guess for V R at the start of the microiterations: ,0 R 1 (and zero for the first microiteration embedded in the first macroiteration).
2. The convergence threshold for the microiterations: (C) fixed threshold δ, or (D) dynamic threshold δ [n] = |Δρ [n] |. These lead to four possible convergence regimes: AC, BC, AD, BD; the latter being our default. Figure 3 illustrates how the number of microiterations evolves. A dynamic precision threshold D reduces the number of microiterations in the beginning of the SCF procedure, simply because the threshold for convergence is looser. Using the converged V R from the previous macroiteration B helps close to SCF convergence, because the orbitals do not change much and the starting guess for the microiterations is also better. Combining those two choices results in the optimal convergence pattern: the convergence threshold is progressively tighter, while at the same time the starting guess for the reaction potential improves. The opposite choice (AC instead of BD) requires a large number of microiterations throughout, whereas the intermediate choices (AD and BC) result in a large number of iterations at the beginning (BC) or at the end (AD). We underline that all four choices converge to the same result for the example in Figure 3, but we can envisage cases  initial guess (C) or guess from previous macroiteration (D). A dynamic threshold (green and red curves) reduces the number of microiterations at the beginning of the SCF procedure. A starting guess from the previous SCF macroiteration (green and blue curves) is effective close to convergence. Combining the two (green curve) is the optimal strategy. The dip observed for the blue and green curves at macroiteration 1 is due to the fact that the macroiteration 0 is a preliminary step and the orbital are not changed progressing from macroiteration 0 to macroiteration 1, but the convergence threshold is tightened. This results in an almost converged reaction potential as a starting guess for the microiterations nested in macroiteration 1.
where convergence could potentially be prevented by choices A and C.

RESULTS
For all systems, the solvation energies have been computed with both Gaussian16 57  Two sets of calculations have been performed. The aim of the first set was to determine a good parametrization for the cavity surface in terms of the atomic radii and the cavity surface thickness. Once a satisfactory parametrization was achieved, an extensive benchmark of solvation energies was performed, by considering the Minnesota Solvent Descriptor Database (MSDD) of Marenich et al. 59 All calculations reported are KS-DFT using the PBE0 functional. 60 Gaussian16 results are obtained with the Def2-TZVP, 61−63 basis set, except where otherwise stated. MRChem results are obtained setting the global precision parameter to 10 −5 . In other words, the obtained absolute energy is correct with at least five digits with respect to the complete basis set (CBS) limit. 47 This is not to be confused with the convergence threshold of a SCF calculation performed with an atomic basis set, which will guarantee the "exact" result within the chosen basis, but where the precision compared to the CBS is limited by the choice of basis.

Cavity Parametrization.
For the parametrization calculations, 4 molecules of different levels of polarity were chosen: water, ethanol, formaldehyde, and ethyne (geometries taken from the MSDD, 59 file names 0217wat, 0045eth, 0069met, and 0030eth). No geometry optimization was performed. They were chosen to give a minimal set of neutral (polar and apolar) systems, to allow for a reliable yet simple data set to identify a good choice of the parameters defining the cavity.
In Gaussian16, the external iteration procedure 64,65 was used to extract the reaction energy from the total energy. a The spheres used for the cavities were atom-centered and used the atoms' Bondi radii 66 scaled by a factor of 1.1, as is standard for Gaussian16. Three different permittivities have been employed: 2.0, 4.0, and 80.0.
In MRChem, the cavity is also built from atom-centered spheres, with each radius R i parametrized as where R i vdW is the Bondi radius 66,67 of the i-th atom, σ i is the width of the cavity boundary, and α i and β i are adjustable parameters. We allowed for granular, sphere-by-sphere flexibility in our implementation of the cavity function. By default, one value is used for each parameter (α, β, σ) for all spheres. The combination α = 1.1 and β = 0.0 would yield matching radii between MRChem and Gaussian16. In the following, we explored results when α values were 1.0, 1.1, 1. The aim of the parametrization is to see how the cavity width σ affects the results of our calculations, compared to a sharp-boundary method, and to choose the combination of α and β coefficients that provides a good correlation between our method and a sharp boundary implementation. The goal is not to replicate results from Gaussian16 implementation: our Outlier species are marked in blue and red when containing bromine and iodine, respectively. The labels refer to A. 5-bromouracil, H 3 C 4 N 2 O 2 Br (n203); B. 5-bromo-3-s-butyl-6-methyl-uracil, H 13 C 9 N 2 O 2 Br (test1013); C. 2-bromoanisole, H 7 C 7 OBr (test5008); D. Bromobenzene, H 5 C 6 Br (0186bro); E. 4-bromopyridine, H 4 C 5 NBr (0573bro); F. 1-bromo-2-chloroethane, H 4 C 2 ClBr (0202bro); G. 5-iodouracil, H 3 C 4 N 2 O 2 I (test2018); H. iodomethane, H 3 CI (test4003); I. iodobenzene, H 5 C 6 I (test4001); J. 1,4-dichlorobenzene, H 4 C 6 Cl 2 (0176pdi); K. 3-bromoanisole, H 7 C 7 OBr (test5009). method has a diffuse cavity layer, whereas the cavity of IEFPCM is a 2-dimensional boundary. This will lead to contributions and errors that are not equivalent. Figure 4 shows the results for the cavity parametrization for α = 1.1 and α = 1.2. Results for α = 1.0 and α = 1.3 are not shown, because they largely overestimate (α = 1.0) or underestimate (α = 1.3) solvation energies, but they are available in the data package available online on DataVerse. 68 We conclude that a cavity parametrization with α = 1.1 and β = 0.5 provides a good correlation with sharp-boundary IEFPCM for all reasonable values of the permittivity and default value of cavity width. This choice of α and β with σ = 0.2 au is the current default in MRChem. For low permittivity (ε = 2.0), Figure 5.1 shows that for neutral species our data is quite close to the main diagonal for small energies, but has a slight systematic deviation for more negative reaction energies (bottom left corner). For ions, Figure 6.1 shows a systematic overestimation with respect to Gaussian16, and a clear distinction between cations and anions. For ε = 4.0 ( Figure 5.2 and 6.2), we see a similar trend, although most data points appear to be closer to the diagonal. For high permittivity ε = 80, Figure 5.3 for neutral species and Figure 6.3 for ionic ones show that the values are now mostly below the diagonal, that is, solvation energies are underestimated compared to Gaussian16. In Figure 5, we can see a set of outlying point with respect to the rest of the data. These points have been identified as species containing bromine b or iodine, c with only one outlier containing chlorine instead. d59 There may be multiple, concomitant reasons for these discrepancies: (a) Bromine and iodine are the only atoms from the fourth and fifth period of the periodic table present in the set; (b) the radii used in the definition of the cavities for these elements might not be appropriate; (c) the different treatment of volume polarization in the two implementations (full account in our model and implicit first-order correction in

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article the IEFPCM model 8,9 ) might affect the description of these molecules, where a more delocalized electronic density is expected. It would be interesting to disentangle the effects of surface and volume polarization, but it is not straightforward to do so and it goes beyond the scope of the present work. The fact that molecules are quite close to the line, especially as the reaction energy becomes small (top right corner), is not surprising. We chose the cavity parameters from a limited set of small molecules. On the other hand, the observed deviations for larger solvation energy are to a large degree systematic, which shows that they could be accounted for, with a more refined parametrization.
Cations tend to have less diffuse density than anions. Therefore, the size of the cavity with respect to the spatial extent of the electronic density is larger for cations than for anions. According to the simple Born model, solvation energy of ions is inversely proportional to the radius of the cavity, which explains the better correlation observed for cations: when the charge distribution is better confined inside the cavity, the difference between a sharp interface formally not accounting for volume polarization and a diffuse one including it, becomes smaller.

Performance.
The current code is a prototype, and we have therefore not yet dedicated attention to improving its performance in terms of computational time and memory footprint. A few general considerations can however be made. The solution of the GPE is technically similar to that of the Helmholtz equation, which we employ to solve the SCF equations. 40,69 It should therefore be possible to achieve linear scaling with respect to the system size once the code is fully optimized. 70 This is a feature of MRA, 37 which is designed to decouple the long-and short-range interactions automatically thanks to the adaptive refinement scheme coupled with the use of the nonstandard form of operators. 38 In this sense, the algorithm should be competitive with implementations of sharp-cavity models that employ the fast multipole method (FMM) to accelerate the matrix-free solution of the PCM equations. 71 A qualitative comparison with the domain decomposition (DD) family of algorithms 9,72,73 is also in order. DD approaches to implicit solvation are, by construction, linear scaling. Furthermore, they are easily recast in a matrix-free form that both reduces the memory footprint and lends itself to further performance boosting via the FMM. 74 However, in our understanding of the algorithm, these advantages of the method are not straightforwardly extended to cavities with diffuse boundaries. Furthermore, when dealing with quantum mechanical source densities, the quantum-classical coupling must rely on volume integrations, e.g., using a DFT grid, to correctly represent the escaped charge. 75 Our algorithm achieves formal simplicity and, in principle, algorithmic efficiency. Real-space methods for the reaction potential can be coupled with GTO methods for the electronic-structure problem, 76 thus making our method of interest beyond multiwavelet-based quantum chemistry. Currently the main bottleneck is constituted by the memory footprint of the functions describing the cavity and the solvent reaction potential, since they extend throughout the whole computational domain. Work is currently in progress to deal with such functions in an efficient way.

CONCLUSIONS
We have implemented, parametrized, and benchmarked a continuum solvation model based on a position dependent permittivity ε(r). 27 Our algorithm performs microiterations, nested within each SCF cycle, to obtain the solvent reaction potential. We overcome convergence issues using KAIN convergence acceleration and an adaptive convergence threshold. Our implementation is robust and introduces only a modest computational overhead.
With a simple parametrization, we have obtained a good correlation with respect to the IEFPCM implemented in Gaussian16, for an extensive library of geometries and a wide range of permittivities. Some systematic deviations have been observed, suggesting that a more sophisticated cavity parametrization could yield even better agreement. An alternative option, which is often challenging for standard solvation models, is to parametrize the permittivity by making use of an isodensity cavity as support. This choice would forego the radius parametrization altogether, but it might pose other challenges, because the cavity gradient must be computed numerically, and the coupling with the density functional must be taken into account.
The performance and stability might be further improved, by considering a different approach to the SCRF microiterations: a square-root parametrization of the electrostatic potential, as suggested by Fisicaro et al., might prove useful. 30 The flexibility of the method will allow for several additional developments, such as the inclusion of charged particles outside the cavity, as well as other contributions to the solvation energy, such as cavitation, dispersion, and repulsion.  (28) which only requires to compute the analytical gradient of the interlocking sphere cavity function C(r). The analytical gradient of the interlocking sphere cavity is as defined by Fosso-Tande and Harrison: 27 The gradient of a single sphere cavity function C α is i k j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z In the implementation, we use a cutoff of 10 −12 for the denominator, in order to avoid numerical discontinuities.

■ ASSOCIATED CONTENT Data Availability Statement
Input and output files for the Gaussian16 and MRChem calculations reported in this work are available on the Norwegian instance of the Dataverse data repository: 10. 18710/TFSWLC. The data package also includes the Jupyter notebooks used to produce the graphs in this paper.