Field-Aware Interfaces in Continuum Solvation - The Journal of

Subscriber access provided by UNIV OF LOUISIANA

Article

Field-Aware Interfaces in Continuum Solvation Matthew Truscott, and Oliviero Andreussi J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b01363 • Publication Date (Web): 03 Apr 2019 Downloaded from http://pubs.acs.org on April 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Field-Aware Interfaces in Continuum Solvation Matthew Truscott and Oliviero Andreussi∗ Department of Physics, University of North Texas, Denton, TX 76207, USA E-mail: [email protected] Abstract Continuum models of solvation are widespread tools for the prediction of solvation free energies of small molecular compounds from first principles. However, the continuum approximation at the core of these approaches limits their accuracy for the modelling of the aqueous solvation of compounds with highly polarized residue or of charged species. This is due to the fact that straightforward definitions of the continuum interface do not account for the reorganization effects induced by these solutes on the positions of surrounding solvent molecules. This kind of problem is usually overcome by stretching the definition of the continuum model, i.e. by using chemical intuition to adjust (usually shrinking) the size of the solvation interface close to the polarized/charged atoms. Nonetheless, this strategy introduces a significant number of additional parameters that need to be tuned and, at the same time, deter the model’s transferability. A transferable solution is instead represented by an improved definition of the continuum interface, able to automatically account for the polarization/charge state of the embedded system. Following recent approaches in the literature, the component of the solute’s electric field normal to the interface can be used as an effective proxy for the net charge of the embedded system. Here we show a simple definition of this field-aware approach as applied to the recently proposed soft-sphere continuum solvation (SSCS) method. In this model, each soft sphere composing the interface is allowed to readjust as a function of the value of the field flux through its surface. This

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

effect introduces a complex dependence of the interface function on both the electronic and the ionic degrees of freedom of the solute. To account for this dependence during optimization procedures (e.g. the SCF loop and geometry optimization algorithms), the analytic derivatives of the new interface are reported and validated with their numerical counterparts. Application of the field-aware procedure to molecular compounds showing pathological behaviours with the standard SSCS approach show that significant improvements can be achieved by specifically tuning the newly introduced parameters.

1

Introduction

In recent years, continuum solvation models 1–3 have been a major point of interest in the field of computational condensed matter. 4 Via these methods, it is possible to capture the complexity of a solvent environment overcoming the need for an explicit description. Such a task can be unnecessarily expensive and simply unfeasible when running exploratory studies on large sets of molecules for desirable properties, and when performing investigations into interfacial processes that require appreciable simulation length or time scales. The fundamental idea of continuum models is instead to separate a solvated system into an atomistic quantum-mechanical description of the solute, coupled with a statistical averaging of the embedding solvent. Comprehensive work has been devoted towards finding a definition of the interface between these two regions that is both physically sound and chemically accurate. 1,2,5?

? –13

The problem of separating the continuum and quantum-mechanical regions of

space is recast as the expression of an interface function in terms of a minimal set of parameters, whilst maintaining the ability to correctly describe the solvent spatial arrangement in terms of just the solute’s degrees of freedom. While most state-of-the-art continuum solvation approaches rely on sharp two-dimensional interfaces, in condensed-matter frameworks a continuous smoothly-defined three-dimensional interface function is often preferred. 7,8,14–19 Local definitions in terms of the solute’s degrees of freedom, either ionic positions or electronic density, have consistently shown good accuracy for small and compact systems, e.g. 2


Page 2 of 41

Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in reproducing solvation free energies of small organic molecules 8,9,20,21 or interfacial energies of close-packed metals. 12,22?

,23

More recently, non-local definitions of the interface

function have been proposed, 10,11 aiming to overcome spurious continuum regions in larger, more complex, open systems, such as molecular aggregates, biomolecules, and open crystal structures. A challenging aspect of the definition of the solute/solvent boundary is represented by the fact that the interface function ought to depend on the local charge state of the solute and, for neutral solutes, on the specific charge polarization of a residue or atom. This dependence is required in order to reproduce the fact that charged solutes (or polar residues) will interact more strongly with the molecules of a polar solvent, thus shrinking their solvation boundaries. The occurrence of a stronger interaction can be seen as a break-down of the continuum approximation, i.e. the assumption that solvent degrees of freedom can be treated as a homogeneous distribution around the solute. As a result of this problem, minimally parametrized continuum models were shown to provide less reliable results for molecules containing acidic and basic functional groups solvated in water, with mean absolute errors much larger than less polar solutes and above chemical accuracy. 8 Moreover, solvation of charged molecules was shown to require a separate parametrization of the interface function, especially for negatively charged solutes. 21 One sensible path that overcomes this problem is in the use of hybrid continuum/explicit approaches, in which part of the solvation shells of the solute is treated with full atomistic details. 24 Nonetheless, several efforts have been devoted in the literature to propose an implicit solution to this problem that would minimize or remove the need for an explicit treatment of the solvent. Two possible strategies have been followed in the literature to address this problem. In the first strategy, the interface function is adjusted ad hoc, either using chemical intuition or as a dependence on the charge state of the solute. For example, the atomic radius of an element entering the definition of the interface is adjusted when the atom is involved in a

3



Page 4 of 41

specific functional group or to improve the overall agreement with experimental results. 25 Alternatively, a retuning of the model is performed separately for anions and cations. 9,21 As a second strategy, the continuum interface is untouched, but chemical intuition is used to perform a different tuning of solute-solvent interaction terms (e.g. non-electrostatic interactions) according to the type of functional group or residue to which the atom belongs. When combined, these strategies have provided continuum models with impressive accuracies on small organic molecules and ions. However, these approaches rely on chemical intuition, which can be less relevant in non-organic solutes, and on a significant increase of tunable parameters, which can lead to a poor transferability of the models. In the ideal scenario, a continuum approach that exploits the knowledge of the charge/polarization state of the embedded quantum-mechanical system would make the use of chemical intuition or the need for independent parameterizations redundant. Recently, extensions of continuum approaches have been proposed that rely on the electrostatic field at the continuum interface as a proxy for the local charge of the embedded solute. Pomogaeva and Chipman proposed a field-extremum model, 26,27 that calculates the electric field normal to the interface boundary, and selects the minimum and maximum values as potential hydrogen bond locations. These field values are then used to adjust the free energy by means of a post processing term with four tunable parameters. The results showed an average mean unsigned error of 0.9kcal/mol for a set of 264 neutral molecules, and 2.4kcal/mol for ions. Similarly, Rahimi et al. introduced a function 28 that serves the same purpose, acting as a correction that varies depending on the electric field normal to the continuum boundary. This model is also able to capture the steric asymmetry of water and consequently similar asymmetry effects in the majority of solvents. Both these approaches use the electrostatic field through the interface to modulate the strength of the interaction between solute and solvent. Alternatively, Sundararaman and Goddard proposed a more challenging approach, 13 where the continuum interface itself is allowed to adjust as a function of the normal electric field. The reported results show impressive accuracy on both neutral and charged compounds. However, the adopted definition

4


Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


introduces a subtle dependence of the continuum interface on both the electronic and ionic degrees of freedom of the solute, which makes the optimization of the electronic density and of the atomic structure more challenging and was not addressed explicitly in Ref. 13 Following these recent developments, we introduce here a field-aware approach for continuum interfaces, complementing its formulation with the detailed description of all its derived properties. Although the proposed formalism can be generally applied to different continuum interfaces, here we focus on the soft-sphere continuum solvation (SSCS) model of Fisicaro et al., 9 featuring an interface function defined from a set of smooth interlocking spherical functions. While the original radii of the spheres were defined as a uniform scaling of van der Waals radii, 29 here an additional sphere-specific scaling factor is introduced, proportional to the magnitude of the field flux passing through each sphere. Since charged solutes typically contain atoms with larger flux densities, the proposed model should primarily affect them over neutral solutes. The extra scaling function captures the asymmetry of charge seen in many investigations and, like the model of Sundararaman and Goddard, 13 automatically adjusts the continuum boundary to fit charged species as appropriate. The article is organized as follows: in Section 2 we review the main definitions of the soft-sphere interface and of the continuum embedding effects defined on it; in Section 3 we introduce the field-aware extension of the soft-sphere model, including the derivation of the derivatives of the field-aware interface with respect to the solute’s electronic density and ionic positions; in Section 4 we present a sensible definition of the field-aware function and of the parameters involved in it; following a short description of the numerical details involved in the simulations, in Section 5 we present results of the approach on specific neutral and charged molecules, as a function of the tunable parameters involved in the field-aware definition; eventually in the conclusions we review the performances of the approach and outline sensible parametrization strategies.

5



2

Page 6 of 41

Methods

2.1

Soft-Sphere Interfaces

Many wide-spread continuum solvation models rely on interface functions that are generated as the union of atom-centered spheres. The polarizable continuum model of Tomasi et al., 1,2 as well as its analogous formulations 30–32 and following extensions, 33,34 are based on sharp interfaces. These interfaces are generated from interlocking atomic spheres whose radii are given by the van der Waals radii scaled by a single homogeneous scaling factor, tuned to reproduce the experimental values of solvation free energies. The two-dimensional interface is then discretized into geometrical elements (tesserae), which provide the base for the application of boundary element methods as well as for the definition of non-electrostatic intreractions 3,35 and more advanced implicit statistical approaches. 36 Defining the interface on ionic positions has the key advantage that for a given electronic optimization procedure (i.e. the SCF loop of Hartree-Fock and Density Functional Theory calculations) the numerical domain and its associated numerical problems do not change. A secondary advantage, as discussed above, is that specific radii can be modified ad hoc to improve the accuracy of the results. Following the same strategy of PCM, but in the context of condensed-matter simulation packages, Fisicaro et al. introduced an interface function defined in terms of the product of atom-centered spherical functions. In contrast with PCM and its sharp two-dimensional boundary, a smooth definition of the interface function was proposed, allowing to exploit for the continuum embedding the same structured grid used for the electronic structure calculation. In particular, given the set of ionic positions {Ra }, the SSCS interface function s({ξ} , {R} ; r) is expressed as

s({ξ} , {Ra } ; r) = 1 −

Y

h0 ({ξa } ; |r − Ra |),

a

6


(1)

Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where {ξ} is the tunable set of parameters for the SSCS function:

{ξ} =

 vdW    R sphere radii     αξ       ∆ξ

homogeneous scaling factor

(2)

softness

Here, αξ and ∆ξ are global tunable parameters and RavdW is the van der Waals radius of atomic species a. For following equations, for brevity, we will use ra = |r − Ra | and ra = r − Ra . The spherical functions h0a (r) are chosen to be scaled error functions centered on each of the ions, 1 r − αξ RavdW h ({ξa } ; r) = 1 + erf . 2 ∆ξ 0

(3)

The error function is an ideal choice of scaling function here due to its well behaved continuous derivatives. By the above definition, the spherical function returns zero when within the solute and one when outside, with a continuous transition between the regimes that is approximately 4∆ wide. This function is typically defined for each atom, as a function of r, the relative position of a point in space with respect to the atom of interest. Similarly to standard PCM implementations, the Universal Force Field parameters are exploited for the atomic radii. 29 In SSCS a different choice was made for the radius of Nitrogen, 9 which was chosen to be the Bondi radius 37 to improve the agreement with experimental results. We would like to stress here that this kind of adjustment goes in the same direction as the similar procedures followed in the quantum-chemistry literature. However, as nitrogen atoms tend to have a similar positive polarization in organic molecules, the adoption of a field-aware approach may circumvent the need for such a redefinition, and one could therefore simply adopt a full internally-consistent set of radii. Following the strategy of the self-consistent continuum solvation (SCCS) model of Andreussi et al., 8 a multitude of embedding effects can be defined on top of the soft-sphere interface function s(r), including both electrostatic (dielectric, electrolyte) and non-electrostatic 7



Page 8 of 41

(cavitation, dispersion, repulsion, etc.) interactions between the solute and its environment. In particular, the incorporation of this continuum regime into density functional theory involves a modification of the total energy functional to include an additional energy of interaction with the continuum that can be recast into a functional of the continuum interface

E tot ρel , {Ra } = E system ρel , {Ra } + E interf ace s ρel , {Ra }

(4)

where ρel is the electronic density of the system. As reviewed in previous works, 4,10 this functional can be expressed as a sum of the electrostatic generalized Poisson equation term, and terms that reflect the geometrical properties (i.e. the quantum surface S and quantum volume V ) of the solvated system,

E

interf ace

Z [s(r)] = −

1 ([s(r)] − 1)|∇φ(r)|2 dr + αS + βV 8π

(5)

where the dielectric function (r) is defined in SSCS in terms of the interface function as

(r) = 0 + (1 − 0 )s(r).

(6)

and the quantum surface and quantum volume are the following functionals of the interface function 9,38 Z |∇s(r)|dr

S [s (r)] =

(7)

Z V [s (r)] =

s(r)dr.

(8)

In order to optimize the electronic and ionic degrees of freedom of the solute, e.g. through the SCF cycle or a geometry optimization algorithm, we are interested in extremizing the energy functional with respect to the electronic density and the atomic positions. The system contribution to the energy can be computed with standard approaches in vacuum. However, the value and the derivatives of the interface term require ad-hoc procedures.

8


Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It is possible to recast the calculation of the interface contributions to the Kohn-Sham potential and inter-atomic forces into a unified framework, by expressing these in terms of the functional derivative of the interface energy with respect to the interface function δE interf ace [s] 1 d δS δV (r) = − (r)|∇φ(r)|2 + α (r) + β (r) δs 8π ds δs δs 1 d ∇s(r) =− (r)|∇φ(r)|2 − α∇ · +β 8π ds |∇s(r)|

(9)

In particular, one can derive the interface contribution to the Kohn-Sham potential and inter-atomic forces as

δE interf ace [s] (r) δρel Z interf ace [s] 0 0 δs 0 δE = (r, r ) (r )dr . el δρ δs

interf ace VKS (r) =

fainterf ace = −∇Ra E interf ace [s] Z δE interf ace [s] (r)dr = − ∇Ra s(r) δs

(10)

(11)

where the ∇Ra notation is adopted to indicate the partial derivative with respect to the position of atom a. The derivatives of the interface function with respect to the electronic density (δs/δρ) and the ionic positions (∇Ra s) are analytically accessible from the constituent definition of the interface. In particular, if the former value is zero (i.e. no dependence of the interface on the electronic density, as in SSCS), there is no interface contribution to the Kohn-Sham potential.

2.2

Field-Aware Interfaces

As discussed in the introduction, standard interface functions show significantly worse performance when applied to specific functional groups, as well as to charged systems. This 9



is due to the inability of these interface functions to account for the solute’s influence on the structure of the environment. For instance, water molecules will cluster more closely to solute molecules that possess some charge imbalance, which in turn results in an inaccurate positioning of the interface function, which ought to reflect this effect by enveloping the solute molecules more tightly. One can account for this effect by using system-specific strategies, for example by using a different set of tuned parameters for anions. 9,21 The mean absolute error in aqueous solvation free energy, as compared to experimental values, is a popular test to determine how well various competing interface models perform. The SCCS model performs well for neutral solutes, with an error close to chemical accuracy of 1 kcal/mol. 8 However further analysis shows that this error is highly dependent on the functional groups possessed by certain molecules, as expected from the above discussion. In particular, mean errors for alkenes, ketones, and aldehydes are shown to be well below 1 kcal/mol in many cases, 8 whereas amines, ethers and acids show errors of up to 4-5 kcal/mol. When applied to cations and anions, this same model has errors of around 2 kcal/mol for cations and 17 kcal/mol for anions, an order of magnitude higher. In contrast, the SSCS model performs comparably for neutral molecules, and again worse for cations (6 kcal/mol) and anions (10 kcal/mol). To compensate for these higher errors, one can re-tune the models specifically for cations and anions, these parameters are then fed into the simulation according to the a priory knowledge of the charge state of the solute. The effect of this is to shrink the interface, thus reducing the volume of the quantum mechanical region. This captures the physical effect of the solvent interacting more closely to the solute molecule. However, it has to be noted that such a strategy may have limited application for simulations involving different species, in which only the global charge state of the system is assigned a priori, while the specific charge state of the each constituents is unknown and, possibly, depends on the embedding environment. For this reason, a new definition of the interface function, able to automatically shrink

10


Page 10 of 41

Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the presence of more polar or charged systems, is required. The use of the electrostatic field at the interface as a local proxy for the charge and hydrogen-bond propensity of the embedded solute was proposed by Chipmann and collaborators 26,27 and, later, by Bardhan and coworkers. 28,39,40 The same strategy was extended to the definition of the interface function in the context of joint density functional theory by Sundararaman and Goddard, 13 but the detailed derivation of the derivatives connected to the modified interface was not reported. Starting from these developments, here we present the complete derivation of a fieldaware extension of the SSCS model. Following a similar approach to the one exploited to define the quantum-surface of system, the component of the vacuum electric field normal to the interface can be expressed as

E n (r) = E(r) · ∇s(r)

(12)

where the electric field in vacuum E(r) can be computed from the knowledge of the solute charge density ρsolute (r) = ρel (r) +

X

za δ (r − Ra ))

(13)

a

as Z E(r) =

G (r − r0 ) ρsolute (r) dr0

(14)

where za is the valence of atom a and the appropriate kernel G (r − r0 ) depends on the type of boundary conditions of the problem. Using the definition of the soft-sphere interface in Eq. (1), the gradient on the right-hand side of Eq. (12) can be expressed as the sum of individual atomic terms, corresponding to the interfaces of the different soft-spheres " ∇s(r) = ∇ 1 −

# Y a

h0a (ra ) = −

X

∇h0a (ra )

a

Y

h0b (rb ) = (s(r) − 1)

b6=a

11


X ∇h0 (ra ) a

a

h0a (ra )

(15)


Page 12 of 41

where ∇h0a (ra ) =

dh0a ∇ra dra

1 (ra − αξ RavdW )2 =√ exp − ∆2ξ π∆

!

ra ra

(16)

To account for a shrinking of a soft-sphere when embedding a charged or polar part of the solute, an additional scaling factor depending on the electric field in vacuum is introduced in the definition of the field-aware interface. As we want an homogeneous scaling of the whole soft-sphere, the scaling factor needs to account for the global field passing throught the individual atomic portions of the interface. Among the possible alternatives, a natural definition is provided by the electric flux through the exposed surface of each sphere, expressed as Z Φa =

Ean (r)dr

Z =−

E(r) · ∇h0a (ra )

Y

h0b (rb )dr

(17)

b6=a

By exploiting the flux of the electric field in vacuum, the integral over the exposed surface of the sphere is a clear proxy of the total net charge (electronic plus ionic) contained in such a spherical region, as per Gauss’ Law. This is different from the CANDLE model, 13 where the scaling depends on the local component of the normal electric field, as applied to a interface function defined on the solute’s electronic density. The product of soft-sphere functions on the right-hand side of the above equation ensures that the flux is not computed for the portions of the sphere contained inside neighboring spheres. This choice makes the connection between flux and contained charge less direct, especially for soft-spheres that have only a small fraction of their surface exposed to the continuum. However, if the overall surface of the soft-sphere is small, the effect of a field-aware scaling of the interface is bound to be negligible. Moreover, neglecting the flux through the intersection of different spheres allow us to avoid the regions of space of covalent bonds, in which small changes in the electronic structure of the system may produce significant differences in electric flux.

12


Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The new field-aware interface function is written in terms of the field flux on the spheres as sˆ(r) = 1 −

Y

a hΦ a (ra ),

(18)

a

with the field-aware soft-spheres defined as 1 r − f (Φa )αξ RavdW h ({ξa } ; r) = 1 + erf . 2 ∆ξ Φa

(19)

The field-aware function f (Φ) introduced above is chosen to ensure the correct scaling of the sphere radius as a function of net charge contained inside the sphere. Its definition and parametrization is described in the following section. Note that, since the field fluxes are integrated over the soft spheres, there now exists a cyclic dependence between the fluxes and the interface function. This can be resolved by a self consistent loop, but it would add unnecessary computational complexity to the problem. Instead, the unscaled soft spheres h0a are adopted for the calculation of the fluxes. This should produce an identical result to any reasonable margin of error: the electronic density does not dramatically vary in the vicinity of the sphere surfaces and therefore the change in flux by modifying the size of the sphere by some incremental amount should be negligible. In order to be able to exploit the new interface function in the Kohn-Sham SCF procedure and for geometry optimization, we need to characterize the derivatives of the field-aware softsphere interface with respect to the electronic density and ionic positions. Even though our starting interface is only defined on the solute’s ionic positions, as the field fluxes depends non-locally on both the electronic and ionic charge densities, both of these derivatives are now different from zero. In particular, we can express the functional derivative of the interface function with respect to the electronic density in terms of the individual field fluxes and their functional derivatives as X ∂ˆ δˆ s s 0 δΦa 0 (r, r ) = (r ) el (r). δρel ∂Φ δρ a a 13


(20)


Page 14 of 41

By applying this result to Eq. (10), we get for the interface contribution to the Kohn-Sham potential interf ace VKS (r)

Z X ∂ˆ s 0 δΦa δE interf ace [s] 0 0 = (r ) el (r) (r )dr ∂Φa δρ δs a Z X δΦa ∂ˆ s 0 δE interf ace [s] 0 0 (r ) (r )dr = (r) δρel ∂Φa δs a

(21)

Note that although the functional derivative of the interface is defined in r for a change in the charge distribution given in r0 , these coordinate systems can be separated cleanly. This is particularly useful for implementation, since it would be numerically cumbersome and expensive to manage a quantity defined on six dimensions. By separating the expression as we have shown, one can now evaluate a three-dimensional integral and multiply by the final functional derivative term. The partial derivative of the interface with respect to the field fluxes can be decomposed as a ∂ˆ s 0 ∂ˆ s dhΦ df a 0 (r ) = (r ) a a ∂Φa ∂hΦ df dΦa a dhΦa df Y Φb 0 = αξ RavdW a0 h (r ) dra dΦa b6=a b b

(22)

where the total derivative of ha is given in Eq. (16), while the derivative of the field-aware function

df dΦa

will depend on its definition (as reported in the following Section).

The partial derivatives of the interface function with respect to the ionic positions is given by ∇Rc sˆ(r) = −

X

a ∇Rc hΦ a (ra )

a

=

X a

Y

b hΦ b (rb )

b6=a

1 dha √ π∆ dra

= ∇Rc s(r) +

X a

Y ra vdW df b δca − αRa ∇rc Φa hΦ b (rb ) ra dΦa b6=a

∇Rc Φa

∂ˆ s (r) ∂Φa

14


(23)

Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, the result can be expressed as a correction to the derivative of the original SSCS interface function. When combined with Eq. (11), the field-aware contribution to interatomic forces can be written as

fcinterf ace

=

fcinterf ace,0

−

XZ

∇Rc Φa

a

δE interf ace [s] ∂ˆ s (r) (r)dr ∂Φa δs

(24)

Comparing Eq.s (24) and (21) allows to appreciate the similarities between the two expressions. The last ingredients missing in the above expressions are the derivatives of the field fluxes with respect to the solute degrees of freedom. In particular, the functional derivative of the field fluxes with respect to the electronic density can be expressed as δ δΦa (r) = δρel δρel (r)

Z

E(r0 ) · ∇h0a (r0 )

Y

h0b (rb0 )dr0

b6=a

! # Y X δ h0b (rb0 )dr0 = el zc δ(rc00 ) dr00 · ∇h0a (r0 ) G(r0 − r00 ) ρel (r00 ) + δρ (r) c b6=a Z Y = G(r0 − r) · ∇h0a (r0 ) h0b (rb0 )dr0 Z "Z

b6=a

(25) It is important to note that the final result above is free from any dependence on the electronic density. This allows us to compute this factor only once, before the SCF loop, for a given set of atomic positions. Obtaining the derivative of the flux with respect to the ionic positions is more long-winded and so we leave the steps to the supporting information, showing only

15



Page 16 of 41

the initial expansion and the final result Z ∇Rc Φa = ∇Rc

E(r) · ∇ha (ra )

Y

hb (rb )dr

b6=a

Z =−

(H [φ (zc δ(rc ))] ∇ha (ra ) + δac H [ha (ra )] E(r))

Y

hb (rb )dr

b6=a,c

Z + (δac − 1)

[E(r) · ha (ra )] ∇hc (rc )

Y

hb (rb )dr (26)

b6=a,c

expressed in terms of the Hessian H [φ (zc δ(rc ))] of the electrostatic potential of atom c and of the Hessian H [ha (ra )] of the soft sphere on atom a. All of the derivatives reported above were tested by comparing their analytic expressions with numerical finite-differences estimates.

2.3

Field-Aware Function

The main ingredient of the field-aware approach is represented by the scaling function that is responsible to shrink the solvation radii as the field flux increases. This function should assume a value of 1 for low positive or negative electric fluxes, so as to not rescale the continuum interface for non-polar atoms or non-charged systems. At the opposite, in order to reproduce the shrinking of the interface close to polar/charged residues, the field-aware function needs to decrease for high positive or negative electric fluxes. Moreover, as demonstrated in multiple preceding publications (e.g. in Ref. 21 ), the mean average error from anions is greater than cations, suggesting an asymmetric relationship between the flux and the size of the atom in question. 39,40 For the above reasons, we chose to introduce a field-aware function using a piece-wise definition similar to the one adopted for the SCCS model, 8 namely

f (Φ) = 1 − f0 t(Φ)(κ − sgn(Φ))2

16


(27)

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where t(Φ) is a piece-wise function defined by

t(Φ) =

    0    

Φ ≤ Φmin

1 (x − sin x) Φmin < Φ < Φmax 2π       1 Φ ≥ Φmax

x = 2π

Φ − Φmin Φmax − Φmin

(28)

(29)

There are four field-specific parameters to be tuned. The maximum scaling of soft spheres corresponding to positive and negative residues is controlled by the combination of a global field-aware intensity, f0 , and an asymmetry parameter, κ. The transition region between a non field-aware regime and the maximum scaling factor is controlled by two thresholds, Φmin and Φmax , of the electrostatic fluxes. A typical example of the behavior of the field-aware function is reported in Figure 1. In principles, since this function directly affects the scaling of the soft-spheres, which was previously globally tuned for all species of a specific charge, a retuning of the SSCS αξ scaling factor is also necessary. In the subsequent analysis, we find it convenient to recast the field cutoff parameters in terms of a mean flux, and a spread parameter. These are related by Φavg = (Φmax − Φmin )/2 and Φspread = Φmax − Φavg . It should also be noted that the flux values are given in Atomic Rydberg units.

3

Computational Details

Pseudopotentials were chosen according to the SSSP (standard solid state pseudopotentials) efficiency set. 41 Kohn-Sham wavefunctions and charge densities were expanded in plane waves up to kinetic energy cutoffs of 30 and 300 Ry respectively. Simulation cells were chosen to be 30.0 a.u. to reduce periodic boundary artefacts in charged systems. Note that when displaying quantities of flux, these are given in Atomic Rydberg units unless 17



field-aware function

1.00 0.98 scaling factor

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 41

0.96 0.94 0.92 0.90 4

2 0 2 electric flux through sphere

4

Figure 1: Example of the field-aware function (f (Φ), labeled scaling factor on the graph) for some typical values of parameters: κ = −0.5, f0 = 0.05, Φmin = 1.0, and Φmax = 3.0.

18


Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


specified. Hence if a sphere has a flux of 8π, the corresponding charge housed within is equivalent to one electron charge. A few selected neutral and charged molecules were studied as a function of the field-aware paramters and the soft-sphere scaling factor. The studied systems were chosen so as to be representative of the different pathological cases encountered in standard continuum models. In particular, water and its corresponding ionic species (hydroxyl anion and hydronium cation) were chosen due to the fact that the solvation free energies of both charged species show large MAE values with the soft-spheres model, even after re-parameterization for charged molecules was performed. In addition, one neutral molecule, one cation, and one anion are also picked to show the potential of this model. These are chosen from a set of 13 neutral molecules, 15 cations and 15 anions, as the worst performing species when calculating the solvation energy by the SSCS. Following this criterium, propanoic acid, the diethylether cation, and the benzyl alcohol anion were selected. Implementation and simulations were performed using a development version of the Environ-1.0 42 plugin for multiscale embeddings in condensed-matter simulations, as coupled with the Quantum Espresso package. 43,44 As discussed previously, the field aware function has four tunable parameters. In addition, the soft-sphere parameters, which for the neutral set were chosen as a best fit including species with internal polarization, need to be reconsidered. For simplicity we begin by just considering how the scaling factor changes. Hence, for preliminary analysis of neutral molecules, we consider values of αξ ranging from 1.12 to 1.20. This is under the assumption that the optimal scaling factor from the soft-spheres is a result of a compromise between molecules without internal polarization effects, which should thus require a larger spherical functions, and molecules with internal polarization effects, which should require smaller spherical functions for specific atoms. For cations and anions, the asymmetry is set to -1 and 1, respectively, to minimize the computational costs of gathering results, since the flux values for the species we test for all have the same sign, negative for anions and positive for cations. A range of electric flux

19



Page 20 of 41

parameters compatible with the values observed for the selected molecules was considered (see fluxes reported in Figure 2), again to minimize computational time. This range may need to be extended for a more comprehensive parameterization of the model that includes multivalence charged systems. The field awareness f0 is chosen to be the upper limit on the reduction of the radii of the spherical functions. This was originally chosen to match the reduction as seen in the soft-spheres model between the neutral set and the anion set, but considering the reparameterization of this scaling as part of the field-aware model, it was found that increasing this value yielded better results. Note that as defined, both the asymmetry factor and the field awareness control the global intensity of the field-aware procedure.

4

Results

Figure 2: Comparison of electric flux values for each individual atom for the six molecules. A color map proportional to the absolute value of the flux is adopted. Atoms showing the largest absolute electric fluxes are identified in Figure 3.

Table 1: Best Parameter Set Molecule Water Propanoic Acid Water (+) Diethyl Ether (+) Water (-) Benzyl Alcohol (-)

Abs Error Asymmetry 0.0163 0.1035 0.2950 1.2634 0.4813 10.9610

0.0 1.0 1.0 1.0 -1.0 -1.0 20

α

f0

1.16 1.20 1.18 1.12 1.20 1.12

0.06 0.05 0.05 0.06 0.07 0.07


Flux Mean Flux Spread 1.5 2.0 1.5 1.5 1.5 2.0

0.6 0.4 0.6 0.6 0.6 1.5

Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3: Set of molecules tested. The atom possessing the maximum absolute electric flux is highlighted by a red circle, and its flux value is reported

21



The field-aware model is able to converge for the majority of parameter combinations, although numerical issues arise with the SCF procedure when specific ranges of parameter values are selected. This is linked to the sharpness of the field-aware potential when the fieldaware function has a transition that is too sharp with respect to the flux. Nevertheless, for the ranges of parameters tested in this study, involving about 3000 different combinations, 99% of the SCF and geometry optimization calculations converged successfully, validating the analytic differentiation of field-aware interfaces reported above. Due to the number of parameters, it is challenging to optimize this model from scratch. Our main aim here is thus to display insights from a broad sweep, as well as observations that may simplify the parameterization process. As this process will benefit from a significant number of experimental and theoretical results on solvation of charged molecules and interfaces (e.g. results of differential capacitance of metal substrates), its completion lies outside of the scope of this work. We start by showing in Figure 2 the final values of the electric flux through each soft sphere for the selected molecules. For full numerical values, we refer to the Supporting Information. The values attained for the molecules chosen match expectations, neutral molecules center around zero flux, with the oxygen atom being where most of the negative charge imbalance lies and anions show large relative charge differences. One can see how the field-aware function acts on individual atoms by comparing Figure 2 to the field-aware function, as in Figure 1. Any flux with an absolute value over the cutoff will result in a minimum value for the scaling factor, which varies from 1 (no scaling) to a value determined by the field awareness and asymmetry parameters. In particular, due to the definition of the asymmetry factor, a deviation of 0.5 away from symmetry produces a heavy bias towards spheres of negative flux. The choice of flux cutoff parameters can be simplified by considering typical flux values for a training set. For this investigation we limit our analysis to consider how an adjustment of parameters for individual species can improve on absolute solvation error values. In Table 1 we report the best performing set of parameters for each

22


Page 22 of 41

Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


considered compound. Despite most of these parameters fall in a similar range of values, the corresponding field-aware functions are clearly not representative of the values of a fully optimized scaling function. For example, in the case of the chosen anions, the absolute flux through the oxygen atoms is much greater than the other atoms and as such, only the sphere around this atom is appreciably altered. The kind of parameterization required will also depend on the broadness of the task for which the field-aware function is fitted. A conservative approach to the field-aware function would be to tune it so that it matches the shrinking of the interfaces observed for charged species. In particular, Fisicaro et al. reported a optimized value of the homogeneous softsphere scaling parameter αξ that goes from a value of 1.12 for neutral compounds to a value of 1.1 for cations, reaching 0.98 for anions. To reproduce this trends, a field-aware function with sharp transitions between low and high flux values would be sufficient. However, if the field awareness is designed to also improve the description of polar groups in neutral molecules, a more smoothly varying field-aware function needs to be considered.

4.1

Propanoic Acid

Propanoic acid (figure 4) is the worst performing neutral molecule for the standard SSCS model, out of the 13 testing neutral molecules. This acidic compound is expected to produce larger errors, due to its polar nature. The SSCS model underestimates the solvation energy by 2.53kcal/mol, suggesting a scaling factor that is too small. Having a function that reduces the size of the spherical functions would not be able to improve on the SSCS result and would only worsen the result. The problem is that the homogeneous scaling factor optimized in the original SSCS model represents a trade-off between apolar and polar groups, and is thus on average smaller the needed for apolar groups, while being larger than needed for polar regions. As a result, when the global scaling factor is allowed to vary above the original soft-sphere value, improved results are observed. In particular, in the region we considered in our five dimensional parameter search, the global scaling factor and the asymmetry had 23



Absolute Solvation Error: Propanoic Acid 1.00

4.8

0.75

4.2

0.50

3.6

0.25

3.0

0.00

2.4

1

asymmetry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 41

0.25

1.8

0.50

1.2

0.75

0.6

1.00 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 mean flux

0.0

Figure 4: Absolute solvation free energy error for propanoic acid in kcal/mol as a function of the mean flux threshold and the asymmetry parameter. The flux spread is fixed at 0.4, the field awareness is fixed at 0.05, while an homogenous scaling parameter αξ = 1.20 is considered

Figure 5: Propanoic Acid. Unscaled interface (left) represents the spherical functions as calculated by standard SSCS. Scaled interface (right) is based off an optimized set of parameters, as reported in Table 1.

24


Page 25 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the greatest influence on the solvation energy for this acid. In fact, a reparameterization with α = 1.20 brought the calculated solvation energy remarkably close to the experimental value (see results in Table 1), a result that is consistent with the effect of interface size on the solvation energy. In addition, promising trends are shown in Figure 5, featuring the behavior of the absolute error on solvation free energies against the asymmetry and the average flux values. This figure shows that for the chosen fixed parameters, putting less weight on negative fluxes provides particularly small errors. Indeed, the region where the asymmetry is negative appears to be more stable with respect to the flux average parameter. This is consistent with the fact that there is a single oxygen atom with appreciable negative flux, which is always accounted for in the range chosen for the maximum cutoff flux.

4.2

Diethyl Ether

As seen in Figure 2, this cation has little variance in flux among its atomic constituents. In the standard SSCS model, the solvation free energy of diethyl ether is overestimated by 16.57kcal/mol even after scaling. Despite the good performance of the SSCS model after reparameterization for the cations, there remain a number of species that have significant solvation errors. In Figure 6 we plot the field awareness against the flux average to show how the absolute error for this species can be drastically improved without a modification of αξ to a lower value, which is the standard strategy when accounting for cations. The region to the right is one where the flux values are too low to be considered by the field-aware function and thus the error is equal to the standard SSCS result for αξ = 1.12. Note that since the minimum and maximum possible scaling is controlled by the choice of parameters under investigation, as long as the model is tuned sensibly, the worst case scenario ought to be as one can see in the region to the right of this plot, i.e. no improvement. On the opposite end, when parameters such as the ones varied in the plot are properly chosen, the absolute solvation error can improve from 18.07kcal/mol to 1.35kcal/mol. Despite the large change in solvation free energy, the overall effect on the shape of the interface is mostly noticeable 25



Absolute Solvation Error: Diethyl Ether Cation 20.0

field awareness

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 41

0.065

17.5

0.060

15.0

0.055

12.5

0.050

10.0

0.045

7.5

0.040

5.0

0.035

2.5

0.030 1.5

2.0

2.5

3.0 3.5 mean flux

4.0

4.5

0.0

Figure 6: Absolute solvation free energy error for diethyl ether in kcal/mol, as a function of mean flux threshold and field awareness parameter. The flux spread is fixed at 0.6, and an homogenouse scaling factor of αξ = 1.20 is adopted. An asymmetry κ = 1 is selected, as in this molecule there are no atoms with negative fluxes.

26


Page 27 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Diethyl Ether. Unscaled interface (visualized on the left half of the molecule) represents the spherical functions as calculated by standard SSCS. Scaled interface (represented on the right half of the molecule) is based off an optimized set of parameters, as reported in Table 1. The effects of the field-aware scaling are small, mostly noticeable for the acidic hydrogen atom bonded to the oxygen.

27



for the acidic hydrogen atom bonded to the oxygen, as reported in Figure 7.

4.3

Benzyl Alcohol

Absolute Solvation Error: Benzyl Alcohol Anion 1.20

19

1.19

18

1.18

17 16

1.17 alpha

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 41

15

1.16

14

1.15

13

1.14

12

1.13 1.12 0.03

11 0.04

0.05 0.06 field awareness

0.07

10

Figure 8: Absolute solvation free energy error for benzyl alcohol in kcal/mol as a function of homogenous scaling αξ and field awareness parameter. The flux average and spread parameters are set to values of 2.0 and 0.5, respectively, so that only a couple of atoms are affected by the field-aware procedure. For the range of parameters we currently have tested for, a potential optimal value for this particular species is somewhat far off. Nonetheless, like in the diethyl ether case, a significant improvement of the results can be obtained by properly selecting the field-aware parameters. In particular, as reported in Figure 8, simply by increasing the field-awareness of the model, one can achieve the same result as obtained by retuning of the SSCS model. The free energy of solvation of benzyl alcohol is overestimated by 19.66kcal/mol before reparameterization, reducing to an error of 10.11kcal/mol with αξ = 0.98. The figure demonstrates how a similar solvation free energy can be achieved without decreasing αξ and instead with a field-aware 28


Page 29 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9: Benzyl Alcohol. Unscaled interface (left) represents the spherical functions as calculated by standard SSCS. Scaled interface (right) is based off an optimized set of parameters as reported in Table 1. function tuned properly. In fact the figure seems to suggest that we could use a larger field awareness parameter. However, by visual inspection of the interface function for the best set of parameters (Figure 9, it appears that the improvement in solvation free energy for this compound is obtained by significantly reducing the size of the soft-sphere around the oxygen atom. Increasing even further this effect would lead to an unphysically small solvation radius. The remaining discrepancy in solvation free energies should thus be ascribed to the aromatic portion of the molecule and to the assumptions of the model used to describe non-electrostatic interactions.

4.4

Water

As in the previous systems, we consider the prototypical cases of the aqueous solvation of water and its derived ionic species, hydroxyl anion and hydronium cation. Although free energy of solvation of a solvent molecule in its own solution represents a pathological application of continuum solvation models, water is an interesting test case, which performs well in its neutral state, but remarkably poorly in its singly charged ionic state due to its small size. Even after the standard scaling that the SSCS performs, while the absolute solvation error from neutral water is only 0.81kcal/mol, the errors from the cation and anion 29



alpha

Absolute Solvation Error: Water Cation 1.20

48

1.19

42

1.18

36

1.17

30

1.16

24

1.15

18

1.14

12

1.13

6

1.12 0.03

0.04


0.07

0

Absolute Solvation Error: Water 1.20

2.00

1.19

1.75 1

alpha

1.18

1.50

1.17

1.25

1.16

1.00

1.15

0.75

1.14

0.50

1.13 1.12 0.03

0.25

1

0.04


0.07

0.00

Absolute Solvation Error: Water Anion 1.20

16

1.19

14

1.18

12

1.17

10

1.16

8

1.15

6

1.14

4

1.13

2 1

alpha

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 41

1.12 0.03

0.04


0.07

0

Figure 10: Absolute solvation free energy errors of neutral water (central panel), hydronium ion (top panel), and hydroxil anion (bottom panel) in kcal/mol as functions of field awareness parameter and homogenous scaling parameter αξ . The asymmetry parameter is set to zero, +1, and -1 for the three systems, respectively. Mean and spread values of the flux thresholds are set to 1.5 and 0.6. 30


Page 31 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are 18.56kcal/mol and 12.13kcal/mol respectively. The maximum flux cutoff is chosen to be close to the minimum flux of water, resulting in a general scaling of all atoms. In the different panels of Figures 10, we can see that a linear relation exists between the increase in the αξ scaling factor and the corresponding optimal value of field awareness f0 . While the slope of this relationship is the same for cation and anion, the slope for neutral molecule is divided by factor of 4 due to the chosen values of the asymmetry parameter (κ = 0 for neutral vs. κ = ±1 for cation/anion) that introduce and additional scaling in the fieldaware strength. Increasing the field-awareness from zero improves the solvation energy with respect to the experimental value, and as expected, the optimal field awareness parameter is different for all three cases. This disagreement in optimal field awareness can be overcome by an appropriate choice of an asymmetry parameter common to the three systems. One possible method for choosing the asymmetry parameter is by combining the optimal field awareness parameter individually obtained for cation (fc ) and anion (fa ) for an asymmetry of +1 and -1, respectively. The maximum shrinking effect obtained in the two systems is thus max fa t (Φ) (−1 − sgn (Φ))2 = 4fa ,

(30)

max fc t (Φ) (1 − sgn (Φ))2 = 4fc .

(31)

Φ

Φ

By imposing that the optimal global values of κ and f0 reproduce both effects, we obtain a system of two equations in two unknowns,

f0 (κ − 1)2 = 4fa ,

(32)

f0 (κ + 1)2 = 4fc .

(33)

Given the above conditions, the determination of κ is trivial

κ=

√ 1 (R + 1 ± 2 R), R−1 31


(34)


where R = fc /fa , and consequently f0 can be found. For the hydroxyl anion and hydronium cation, parameters that produced reliable solvation free energy estimates for both would be κ = −0.16 and f0 = 0.24. Simulations exploiting a common set of parameters, including the two derived above, produce a systematic improvement in the accuracy of the calculations for all water species (Φmin = 0.1, Φmax = 2.9). Solvation free energy errors of 0.55kcal/mol, 2.67kcal/mol, and 3.26kcal/mol are obtained for neutral, cationic, and anionic water, respectively.

5

Conclusions

In conclusion, a new field-aware approach for continuum interfaces is proposed and analyzed. The field-aware interface has analytical, smoothly varying derivatives, allowing for its use in accurate geometry optimization or in molecular dynamics simulations. By exploiting reciprocal-space techniques to compute the electric field and derived properties, the overhead linked to this non-local definition is substantially reduced. The field-aware interface requires a number of new parameters, which can tuned to reproduce experimental data on solvation free energies of neutral and charged compounds. An analysis of the main parameters involved in the model and of their effects on solvation free energies of a significant set of test molecules was also presented. Significant improvements over the standard approaches are demonstrated for a specifically tuned set of parameters. Nonetheless, an optimal parametrization aimed at improving the description of both neutral and charged compounds will require a more systematic search over the parameter space. In these regards, it will be important to focus the search of the individual parameters on specifically selected results, in order to isolate the physical effects introduced by this model. While experimental data on solvation free energies of neutral and charged molecules ? will be the main ingredient of the parameterization, first-principles molecular dynamics simulations of salts dissolution can represent a second unbiased source of determination of some of the parameters of the model. 45 Eventually,

32


Page 32 of 41

Page 33 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parametrization on the physical properties of electrochemical interfaces, e.g. the differential capacitance of noble metal substrates, 46,47 will also represent a potential strategy for the tuning of the presented approach.

Supporting Information Available Derivatives of the field fluxes with respect to the quantum-mechanical degrees of freedom (electronic density and ionic positions). Detailed values of the field fluxes computed for each atom of the tested molecules. Best performing sets of parameters for each tested molecule and ion.

Acknowledgement We acknowledge start-up funding support from the University of North Texas (UNT). We acknowledge computational support from the computing facilities at UNT and through the Center for Advanced Scientific Computing and Modeling (CASCaM).

References (1) Tomasi, J.; Persico, M. Molecular-Interactions in Solution - an Overview of Methods Based on Continuous Distributions of the Solvent. Chemical Reviews 1994, 94, 2027– 2094. (2) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chemical Reviews 2005, 105, 2999–3093. (3) Cramer, C. J.; Truhlar, D. G. Implicit solvation models: Equilibria, structure, spectra, and dynamics. Chemical Reviews 1999, 99, 2161–2200.

33



(4) Andreussi, O.; Fisicaro, G. Continuum embeddings in condensed-matter simulations. International Journal of Quantum Chemistry 2019, 119, e25725. (5) Pascual-Ahuir, J. L.; Silla, E. GEPOL: An improved description of molecular surfaces. I. Building the spherical surface set. Journal of Computational Chemistry 1990, 11, 1047–1060. (6) Pascualahuir, J. L.; Silla, E.; Tunon, I. Gepol - an Improved Description of MolecularSurfaces .3. A New Algorithm for the Computation of a Solvent-Excluding Surface. Journal of Computational Chemistry 1994, 15, 1127–1138. (7) Fattebert, J. L.; Gygi, F. Density functional theory for efficient ab initio molecular dynamics simulations in solution. Journal of Computational Chemistry 2002, 23, 662– 666. (8) Andreussi, O.; Dabo, I.; Marzari, N. Revised self-consistent continuum solvation in electronic-structure calculations. The Journal of Chemical Physics 2012, 136, 064102. (9) Fisicaro, G.; Genovese, L.; Andreussi, O.; Mandal, S.; Nair, N. N.; Marzari, N.; Goedecker, S. Soft-Sphere Continuum Solvation in Electronic-Structure Calculations. Journal of Chemical Theory and Computation 2017, 13, 3829–3845. (10) Andreussi, O.; Hormann, N. G.; Nattino, F.; Fisicaro, G.; Goedecker, S.; Marzari, N. Solvent-aware Interfaces in Continuum Solvation. Journal of Chemical Theory and Computation 2019, acs.jctc.8b01174. (11) Quan, C.; Stamm, B. Mathematical analysis and calculation of molecular surfaces. Journal of Computational Physics 2016, 322, 760–782. (12) Sundararaman, R.; Schwarz, K. A.; Letchworth-Weaver, K.; Arias, T. A. Spicing up continuum solvation models with SaLSA: The spherically averaged liquid susceptibility ansatz. The Journal of Chemical Physics 2015, 142, 054102. 34


Page 34 of 41

Page 35 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(13) Sundararaman, R.; Goddard, W. A. The charge-asymmetric nonlocally determined local-electric (CANDLE) solvation model. The Journal of Chemical Physics 2015, 142, 064107. (14) Fattebert, J.-L.; Gygi, F. First-principles molecular dynamics simulations in a continuum solvent. International Journal of Quantum Chemistry 2003, 93, 139–147. (15) Jinnouchi, R.; Anderson, A. B. Aqueous and Surface Redox Potentials from SelfConsistently Determined Gibbs Energies. The Journal of Physical Chemistry C 2008, 112, 8747–8750. (16) Scherlis, D. A.; Fattebert, J.-L.; Gygi, F.; Cococcioni, M.; Marzari, N. A unified electrostatic and cavitation model for first-principles molecular dynamics in solution. The Journal of Chemical Physics 2006, 124, 074103. (17) Sánchez, V. M.; Sued, M.; Scherlis, D. A. First-principles molecular dynamics simulations at solid-liquid interfaces with a continuum solvent. The Journal of Chemical Physics 2009, 131, 174108. (18) Dziedzic, J.; Helal, H. H.; Skylaris, C.-K.; Mostofi, A. A.; Payne, M. C. Minimal parameter implicit solvent model for ab initio electronic-structure calculations. EPL (Europhysics Letters) 2011, 95, 43001. (19) Ringe, S.; Oberhofer, H.; Hille, C.; Matera, S.; Reuter, K. Function-Space-Based Solution Scheme for the Size-Modified PoissonâĂŞBoltzmann Equation in Full-Potential DFT. Journal of Chemical Theory and Computation 2016, 12, 4052–4066. (20) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Accounts of Chemical Research 2008, 41, 760–768. (21) Dupont, C.; Andreussi, O.; Marzari, N. Self-consistent continuum solvation (SCCS): The case of charged systems. The Journal of Chemical Physics 2013, 139, 214110. 35



(22) Hörmann, N. G.; Andreussi, O.; Marzari, N. Grand canonical simulations of electrochemical interfaces in implicit solvation models. The Journal of Chemical Physics 2019, 150, 041730. (23) Mathew, K.; Sundararaman, R.; Letchworth-Weaver, K.; Arias, T. A.; Hennig, R. G. Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways. The Journal of Chemical Physics 2014, 140, 084106. (24) Hörmann, N. G.; Guo, Z.; Ambrosio, F.; Andreussi, O.; Pasquarello, A.; Marzari, N. Band alignment at semiconductor-water interfaces using explicit and implicit descriptions for liquid water. In preparation 2019, (25) Barone, V.; Cossi, M.; Tomasi, J. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. The Journal of Chemical Physics 1997, 107, 3210–3221. (26) Pomogaeva, A.; Chipman, D. M. Field-Extremum Model for Short-Range Contributions to Hydration Free Energy. Journal of Chemical Theory and Computation 2011, 7, 3952–3960. (27) Pomogaeva, A.; Thompson, D. W.; Chipman, D. M. Modeling short-range contributions to hydration energies with minimal parameterization. Chemical Physics Letters 2011, 511, 161–165. (28) Molavi Tabrizi, A.; Goossens, S.; Mehdizadeh Rahimi, A.; Knepley, M.; Bardhan, J. P. Predicting solvation free energies and thermodynamics in polar solvents and mixtures using a solvation-layer interface condition. The Journal of Chemical Physics 2017, 146, 094103. (29) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. Uff, a Full Periodic-Table Force-Field for Molecular Mechanics and Molecular-Dynamics Simulations. Journal of the American Chemical Society 1992, 114, 10024–10035. 36


Page 36 of 41

Page 37 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Klamt, A.; Schuurmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. Journal of the Chemical Society, Perkin Transactions 2 1993, 799–805. (31) Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. The Journal of Physical Chemistry 1995, 99, 2224–2235. (32) Chipman, D. M. New formulation and implementation for volume polarization in dielectric continuum theory. The Journal of Chemical Physics 2006, 124, 224111. (33) Cancès, E.; Mennucci, B.; Tomasi, J. A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics. The Journal of Chemical Physics 1997, 107, 3032–3041. (34) Mennucci, B.; Cancès, E.; Tomasi, J. Evaluation of Solvent Effects in Isotropic and Anisotropic Dielectrics and in Ionic Solutions with a Unified Integral Equation Method: Theoretical Bases, Computational Implementation, and Numerical Applications. The Journal of Physical Chemistry B 1997, 101, 10506–10517. (35) Amovilli, C.; Mennucci, B. Self-Consistent-Field Calculation of Pauli Repulsion and Dispersion Contributions to the Solvation Free Energy in the Polarizable Continuum Model. The Journal of Physical Chemistry B 1997, 101, 1051–1057. (36) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. The Journal of Physical Chemistry A 1998, 102, 5074–5085. (37) Bondi, A. van der Waals Volumes and Radii. The Journal of Physical Chemistry 1964, 68, 441–451. (38) Cococcioni, M.; Mauri, F.; Ceder, G.; Marzari, N. Electronic-Enthalpy Functional for Finite Systems Under Pressure. Physical Review Letters 2005, 94, 145501. 37



(39) Bardhan, J. P.; Jungwirth, P.; Makowski, L. Affine-response model of molecular solvation of ions: Accurate predictions of asymmetric charging free energies. The Journal of Chemical Physics 2012, 137, 124101. (40) Bardhan, J. P.; Knepley, M. G. Communication: Modeling charge-sign asymmetric solvation free energies with nonlinear boundary conditions. The Journal of Chemical Physics 2014, 141, 131103. (41) Prandini, G.; Marrazzo, A.; Castelli, I. E.; Mounet, N.; Marzari, N. A Standard Solid State Pseudopotentials (SSSP) library optimized for accuracy and efficiency (Version 1.0, data download). 2018. (42) Andreussi, O.; Nattino, F.; Dabo, I.; Timrov, I.; Fisicaro, G.; Goedecker, S.; Marzari, N. Environ 1.0: a library for environment effect in first-principles simulations of materials. Environ_1.0 2018, www.quantum\–environ.org. (43) Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 2009, 21, 395502. (44) Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. Journal of Physics: Condensed Matter 2017, 29, 465901. (45) Ringe, S.; Oberhofer, H.; Reuter, K. Transferable ionic parameters for first-principles Poisson-Boltzmann solvation calculations: Neutral solutes in aqueous monovalent salt solutions. The Journal of Chemical Physics 2017, 146, 134103. (46) Nattino, F.; Truscott, M.; Marzari, N.; Andreussi, O. Continuum models of the electrochemical diffuse layer in electronic-structure calculations. The Journal of Chemical Physics 2019, 150, 041722.

38


Page 38 of 41

Page 39 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(47) Sundararaman, R.; Letchworth-Weaver, K.; Schwarz, K. A. Improving accuracy of electrochemical capacitance and solvation energetics in first-principles calculations. The Journal of Chemical Physics 2018, 148, 144105.

39



TOC Graphic

40


Page 40 of 41

Page 41 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Biography Dr. Andreussi studied Physical Chemistry at Scuola Normale Superiore of Pisa and at the Chemistry Department of the University of Pisa, graduating in 2003 with a thesis on modelling the effect of metal nanoparticles on molecular excited states using continuum approaches, under the supervision of Prof. Tomasi. During his Ph.D. in the group of Prof. Parrinello, Dr. Andreussi familiarized with the field of condensed-matter and statistical mechanics simulations, using first-principles and classical molecular dynamics simulations to study water overlayers on different substrates. After graduation, in 2008 he joined as a postdoctoral associate the research group of Prof. Marzari at MIT where he kept investigating solvents and solvation effects, by means of classical and multiscale simulations. Following additional research and teaching experiences, in 2018 he joined the University of North Texas as an Assistant Professor in the Department of Physics. His main interests are in the development and implementation of continuum and multiscale embedding approaches, with applications in a significant range of topics, including the study of room-temperature ionic liquids, the modeling of devices composed by metal nanoparticles and light-harvesting proteins, and the characterization of heterogeneous catalysis and electrochemistry.

41


Field-Aware Interfaces in Continuum Solvation - The Journal of

Recommend Documents