Soft-Sphere Continuum Solvation in Electronic-Structure Calculations

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Soft-sphere continuum solvation in electronic-structure calculations Giuseppe Fisicaro, Luigi Genovese, Oliviero Andreussi, Sagarmoy Mandal, Nisanth N. Nair, Nicola Marzari, and Stefan Goedecker J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 19 Jun 2017 Downloaded from http://pubs.acs.org on June 20, 2017

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 56

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

Journal of Chemical Theory and Computation

Soft-sphere continuum solvation in electronic-structure calculations G. Fisicaro,∗,† L. Genovese,‡ O. Andreussi,¶,§ S. Mandal,k N. N. Nair,k N. Marzari,§ and S. Goedecker† †Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland ‡Laboratoire de simulation atomistique (L_Sim), SP2M, INAC, CEA-UJF, F-38054 Grenoble, France ¶Institute of Computational Science, Universita’ della Svizzera Italiana,Via Giuseppe Buffi 13, CH-6904 Lugano, Switzerland §Theory and Simulations of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, Station 12, CH-1015 Lausanne, Switzerland kDepartment of Chemistry, Indian Institute of Technology Kanpur, Kanpur 208016, India E-mail: [email protected]

Abstract We present an implicit solvation approach where the interface between the quantummechanical solute and the surrounding environment is described by a fully continuous permittivity built up with atomic-centered “soft” spheres. This approach combines many of the advantages of the self-consistent continuum solvation model in handling solutes and surfaces in contact with complex dielectric environments or electrolytes in electronic-structure calculations. In addition it is able to describe accurately both

1

ACS Paragon Plus Environment


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

neutral and charged systems. The continuous function, describing the variation of the permittivity, allows to compute analytically the non-electrostatic contributions to the solvation free energy that are described in terms of the quantum surface. The whole methodology is computationally stable, provides consistent energies and forces, and keeps the computational efforts and runtimes comparable to those of standard vacuum calculations. The capabilitiy to treat arbitrary molecular or slab-like geometries as well as charged molecules is key to tackle electrolytes within mixed explicit/implicit frameworks. We show that, with given, fixed atomic radii, two parameters are sufficient to give a mean absolute error of only 1.12 kca/mol with respect to the experimental aqueous solvation energies for a set of 274 neutral solutes. For charged systems, the same set of parameters provides solvation energies for a set of 60 anions and 52 cations with an error of 2.96 and 2.13 kcal/mol, respectively, improving upon previous literature values. To tackle elements not present in most solvation databases, a new benchmark scheme on wettability and contact angles is proposed for solid-liquid interfaces, and applied to the investigation of the stable terminations of a CdS (11¯20) surface in an electrochemical medium.

1

Introduction

The computational study of matter in various environments is a continuously growing field in solid state physics and chemistry. Systems of interest are for instance molecules, clusters or surfaces in contact with solvents. 1 A possibility is to include explicitly in the simulation domain all the atoms and molecules of the solution. Such an explicit treatment is in principle the natural way to account for solvent effects in a first-principles scheme. However, this approach enormously increases the computational cost and limits at the same time the size of the system contained in the explicit dielectric medium. 2 As a consequence, the study of the solute-solvent interactions at length scales larger than molecular sizes, or the generation of long molecular dynamics trajectories, require big computational efforts.

2


Page 2 of 56

Page 3 of 56

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


Numerically expensive investigations in material science such as structure predictions 3 or reaction path determinations 4 would become very time-consuming if not altogether unaffordable. During an exploration of the potential energy surface, most of the minima would arise from rearrangements of solvent molecule. For an hydrated surface, a lot of time would consequently be spent moving water molecules rather than sampling the solid-water interface that is of interest. The explicit description of water also suffers from some well-known limitations of current state-of-the-art ab-initio methods, in particular regarding its structural and dielectric properties. 5–8 The alternative is an implicit description of the solvent effects, while still treating the other parts of the system explicitly on an atomic level. Such an explicit/implicit treatment requires three main ingredients: • A dielectric cavity represented by a proper function ǫ(r) mimicking the surrounding solvent of a solute as a continuum dielectric; • A solver for the generalized Poisson equation 9 ∇ · ǫ(r)∇φ(r) = −4πρ(r),

(1)

where φ(r) is the potential generated by a given charge density ρ(r); • A model for the non-electrostatic terms to the total free energy of solvation. The dielectric function ǫ(r) has to take the value of one where the solute is placed to solve a vacuum-like quantum problem, and the bulk dielectric constant ǫ0 outside. Several implicit solvation models have been reported in the quantum chemistry literature, starting from the earliest work of Onsager 10 to the widespread polarizable continuum model (PCM) developed by Tomasi and co-workers. 11–13 In the PCM formulation, the cavity surrounding the solute is described by a sharp and discontinuous dielectric function ǫ(r), and a polarization charge density is exactly local3



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

ized at the interface between the vacuum and the dielectric. In this way, the dielectric environment is represented by an effective surface polarization charge, reducing the threedimensional dielectric problem into a two-dimensional one. In its original formulation this approaches suffers from discontinuities in the atomic forces and numerical singularities, that are a consequence of sharp cavities. This drawback hindered a PCM extension to ab-initio molecular dynamics (MD). Modified versions have been proposed in the literature, 14,15 as well as advanced definitions of the dielectric cavity in terms of an isodensity of the solute charge density. 16,17 Another possibility is to represent the interface between the region inside the cavity and the implicit solvent region outside the cavity by a smooth function of the solute electronic density ρ(r). 18–20 Non-electrostatic contributions to the solvation free energy can be expressed as function of the cavity volume and surface area. 19,21,22 For a charge-dependent approach ǫ[ρ](r) the functional derivative with respect to ρ(r) of the Hamiltonian delivers new additional terms to the Kohn-Sham (KS) potential, both from the electrostatic energy and the non-electrostatic contributions. 22 These formulations are elegant and easy to be parametrized. The cavity changes self-consistently during the wave function optimization and forces do not depend explicitly on the shape of the cavity. However, there are a number of drawbacks. The smooth transition between the quantum and continuum regions gives rise to an inherently three-dimensional problem, while the self-consistent nature of the approach requires solving a new problem at each step of electronic optimization: thus the approach is intrinsically computationally more demanding than the rigid 2D alternatives such as PCM. Moreover, accuracy might be affected due to the low number of parameters; neutral and charged molecules cannot be treated simultaneously, as these require different parametrizations, 23 which lead to different Hamiltonians; furthermore, it is not possible to locally modify cavity sizes looking only to the electronic density. Greater flexibility can be added to these smooth functional approaches by defining them in terms of a fictitious electronic density, which is kept rigid and defined in terms of pa-

4


Page 4 of 56

Page 5 of 56

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


rameterized functions centered at atomic positions. 24 Although this approach may overcome some of the computational limitations of self-consistent alternatives, its original formulation only focused on the electrostatic terms of the solvation free energy. To make an mixed explicit/implicit treatment attractive to material science, it is necessary to fulfill several requirements: • Accurate forces and a numerical cost comparable to standard vacuum calculations to make molecular dynamics or extensive potential energy surface (PES) explorations feasible; • A small number of cavity parameters, in order to easily fit a new environment or find the best solvent conditions in a given range; • Exact treatment of molecular or slab-like geometries, typical of systems plunged in a solvent; • Ability to treat neutral and charged molecules simultaneusly, in order to tackle complex interfaces (e.g. a double layer). The self-consistent continuum solvation (sccs) model developed by Andreussi et al. 22 goes in this direction, satisfying the first three conditions. In the present work an improved implicit solvation model is proposed. It preserves the flexibility of PCM to locally vary the cavity’s size around each atom as well as the straightforward definition of the nonelectrostatic contributions to the total energy of the charge-dependent solvation models, function of the cavity volume and surface. 22 The dielectric cavity is based on analytic smooth spheres centered on each atom, and it is fully continuous from the vacuum-like inner regions to the external bulk solvent. Because of these features, we named our approach “soft-sphere” model. The proposed formulation satisfies all the requirements listed above. It gives the same accuracy for the ionic forces as an ordinary vacuum calculation. Compared to the similar approache of Sanchez et al., 24 that is based on functionals of a fictitious solute density, our 5



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

definition of the cavity in terms of atom-centered smooth functions simplifies significantly the calculation of derivatives. In addition, the continuous function, describing the variation of the permittivity, allows us to compute analytically all the non-electrostatic contributions to the free energy. Once we fix the atomic radii, a given environment can be fitted with two parameters in a straightforward way. Details of the model will be presented on Sec. 2. Our implementation has a small numerical overhead, thanks to an efficient solver for the generalized Poisson equation, 9 based on a preconditioned conjugate gradient (PCG) algorithm. The method has been tested within Kohn-Sham density functional theory (DFT) as implemented in the BigDFT package. 25,26 Thanks to the underlining electrostatic solver of this package, 27 the solvation problem for slab or isolated systems is solved with the correct boundary conditions. An extensive parametrization study has been carried out both for neutral molecules and ions in water and two non-aqueous solvents. Results are reported in Sec. 3. Our tests suggest that our approach can well handle slab-like boundary conditions. In Sec. 4 we propose a benchmarking database for slab like geometries, that is based on the wetting properties and solvent contact angles of solvated surfaces. Standard protocols parametrize implicit solvation models on a database of experimental free solvation energies for neutral and ionic solutes in several solvents. 28 However, if applied to solid-liquid interfaces in material science, several issues arise. The molecules of the standard benchmark sets have a limited representation of some elements of the periodic table. Many of these underrepresented elements, such as transition metal atoms, play however an important role in many materials science problems. Solar-energy harvesting in dye-sensitized 29 and hybrid cells 30,31 or electro-catalytic water splitting 32,33 are just two simple examples where such materials are extensively exploited. Finally, as an illustrative application of the approach we present a test case of relevance in material science, namely the prediction of stable terminations of the CdS (11¯20) surface

6


Page 6 of 56

Page 7 of 56

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 an electrochemical medium.

2

Soft-sphere model

The electrostatic solvation energy is defined as the difference between the total energy of a given atomic system in the presence of the dielectric environment Gel and the energy of the same system in vacuum G0 ∆Gel = Gel − G0 .

(2)

A full comparison with experimental solvation energies needs the inclusion of non-electrostatic contributions to obtain the total free energy of solvation ∆Gsol . In this case the main terms in the solute Hamiltonian, as introduced by PCM, 12 are

∆Gsol = ∆Gel + Gcav + Grep + Gdis + ∆Gtm + P ∆V,

(3)

where ∆Gel is the electrostatic contribution, Gcav the cavitation energy, i.e. the energy necessary to build up the solute cavity inside the solvent medium. Grep is a repulsion term representing the continuum counterpart of the short-range interactions induced by the Pauli exclusion principle, whilst the dissociation energy Gdis reflects van der Waals interactions. The thermal term ∆Gtm accounts for the vibrational and rotational changes and, finally, P ∆V includes volume changes in the solute Hamiltonian. In the case of a charge density dependent cavity, it has been shown that the nonelectrostatic terms Gcav , Grep , and Gdis can be expressed as linear functions of the “quantum surface” S and “quantum volume” V of the dielectric cavity. 19,22 In particular, Scherlis et al. 19 extended the work of Fattebert and Gygi 18 including a cavitation term expressed as a product of the experimental surface tension of the solvent γ and the quantum surface S of the solute cavity Gcav = γS.

7


(4)


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 8 of 56

Both the quantum surface and volume can be expressed in terms of the dielectric function ǫ(r). Neglecting as in other models the thermal contribution ∆Gtm to the total energy and considering that for simulations at standard pressure the term P ∆V does not play any role, the total solvation free energy can be modeled by the following expression 22

∆Gsol = ∆Gel + (α + γ)S + βV,

(5)

where α and β are solvent-specific parameters which can be fitted to experimental data such as solvation energies of neutral molecules and ions. The linear relation between the repulsion energy and the quantum surface follows the idea of traditional PCM approaches, since in PCM models Grep is proportional to the solute electronic density lying outside the dielectric cavity 34 and this last quantity can be linearly correlated to the cavity surface. The linear model for the non-electrostatic terms follows the same lines as the solvation models developed by Cramer and Truhlar, 13,35 i.e. the SMx family, where these contributions are evaluated as a sum over all the atoms, each one proportional to the solvent-accessible surface area, albeit they use more variables. Although Eq. (5) has been introduced within a chargedependent model, 22 we will use the same expression for the total solvation free energy for the soft-sphere model after a suitable generalization of the quantum surface and volume. These last non-electrostatic contributions do not change during the wavefunction optimization as they do not depend anymore on the electronic charge density (now the dielectric cavity is a function of the atomic coordinates). As a consequence only the electrostatic contribution Gel needs to be computed via an optimization of the electronic wavefunction at the ab-initio level. Appendix A outlines the main equations that have to be implemented in a DFT code for the wavefunction optimization procedure in presence of an implicit environment. The spatially varying dielectric function ǫ[ρ](r) or ǫ(r, {Ri }) has to meet several conditions: • Go monotonically from a value of 1 inside the cavity (mimicking vacuum) to the bulk 8


Page 9 of 56

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


dielectric constant ǫ0 ; • Be smooth in the transition region to guarantee a proper discretization on a three dimensional grid; • Have as many continuous derivatives as possible to avoid convergence and numerical issues during the iterative wavefunction optimization loop. The soft-sphere cavity, which satisfies all these requirements, will be presented in Sec. 2.1. Definitions of the quantum surface and volume will be given in Sec. 2.2. For a dielectric cavity where ǫ(r, {Ri }) is an explicit function of the atomic coordinates Ri , additional contributions arise to the atomic forces with respect to their gas-phase formulation. Formulas for the forces are presented in Appendix A.3. A test for the energy and force accuracy is reported on Sec. 2.3. The methodology described in this paper as well as a complementary method based on a charge density dependent continuum solvation model, developed by Andreussi et al., 22 have been implemented in the BigDFT software package.

2.1

Definition of the dielectric cavity

We define the soft-sphere cavity by means of atomic-centered interlocking spheres described by continuous and differentiable h functions which smoothely switch from 0 inside the sphere to the value 1 outside. Their product, which is still continuous and differentiable in the whole domain, reproduces the analytic cavity

ǫ(r, {Ri }) = (ǫ0 − 1)

(

Y i

h({ξ}; kr − Ri k)

)

+ 1,

(6)

where ǫ0 is the dielectric constant of the surrounding medium, { ξ} a set of parameters describing the spheres and kr − Ri k the distance from the sphere center Ri . The derivative

9



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 10 of 56

of the dielectric function with respect to Ri is ∂ǫ(r, {Ri }) ǫ(r, {Ri }) − 1 ∂h({ξ}; kr − Ri k) = ∂Ri h({ξ}; kr − Ri k) ∂Ri

(7)

for a given function h. This function can be expressed in terms of the error function kr − Ri k − ri 1 1 + erf , h(ri , ∆; kr − Ri k) = 2 ∆

(8)

where ri are the sphere radii which depend on the particular atomic species, and ∆ a parameter (fixed for all atoms) which controls the transition region (≈ 4∆ wide) where the polarization charge is located. This choice for h guarantees a “soft” dielectric function in Eq. (6) that is analytic everywhere except at the point r = Ri . With our choice of ∆ the discontinuities in the derivatives are however exponentially small and do not have any negative effect on our calculations. Following the lines of all PCM implementations, 11,12,35,36 ri are taken to be equal to the van der Waals radii rivdW of the corresponding elements scaled by a constant multiplying factor f (i.e. ri = f rivdW ). Finally, the soft-sphere cavity for a system of atoms with coordinates Ri is uniquely defined by the set of parameters {ǫ0 , rivdW , ∆, f }. The solvation model needs three additional parameters α, β and γ for the non-electrostatic terms to the total solvation energy [Eq. (5)]. The parametrization of all parameters will be reported on Sec. 3.

2.2

Geometrical contributions

Eq. (3) allows us to compute the total free energy of solvation. Thanks to Eq. (5) the integration at DFT level of the non-electrostatic terms becomes straightforward. The evaluation of the free energy and its functional derivative can be expressed in term of the quantum

10


Page 11 of 56

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


surface and volume. Introducing the function

θ[ǫ](r) =

ǫ0 − ǫ(r) , ǫ0 − 1

(9)

which takes a value of one inside the dielectric cavity and zero externally, the quantum volume V of the solute can be defined as

V [ǫ] =

Z

(10)

drθ[ǫ](r).

The integral has to be evaluated in the whole simulation domain. The associated quantum surface S can be defined from the gradient of θ[ǫ]: 21

S[ǫ] =

Z

1 dr|∇θ[ǫ](r)| = ǫ0 − 1

Z

dr|∇ǫ(r)|.

(11)

The functional derivatives of the quantum volume and surface with respect to the dielectric function are: δV [ǫ] 1 =− , δǫ(r) ǫ0 − 1 Z δS[ǫ] 1 δ∂i ǫ(r′ ) ∂i ǫ(r′ ) = dr′ δǫ(r) ǫ0 − 1 δǫ(r) |∇ǫ(r′ )| ∇2 ǫ(r) ∂i ǫ(r)∂j ǫ(r)∂i ∂j ǫ(r) 1 , − = ǫ0 − 1 |∇ǫ(r)|3 |∇ǫ(r)|

(12)

(13)

where repeated indices are considered as summed. Thanks to the continuous and differentiable permittivity ǫ(r), the function θ[ǫ](r) as well as the quantum volume and surface can be analytically computed. This holds also for their derivatives of Eq. (11-13).

11



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

2.3

Energy and force accuracy

Accurate energies and forces are fundamental, especially for molecular dynamics and structure prediction calculations. Smooth and analytical cavities can guarantee accuracy levels comparable to the ones in vacuum. Appendix A.3 shows the additional terms to the forces for the soft-sphere model. These can be analytically computed, [see Eq. (48)]. In addition, it is worth noting that all these analytic contributions to the forces can be computed nonce once during a wave function optimization. The extent of the transition region (≈ 4∆ wide) can affect the energy and force accuracy. A large ∆ can lead to a not-physical superposition of large dielectric cavities, while a small value can require a high resolution mesh which becomes computationally expensive. As a consequence a compromise is needed, considering that the more mesh points fall in the transition region, the more accurate is the description of solvent effects over the solute and the polarization charge there located. In the following we report a test of the energy and force accuracy for a water molecule in implicit water. Details of the model’s parameter setting will be given in Sec. 3. In the present test we used the parametrization reported in row 3 of Table 1. Forces can be tested by means of the test-forces tool of BigDFT, an internal program which compares the variation of the total energy ∆E with the path integral of the forces P −Fa · dra . As consequence of Eq. (46), their difference δ = ∆E + a Fa · dra should be

zero. Due to the finite discretization and the numerical integration, such value depends on the 3 dimensional spatial grid [hx , hy , hz ] of the simulation domain (where the KS problem is solved). It reaches δ ∼ 10−6 Ha in gas phase for very accurate calculations. Typical values of hi (with i = x, y, z) are in most cases between 0.3 and 0.6 bohr. Present test runs make use of an uniform grid spacing in all three directions, i.e. hgrid := hx = hy = hz . Fig. 1 shows the accuracy δ of the force computation as function of ∆. Solid lines represent the value of δ for a test-forces run of the H2 O molecule in vacuum for a given hgrid . The accuracy of the forces is comparable to that of vacuum runs, for ∆ values between 0.2 and 0.8 a.u.; 12


Page 12 of 56

Page 13 of 56

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


deviations start to take place at ∆ = 0.1 a.u., while smaller values affect the convergence of the generalized Poisson solver due to the increased sharpness of the transition region. In Fig. 1b we reported the accuracy of the solvation free energy ∆Gsol , i.e. Eq. (5), for the same set of calculations. The accuracy is defined as the difference between values extracted for a given hgrid and the most accurate grid we tested at hgrid = 0.2 a.u.. Varying ∆ the same accuracy is recovered for hgrid ≥ 0.4 a.u.. An hgrid = 0.4 a.u. guarantees an accuracy of 10−2 kcal/mol. Similar results obtained with other small organic molecules lead to similar conclusions. Performance of the soft-sphere model is discussed in Appendix B for the same system.

3

Parametrization

The soft-sphere solvation model is defined by the parameter set {ζ} = {ǫ0 , rivdW , ∆, f ; α, β, γ}, i.e. six parameters plus an atomic radius rivdW for each distinct atomic species of the quantum system in contact with the implicit solvent. Only the first four parameters are related to the dielectric cavity and the electrostatic contribution ∆Gel . To parametrize the solvation model we need to map the {ζ} = {ǫ0 , rivdW , ∆, f ; α, β, γ} parameter space, and minimize the mean absolute error (MAE) over a given collection of experimental free solvation energies. In the following we will report the parametrization for a water environment, both for neutral molecules and ions. In addition, we have parametrized two non-aqueous solvents, namely mesitylene and ethanol. Structures and experimental data have been retrieved from the Minnesota Solvation Database, version 2012. 28,37–39 Solutes used as benchmark contain at most the following elements: H, C, N, O, F, Si, P, S, Cl, Br, and I. Soft norm-conserving pseudopotentials including a non-linear core correction 40,41 along with PBE functional were used to describe the core electrons and exchange-correlation for all calculations of the present section. The Libxc 42 library was exploited for the calculation of the functionals. Solvation free energies were computed via Eq. (5) using atomic structures

13



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 1: H2 O molecule in implicit water described by the soft-sphere model: (a) force accuracy δ obtained by means of the test-forces tool of BigDFT (see text) as function of the cavity parameter ∆ for different grid spacing hgrid (colored solid lines represent the accuracy in vacuum runs for each hgrid ); (b) accuracy of the solvation free energy ∆Gsol with respect to the spacial grid hgrid for different ∆.

14


Page 14 of 56

Page 15 of 56

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


optimized both in vacuum and in the presence of the implicit solvent. BigDFT represents the wavefunction in a wavelet basis set. These basis functions are centered on a cartesian mesh of resolution [hx , hy , hz ]. 25 The generalized Poisson equation, i.e. Eq. (1), has been solved via the recently developed generalized Poisson solver based on Interpolating Scaling Functions. 9 Its preconditioner, based on the BigDFT vacuum Poisson solver, allows to solve the electrostatic problem in some ten iterations with an exact treatment of free, surface and periodic boundary conditions. 27,43 A new version of the algorithm which includes a proper use of an input guess is reported in Appendix B. All simulations reported in the present study use an uniform grid of hx = hy = hz = 0.40 a.u. and free boundary conditions.

3.1

Neutral molecules

The soft-sphere model has been benchmarked on a collection of 274 aqueous solvation free energies to reproduce a water environment. At first, to speed up the exploration of the parameter space, we trained the model on a representative subset, which consists of 13 molecules spanning the main functional groups of the entire set. 44 3.1.1

Cavity parameters

In order to reduce the size of the parameter space (in principle a MAE should be extracted for each {ζ} setup), we proceeded fixing the largest possible number of parameters. Since we are benchmarking the model on aqueous solvation free energies, the experimental value at low frequency and ambient conditions of ǫ0 = 78.36 has been used. A critical point in PCM approaches is the choice of shape and size of the cavity, which should reflect the charge distribution of the solute. In the first formulation of this model, 45 PCM atomic radii ri were chosen to be proportional to the van der Waals radii. 46,47 Importing a consistent set of van der Waals radii rivdW which spans the whole periodic table is more meaningful than fitting radii for all the elements which would presumably lead to some overfitting. We decided to import rivdW from the Unified Force Field (UFF) parameterization 48 for all elements of the 15



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

periodic table. Hydrogen atoms have been considered explicit in the present work. ∆ has been fixed to 0.5 a.u. for all atoms, generating a transition region of ≈ 2 a.u.. Different choices of ∆ will be also investigated and discussed at the end of the paragraph. Finally, only the multiplying factor f for the atomic radii are left free to vary, reducing the four-parameter problem for the cavity definition to a single-parameter optimization. The flexibility to locally change atomic radii has also another advantage. Considering a wet surface, the radii for atoms which are in the bulk region may be further enlarged to ensure that the dielectric permittivity is exactly 1 for all the domain which is inaccessible by the solvent. This condition is not always satisfied by an implicit solvation model and has to be checked. Since we imported the non-electrostatic contributions from the sccs approach, we started the present parametrization from the same non-electrostatic proportional factors (α + γ = 50.0 dyn/cm, β = −0.35 GPa). Finally, we ended up with a single-parameter optimization for f to define the model. A value of f = 1.12 minimizes the prediction error of the experimental aqueous solvation energies over the representative set of 13 molecules, giving a MAE of 1.39 kcal/mol. Exporting the same {ζ} setup (with the optimized f ) to the whole set of 274 neutral molecules, we obtained a MAE of 1.10 kcal/mol. This demonstrates that the soft-sphere model gives about the same accuracy as the charge density dependent sccs model (MAE of 1.14 kcal/mol on the same set of 274 neutrals, see Table 1) and that it is possible to export the linear approach for the non-electrostatic terms (implemented and parametrized for the sccs model). A careful analysis of the error distribution over all functional groups, allowed us to fine tuning some van der Waals atomic radii rivdW . We found that a slight variation of the Nitrogen radius from the UFF 48 value 1.83 Å to the one of Bondi 46 1.55 Å improves the whole parametrization. This trend reflects the interaction of the Nitrogen lone pair with Hydrogen atoms of water, making their reciprocal distances smaller. The parametrization curve varying f over the 13 molecules has been reported in Fig. (2) (black dots, full line).

16


Page 16 of 56

Page 17 of 56

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


Table 1:

Mean Absolute Errors in aqueous solvation free energies (kcal/mol) for a set of 274 neutral molecules obtained with the soft-sphere solvation model. atomic radii rivdW

cavity multiplying factor f

α + γ [dyn/cm]

β [GPa]

MAE [kcal/mol]

UFF 48

1.12

50.0

-0.35

1.10

UFFa

1.16

11.5

0.00

1.12

UFFa

1.12

50.0

-0.35

1.03

Pauling 47

1.36

50.0

-0.35

1.17

Bondi 46

1.32

50.0

-0.35

1.06

a

Van der Waals atomic radii from the Unified Force Field setup 48 with Bondi’s radius 46 for Nitrogen.

The MAE minimum is located again at f = 1.12 with a MAE of 1.15 kcal/mol for the 13 molecules. Applying the new setup to the full set of 274 molecules gives a MAE of 1.03 kcal/mol (see Table 1). This remarkable result has been achieved with the crude import of the linear non-electrostatic model from sccs. The inset of Fig. 2 reports a global comparison between experimental and ab-initio aqueous solvation free energies for the full set of 274 neutrals. All dots stay close to the diagonal, reflecting the accuracy of the soft-sphere model over all functional groups. PCM formulations in the literature mainly use van der Waals radii from Pauling 47 or Bondi 46 ’s work. Even though these sets work well for organic molecules, values for the full periodic table are not available, making the UFF setup more tempting for general ab-initio calculations. Nevertheles we explored the performance of our solvation model with the van der Wall radii of Pauling and Bondi’s. Following the same procediure as in our the previous optimization, we first performed an optimization of the cavity multiplying factor f over the representative set of 13 molecules, and only afterwards for the entire {ζ} setup of 274 neutral molecules. Both f and MAEs for the full set are reported in Table 1. The optimized f for both sets is larger than the one obtained with the UFF implementation. This is related to Hydrogen radius, being 1.20 Å in the Pauling and Bondi’s work, and 1.443 Å in the UFF. The overall accuracy does not improve with respect to the UFF parametrization, where MAEs of 1.17 and 1.06 kcal/mol have been found for Pauling and Bondi’s radii, respectively. 17



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 2: Mean Absolute Error in aqueous solvation free energies for the representative set of 13 neutral molecules as function of the cavity multiplying factor f (UFF radii rivdW 48 with Bondi’s radius 46 for Nitrogen, α + γ = 50.0 dyn/cm and β = −0.35 GPa). The inset reports experimental and ab-initio aqueous solvation free energies for the full set of 274 neutral molecules with the same parameters setup (optimized f = 1.12).

So far, all parameters related to the dielectric cavity have been optimized, except ∆. A window of 0.5 < ∆ < 0.65 does not change the results in term of accuracy (increasing ∆ request a proportional increase of f to recover the same MAE). 3.1.2

Non-electrostatic parameters

The last optimization step treats the non-electrostatic parameters α, β and γ. Assuming a surface tension of γ = 72 dyn/cm for water at room temperature, only α and β need to be optimized. The sum α + γ can be considered as a single parameter since both variables are multiplying factors of the quantum surface S in Eq. (5). We extrapolated all the β and α +γ couples which minimize the MAE over the representative set of 13 molecules (determining also the optimized f for each fixed couple). As noticed in Ref. [ 22], we found an almost perfectly linear behavior between the error, i.e. the MAE, as function of β and α + γ. This was probably due to the small variation in size of the molecules considered in the training 18


Page 18 of 56

Page 19 of 56

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


set, which makes the volume and the surface terms to be linearly related. All MAE minima lie approximately on the same line β = c(α + γ) + d, where c = −9.1 nm−1 and d = 0.1045 GPa. The best parametrization still correspond to the sccs optimal setup with α + γ = 50.0 dyn/cm and β = −0.35 GPa. However, given the strong correlation between the two tunable parameters, it is important to stress that the optimal parametrization reported above may be strongly linked to the type of systems (small isolated molecules) under study. Whether these parameters could be used in systems with very different sizes, e.g. large molecular and supramolecular systems or solids described using a slab geometry, is not obvious. For this reason, unless a more refined optimization that takes into account experimental data on large systems is performed, a parametrization only relying on the surface term, thus putting β = 0 GPa, should be more transferable and preferable for the above applications. For surface calculations, β = 0 GPa is mandatory to avoid an unphysical dependence of the energy on the slab volume. 3.1.3

Results

The final MAEs for both surface (row 2) and small cluster (row 3) parametrizations over the full set of 274 experimental solvation energies are 1.12 and 1.03 kcal/mol, respectively (see Table 1). Aqueous solvents are widely used and of great relevance in wet chemistry. Investigations may require a mixed explicit/implicit treatment of a water environment. The soft-sphere model predicts a solvation free energy for H2 O and (H2 O)2 in implicit water of −5.74 and −10.19 kcal/mol for the surface parametrization, and −5.60 and −10.34 kcal/mol for the small clusters one. Experimental values are −6.31 and −11.27 kcal/mol, respectively. 28 Predictions are within 1.1 kcal/mol and both parametrizations perform pretty well. Marenich et al. in Ref. [ 38] explored performances of several implicit solvation models. These approaches hold as common feature the use of dozen of parameters to model all contributions to the solvation free energy. Among others, they tested the SM8 38 model of the SMx family, the Integral Equation Formalism Polarizable Continuum Model 45 of Gaussian

19



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

03 (IEF-PCM), 49,50 the Conductor-like PCM model 51 in GAMESS 52 (C-PCM/GAMESS), Jaguar’s Poisson-Boltzmann self-consistent reaction field solver 53,54 (PB/Jaguar) and the Generalized Conductor-like Screening Model in NWChem 55 (GCOSMO/NWChem). MAEs of these methods for aqueous solvation free energies are been reported in Table 2. for the set of 274 neutrals Data from the recent SM12 implementation 39 have been also included. In addition we tested the charge-dependent sccs model implemented within the BigDFT package. Using only two parameters (surface parametrization), the soft-sphere model applied to the same dataset lies in the same range of accuracy. A carefully tuning of atomic radii and/or a functional dependence of delta over the periodic table could further lower the MAE and improve the whole accuracy. Table 2: Mean Absolute Errors in aqueous solvation free energies (kcal/mol) for several solvation models (MAEs from Ref. [ 38]). Model benchmarks refer to same set of 274 neutrals, 60 anions and 52 cations of the Minnesota Solvation Database, version 2012. 28 Method

neutrals

cations

anions

soft-spherea

1.12

2.13

2.96

sccs 23

1.14b

2.27c

5.54c

SM8 38

0.55

2.70

3.70

SM12 39

0.59

2.90

2.90

PB/Jaguar 38

0.86

3.10

4.80

IEF-PCM 38

1.18

3.70

5.50

C-PCM/GAMESS 38

1.57

7.70

8.90

8.17

11.00

7.00

GCOSMO/NWChem 38 a

3.2

b

parametrization of row 2 Table 1; sccs implemented in BigDFT; c sccs for ions corresponds to a reduced set of 55 anions and 51 cations 23 of the same Minnesota data set.

Anions and cations

In order to strictly validate the soft-sphere model, we benchmarked it on a set of aqueous free energies of solvation for 112 singly-charged ions (60 anions and 52 cations). The Minnesota Solvation Database, version 2012, 28,37,39 contains two sets of them, i.e. the unclustered and 20


Page 20 of 56

Page 21 of 56

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 selectively clustered data set. The last is identical to the first, except for 31 ions which are clustered with an explicit water molecule. More details can be found in the article discussing the SM6 parametrization. 37 As done in the last articles of the SMx family, 38,39 and in order to make comparisons with others solvation models, we benchmarked the soft-sphere model on the selectively clustered data set. The two main model parametrizations have been explored in row 2 and 3 of Table 1, except for the radii multiplying factor f which was left free to be optimized. Due to the nonlinear effects in the polarization of the surrounding dielectric with high fields, the factor f , or alternatively all the radii of atoms bearing the ionic charge, need to be reconsidered. 11,12 Furthermore, ab-initio gas-phase optimizations and Monte Carlo simulations of systems in aqueous solution suggested that solvent molecules of the first solvation shell come closer to a charged molecule than to a neutral. 56 Although most PCM models implement a charge dependency on the atomic radii, 11,57 we decided to reproduce experimental solvation energies optimizing the global cavity factor f both for anions and cations. A similar approach has been also investigated by others authors. 56 This procedure is less accurate than modifying single atomic radii, especially when the solute molecule sizes become larger or for extended interfaces. However, it is still a good approximation for small size ions. A more detailed approach with a fine tuning of the individual van der Waals radii and their dependence on the local atomic charges is beyond the scope of the present paper. Table 3 reports MAEs both for anions and cations. Fixing the non-electrostatic parameters to α + γ = 50.0 dyn/cm and β = −0.35 GPa, the optimization procedure on f suggested an f = 0.98 as optimal multiplying factor for anions and f = 1.10 for cations, with MAEs of 3.05 and 2.07 kcal/mol, respectively. These results are in perfect agreement with similar findings obtained in the parametrization of self-consistent cavities. 23 Fig. 3 reports experimental and ab-initio aqueous solvation free energies both for anions (empty blue circles) and cations (full black circles). Concerning cations, dots belonging to the diagonal are mainly cations containing Nitro-

21



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

gen. On the other side, the model tends to underestimate free solvation energies for four different types of alcohols, two types of ketons and the water cation. This trend suggests that smaller cavities are needed for these cations due to the stronger electrostatic interaction between the solvent and the solute. Excluding these last seven solutes, the overall MAE for the remaining 45 cations decreases from 2.07 to 1.33 kcal/mol. A factor of f = 1.04 allows to improve the prediction for the 7 units (alcohols, ketons and the water cation) where their MAE is 2.07 kcal/mol. This fine analysis has been conducted to identify margins for improvements. Similar arguments are valid for anions, but a systematic study of functional groups is less straightforward since oxygen comes into play in the majority of the units. An optimization of f keeping fixed α + γ = 11.5 dyn/cm and β = 0 GPa, leads to a similar accuracy for the prediction of solvation free energies (see Table 3). MAEs are 2.96 and 2.13 kcal/mol for anions and cations, respectively. We see that the soft-sphere model performs very well for ions in water, considering that the estimated experimental average uncertainty for solvation free energies of ionic solutes is 3 kcal/mol. 37 It is worth noting the high accuracy has been reached, by only optimizing one parameter, i.e. f , with respect to the parametrization for neutral solutes. Table 2 reports MAEs for the same set of ions in water obtained with others solvation models. All data refer to the same selectively clustered data set of 60 anions and 52 cations of the Minnesota Solvation Database, version 2012, except sccs which corresponds to a reduced set of 55 anions and 51 cations 23 of the same database. The soft-sphere model outperforms all current implicit solvation models, with a MAE of 2.57 kcal/mol over the full ion dataset. The possibility to use the same α, β and γ and to change only the radii upon ionization allows to tackle systems which contains simultaneously neutral and charged complexes. This is not possible in charge density dependent approaches which can not distinguish between atoms and need different non-electrostatic parametrizations (and then different Hamiltonians) for neutrals and ions. Since both parametrizations of Table 1, i.e. row 2 and 3, are able to accurately reproduce

22


Page 22 of 56

Page 23 of 56

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


solvation energies of neutrals and ions, and taking into account that surface calculations require β = 0 GPa and that reducing the parameter space is always advantageous, we can just use the quantum surface S in the model of non-electrostatic contributions in Eq. (5). The whole solvation model is thus uniquely defined by two parameters, i.e. f and α.

Figure 3: Experimental and ab-initio aqueous solvation free energies for a set of 112 single-charged ions, of which 60 anions (empty blue circles) and 52 cations (full black circles), obtained with the soft-sphere model (UFF radii rivdW 48 with Bondi’s radius 46 for Nitrogen, α + γ = 50.0 dyn/cm and β = −0.35 GPa, f = 0.98 for anions and 1.10 for cations). Table 3: Mean Absolute Errors in aqueous solvation free energies (kcal/mol) for a set of 112 ions (60 anions and 52 cations) obtained with the soft-sphere modela . cavity multiplying factor f

α + γ [dyn/cm]

β [GPa]

MAE [kcal/mol]

Anions

0.98 1.00

50.0 11.5

-0.35 0.00

3.05 2.96

Alcohols and ketons cations

1.04 1.04

50.0 11.5

-0.35 0.00

2.07 1.88

Cations without alcohols and ketons

1.10 1.11

50.0 11.5

-0.35 0.00

1.33 1.46

All Cations

1.10 1.10

50.0 11.5

-0.35 0.00

2.07 2.13

-

50.0 11.5

-0.35 0.00

2.60 2.57

Total ions UFF radii rivdW

48

with Bondi’s radius for Nitrogen.

23


a


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

3.3

Non-aqueous solvents

The study of chemical reactions within non-aqueous solvents is as also important. Therefore it is worth investigating how the parametrization has to be modified for a non-aqueous solvent. As test cases, we took two solvents, namely mesitylene and ethanol, which have low dielectric constants ǫ0 of 2.26 58 and 24.85, respectively at 20◦ C. To optimize the soft-sphere parameters {ζ} = {ǫ0 , rivdW , ∆, f ; α, β, γ}, experimental solvation free energies for a set of seven (mesitylene) and eight (ethanol) organic molecules have been taken from the Minnesota Solvation Database. 28 Changing only ǫ0 to the experimental value of the non-aqueous solvent and taking all others parameters from the water parametrization (row 2 or 3 of Table 1) does not reproduce the experimental data, underestimating all solvation free energies with an overall MAE of 3.42 and 2.84 kcal/mol for mesitylene and ethanol, respectively. Turning off the non-electrostatic contributions (α + γ = 0 dyn/cm, β = 0 GPa), and using the proper ǫ0 has the same effect with MAEs of 3.42 and 1.74 kcal/mol. As a consequence each solvent needs its own parametrization. PCM approaches usually have solvent-dependent parameters, like the multiplying factor for the atomic radii or several collections of parameters for the non-electrostatic energetic. We therefore set in addition the values rivdW for the radii to the UFF values 48 with the exception of Nitrogen where we use the Bondi’s radius and used again ∆ = 0.5 as we did for water. The surface tension at 20◦ C of mesitylene is γ = 28.80 dyn/cm, 59 whilst for ethanol is γ = 22.10 dyn/cm. Since solvation energies are well reproduced with only the surface term in Eq. (5), we end up with a two-parameters optimization, i.e. f and α. First we optimized α, keeping f to the water value of f = 1.16. This provides the major correction to the overall MAE. Then f and α have been fine tuned concurrently. A value of f = 1.22 and α + γ = −12.0 dyn/cm provide a MAE of 0.71 kcal/mol for mesitylene. Doing the same for ethanol, a f = 1.22 and α + γ = −4.0 dyn/cm allow to reproduce solvation free energies with a MAE of 1.28 kcal/mol. Accuracies are comparable to the SM8 model, which gives a MAE of 0.40 (mesitylene) and 1.53 (ethanol) kcal/mol for the same data set of neutrals. 38 24


Page 24 of 56

Page 25 of 56

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


Hence our approach has an accuracy comparable to similar state-of-the-art methods, despite the reduced number of parameters. Nonetheless, it is important to stress that the reported parameterizations may be affected by the reduced size of the fitting set: the limited numbers and the narrow ranges of experimental solvation free energies used for the fit may give rise to a parameterization which is not transferable to a broader set of compounds. This is particularly striking in a model which only relies on two tunable parameters, for which the optimized non-electrostatic terms end up being always negative. Thus, similarly to the electrostatic term, these terms contribute to a stabilization of the solute, excluding the possibility of modeling insoluble compounds with positive solvation free energies. On the other hand, the optimization procedure gave a prefactor for the atomic radii f greater than the corresponding value in water. This trend is in agreement with the larger sizes of mesitylene and ethanol molecules with respect to H2 O, and meets previous studies of solvated molecules where the cavity is defined as a convolution of both solute and solvent electronic densities. 60 Including an optimization on β or additional dependences for the nonelectrostatic terms on specific solvent descriptors like the refractive index and the Abraham’s hydrogen bond acidity and basicity parameters, 38 can further improve the prediction power on non-aqueous solvents and make the soft-sphere model universal, where “universal” means applicable to all solvents.

4

Solid-liquid interfaces and contact angle

As a further benchmark system for implicit solvation models we consider the contact angle θC that a liquid drop on a surface forms with the surface plane. It is experimentally accessible quantity and a collection on several substrates can be used as test set. Sessile-drop and captive-bubble techniques are common experimental reference methods. 61 This benchmarking approach is alos useful for the parametrization of elements of the periodic table that are not represented in standard databases of experimental solvation energies.

25



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 56

Although the idea of a direct comparison of predicted and experimental contact angles seems straightforward, its is not technically challenging. Clean surfaces do not exist in nature where reconstructions, adsorbed molecules or radicals, and defects usually take place. An implicit solvation approach gives results for ideal surfaces whereas experimental measurements may be influenced by all these deviations from the ideal behaviour. Other simulation techniques could be applied where both the surface and the liquid are treated explicitly (i.e. molecular dynamics, Monte Carlo methods or others 62 ). The proposed benchmark provides necessary conditions that an implicit scheme has to fulfill if applied to a solid-liquid interface. Simple conditions have to be verified like partial or complete wetting, and the degree of the solid-liquid interaction by approximately comparing θC . The contact angle quantifies the wettability of a surface and reflects the equilibrium between liquid (drop), solid (surface) and vapor-phase interactions. Looking to the contact line at the surface and imposing the thermodynamic equilibrium between surface tensions of the three phases, a relation can be deduced for the contact angle in terms of the interfacial energies solid-vapor γSG , solid-liquid γSL , and liquid-vapor γLG (i.e. the surface tension)

cos θC =

γSG − γSL . γLG

(14)

This relation is known as the Young’s equation. 63,64 Since the three surface energies at equilibrium form the side of a triangle, partial wetting occurs only when the triangle inequalities γij < γjk + γik is satisfied. θC can easily be related to the work of adhesion of the solid and liquid phases in contact, i.e. WLS = γSG + γLG − γSL , leading to the Young-Dupré equation 65 WLS = γLG (cos θC + 1). The work of adhesion is the work which must be done to separate two adjacent materials at their phase boundary (liquid-liquid or solid-liquid) from one another. Conversely, it is the energy which is released in the process of wetting and is a measure of the strength of the contact between two phases. γSG and γLG are the energies required to create a new surface

26


Page 27 of 56

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


whilst γSL represents the work that has to be done in order to form the interface. These quantities are difficult to determine experimentally and only θC and γLG can be directly determined. Contact angle measurements can give access to γSG and γLG indirectly via the Young equation. In this case γSL is assumed to be a function of γSG and γLG . Several functional relations are present in the literature, 66–69 taking into account the polar or nonpolar character of the solvent as well as its acid or basic nature. 70 In order to classify hydrophobic and hydrophilic surfaces, van der Oss exploited the concept of free energy of solvation ∆GSL both for molecules and condensed matter systems. 71 It can be related to the work of adhesion being ∆GSL = −WLS Another parameter, especially useful to characterize hydrophilic surfaces, is the spreading coefficient S defined as S = γSG − γLG − γSL . It represents the work performed to spread a liquid over a unit surface area of a clean and non-reactive solid (or another liquid) at constant temperature and pressure and in equilibrium with liquid vapor. Coupled with the Young’s equation it leads to S = γLG (cos θC − 1), which means that a liquid drop partial wets the surface only for S < 0, while a total wetting occurs for S > 0. Zero wetting takes place when WLS < 0. These relations can be made explicit by rewriting cos θC as

cos θC =

S WLS −1= + 1. γLG γLG

(15)

It is also clear that WLS > γLG implies that θC < 90◦ , i.e. the surface has an hydrophilic character, whereas it is hydrophobic in the opposite case. Furthermore, both WLS and S should be similar in magnitude to γLG to achieve partial-wetting conditions (i.e. 0 ≤ WLS /γLG ≤ 2, −2 ≤ S/γLG ≤ 0). Once θC is experimentally determined, these conditions need to be fulfilled, providing useful fixed points to benchmark a solvation model. For a slab, the DFT surface energy γSG is calculated according to

γSG =

1 N − N Ebulk lim Eslab-SG 2A N →+∞ 27


(16)


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 56

N where N is the number of atoms within the slab, A is the slab’s cross section area, Eslab-SG

is the total energy of the N -atoms relaxed slab, and Ebulk the energy of the bulk system per atom. The solid-liquid interface energy γSL can be easily computed by means of a similar equation γSL =

1 N − N Ebulk , lim Eslab-SL 2A N →+∞

(17)

N considering Eslab-SL as the total energy of the slab in contact with the implicit solvent. This

approach gives direct access to γSL and it is not necessary to specify an analytical functional dependence of γSL in terms of γSG and γLG in order to estimate the contact angle via Eq. (14). We determined spreading the coefficient S, the work of adhesion WLS and the contact angles θC for a collection of surfaces in contact with water. Table 4 reports computed data for a silver (001) surface, a cleaved and reconstructed (001) surface of SiO2 α-quartz 72 and two carbon-based surfaces, i.e. diamond (001) and graphene. Thanks to the implemented UFF tabulation of van der Waals radii, our soft-sphere model is able to handle such materials. We used the parametrization reported in the second row of Table 1, valid for surface calculations (β = 0 GPa). It is worth noting that additional benchmark procedures have to be carried out for the cavity parameters {ζ} = {ǫ0 , rivdW , ∆, f ; α, β, γ} to guarantee a consistent physical description of the solid-liquid interface. In particular van der Waals radii rivdW for elements not present in the regular set of organic molecules have to be monitored. Small values can lead to no-physical dissolution of surface atoms. The computation of Ebulk is not necessary because θC , WLS and S are all functions of the difference γSG − γSL , and Ebulk of Eq. (16,17) is cleared. From this side, cos θC can be N N seen as the solvation energy per unit area ∆Gsol slab = (Eslab-SL − Eslab-SG )/2A divided by γLG ,

namely cos θC = −

Gsol slab . γLG

(18)

For water at room temperature γLG has the value 72 dyn/cm. The parameters for the DFT

28


Page 29 of 56

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


calculations (grid spacings, Monkhorst-Pack k -point mesh) have been chosen in BigDFT to reach a surface energy accuracy of 0.01 J/m2 . All systems have been relaxed both in vacuum and in the presence of the implicit aqueous environment until all the forces were less than 5 meV/Å. BigDFT allows to use exact surface boundary conditions, avoiding spurious interactions in the direction orthogonal to the surface. As an additional check, we compared vacuum surface energies [Eq. (16)] with DFT literature data, recovering a good agreement for all systems. Table 4: Simulated work of adhesion WLS , spreading coefficient S and contact angle θC for several surfaces in contact with implicit water described by the soft-sphere model. slab layers

WLS

S

θC

[mJ/m2 ]

[mJ/m2 ]

[degree]

silver (001)

8

479

33

-

SiO2 α-quartz cleaved (001)

18

181

35

-

SiO2 α-quartz reconstructed (001)

27

76

-69

87

Diamond (001)

12

69

-76

92

Graphene

1

62

-82

97

Unless an oxide layer is present, clean metal surfaces are hydrophilic meaning that water completely or nearly completely wets them (gold, silver, copper and others), 61,73 The positive value of S for clean (001) silver surface (see Table 4) predicted by the soft-sphere model agrees with its hydrophilic character, allowing for a complete wetting of its surface. A similar validation can be extended to other metallic materials and represents a preliminary checkpoint if implicit solvation approaches need to be integrated in such solid-liquid investigations. Furthermore, as a consequence of their bulk bonding strength, surfaces can also be divided into high- and low-energy surfaces, The stronger the chemical bonds are, i.e. the higher the surface energy is, the more easily complete wetting is reached. In our set of surfaces, the cleaved SiO2 α-quartz surface falls in this last class. Among the different SiO2 highly-reactive surfaces, the (001) was chosen because it is the most stable. 72 Applying Eq. (16) we recovered 29



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

a surface energy of γSG = 2.22 J/m2 . This value is consistent to the one reported in Ref. [ 72], from where the cleaved and reconstructed SiO2 structures come from. Implicit solvation correctly predicts its hydrophilic character with a spreading coefficient S > 0 as reported in table 4. Conversely, the reconstructed SiO2 α-quartz (001) has a γSG o 0.39 J/m2 . Because of its negative spreading coefficient, implicit solvation predicts a partial wetting for this surface with a contact angle of θC = 87◦ and a work of adhesion of WLS = 76 mJ/m2 . Therefore the implicit procedure is able to predict different θC for high and low energy surfaces of the same material. A direct comparison with experimental contact angles is difficult since real surfaces undergo reconstruction as well as hydroxylation, contamination and patterning concurrently. 74 Similar considerations hold for the clean (001) diamond surface. Molecular dynamics and Monte Carlo simulations predicted an hydrophobic character for an hydrogen-terminated diamond slab. 62,75 In our case a work of adhesion of 69 mJ/m2 has been obtained, together with a θC = 92◦ degree. Remaining within the carbon-material class, we applied the implicit solvation scheme to a graphene sheet. Graphene is a two dimensional arrangement of carbon atoms, bonded covalently to form an honeycomb lattice due to the sp2 orbital hybridization. Weak van der Waals interactions take place with water molecules once graphene is in contact with an aqueous solution, allowing for a partial wetting of the surface. Several experimental and simulation studies investigated the wetting features of such carbon allotropes. For graphene the measured a contact angle is 127◦ , 76 although a theoretical study demonstrated a discrepency with the contact angle of graphite, 77,78 suggesting a θC of the order of 95-100◦ . 79 For the partial wetting behavior, data extracted with our soft-sphere solvation scheme is in agreement with these experimental and theoretical studies, predicting a negative spreading coefficient, a contact angle of 97◦ and a work of adhesion of 62 mJ/m2 . A lowering of the cavity prefactor to f = 1.12 does not affect the whole wetting predictions for all studied solid-liquid interfaces in terms of the hydrophobic/hydrophilic character.

30


Page 30 of 56

Page 31 of 56

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


5

Structure of CdS (11¯20) surface in electrochemical media

As an application, we investigated the structure of the CdS (11¯20) surface in contact with water at different pH and bias potentials, mimicking electrochemical conditions. CdS is among the few promising photocatalytic materials for the water splitting reaction, due to its appropriate band gap and band positions. 80–82 Recently, it was reported that nanorods and other nanoporous materials based on CdS show very high photon-to-hydrogen conversion efficiency when used together with hole-scavenger molecules. 83–85 It was observed experimentally that a high pH is necessary for the efficient photo-catalytic activity of CdS. 84 However, photocorrosion was seen after prolonged photo-exposure at such pH levels. 84,85 Understanding the structure and reactivity of the CdS surfaces at electrochemical conditions can provide valuable informations on improving the CdS-based photocatalytic systems. We employed here the ab-initio thermodynamics based “computational hydrogen electrode” (CHE) approach 86,87 to study the structure of the most exposed CdS (11¯20) surface at different pH and bias potentials. Solvent effects have been included using the implicit soft-sphere solvation model. Most of the CHE computations in literature were carried out under in vacuo conditions, 87–92 and a qualitative understanding of solvent effects on the results of such computations remain unclear. 93–96 Recently the CHE approach was used to study the structure of CdS surfaces in contact with vacuum. 97 Here we have computed the formation free energies of various surface terminations (O, OH, OOH), varying pH and bias potentials. The following reactions were considered for the formation of these terminations: 2H2 O∗ → 2O∗ + 4H+ + 4e− ,

(19)

2H2 O∗ → 2HO∗ + 2H+ + 2e− ,

(20)

31



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

2H2 O∗ + 2H2 O → 2HOO∗ + 6H+ + 6e− .

Page 32 of 56

(21)

The above equations include two moles of H2 O∗ as there are two adsorption sites in the surface model that we have considered. Here O∗ , HO∗ , HOO∗ and H2 O∗ represent O, HO, HOO and H2 O terminated surfaces, respectively. Following Ref. [ 86], the free energy change ∆G of the above reactions with reference to a standard hydrogen electrode (SHE) at pH2 = 1 bar and T = 298 K as a function of bias potential U and pH can be calculated, respectively, as ZPE ∆G2O∗ = ∆E2O∗ + ∆E2O ∗ − T ∆S2O∗

+4∆G(pH) − 4eU,

(22)

ZPE ∆G2HO∗ = ∆E2HO∗ + ∆E2HO ∗ − T ∆S2HO∗

+2∆G(pH) − 2eU,

(23)

ZPE ∆G2HOO∗ = ∆E2HOO∗ + ∆E2HOO ∗ − T ∆S2HOO∗

+6∆G(pH) − 6eU.

(24)

Here ∆E ZPE and T ∆S are the zero-point energy correction and the vibrational entropic contribution to the free energy, respectively, and ∆G(pH) = −kB T pH ln 10. ∆E2O∗ , ∆E2HO∗ and ∆E2HOO∗ are calculated from DFT calculations as follows:

∆E2O∗ = E2O∗ + 2EH2 − E2H2 O∗ ,

(25)

∆E2HO∗ = E2HO∗ + EH2 − E2H2 O∗ ,

(26)

∆E2HOO∗ = E2HOO∗ + 3EH2 − E2H2 O∗ − 2EH2 O .

(27)

32


Page 33 of 56

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 energies of the surface-adsorbed structures were calculated with and without the implicit solvation model. Preliminary calculations showed that a slab built up from an 1×1 supercell with 8 atomic layers, i.e. Cd16 S16 , and a k-point sampling 5×1×5 adequately converged the energy (see the Supporting Information). Zero-point energy and vibrational entropic contributions were computed from the harmonic frequencies, extracted by means of a finite difference method. For these last two contributions, a 2×1×2 k-points mesh guarantees reasonable convergence. Further, we found that both ∆E ZPE and T ∆S were nearly the same with and without continuum solvent (see the Supporting Information). Thus, we have not included solvent effects while computing these quantities. All computations were performed with spin-polarized DFT in a wavelet basis uising the PBE 98 exchange-correlation functional as implemented in the BigDFT program. 25,26,42 Surface boundary conditions were employed for solving the Poisson equations. 43 The real space grids along all the three directions were set to 0.45 Bohr. Soft norm-conserving pseudopotentials including non-linear core correction 40,41 were used to describe the core electrons. For the dielectric cavity, a van der Waals radius for Cd of 1.58 46 Å has been used, since the UFF value 1.424 Å leads to the surface dissolution. The implicit model predicts a spreading coefficient S of 28 mJ/m2 for the CdS (11¯20) surface in contact with water, confirming its hydrophilic character (S > 0). Different contributions to ∆G at standard conditions (with reference to SHE) are given in Table 5. The zero-point energy and the entropy of H2 (g) and H2 O(g) were taken from an experimental database. 99 GH2 O(l) was computed as GH2 O(g) (298 K, 1 bar)−∆Gv (298 K, 1 bar), 100 where ∆Gv is the free energy of vaporization, equal to 0.09 eV computed based on an experimental database 99 (see also the Supporting Information). We have considered adsorption at both surface Cd and S sites for all the adsorbates. It was found that except for oxygen, all the other adsorbates prefer to bind on Cd sites, consistent with previous works. 97 In Fig. 4 we show the free energy change ∆G for different surface terminations as a function of bias potential U at pH = 0 and pH = 14. We report results both for pure gas-

33



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 34 of 56

phase conditions (solid lines) and in the presence of an implicit water environment (dotted lines). For small bias potentials, the most stable surface is H2 O terminated, where molecular water is adsorbed on surface Cd atoms at both pH values. The HOO terminated surface is stable above 3.12 V (2.53 V) bias potential with (without) implicit water at pH = 0, whilst in the range from 0.93 V (1.24 V) to 3.12 V (2.53 V) the O terminated surface is thermodynamically preferred under these conditions. Similar conclusions hold for pH = 14. The inclusion of the solvent interaction results in a substantial shift of the relative stability with respect to the bias potentials. For e.g., the threshold bias for observing HOO terminated surface has increased by nearly 0.6 V by including the implicit solvent in the calculations. Table 5: Contributions to the free energy change ∆G computed at 298 K and 1 bar with reference to the standard hydrogen electrode for O, HO and HOO terminated CdS (11¯20) surfaces. ∆E Solvent and ∆E are the change in energy calculated with and without implicit solvent, respectively [Eqs. (25-27)].

6

∆E [eV]

∆E Solvent [eV]

∆E ZPE [eV]

T ∆S [eV]

2O∗ (S)

6.39

5.14

-0.63

0.79

2HO∗ (Cd)

5.38

4.73

-0.40

0.56

2HOO∗ (Cd)

10.96

10.87

-0.82

0.10

Conclusions

A soft-sphere solvation model for ab-initio electronic-structure calculations has been developed, parametrized and tested. The transition from the dielectric cavity to outside is continuous and differentiable, allowing for the analytical calculation of the additional terms to the forces as well as the cavity-dependent non-electrostatic contributions to the total energy that are described in terms of the quantum surface. As consequence we could obtain very high accuracy on energies and forces. The computational cost does not increase significantly with respect to standard gas-phase calculations. We expect that these features make our methods useful for numerous applications in material science. Fixed the atomic radii, only two fitting 34


Page 35 of 56

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 4: Free energy change ∆G as a function of the bias potential U for O, HO, HOO and H2 O terminated CdS (11¯ 20) surfaces at pH = 0 and pH = 14. Results obtained in pure gas-phase conditions (solid lines) and with the presence of an implicit water environment (dotted lines) are shown. Various threshold bias potentials are marked.

35



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

parameters are needed to define the model, independently of the cluster or slab geometry. Mean absolute errors with respect to experimental solvation energies for neutral molecules and ions are close to the performances of others solvation approaches with a larger number of parameters. In particular for ions the implemented approach outperforms in accuracy all existing implicit solvation models. The flexibility of the model to locally change atomic radii allows to simultaneously simulate neutrals and ions, conditions typical of a double layer of an electrolyte in contact with a surface. A new benchmark protocol for solid-liquid interfaces is proposed in terms of wettability and contact angles. This allows to find parameters for elements not represented in current solvation databases. Results are shown for several metallic and insulating surfaces in contact with water. As a showcase application we studied the terminations of CdS (11¯20) surface in an electrochemical medium as a function of bias potential and pH. We observed that the inclusion of solvent effects introduces a substantial shift in the bias potentials at which surface terminations change. The relative stabilities of several adsorbed species are affected by the presence of solvent. As a consequence, its inclusion in computations which are based on the CHE approach is crucial.

Acknowledgement G. Fisicaro thanks Dr. D. Tomerini for continuous useful discussions. This work was done within the PASC and NCCR MARVEL projects. Computer resources were provided by the Swiss National Supercomputing Centre (CSCS) under Project ID s707. L. Genovese acknowledges also support from the EXTMOS EU project. Supporting Information Convergence tests are reported varying the k -point sampling for the total energy as well as the zero-point energy E ZPE and the vibrational entropic contribution T S of the clean CdS surface. Formulas are presented for the free energy of liquid water GH2 O(l) (298 K, 1 bar). This information is available free of charge via the Internet at http://pubs.acs.org. 36


Page 36 of 56

Page 37 of 56

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


A

Set up of KS Self-Consistent cycle under implicit environments

In this section we will outline the main equations that have to be implemented in a DFT code for the wavefunction optimization procedure in presence of implicit environments. Formulas reported are general and do not depend on the particular functional form of the dependence on the atomic coordinates Ri or on the electronic charge density ρ(r) for the cavity.

A.1

Electrostatic Energy

For an implicit environment the Hartree Potential might be seen as a functional of both the electronic density and the dielectric function, i.e. φ = φ[ρ, ǫ; Ω]. Here Ω represents the real-space domain (together with the appropriate boundary conditions) choosen to solve Eq. (1). Clearly, linearity is preserved with respect to ρ, namely φ[ρ1 +ρ2 , ǫ; Ω] = φ[ρ1 , ǫ; Ω]+ φ[ρ2 , ǫ; Ω]. A given electrostatic potential is related to an electrostatic energy 1 Eel [ρ, ǫ; Ω] ≡ 2

Z

drρ(r)φ[ρ, ǫ; Ω](r) ,

(28)

Ω

that by construction is quadratic in the charge density ρ, which implies

Eel [ρ1 + ρ2 , ǫ; Ω] = Eel [ρ1 , ǫ; Ω] + Eel [ρ2 , ǫ; Ω] +

Z

drφ[ρ1 , ǫ; Ω](r)ρ2 (r) .

(29)

Ω

It can be easily demonstrated that the functional derivative of the electrostatic energy with respect to the dielectric function is: 1 δEel [ρ, ǫ; Ω] = − |∇φ[ρ, ǫ; Ω](r)|2 . δǫ(r) 8π

(30)

As the electrostatic environment is modified, also the terms describing the ionic energy and potential must be defined accordingly. In the definition of the external potential con37



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 38 of 56

tribution, we have to treat consistently the contribution of v(r) coming from the solution of the generalized Poisson equation where ionic charges ρaion can be given as input. In other terms, we can define v(r) = vion (r) + v∆ (r) where vion = φ[ρion , ǫ; Ω], and ρion =

P

a

(31)

ρaion . Here v∆ contains extra terms, localised around

the atoms, coming for instance from the definition of the pseudopotentials. In this framework, the ionic energy of the system can equivalently be defined as the electrostatic energy of the above generalized Poisson equation minus the self-energy contribution of each of the ions in vacuum. Eion [ǫ] = Eel [ρion , ǫ; Ω] − Eself ,

(32)

where Eself =

X

Eel [ρaion , 1; R3 ]

(33)

a

A.2

Calculating the Ground-State energy in the Kohn-Sham formulation

The total energy within the Kohn-Sham formalism and the soft-sphere solvation model might be defined as:

E[ρ, ǫ] = Eion [ǫ] + Ts [ρ] +

Z

v(r)ρ(r)dr + Eel [ρ, ǫ; Ω] + Exc [ρ] + (α + γ)S[ǫ] + βV [ǫ]

= Ts [ρ] + E∆ [ρ] + Exc [ρ] + Eel [ρ + ρion , ǫ; Ω] − Eself + (α + γ)S[ǫ] + βV [ǫ] .

(34)

where Ts [ρ] is the electron kinetic energy, Eel [ρ + ρion , ǫ; Ω] the total electrostatic energy related to the ionic charges ρion and electronic charge density ρ, Exc [ρ] the exchange-correlation term. The additional terms E∆ [ρ] and Eself are both related to the ionic charges. In parR ticular E∆ = drv∆ (r)ρ(r). The non-electrostatic terms to the total energy come from Eq.

(5).

38


Page 39 of 56

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


Within this formula, the Kohn-Sham hamiltonian can be defined as HKS = −1/2∇2 +VKS , where

VKS [ρ](r) = Vxc (r) + v∆ (r) + φtot (r) + Vextra (r) , ( ) Z ′ 1 δS[ǫ] β δǫ(r) δǫ(r ) Vextra (r) = − dr′ |∇φtot (r′ )|2 + (α + γ) − . 8π δρ(r) δǫ(r′ ) ǫ0 − 1 δρ(r)

(35) (36)

Here we called the total solute charge ρtot = ρ + ρion and the corresponding potential φtot = φ[ρtot , ǫ; Ω]. In Eqs. (35) and (36) an explicit dependence of the dielectric function from the charge density has been considered, namely ǫ = ǫ[ρ], and δS[ǫ]/δǫ(r′ ) is explicited in Eq. (13). When ǫ is a local function of the density, 18,22 we have δǫ(r′ )/δρ(r) = ǫ(r)δ(r ˙ − r′ ) and ∂i ǫ = ǫ∂ ˙ i ρ. Therefore ǫ˙ ∂i ρ∂j ρ∂i ∂j ρ ∇2 ρ ǫ˙ δS − = δρ ǫ0 − 1 |ǫ| ˙ |∇ρ|3 |∇ρ| |ǫ| ˙ ∂i ρ =− ∂i ǫ0 − 1 |∇ρ| α+γ 1 |∇φtot (r)|2 + Vextra (r) = −ǫ(r) ˙ 8π ǫ0 − 1 ∇2 ρ(r) β ∂i ρ(r)∂j ρ(r)∂i ∂j ρ(r) − + |∇ρ(r)| ǫ0 − 1 |∇ρ(r)|3

(37) (38) (39)

As the dielectric function ǫ might evolve with the charge density, it is important to reevaluate the ionic contribution to the total energy for each of the density optimisation steps. In the above equation we have also assumed that ǫ/| ˙ ǫ| ˙ = −1 as the dielectric function should increase with the decreasing of the density. The DFT energy can be alternatively obtained from the band structure term of the KS Hamiltonian: EBS =

X j

fj hψj |HKS |ψj i = Ts [ρ] + Epot

39


(40)


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 40 of 56

where the potential energy is defined as the expectation value of the local potential:

Epot =

Z

drVKS (r)ρ(r) = Eeltot − Eion − Eself + Eel + E∆ + EVxc + EVex .

(41)

Here we called Eeltot = Eel [ρtot , ǫ; Ω], and EVex = EVxc =

Z

Z

drVextra (r)ρ(r) ,

(42)

drVxc (r)ρ(r) ,

(43)

By defining the Hartree energy of the cavity as c EH ≡ Epot − Eeltot − E∆ − EVxc ,

(44)

we can recover the formal expression for the KS energy:

c E = EBS − EH − EVxc − Eself .

(45)

The above expression might be interesting as it does not require an explicit evaluation of the double-counting term EVex .

A.3

Forces for the soft-sphere model

Forces are defined as minus the derivative of the total energy E[ρ, ǫ] [Eq. (34)] with respect to the atomic positions: Fa = −

dE . dRa

(46)

Here the subscript a points to a particular atom. The Hellmann-Feynman theorem states that, at self-consistency, the total derivative becomes the derivative with respect to the explicit dependence on the atomic positions. Since the electronic kinetic energy Ts [ρ] and the exchange-correlation contribution Exc [ρ] do not 40


Page 41 of 56

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


depend explicitly on the atomic coordinates, the only terms contributing to the forces are related to the total electrostatic energy and to the non-electrostatic contributions of the surrounding dielectric related to the quantum surface and volume. Contributions also come from the additional terms E∆ [ρ] and Eself . Being φtot = φ[ρtot , ǫ; Ω] the electrostatic potential generated by ρ + ρion , the contribution coming from the total electrostatic energy Eel [ρ + ρion , ǫ; Ω] can be split in different parts: ∂Eel = ∂Ra

Z

∂ρion (r) 1 drφtot (r) − ∂Ra 8π

Z

dr

∂ǫ(r) |∇φtot (r)|2 , ∂Ra

(47)

which gives rise to the following expression ∂E∆ − Fa = − ∂Ra

Z

∂ρion (r) ∂R Z a ∂ǫ(r) δS[ǫ] δV [ǫ] 1 2 − dr (α + γ) . (48) +β − |∇φtot (r)| ∂Ra δǫ(r) δǫ(r) 8π

drφtot (r)

By symmetry ∂Eself /∂Ra = 0. The first two terms can be calculated as in vacuum. Considering the particular case of our soft-sphere cavity defined by Eq. (6), all non-electrostatic additional terms can be analytically evaluated. To compute the nabla of φtot , central, forward and backward finite difference filters of order 16 have been used, which match the accuracy of the underlying vacuum Poisson solver. 27

B

Input guess implementation for the PCG algorithm

In our previous work 9 we reported a preconditioned conjugate gradient (PCG) algorithm to solve the generalized Poisson equation, i.e. Eq. (1), to compute the electrostatic potential φ(r) generated by a given charge density ρ(r). We demonstrated that it is able to solve the electrostatic problem with some ten iterations reaching an accuracy of ∼ 10−10 arb. units with respect to analytic functions. During a self-consistent field (SCF) loop, the generalized Poisson equation has to be 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

Page 42 of 56

solved several times as much as the number of iterations needed to reach the wavefunction convergence. Solving the problem repeatedly from scratch is computationally inefficient, especially when the wavefunction is approaching its final configuration. A proper use of an input guess for the potential φ(r) from previous calculations can considerably speed up the solution of the electrostatic problem and, by consequence, the whole SCF loop making implicit solvation and gas phase runtime closer. In our previous paper we implemented a standard input guess approach for the electrostatic potential, where the generalized Poisson operator A=−

1 ∇ · ǫ(r)∇ 4π

(49)

was applied to the potential input guess φ0 and the minimization procedure started from the initial gradient r0 = ρ − Aφ0 . Here ǫ(r) is the continuum dielectric function describing the dielectric cavity. This approach makes use of a finite difference filter for the nabla operator ∇, which could affect the overall converge especially when the mesh of the simulation box is too coarse or numbers become small with respect to the numerical noise. Such an effect could definitely corrupt the overall PCG performances. In order to prevent that, a solution has been found and reported in Algorithm 1. The implemented idea follows the same mathematical argument which allowed us to avoid the use of finite difference filters within the PCG scheme (see Ref. [ 9] for details). The minimization procedure starts from the initial gradient r0 . P is the preconditioner which inverse has to be applied to the residual vector rk returning the preconditioned residual vk , and, finally, φk is the solution of the generalized Poisson equation. The convergence criterion is imposed on the Euclidean norm of the residual vector rk . Both performances and accuracy of the suggested PCG scheme (with input guess) has been tested for a single-point DFT energy evaluation of a water molecule in implicit water described by the soft-sphere solvation model. Fig. 5 report the gradient norm during the wave function optimization. The SCF loop converged in 15 iterations. Red vertical dots represent the number of PCG iterations needed to solve the generalized Poisson equation on 42


Page 43 of 56

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


Algorithm 1 Preconditioned conjugate gradient (PCG) algorithm with input guess 1 1: compute q = 4π

√

√ ǫ ∇2 ǫ

2: r0 = ρ − φ0 q 3: v0 = P −1 r0 4: r1 = (φ0 − v0 )q 5: set p0 = s0 = 0, φ1 = v0 6: for k = 1, ... do 7:

if krk k2 ≤ τ exit

8:

vk = P −1 rk

9:

pk = vk + βk pk−1 (where βk =

(vk ,rk ) (vk−1 ,rk−1 ) , k

6= 0 )

10:

sk = Apk = vk q + rk + βk sk−1 (since Avk = vk q + rk )

11:

αk =

12:

φk+1 = φk + αk pk

13:

rk+1 = rk − αk sk

(vk ,rk ) (pk ,sk )

14: end for

43



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

each SCF step. Beside the input guess, a dynamic exit from the PCG loop has been also used as accurate potentials are not needed during the first stages of SCF optimization. Hence the exit threshold τ from the PCG loop on the residual norm krk k2 is made dynamic, set equal to the norm of the Kohn-Sham energy gradient from the current wave function optimization and allowed to vary between τmin and τmax (setting τmin = 10−4 and τmin = 10−6 guarantees a correct SCF convergence while lowering the overall electrostatic calculation time for the water molecule in implicit water). Fig. 5 shows that a mean of 4 PCG iterations are needed for each SCF iteration, which ultimately consists of a standard Poisson application and a mere vector multiplication. Tests on a larger system have been reported in our previous work 9 for a protein PDB ID: 1y49 (122 atoms) and compared to an optimized version of the COSMO 51 Poisson Solver for dielectric environments. 101 Applying the soft-sphere cavity with the new input guess implementation (Algorithm 1), we recovered a ratio of the wave function optimization runtime in solvent and gas phase of 1.16.

Figure 5: SCF convergence during a single-point DFT energy evaluation for water molecules in implicit water. Red vertical dots represent the number of PCG iterations to solve the generalized Poisson equation. Input guess and dynamic exit from the PCG loop have been used.

44


Page 44 of 56

Page 45 of 56

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


References (1) Dabo, I.; Li, Y.; Bonnet, N.; Marzari, N. In Fuel Cell Science: Theory, Fundamentals and Bio-Catalysis; Wieckowski, A., Nørskov, J., Eds.; John Wiley & Sons, Inc., 2010; pp 415–431. (2) White, J. A.; Schwegler, E.; Galli, G.; Gygi, F. The solvation of Na+ in water: Firstprinciples simulations. J. Chem. Phys. 2000, 113, 4668–4673. (3) Goedecker, S. Minima hopping: An efficient search method for the global minimum of the potential energy surface of complex molecular systems. J. Chem. Phys. 2004, 120, 9911–9917. (4) Schaefer, B.; Mohr, S.; Amsler, M.; Goedecker, S. Minima hopping guided path search: An efficient method for finding complex chemical reaction pathways. J. Chem. Phys. 2014, 140, –. (5) Sprik, M.; Hutter, J.; Parrinello, M. Ab initio molecular dynamics simulation of liquid water: Comparison of three gradient-corrected density functionals. J. Chem. Phys. 1996, 105, 1142–1152. (6) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. J. Chem. Phys. 2004, 120, 300–311. (7) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. The influence of temperature and density functional models in ab initio molecular dynamics simulation of liquid water. J. Chem. Phys. 2005, 122, 014515. (8) Sit, P. H.-L.; Marzari, N. Static and dynamical properties of heavy water at ambient conditions from first-principles molecular dynamics. J. Chem. Phys. 2005, 122, 204510. 45



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

(9) Fisicaro, G.; Genovese, L.; Andreussi, O.; Marzari, N.; Goedecker, S. A generalized Poisson and Poisson-Boltzmann solver for electrostatic environments. J. Chem. Phys. 2016, 144 . (10) Onsager, L. Electric Moments of Molecules in Liquids. J. Am. Chem. Soc. 1936, 58, 1486–1493. (11) Tomasi, J.; Persico, M. Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent. Chem. Rev. 1994, 94, 2027–2094. (12) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (13) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161–2200. (14) Senn, H. M.; Margl, P. M.; Schmid, R.; Ziegler, T.; Blöchl, P. E. Ab initio molecular dynamics with a continuum solvation model. J. Chem. Phys. 2003, 118 . (15) Lange, A. W.; Herbert, J. M. A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach. J. Chem. Phys. 2010, 133 . (16) Foresman, J. B.; Keith, T. A.; Wiberg, K. B.; Snoonian, J.; Frisch, M. J. Solvent Effects. 5. Influence of Cavity Shape, Truncation of Electrostatics, and Electron Correlation on ab Initio Reaction Field Calculations. J. Phys. Chem. 1996, 100, 16098– 16104. (17) Wiberg, K. B.; Keith, T. A.; Frisch, M. J.; Murcko, M. Solvent Effects on 1,2Dihaloethane Gauche/Trans Ratios. J. Phys. Chem. 1995, 99, 9072–9079. (18) Fattebert, J.-L.; Gygi, F. Density functional theory for efficient ab initio molecular dynamics simulations in solution. J. Comput. Chem. 2002, 23, 662–666. 46


Page 46 of 56

Page 47 of 56

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


(19) 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. J. Chem. Phys. 2006, 124, –. (20) 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 (Europhys. Lett.) 2011, 95, 43001. (21) Cococcioni, M.; Mauri, F.; Ceder, G.; Marzari, N. Electronic-Enthalpy Functional for Finite Systems Under Pressure. Phys. Rev. Lett. 2005, 94, 145501. (22) Andreussi, O.; Dabo, I.; Marzari, N. Revised self-consistent continuum solvation in electronic-structure calculations. J. Chem. Phys. 2012, 136, –. (23) Dupont, C.; Andreussi, O.; Marzari, N. Self-consistent continuum solvation (SCCS): The case of charged systems. J. Chem. Phys. 2013, 139 . (24) Sánchez, V. M.; Sued, M.; Scherlis, D. A. First-principles molecular dynamics simulations at solid-liquid interfaces with a continuum solvent. J. Chem. Phys. 2009, 131, 174108. (25) Genovese, L.; Neelov, A.; Goedecker, S.; Deutsch, T.; Ghasemi, S. A.; Willand, A.; Caliste, D.; Zilberberg, O.; Rayson, M.; Bergman, A.; Schneider, R. Daubechies wavelets as a basis set for density functional pseudopotential calculations. J. Chem. Phys. 2008, 129 . (26) http://www.bigdft.org (27) Genovese, L.; Deutsch, T.; Neelov, A.; Goedecker, S.; Beylkin, G. Efficient solution of Poisson’s equation with free boundary conditions. J. Chem. Phys. 2006, 125, –. (28) Marenich, A. V.; Kelly, C. P.; Thompson, J. D.; Hawkins, G. D.; Chambers, C. C.;

47



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

Giesen, D. J.; Winget, P.; Cramer, C. J.; Truhlar, D. G. Minnesota Solvation Database - version 2012 ; University of Minnesota, Minneapolis, 2012; p 3. (29) O’Regan, B.; Gratzel, M. A low-cost, high-efficiency solar cell based on dye-sensitized colloidal TiO2 films. Nature 1991, 353, 737–740. (30) Kojima, A.; Teshima, K.; Shirai, Y.; Miyasaka, T. Organometal Halide Perovskites as Visible-Light Sensitizers for Photovoltaic Cells. J. Am. Chem. Soc. 2009, 131, 6050–6051. (31) Lee, M. M.; Teuscher, J.; Miyasaka, T.; Murakami, T. N.; Snaith, H. J. Efficient Hybrid Solar Cells Based on Meso-Superstructured Organometal Halide Perovskites. Science 2012, 338, 643–647. (32) Kudo, A.; Miseki, Y. Heterogeneous photocatalyst materials for water splitting. Chem. Soc. Rev. 2009, 38, 253–278. (33) Nocera, D. G. The Artificial Leaf. Acc. Chem. Res. 2012, 45, 767–776. (34) Amovilli, C.; Mennucci, B. Self-Consistent-Field Calculation of Pauli Repulsion and Dispersion Contributions to the Solvation Free Energy in the Polarizable Continuum Model. J. Phys. Chem. B 1997, 101, 1051–1057. (35) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Acc. Chem. Res. 2008, 41, 760–768. (36) Barone, V.; Cossi, M.; Tomasi, J. Geometry optimization of molecular structures in solution by the polarizable continuum model. J. Comput. Chem. 1998, 19, 404–417. (37) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SM6:âĂĽ A Density Functional Theory Continuum Solvation Model for Calculating Aqueous Solvation Free Energies of Neutrals, Ions, and Solute-Water Clusters. J. Chem. Theory Comput. 2005, 1, 1133–1152.

48


Page 48 of 56

Page 49 of 56

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


(38) Marenich, A. V.; Olson, R. M.; Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SelfConsistent Reaction Field Model for Aqueous and Nonaqueous Solutions Based on Accurate Polarized Partial Charges. J. Chem. Theory Comput. 2007, 3, 2011–2033. (39) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Generalized Born Solvation Model SM12. J. Chem. Theory Comput. 2013, 9, 609–620. (40) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. (41) Willand, A.; Kvashnin, Y. O.; Genovese, L.; Vázquez-Mayagoitia, A.; Deb, A. K.; Sadeghi, A.; Deutsch, T.; Goedecker, S. Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations. J. Chem. Phys. 2013, 138, –. (42) Marques, M. A.; Oliveira, M. J.; Burnus, T. Libxc: A library of exchange and correlation functionals for density functional theory. Comput. Phys. Comm. 2012, 183, 2272 – 2281. (43) Genovese, L.; Deutsch, T.; Goedecker, S. Efficient and accurate three-dimensional Poisson solver for surface problems. J. Chem. Phys. 2007, 127, –. (44) Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field. J. Chem. Theory Comput. 2010, 6, 1509–1519. (45) Miertuˇs, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilization of {AB} initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117 – 129. (46) Bondi, A. van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441–451. (47) Pauling, L. The Nature of the Chemical Bond, 3rd ed.; Cornell University Press:Ithaca, NY, 1960; p 257. 49



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

(48) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard III, W. A.; Skiff, W. M. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035. (49) 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. J. Chem. Phys. 1997, 107, 3032–3041. (50) 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. J. Phys. Chem. B 1997, 101, 10506–10517. (51) Klamt, A.; Schuurmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799–805. (52) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347–1363. (53) Tannor, D. J.; Marten, B.; Murphy, R.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Honig, B.; Ringnalda, M.; Goddard, W. A. Accurate First Principles Calculation of Molecular Charge Distributions and Solvation Energies from Ab Initio Quantum Mechanics and Continuum Dielectric Theory. J. Am. Chem. Soc. 1994, 116, 11875– 11882. (54) Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. New Model for Calculation of Solvation Free Energies: Correc-

50


Page 50 of 56

Page 51 of 56

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


tion of Self-Consistent Reaction Field Continuum Dielectric Theory for Short-Range Hydrogen-Bonding Effects. J. Phys. Chem. 1996, 100, 11775–11788. (55) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; Dam, H. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; de Jong, W. NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations. Comput. Phys. Commun. 2010, 181, 1477 – 1489. (56) Orozco, M.; Luque, F. Optimization of the cavity size for ab initio MST-SCRF calculations of monovalent ions. Chem. Phys. 1994, 182, 237 – 248. (57) Miertuš, S.; Bartoš, J.; Trebatická, M. Dependance of atomic radii and volumes on the electron distribution in solute molecule and on solute-solvent interaction. J. Mol. Liq. 1987, 33, 139 – 156. (58) Haynes, W. M. CRC Handbook of Chemistry and Physics; CPC Press, Boca Raton, 2012; Chapter 6, p 187. (59) Earhart, H. W.; Komin, A. P. Kirk-Othmer Encyclopedia of Chemical Technology; John Wiley & Sons, Inc., 2000; p 1. (60) Sundararaman, R.; Schwarz, K. A.; Letchworth-Weaver, K.; Arias, T. A. Spicing up continuum solvation models with SaLSA: The spherically averaged liquid susceptibility ansatz. J. Chem. Phys. 2015, 142, 054102. (61) Drelich, J.; Chibowski, E.; Meng, D. D.; Terpilowski, K. Hydrophilic and superhydrophilic surfaces and materials. Soft Matter 2011, 7, 9804–9828. (62) Sedlmeier, F.; Janecek, J.; Sendner, C.; Bocquet, L.; Netz, R. R.; Horinek, D. Water at polar and nonpolar solid walls (Review). Biointerphases 2008, 3 . (63) Young, T. An Essay on the Cohesion of Fluids. Phil. Trans. R. Soc. Lond. 1805, 95, 65–87. 51



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

(64) Tadmor, R. Line Energy and the Relation between Advancing, Receding, and Young Contact Angles. Langmuir 2004, 20, 7659–7664. (65) Dupré, A. M.; Dupré, P. Théorie mécanique de la chaleur ; Gauthier-Villars Paris, 1869; p 369. (66) Girifalco, L.; Good, R. A theory for the estimation of surface and interfacial energies. I. Derivation and application to interfacial tension. J. Phys. Chem. 1957, 61, 904–909. (67) Good, R.; Girifalco, L. A theory for estimation of surface and interfacial energies. III. Estimation of surface energies of solids from contact angle data. J. Phys. Chem. 1960, 64, 561–565, cited By 356. (68) Fowkes, F. Additivity of intermolecular forces at interfaces. I. Determination of the contribution to surface and interfacial tensions of dispersion forces in various liquids. J. Phys. Chem. 1963, 67, 2538–2541. (69) Li, D.; Neumann, A. Equation of state for interfacial tensions of solid-liquid systems. Adv. Colloid Interface Sci. 1992, 39, 299 – 345. (70) van Oss, C. J. Acid-base interfacial interactions in aqueous media. Colloids Surf. A: Physicochem. Eng. Aspects 1993, 78, 1–49. (71) C. J. van Oss, Interfacial Forces in Aqueous Media, Second Edition; CRC Press, 2006; pp 131–155. (72) Goumans, T. P. M.; Wander, A.; Brown, W. A.; Catlow, C. R. A. Structure and stability of the (001) α-quartz surface. Phys. Chem. Chem. Phys. 2007, 9, 2146–2152. (73) Schrader, M. E. Wettability of clean metal surfaces. J. Colloid Interface Sci. 1984, 100, 372 – 380.

52


Page 52 of 56

Page 53 of 56

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


(74) Kim, E.-K.; Kim, J.-Y.; Kim, S. S. Significant change in water contact angle of electrospray-synthesized SiO2 films depending on their surface morphology. Surf. Interface Anal. 2013, 45, 656–660. (75) Sendner, C.; Horinek, D.; Bocquet, L.; Netz, R. R. Interfacial Water at Hydrophobic and Hydrophilic Surfaces: Slip, Viscosity, and Diffusion. Langmuir 2009, 25, 10768– 10781. (76) Wang, S.; Zhang, Y.; Abidi, N.; Cabrales, L. Wettability and Surface Free Energy of Graphene Films. Langmuir 2009, 25, 11078–11081. (77) Fowkes, F. M.; Harkins, W. D. The State of Monolayers Adsorbed at the Interface Solid-Aqueous Solution. J. Am. Chem. Soc. 1940, 62, 3377–3386. (78) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. On the Water-Carbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes. J. Phys. Chem. B 2003, 107, 1345–1352. (79) Taherian, F.; Marcon, V.; van der Vegt, N. F. A.; Leroy, F. What Is the Contact Angle of Water on Graphene? Langmuir 2013, 29, 1457–1465. (80) Kudo, A.; Miseki, Y. Heterogeneous photocatalyst materials for water splitting. Chem. Soc. Rev. 2009, 38, 253–278. (81) Walter, M. G.; Warren, E. L.; McKone, J. R.; Boettcher, S. W.; Mi, Q.; Santori, E. A.; Lewis, N. S. Solar Water Splitting Cells. Chem. Rev. 2010, 110, 6446–6473. (82) Ashokkumar, M. An overview on semiconductor particulate systems for photoproduction of hydrogen. Int. J. Hydrogen Energy 1998, 23, 427 – 438. (83) Bao, N.; Shen, L.; Takata, T.; Domen, K. Self-Templated Synthesis of Nanoporous CdS Nanostructures for Highly Efficient Photocatalytic Hydrogen Production under Visible Light. Chem. Mater. 2008, 20, 110–117. 53



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

(84) Simon, T.; Bouchonville, N.; Berr, M. J.; Vaneski, A.; Adrović, A.; Volbers, D.; Wyrwich, R.; Döblinger, M.; Susha, A. S.; Rogach, A. L.; Jäckel, F.; Stolarczyk, J. K.; Feldmann, J. Redox shuttle mechanism enhances photocatalytic H2 generation on Ni-decorated CdS nanorods. Nat. Mater. 2014, 13, 1013–1018. (85) Kalisman, P.; Nakibli, Y.; Amirav, L. Perfect Photon-to-Hydrogen Conversion Efficiency. Nano Lett. 2016, 16, 1776–1781. (86) Nørskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; Jónsson, H. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108, 17886–17892. (87) Valdés, Á.; Qu, Z.-W.; Kroes, G.-J.; Rossmeisl, J.; Nørskov, J. K. Oxidation and Photo-Oxidation of Water on TiO2 Surface. J. Phys. Chem. C 2008, 112, 9872–9879. (88) Bajdich, M.; García-Mota, M.; Vojvodic, A.; Nørskov, J. K.; Bell, A. T. Theoretical Investigation of the Activity of Cobalt Oxides for the Electrochemical Oxidation of Water. J. Am. Chem. Soc. 2013, 135, 13521–13530. (89) Neufeld, O.; Toroker, M. C. Platinum-Doped α-Fe2O3 for Enhanced Water Splitting Efficiency: A DFT+U Study. J. Phys. Chem. C 2015, 119, 5836–5847. (90) Zhang, X.; Klaver, P.; van Santen, R.; van de Sanden, M. C. M.; Bieberle-Hütter, A. Oxygen Evolution at Hematite Surfaces: The Impact of Structure and Oxygen Vacancies on Lowering the Overpotential. J. Phys. Chem. C 2016, 120, 18201–18208. (91) Nguyen, M.-T.; Piccinin, S.; Seriani, N.; Gebauer, R. Photo-Oxidation of Water on Defective Hematite(0001). ACS Catal. 2015, 5, 715–721. (92) Hellman, A.; Iandolo, B.; Wickman, B.; Grönbeck, H.; Baltrusaitis, J. Electrooxidation of water on hematite: Effects of surface termination and oxygen vacancies investigated by first-principles. Surf. Sci. 2015, 640, 45 – 49. 54


Page 54 of 56

Page 55 of 56

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


(93) Sakong, S.; Naderian, M.; Mathew, K.; Hennig, R. G.; Groß, A. Density functional theory study of the electrochemical interface between a Pt electrode and an aqueous electrolyte using an implicit solvent method. J. Chem. Phys. 2015, 142, –. (94) Sakong, S.; Groß, A. The Importance of the Electrochemical Environment in the Electro-Oxidation of Methanol on Pt(111). ACS Catal. 2016, 6, 5575–5586. (95) Lespes, N.; Filhol, J.-S. Using Implicit Solvent in Ab Initio Electrochemical Modeling: Investigating Li+/Li Electrochemistry at a Li/Solvent Interface. J. Chem. Theory Comput. 2015, 11, 3375–3382. (96) Letchworth-Weaver, K.; Arias, T. A. Joint density functional theory of the electrodeelectrolyte interface: Application to fixed electrode potentials, interfacial capacitances, and potentials of zero charge. Phys. Rev. B 2012, 86, 075140. (97) Zhou, Z.; Han, F.; Guo, L.; Prezhdo, O. V. Understanding divergent behaviors in the photocatalytic hydrogen evolution reaction on CdS and ZnS: a DFT based study. Phys. Chem. Chem. Phys. 2016, 18, 16862–16869. (98) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (99) NIST Computational Chemistry Comparison and Benchmark Database; NIST Standard Reference Database Number 101, Release 18, October 2016, Editor: Russell D. Johnson III. (100) Klaumünzer, M.; Mačković, M.; Ferstl, P.; Voigt, M.; Spiecker, E.; Meyer, B.; Peukert, W. Phase Transition Behavior and Oriented Aggregation During Precipitation of In(OH)3 and InOOH Nanocrystals. J. Phys. Chem. C 2012, 116, 24529–24537. (101) Liu, F.; Luehr, N.; Kulik, H. J.; MartÃŋnez, T. J. Quantum Chemistry for Solvated

55



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 ACS Paragon Plus Environment

Page 56 of 56

Soft-Sphere Continuum Solvation in Electronic-Structure Calculations

Recommend Documents