Implicit solvation using a generalized finite-difference approach in

Subscriber access provided by University of Sunderland

Condensed Matter, Interfaces, and Materials

Implicit solvation using a generalized finite-difference approach in CRYSTAL: implementation and results for molecules, polymers and surfaces Frédéric Labat, Bartolomeo Civalleri, and Roberto Dovesi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00762 • Publication Date (Web): 22 Oct 2018 Downloaded from http://pubs.acs.org on October 23, 2018

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

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

Page 1 of 48 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

Implicit solvation using a generalized finite-difference approach in CRYSTAL: implementation and results for molecules, polymers and surfaces Frédéric Labat,∗,† Bartolomeo Civalleri,‡ and Roberto Dovesi‡ PSL Research University, Chimie Paristech-CNRS, Institut de Recherche de Chimie de Paris, 11 rue P. et M. Curie, 75005 Paris, France , and Dipartimento di Chimica IFM, Università di Torino and NIS – Nanostructured Interfaces and Surfaces – Centre of Excellence, Via P. Giuria 7, 10125 Torino, Italy E-mail: [email protected]

Abstract We present the implementation of an implicit solvation model in the CRYSTAL code. The solvation energy is separated into two components: the electrostatic contribution arising from a self-consistent reaction field treatment obtained within a generalized finite-difference Poisson model, augmented by a non-electrostatic contribution proportional to the solvent-accessible surface area of the solute. A discontinuous dielectric boundary is used, along with a solvent-excluded surface built from interlocking atomcentered spheres on which apparent surface point charges are mapped. The procedure is ∗

To whom correspondence should be addressed PSL Research University, Chimie Paristech-CNRS, Institut de Recherche de Chimie de Paris, 11 rue P. et M. Curie, 75005 Paris, France ‡ Dipartimento di Chimica IFM, Università di Torino and NIS – Nanostructured Interfaces and Surfaces – Centre of Excellence, Via P. Giuria 7, 10125 Torino, Italy †

1

ACS Paragon Plus Environment

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

general and can be performed both at the Hartree-Fock and Density Functional Theory levels, with pure or hybrid functionals, for systems periodic in 0, 1 and 2 directions, that is for isolated molecules and extended polymers and surfaces. The Poisson equation resolution and apparent surface charge formalism is first validated on model analytical test cases. The good agreement obtained on solvation free energies is further confirmed by calculations performed on a large test set of 501 neutral molecules, for which a mean unsigned error of 1.3 kcal/mol is obtained when compared to the available experimental data. Importantly, the self-consistent reaction field procedure converges well for all molecules tested. This is further verified for all polymers and surfaces considered. In particular, for periodic systems, results obtained on an infinite glycine chain and on the wettability parameters of SiO2 surfaces are in good agreement with previously-published data. The size-extensivity of the energetic terms involved in the electrostatic contribution to the solvation energy is also well verified. These encouraging results constitute a first step to take into account complex environments in the CRYSTAL code, potentially allowing for a more accurate modeling of complex processes for both periodic and non-periodic systems.

1

Introduction

Taking into account a complex environment with quantum chemistry methods still remains challenging nowadays. The role of the environment surrounding the quantum system can however be crucial to obtain a realistic picture of several molecular properties, including photophysical properties of molecules 1 and extended polymers, 2 or reactivity on surfaces 3,4 for instance. While a full treatment of the whole model with its environment at the quantum level is conceptually simple, it is generally problematic or even impossible since it leads to a vast increase in the number of atoms and degrees of freedom. Two main strategies have thus been proposed to describe the environment at lower computational cost. Embedding techniques such as polarizable embedding 5–7 or self-consistent embedding 8–10 retain the atomistic

2


Page 2 of 48

Page 3 of 48 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


description of the environment, while dielectric continuum models 11,12 consider the environment as a dielectric continuum, that is a structure-less medium. The widespread use and development of these latter models can probably be related to their conceptual simplicity, which relies on defining a cavity in which the model is placed, selecting a dielectric form along with suitable parameters defining certain properties of the environment described as a continuum and accounting for the interaction of the model and its surroundings. In particular, this has been the model of choice to describe solvation effects implicitly at low-computional cost. 11 In such models, the solute-solvent interaction is separated in different contributions summing to the total solvation energy ∆Gtot s according to: elec ∆Gtot + ∆Gnelec . s = ∆Gs s

(1)

Electrostatic interactions (∆Gelec s ) are obtained by solving the Poisson equation, while ) correspond to the energy involved in the creation of non-electrostatic interactions (∆Gnelec s the solute cavity along with possible additional quantistic contributions such as dispersion and exchange. 11,13–15 The most popular methods for solving the Poisson equation discretize the domain of interest into small regions. Those methods include finite difference, 16–23 finite element 24–29 and boundary element methods. 30–34 In the molecular quantum chemistry field, the Polarizable Continuum Model (PCM) and its variants are probably the most popular implicit solvation models, since they were the first to introduce a molecular-shaped cavity that allows in principle to deal with any molecule of arbitrary shape. The interaction of the solute and the solvent is represented by a reaction potential which is incorporated into the solute Hamiltonian. By self-consistently solving for the reaction potential and the solute charge density, the mutual polarization of the solute and the solvent can then be obtained. Since the pioneering work of Miertuš et al., 30 PCM models have been constantly improved, 11,12 allowing in particular to accurately compute gradients and various response properties 35 with the

3



introduction of the continuous surface charge formalism 36 and subsequent variations. 37–39 Efforts have also been directed to non-equilibrium solvation description 40,41 and to the coupling with molecular mechanics approaches to describe large systems for instance. 42 More recently, efficiency for large systems has also been gained by solving the Poisson equation with the domain-decomposition paradigm. 43–45 In its traditional formulation, the PCM is a boundary element method using a discontinuous dielectric, whose conductor-like screening model (COSMO) 46 and conductor-like PCM (C-PCM) 47 variants have also been extended to periodic systems such as polymers and surfaces. 48–50 In materials modeling, alternative schemes based on a continuous dielectric have recently been proposed. All of them are derived from the original works of Fattebert and Gygi, 51,52 who proposed to formulate the dielectric constant of the solvent as a continuous function of the solute electron density, and to self-consistently solve for the reaction potential and the solute charge density. This is a natural formulation for density functional theory (DFT), which has seen various modifications and implementations in solid-state codes such as BigDFT, 53,54 CP2K, 55 FHI-aims, 56 ONETEP, 57,58 Quantum Espresso 53,59–62 or VASP 63,64 for instance. Here, we adopted a different strategy, relying on the original formulation of a discontinuous dielectric boundary. The solvation energy is separated into two components: the electrostatic contribution arising from a self-consistent reaction field treatment obtained within a generalized finite-difference Poisson model, augmented by a non-electrostatic contribution proportional to the solvent-accessible surface area of the solute. The solute cavity is a solvent-excluded surface, on which apparent surface charges are mapped to represent the solvent polarization due to the solute charge distribution, which are mutually equilibrated until self-consistency. This strategy is general and can be applied both at the Hartree-Fock (HF) and DFT levels, considering either pure or hybrid functionals. In addition, in view of a general implementation possibly involving salt effects even under non-linear conditions, a finite-difference approach has been chosen. Such methods have indeed already been successfully applied to solve both the Poisson and the Poisson-Boltzmann equations. 65 On the 4


Page 4 of 48

Page 5 of 48 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


other hand, boundary element methods typically involve computationally-expensive volume integrals when considering this latter, 66 making finite-difference approaches possibly more suitable. This paper is organized as follows. In Section 2, we present the methodology chosen to obtain the total solvation energy. We first detail the finite-difference approach adopted for the electrostatic contribution, its generalization to non-cubic grids and its extension to infinite periodic systems in one and two dimensions. Details concerning the solvent-excluded surface construction are then given, and the non-electrostatic contribution calculation is also commented. Section 3 presents the computational details. Results are then reported and discussed in Section 4, and conclusions are finally drawn in Section 5.

2

Methodology

We first begin by presenting the finite-difference methodology chosen for obtaining the electrostatic contribution to the solvation energy, and we introduce its generalization to non cubic grids and infinite periodic systems in one and two dimensions. The solvent-excluded surface construction along with the apparent surface charge position determination are then described. Finally, details regarding the non-electrostatic contribution calculation are given.

2.1

Electrostatic contribution

In a non-homogeneous dielectric distribution (r), the Poisson equation (PE) relates the total electrostatic potential φ(r) of the system to the solute charge density (ρint (r)) according to:

∇ · [(r) · ∇φ(r)] + 4πρint (r) = 0.

(2)

In a discontinuous dielectric treatment of the solute-solvent interface (see Fig. 1a), only two values for (r) are possible: in in Ωin (within the cavity) and out in Ωout (outside the cavity). In a finite-difference approach, the physical properties of the system are mapped onto 5



Γ

Page 6 of 48

Ωout

Ωin

k j

j

i

i

k

(a)

(b)

Figure 1: (a) Two dimensional, (10 × 10) representation of the electrostatic system for a CH3 CONH2 molecule showing: interior of cavity (Ωin ), solvent region (Ωout ), and grid edges (Γ). Small and large grey circles represent grid points where the PE is solved using an iterative solver. Orange circles at the edges represent grid points where potentials are permanently assigned in order to establish Dirichlet boundary conditions. Small red circles correspond to grid-line centers where the dielectricum data are assigned, either to the interior (in ) or exterior (out ) value in the sharp discontinuous dielectric model at the solute/solvent interface. Large grey circles are boundary grid points (BGP). Adapted from Ref. 17; (b) Elementary volume element for a three dimensional cubic grid, showing corresponding numbering of the six neighboring pairs of dielectric (in red) and grid points (in blue) from 1 to 6 of a central grid point (in black, 0).

6


Page 7 of 48 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 three dimensional lattice of points and derivatives involved in Eq. (2) are approximated as differences between the function at these points. The PE must be satisfied everywhere on the three dimensional grid, and especially at all interior grid points (Ωin ∪ Ωout , see Fig. 1a). Here, the solute charge density ρint (r) in Eq. (2) is represented by atomic charges distributed on the three dimensional grid using a quadratic inverse interpolation algorithm, where each atomic charge is spread onto the nearest 27 grid points. 67 For a cubic grid with N grid points per direction, the potential at any grid point (φo ) depends on the charge (q0 ) at that grid point, the grid spacing h, as well as the potentials and dielectric values of the six neighboring grid points (φn and n , n ∈ [1; 6], 7-point stencil formula, see Fig. 1b) 16 according to: φ0 =

6 X

!

n φn + 4πq0 /h /

n=1

6 X

(3)

n .

n=1

This has been commonly used in previous implementations in DELPHI, 68 UHBD 69 or MEAD 70 for instance. For a lattice-shaped grid with Nx Ny Nz points, a generalization of Eq. 3 can be obtained and the potential at any grid point can be expressed as:

"

2 4 6 hx hz X hx hy X hy hz X n φn + sin β n φn + sin γ n φn + 4πq0 φ0 = sin α hx n=1 hy n=3 hz n=5 # " 2 4 6 hx hz X hx hy X hy hz X n + sin β n + sin γ n , / sin α hx n=1 hy n=3 hz n=5

# (4)

with hx , hy and hz the grid spacings along the x, y and z directions, α, β and γ the inclination angles of the three axes along a, b and c, and V the volume of an elementary grid element. Full derivation of Eq. (4) is given in Appendix A. Application of Eq. (3) or (4) results in a system of Nx Ny Nz linear equations which are iteratively solved with permanently assigned potentials on Γ (see Fig. 1a, Dirichlet boundary conditions) along non-periodic directions until a convergence criteria is reached to get the total electrostatic potential at each grid point. In the current implementation, an optimal

7



Page 8 of 48

successive over-relaxation (O-SOR) 71,72 algorithm is used to solve this set of linear equations according to: (n)

(n−1)

φ0 = (1 − ω) · φ0

+ ω · φ0

(5)

where n and (n − 1) correspond to the current and previous interations, respectively, ω is the optimal over-relaxation parameter which is estimated from the spectral radius of the corresponding Gauss-Seidel approach as proposed in Ref. 19, and φ0 is obtained with Eq. (3) or (4). For molecules, the equivalent dipole to the charge distribution is evaluated, and the Debye-Hückel potentials are computed at all grid points i belonging to Γ according to: 68

φΓi = q+ ·

e(−r−,i /λD ) e(−r+,i /λD ) + q− · out r+,i out r−,i

(6)

where r+,i and r−,i are the distances from grid point i to the center of positive and negative charge, respectively, q+ and q− are the sums of all positive and negative charges of the solute, and λD is the Debye-Hückel length which is considered as infinite for a solvent without salt effects. For periodic systems, Ewald formulas for systems with reduced periodicity are used. 73 More precisely, for a polymer periodic in the x direction with N atoms with charges qn (n ∈ [1; N ]) per unit cell, a lattice parameter a and a translation vector T, the electrostatic potential at site with coordinates xm is given by: 0

φ

1D

(xm ) =

N XX T

−

1 a

N erfc (ξ|xm − xn + T|) 1 X X qn + qn e−ik1 (xm −xn ) K0 k12 /4ξ 2 , ρ2mn ξ 2 |xm − xn + T| a k 6=0 n=1 n=1 1

N X

2ξ qn γ + log ρ2mn ξ 2 + E1 ρ2mn ξ 2 − √ qm π n=1

(7)

where K0 (u, v) is an incomplete modified Bessel function of the second kind defined as:

K0 (u, v) =

Z

1

∞

dt −ut−v/t e t

(8)

and evaluated using an adaptive Romberg method, E1 (v) is the exponential integral defined 8


Page 9 of 48 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


as: E1 (v) =

Z

∞

1

1 −vt e dt = t

Z

∞

v

1 −t e dt, t

(9)

γ is the Euler Mascheroni constant, ξ is a parameter controlling the relative decay of the p real and reciprocal space sums, ρmn = (zm − zn )2 + (ym − yn )2 , k1 forms a discrete set of k points in the reciprocal space, and N 0 indicates that the term (n = m, T = 0) is excluded

from the real space sum. Similarly, for a slab periodic in the x and y directions with unit cell area A and a translation vector T: 0

φ

2D

(xm ) = +

N XX T n=1 N X

π A

qn

qn

n=1

erfc (ξ|xm − xn + T|) |xm − xn + T|

X k6=0

1 e−ik·(rm −rn ) g k, zm − zn , ξ k

√ N 2 πX 2ξ 1 −ξ2 (zm −zn )2 √ − qn e + π (zm − zn ) erf (ξ (zm − zn )) − √ qm A n=1 ξ π where: kz

g(k, z, ξ) = e erfc

k + ξz 2ξ

−kz

+e

erfc

k − ξz . 2ξ

(10)

(11)

and k forms a discrete set of vectors in the reciprocal space. This corresponds to the wellknown Parry formulation. 74 Apparent surface charges (ASC) values are determined from converged electrostatic potentials of boundary grid points (BGP), that is from grid points with at least one dielectric neighboring point in a different dielectric media than the other five (see Fig. 1a). For an ASC obtained from a BGP with distributed charge q0 , potential φ0 and a cubic grid with spacing h, it has already been shown 23 that the ASC value can be obtained from the following equation: q ASC = −q0 +

3h 2π

φ0 −

6 1X

6

!

φn .

n=1

(12)

This can also be generalized to any grid shape from Gauss’ law (see Appendix B for

9



Page 10 of 48

derivation): q ASC = −q0 +

V 2 2 2 2 2 2 (2φ − φ − φ ) . h h h (2φ − φ − φ ) + h h (2φ − φ − φ ) + h 0 5 6 0 1 2 0 3 4 y x y z x z 4πh2x h2y h2z (13)

By including this set of point charges into the solute Hamiltonian, a self-consistent reaction field (SCRF) treatment of solvation is hence possible. 22 Since only a surface polarization is considered in the present implementation, neglecting volume polarization, 75 a renormalization procedure 76 of the values of the ASC is applied at each SCRF cycle to verify the relation existing between the total solute charge Qs and the theoretical sum of ASC (Qth ASC ): Qth ASC = −

out − 1 · Qs . out

The difference th ∆QASC = Qth ASC − QASC = QASC −

(14)

X

qnASC

(15)

n

between the theoretical sum Qth ASC and the computed sum QASC before renormalization can then be used to monitor the quality of the ASC formalism employed. When evaluating the electrostatic contribution to the solvation energy, the fundamental contribution is the reaction field energy (Erf ), which corresponds to one-half of the interaction energy between the solute and the ASC. For a system with N ASC (NASC ), this contribution is obtained as a Coulomb interaction between the solute potentials at the positions of the ASC (φASC ) and their values (qnASC ), according to: 11,23 n

Erf =

NASC 1 f int f 1 X , Ψ V Ψ = qnASC · φASC n 2 2 n=1

(16)

where Ψf is the wavefunction of the solute in solution, and V int is the interaction operator between the solute and the ASC. For periodic systems, ASC belong to the reference unit cell and evaluation of φASC is made with Ewald formulas for systems with reduced periodicity 73 n

10


Page 11 of 48 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


as described above. Since the interaction operator V int describes the interaction between the solute and an averaged distribution of solvent, the electrostatic contribution to the solvation energy has the status of a free energy and is obtained as:

∆Gelec = Ψf H 0 + V int Ψf − Ψ0 H 0 Ψ0 − Erf s = ∆EEN + Erf

(17)

where H 0 is the solute Hamiltonian in gas-phase, Ψ0 is the solute wavefunction in gas-phase,

∆EEN = Ψf H 0 Ψf − hΨ0 |H 0 |Ψ0 i is the change in the electronic energy including nuclear

repulsion of the solute when it goes from the gas phase into solution and Erf is the reaction

field energy defined in Eq. (16).

2.2

Solvent-excluded surface construction

Different choices for the solute cavity are available, and their suitability for impliciy solvation calculations is still an open question. 77,78 Here, the solvent-excluded surface (SES) 79,80 has been chosen since it produces a smooth outer-surface contour by rolling a solvent probe over the solute. This was found especially important to produce a suitable cavity for periodic systems in two dimensions, in order to avoid cavities and voids of the solute to be filled with high-dielectric regions corresponding to the solvent, leading to an unphysical description of the solvation process. Three types of pieces form the SES: 80 sphere patches, tori and spherical triangles, when the solvent probe touches simultaneously one, two and three atoms, respectively. The first correspond to contact regions of the SES, while the other two belong to the reentrant portions. To build the SES, the algorithm chosen is a grid-based approach proposed by Rocchia et al., 23 in which spherical triangles are omitted. Briefly, the Van der Waals surface (VDWS) is first build analytically from the atom-centered spheres of chosen radii, and a first list of BGP is determined. Then, a list of all atom pairs having a circle of intersection (COI) 11



Page 12 of 48

is generated from solvent-expanded atomic radii and a number of equally spaced points are computed on each COI to make these circles almost continuous lines. The points lying inside any atom’s accessible surface are eliminated, leaving a skeleton of solvent-accessible points placed on the different COI found. Using distance calculations between these remaining points and the dielectric neighbors of the set of BGP found above and belonging to Ωout , the reentrant regions of the SES are iteratively constructed. Any of these dielectric points which are farther from the solvent-accessible points than a solvent probe radius are updated to in since they belong to Ωin , and the status of each BGP is checked since some will no longer be BGP while new grid points will become BGP. When the list of BGP is not modified anymore, the procedure is stopped. This typically requires 2 or 3 iterations, even for the largest systems tested in this article. Since many distance calculations are required, the whole procedure is accelerated by using a cubing algorithm to avoid unnecessary checks between distant atom/atom or atom/point pairs. Finally, the final set of BGP is projected onto the resulting SES distinguishing between points belonging to contact and reentrant sections, as described in Ref. 23, leading to the set of ASC included into the solute Hamiltonian. For periodic systems, the same procedure is applied with some slight modifications: (i) supercells are built with 3 and 27 (3×3) unit cells respectively for polymers and slabs; (ii) only BGP belonging to the reference unit cell are considered and projected onto the SES, resulting in an open cavity along the periodic directions. Fig. 2 presents examples of ASC mapped onto the SES for both a non-periodic and a periodic system.

2.3

Non-electrostatic contribution

The non-electrostatic contribution to the total solvation energy (∆Gnelec ) is chosen propors tional to the solvent-accessible surface area (SASA) of the solute, which is the surface traced by the center of a solvent probe rolling onto the solute, according to: ∆Gnelec = 0.860 + 0.005 · SASA. s 12


(18)

Page 13 of 48 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)

(b)

(c)

Figure 2: ASC mapped onto the SES obtained for: (a) a CH3 CONH2 molecule, with 284 point charges; (b) side and (c) top views of a 3 layer slab of MgO (001), with 72 point charges in the unit cell (shown as solid yellow lines). Yellow spheres correspond to ASC in contact regions of the SES, while green ones belong to the reentrant sections. All SES obtained with a grid spacing of about 0.5 Å along all directions and 65 grid points along non periodic directions. Note that atomic radii have been scaled down for clarity. of This formula has been obtained from a least-square fit of the dependence of ∆GExp s hydrocarbons in water to the calculated SASA. 21 In the present implementation, the SASA is numerically estimated using the Shrake and Rupley algorithm, 81 which distributes Nid equally-spaced points on each solvent-expanded atomic sphere i and determines those points which do not lie inside any other atom’s accessible surface Nia (see Fig. 3). The SASA is then estimated as: SASA =

X i

Nia 4π(Ri + Rp ) · d Ni 2

(19)

with Ri and Rp the radii of atom i and of the solvent probe, respectively. For periodic systems, only points belonging to the unit cell are considered. We should note here that Eq. (18) does neither depend on the solute accessible volume, which might be important for solutes with large sizes, 22 nor explicitly on quantistic effects such as dispersion and exchange which are included in other models. 11,13–15 However, the idea of choosing ∆Gnelec proportional to the SASA can be related to that found s

13



Probe

SAS

Figure 3: Shrake and Rupley algorithm illustrated for a system with two atomic spheres, showing van der Waals spheres as dashed line, solvent probe (in blue) rolling over solventexpanded atomic spheres as solid black lines. Green points represent points which do not lie inside any other atom’s accessible surface and are therefore considered as accessible to the solvent, while red ones are buried and thus eliminated. in other approaches such as the Cavitation, Dispersion, and Solvent-structure effects (CDS) approach 82–84 which accounts for these effects by employing empirical atomic surface tension terms and SASA of atoms.

3

Computational details

All calculations have been carried out with a locally modified version of the Crystal code, 85,86 which uses atom centered Gaussian orbitals as a basis set. DFT calculations have been performed using the global hybrid B3LYP functional. 87 In some cases, the PBE functional 88 has also been used for comparison purposes. All molecular and polymer calculations have been performed with 6-311G(d, p) basis sets, while Mg and O atoms were described with 8-61G and 8-51G basis sets for the MgO slab calculations. Suitable k points grids were used to ensure convergence of total energies of the periodic systems considered to within 10−7 Hartree. Numerical DFT integration has been performed considering 75 radial points and 974 angular points, ensuring an error on the integrated electron density to an accuracy of 10−5 e per unit cell or molecule. In all cases, Mulliken atomic charges obtained after the self-consistent field (SCF) procedure were used as an input to iteratively solve the PE using the O-SOR algorithm described above and a convergence criteria of 10−4 kT/e (∼ 2.6 10−6 V at 298.15 K) on the electrostatic 14


Page 14 of 48

Page 15 of 48 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


potential was selected. The SCRF is considered as converged when the difference in total energy between two successive cycles is below 2 · 10−5 Hartree, that is about 10−2 kcal/mol. A grid spacing of about 0.5 Å is used along each grid direction and the number of grid points considered is chosen so that the solute fills at most 50% of the grid along non-periodic directions. These choices were found to be sufficient to converge the presented data to within specified digits. All calculations reported considered water as the solvent, that is: out is taken as 78.3553, in as 1 and the solvent probe radius is 1.385 Å. Atomic radii of the H, C, N, O and S atoms were taken from the PARSE parameter set: 22 1.15, 1.90, 1.60, 1.60 and 1.90 Å, respectively. Unless explicitly specified otherwise, Bondi radii 89 have been used for all the other elements. No scaling factors were applied to improve the agreement with the available experimental data or with previously-published works.

4

Results and discussion

4.1

Analytical cases

Before discussing results obtained on complex systems, that is molecules, polymers and surfaces, we first briefly present data obtained on model systems constituted by a single or multiple charged spheres for which analytical solutions are available. In all cases, both the SCF and SCRF procedures are completely skipped for these preliminary calculations and the electrostatic contribution to the solvation energy is the reaction field component only. This allows us to test the implementation of the Poisson equation solver and apparent surface charge formalism. A cubic grid is chosen throughout. Six selected cases, already proposed in previous finite-difference implementations, 90,91 have been chosen since they cover a large range of solvation energies with either a single or multiple atomic charges.

15



4.1.1

Born ion

A Born ion corresponds to a sphere with radius R and dielectric constant in centered on a ±q charge placed in a medium with dielectric constant out . The solvation energy corresponds to the difference in work to charge the ion in the two environments according to: 92 ∆Gel s = −

q2 (in − out ) . 2R

(20)

Fig. 4 presents reaction field energies computed with a (81×81×81) grid and a grid spacing h = 0.500 Å, with R ranging from 1 to 3 Å and q varying from ±1 to ±3e. In all cases, a very good agreement between the analytical and numerical data is obtained. Smaller grid spacings of h = 0.250 or 0.125 Å lead to differences lower than 10−3 kcal/mol on the computed solvation energies when compared to the h = 0.500 Å case. 0

∆Gel s (kcal/mol)

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 16 of 48

−500

−1000

−1500

q = ±1e q = ±2e q = ±3e 1

1.5

2

2.5

3

R (˚ A)

Figure 4: Electrostatic solvation energies for a Born ion of charge q = ±1, ±2 and ±3e, with radius R ranging from 1 to 3 Å. Computed data (as dots) obtained with a (81×81×81) cubic grid and a grid spacing of h = 0.500 Å. Analytical values obtained with Eq. (20) shown as solid lines.

4.1.2

Kirkwood models

Five additional test cases where multiple charges are distributed within a single 2 Å atomic sphere centered at (0.0,0.0,0.0) Å have been considered: 1. two positive unit charges symmetrically placed at (1.0,0.0,0.0) and (-1.0,0.0,0.0) Å 16


Page 17 of 48 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. two positive unit charges symmetrically placed at (1.0,0.0,0.0) and (-1.0,0.0,0.0) Å, and two negative unit charges symmetrically placed at (0.0,1.0,0.0) and (0.0,-1.0,0.0) Å 3. two positive unit charges symmetrically placed at (1.2,0.0,0.0) and (-1.2,0.0,0.0) Å, and two negative unit charges symmetrically placed at (0.0,1.2,0.0) and (0.0,-1.2,0.0) Å 4. six positive unit charges placed at (0.4,0.0,0.0), (0.0,0.8,0.0), (0.0,0.0,1.2), (0.0,0.0,-0.4), (-0.8,0.0,0.0) and (0.0,-1.2,0.0) Å 5. six positive unit charges placed at (0.2,0.2,0.2), (0.5,0.5,0.5), (0.8,0.8,0.8), (-0.2,0.2,0.2), (0.5,-0.5,0.5) and (-0.8,-0.8,-0.8) Å Table 1 presents computed reaction field energies for these five cases using a cubic grid. Additional data obtained with different grid sizes are given in Table S1 of the supplementary information, the center of the grid being always located at the single sphere center. Analytical Table 1: Reaction field energies of the five Kirkwood models used as test cases and numbered as (1) to (5) (see text for definitions). All energetic data in kcal/mol. The data in parenthesis correspond to the deviation (in %) of the best numerical value (in bold face) from the analytical value. Grid size (1) (2) (3) (4) (5) (81×81×81) -349.82 -64.06 -137.15 -2988.32 -3123.69 (81×81×81) -350.11 -64.12 -137.52 -2990.03 -3126.00 (81×81×81) -350.00 -63.84 -137.20 -2989.94 -3126.50 (81×81×81) -349.97 -63.70 -137.68 -2990.30 -3125.09 (81×81×81) -349.82 -63.43 -136.52 -2989.51 -3126.51 (161×161×161) -349.71 -63.13 -136.01 -2989.40 -3125.05 (-0.00) (0.50) (0.52) (0.00) (±0.02) Analytical -349.73 -62.81 -135.40 -2989.30 -3124.30 h/Å 0.50 0.25 0.20 0.15 0.10 0.05

results are generally well-reproduced for the five cases considered, the largest errors (∼0.50 %) being obtained for cases (2) and (3). A very good agreement with previously-published data obtained using a finite-difference approach 91 is obtained. In addition, as shown in Table S1, the grid sizes weakly affect the data presented in Table 1, except for case (4) with h =0.5 Å for which non-negligible differences are obtained. It should however be noted that 17



cases (3) and (4) correspond to a hypothetical case in which charges are artificially very close to the dielectric boundary. When considering real solutes with atomic radii normally larger than 1.0 Å, such a case should not occur. Although no common trend can be clearly evidenced between grid spacing and relative error when compared to the analytical data for all five cases, very fine grid spacings (h ≤ 0.10 Å) are generally required to reproduce almost quantitatively the analytical data of these difficult test cases. This might be related to the fact that no specific conditions are applied to enforce continuity conditions at the dielectric interface in the current implementation. In particular, the continuity condition on the electric displacement 93 is not imposed. Approaches such as the finite-difference matched interface and boundary 94 which specifically enforces continuity conditions at the dielectric interface have already been shown to more precisely reproduce the analytical data with coarser grids for such difficult test cases. These first six test cases however show that the implemented procedure based on ASC performs well, nicely reproducing analytical data in most cases.

4.2

Molecules

We now turn to the case of non-periodic systems with complex cavities, that is to the case of molecules for which a vast amount of experimental data is available. Although the analytical test cases presented above are useful to test the code, as already shown, reproducing quantitatively the analytical data requires very fine grid spacings in the current implementation, which are unrealistic for application to large scale systems since they lead to a very large number of ASC which have to be included in the solute Hamiltonian at each step of the SCRF procedure. This is exemplified in Fig. 5 for CH3 CONH2 . Although the time needed to build the SES and map ASC onto it is negligible even at a very fine grid spacing of 0.10 Å with only 0.23 s, the number of ASC significantly increases when the grid spacing is below 0.20 Å, reaching 7228 at h =0.10 Å. A more reasonable number of ASC is however found when the grid spacing is around 0.4 or 0.5 Å, with values of 451 and 284, 18


Page 18 of 48

Page 19 of 48

respectively, that is a few hundreds which is the range targeted for routine applications. 1

6 000

0.6 4 000 0.4 2 000

0.2 0

0

0.2

0.4

0.6

0.8

1

Number of ASC

8 000

0.8

Time (s)

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


0

h (˚ A)

Figure 5: Time to build the SES and project ASC (in black) and corresponding number of ASC (in blue) for a CH3 CONH2 molecule with a (81×81×81) grid and varying grid spacings h. All data obtained on a single Intel Xeon CPU E5-2630 @ 2.30GHz. The grid resolution not only influences the number of ASC but also the computed solvation energies. In particular, from Fig. 6a, we can see how the total solvation energy of CH3 CONH2 depends on the grid spacing h. When considering values below 0.5 Å, ∆Gtot s converges to below 0.2 kcal/mol, which is precise enough for the kind of systems targeted. Interestingly, coarse grids lead to large errors on the total surface charge at the dielectric boundary, while negligible errors are obtained for reasonable grid resolutions (h < 0.5 Å, see Eq. (14)). These small errors can be related to the fact that the PE is solved using atomic charges distributed on a grid. Although the amount of solute electronic charge lying outside of the cavity can be expected to be small for the chosen solute molecule, this is an encouraging point since the escaped charge density can be a problem for formalisms based on charge density when the solute is an anion or an excited state or when an external field is acting on the solute for instance. 11,95 In addition, since a grid-based approach is used, the solute centering on the grid might also affect computed solvation energies. Fig. 6b shows the variation of ∆Gtot s for CH3 CONH2 when the solute is shifted continuously along the x, y and z axis of the grid by up to 2 Å. We can note that ∆Gtot s periodically oscillates around a mean value, the period of the oscillation 19



corresponding exactly to the grid spacing used in the calculation. However, differences on the solvation energies computed are small, ranging from 0.02 kcal/mol along x and z to about 0.20 kcal/mol along y, resulting in deviations from the reference value of 0.2 and 1.8 %, respectively. These two preliminary tests therefore show that the computed solvation energy is slightly grid-dependent for molecules but that differences obtained on solvation energies for this first molecular example remain reasonable. However, we can reasonably expect this dependency to be lower in the case of periodic systems which is our main focus, since the solute position on the grid along a periodic direction is determined by the lattice definition. It is also noteworthy that the grid edge positions where Dirichlet boundary conditions are applied have a negligible effect on the computed data. This is shown in Table S2 of the supplementary information, where a difference of about 0.003 kcal/mol is obtained between the highest and the lowest grid fillings considered. Finally, we carried out calculations on the Mobley test set, 96–98 with 501 neutral molecules containing the following chemical elements: H, C, N, O, F, P, S, Cl, Br and I. Most of the molecules considered in this test set are relatively small and rigid, and belong to the following families: alkanes, alkenes and dienes, alkynes, arenes, alcohols and phenols, ethers, ketones and aldehydes, carboxylic acids, esters, amines, amides, nitriles, nitrohydrocarbons, nitrogen heterocycles, thiols, sulfides, as well as chlorinated, fluorinated, iodated and multihalogenated species, and phosphorous compounds. We focus on neutral molecules here since we are mainly interested in applications of the implicit solvation model to infinite periodic models, where a majority of systems have neutral unit cells. All structures have been optimized at the B3LYP/6-311G(d, p) level and verified to be true minima by frequency calculations. Mean unsigned and signed errors (MUE and MSE) of 1.3 and -0.2 kcal/mol are obtained, along with a root-mean-square error (RMSE) of 1.6 kcal/mol, with a SCRF convergence obtained within 4 to 5 cycles on average. These data are in line with results reported for other SCRF models on similar test sets. 54,56,99 For instance, with a test set of 274 similar neutral molecules, 99 MUE and MSE of 1.2 and 0.6 kcal/mol have been reported 20


Page 20 of 48

10−2

−9.2

10−3

−9.4

10−4

−9.6

10−5

−9.8

10−6

∆QASC (e)

−9

−10 0.2

0.4

0.6

0.8

1

10−7

h (˚ A) (a) −9

∆Gtot s (kcal/mol)

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 21 of 48

x y z

−9.2 −9.4 −9.6 −9.8 −10

0

0.5

1

1.5

2

shift (˚ A) (b)

Figure 6: Computed total solvation energy (∆Gtot s ) of the CH3 CONH2 molecule in water, showing influence of (a) the grid resolution; (b) the solute position on the grid. A (65×65×65) grid with a grid spacing of 0.5 Å is used. The dashed black line corresponds to the experimental solvation energy. The mean error on ASC during the whole SCRF before renormalization (see Eq. (15)) is also reported in the top graph. at the mPW1PW/6-31G(d) level with the IEF-PCM/UAHF approach, thus validating the present implementation. Fig. 8 also presents computed total SASA of all solutes, compared to analytical values obtained using a sophisticated enveloping triangulation analytical approach with the CAVE program. 100 An excellent agreement is obtained, with a MUE of 1.43 Å2 , confirming the validity of the implementation of the numerical procedure adopted for SASA calculation.

21



5

Calculated

0

−5 −10 −15 −15

−10

−5

0

5

Exp

Figure 7: Computed total solvation energies in water compared to the experimental data for the Mobley test set with 501 neutral molecules. The black and red dashed lines correspond to errors of ±1.0 and ±1.5 kcal/mol compared to the experimental data, respectively. All data in kcal/mol. 500

400

Numerical

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 22 of 48

300

200

100 100

200

300

400

500

Analytical Figure 8: Computed total SASA for the Mobley test set with 501 neutral molecules compared to the analytical data obtained with CAVE. 100 All data in Å2 , obtained with 500 points per atomic spheres (see § 2.3 for details).

4.3

Periodic systems

We now turn to the case of infinite periodic systems in one and two dimensions, that is to polymers and slabs. Very few solvation energies have been reported for such systems so far, probably because no experimental solvation energies are directly available. We first discuss the size-extensivity properties inherent to any periodic system in one or two dimensions, and 22


Page 23 of 48 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


then we investigate the case of an infinite glycine chain and the wettability properties of SiO2 surfaces. 4.3.1

Size-extensivity

One important property that should be verified for periodic systems is the size-extensivity of the energetic terms involved in Eq. (17). More precisely, supercells of given sizes can be considered, and the energetic terms obtained can be compared to that of the reference unit cell. For a one dimensional periodic system such as a polymer, supercells with n unit cells (Sn ) should display a n dependence. On the other hand, when a (n × n) supercell of a slab is considered (Sn×n ), a n2 dependence of the data is expected. Fig. 9a presents the ratio of the computed data for different supercells and the reference unit cell for a glycine polymer. Corresponding data is shown on Fig. 9b for a MgO (001) surface slab with 3 layers, each layer being composed of a Mg-O dimer. It is clear that the computed data nicely follows the expected n and n2 behaviours. For instance, for the MgO supercells, the largest error obtained on the ∆Gel s term is of about 0.01 % for the S10×10 MgO supercell. In addition, from Fig. 9c and 9d, the error obtained on the computed total surface charge of the open cavities of both systems is very small (see Eq. (14)). This conclusion holds for all tested supercell sizes. Although these very low values can again be related to atomic charges used to solve the PE, it nonetheless indicates that the implementation of the Poisson model with apparent surface charge formalism correctly performs for infinite periodic systems with open cavities. 4.3.2

Glycine infinite chain

As already done in Ref. 48 and 49, oligomers of N-acetyl-Gly-Gly-C-methylamine with increasing lengths are compared to an infinite glycine chain. Since no structural data were available from these previous works, we first optimized the glycine polymer at the PBE level, and constructed oligomers from this optimized structure. A lattice parameter of 7.236 Å was 23



100

8 n ∆EEN Erf ∆Gel s

7

64

Ratio

Ratio

6 5

49

4

36 25

3 2

n2 ∆EEN Erf ∆Gel s

81

2

3

4

5

6

7

16 9 4

8

2

3

4

5

n

(a)

7

8

9

10

7

8

9

10

(b)

10−5

10−5

|∆QASC | (e)

10−6

10−7

6

n

|∆QASC | (e)

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

Page 24 of 48

10−6

2

3

4

5

6

7

10−7

8

2

3

4

5

6

n

n

(c)

(d)

Figure 9: Ratio of the energetic terms involved in Eq. (17) when compared to the data computed for the reference unit cell for supercells of (a) a glycine polymer chain with n unit cells and (b) a 3 layers MgO (001) slab supercell with (n × n) unit cells; (c) and (d) corresponding mean error on ASC during the whole SCRF before renormalization (see Eq. (15)). All data obtained with 65 grid points along the non periodic directions, and a grid spacing of 0.5 Å. UFF radii 101 are used in the MgO case. obtained with a 6-311G(d, p) basis set. Table 2 reports all computed total solvation energies. Overall, for all systems considered, a good agreement with previously-published data is obtained. In particular, when going from n = 1 to n = 3, the solvation energy increase is well 24


Page 25 of 48 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


reproduced and the value obtained for the infinite chain is in line with the C-PCM data. 48 As expected, the SCRF increases the polarization of the oligomers and of the infinite chain, resulting in larger ∆Gtot s when a full SCRF treatment is applied than when only a static response of the solvent to the solute charge distribution is considered (zeroth order reaction field, ZORF). The large differences obtained on ∆Gtot s between the oligomers and the infinite polymer also outline that a periodic approach is mandatory to correctly account for solvation of infinite periodic systems. In fact, only the difference between ∆Gtot s data computed for oligomers of consecutive lengths (data in parenthesis, see Table 2) slowly decreases towards the value obtained for the infinite periodic system. Table 2: Total solvation energies (in kcal/mol) of CH3 CO-(NHCH2 CONHCH2 CO)n -NHCH3 oligomers with increasing chain lengths and corresponding value for the infinite periodic polymer. For the periodic case, the data reported is that of the C4 N2 O2 H6 repeating unit. A 0.5 Å grid spacing is used along all directions, and the grid sizes considered are (65×65×65), (97×97×97), (127×127×127), (158×158×158), (189×189×189) and (14×65×65) for the n = 1, 2, 3, 4, 5 and periodic (+∞) calculations, respectively. SCRF and ZORF correspond to self-consistent and zeroth order reaction field, respectively. Data in parenthesis correspond to the difference in total solvation energies between oligomers of consecutive lengths. n This work SCRF ZORF

1

2

3

a b

5

-18.13

-27.81 -37.11 -46.23 -55.24 (-9.68) (-9.30) (-9.12) (-9.01) -17.01 -26.53 -35.70 -44.70 -53.59 (-9.52) (-9.17) (-9.00) (-8.89)

Previous works ZORFa -20.63 SCRFb

4

-30.03 -38.99 (-9.40) (-8.96) -22.11 – –

– – –

– – –

+∞ -8.38 -8.20 -8.85 -10.28

: C-PCM using Molecular Mechanics and static AMBER atomic charges, from Ref. 48. : COSMO, from Ref. 49.

4.3.3

Wettability of surfaces

A large research effort has been devoted to the determination of the wetting behaviour of surfaces at the experimental level in order to obtain specific wetting properties 102 targeted 25



Page 26 of 48

at various applications both at small and large scales. 103 At the most basic level, static wetting is described by Young’s equation which relates the surface tensions of solid-liquid (γSL ), solid-gas (γSG ) and liquid-gas (γLG ) to the equilibrium contact angle (θC ) the liquid drop makes with a surface: γSL − γSG + γLG cos θC = 0.

(21)

Partial wetting is obtained when γSG < γSL + γLG while complete wetting is obtained when γSG = γSL + γLG . Based on the contact angle, two additional quantities can be derived. The spreading coefficient S characterizes the different wetting states of a surface. More precisely, it represents the surface free energy γSG relative to its value for complete wetting:

S = γSG − (γSL + γLG ) = γLG (cos θC − 1).

(22)

On the other hand, the work of adhesion (WA ) of the liquid-solid is a measure of the strength of the contact between the liquid and solid phases:

WA = γSG − γSL + γLG = γLG (cos θC + 1).

(23)

From a modeling viewpoint, since cos θC can be seen as the solvation energy per unit area (A) of a slab and can be obtained as: 49,54

cos θC = −

1 γLG

∆Gtot s · , 2A

(24)

the quantities S and WA can be directly derived from ∆Gtot s through equations (22) and (23), and the contact angle can be obtained in cases where partial wetting occurs (S ≤ 0). The surface tension of the liquid-gas (γLG ) is taken as 72.80 mN/m for water. Table 3 presents computed WA , S and θC values obtained at the PBE and B3LYP levels for SiO2 α-quartz (001) surfaces. These surfaces have been chosen due to their wellknown hydrophilic or hydrophobic characters depending on the surface termination, which is 26


Page 27 of 48 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


particularly important in the biology field for instance. 104 Despite using gas-phase optimized structures in the present work, values obtained at the PBE level are in good agreement with previously-published data. 54 In particular, the complete (S > 0) and partial wetting (S < 0) of the cleaved and reconstructed surfaces respectively are well recovered, both at the PBE and B3LYP levels. Computed solvation energies converge well with the slab thickness (see Fig. 10a). The higher solvation energy (∼ −6.6 kcal/mol) obtained for the cleaved surface than for the reconstructed one (∼ −1.3 kcal/mol) can be related to the undercoordinated Si and O atoms found at the surface of the former, with protuding O atoms highly exposed to the implicit solvent (see Fig. 10b), 105 resulting in a marked hydrophilic character and complete wetting. On the other hand, in the reconstructed surface, Si and O atoms at the surface are fully-coordinated (see Fig. 10b) showing six-membered rings with siloxane bonds which are known to be hydrophobic, 104,106,107 leading to a smaller solvation energy and only partial wetting. Table 3: Computed works of adhesion (WA , mJ/m2 ), spreading coefficients (S, mJ/m2 ) and contact angles (θC , degrees) for cleaved and reconstructed SiO2 α-quartz (001) surfaces in water, with 18 and 27 atomic layers, respectively. A grid filling of the solute of about 40% along the non-periodic direction is used.

cleaved reconstructedb a b

WA B3LYP PBE PBEa 180 156 181 94 90 76

S B3LYP PBE 35 11 -51 -56

PBE 33 -69

a

θC B3LYP PBE – – 73 77

PBEa – 87

: solvent-optimized structure, from Ref. 54. : from Ref. 105.

5

Conclusions

In this paper, we have presented the implementation of an implicit solvation model in the CRYSTAL code. The solvation energy is separated into two components: the electrostatic contribution arising from a self-consistent reaction field treatment obtained within a generalized finite-difference Poisson model, augmented by a non-electrostatic contribution propor27



0


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 48

−5

−10

−15

9

18

nL

27

36

(a)

(b)

Figure 10: (a) Computed ∆Gtot s for the cleaved (solid) and reconstructed (dashed line) slabs as a function of the number of atomic layers (nL ) obtained at the PBE (blue) and B3LYP (black) levels. Gas-phase optimized structures are considered; (b) side views of the topmost atomic layers of the cleaved (top) and the reconstructed (bottom) SiO2 α-quartz slabs. Beige and red spheres correspond to Si and O atoms, respectively. tional to the solvent-accessible surface area of the solute. A discontinuous dielectric boundary is used, along with a solvent-excluded surface built from interlocking atom-centered spheres on which apparent surface point charges are mapped. The procedure is general and can be performed at the Hartree-Fock and Density Functional Theory levels, both with pure and hybrid functionals, for molecules and systems periodic in 1 (polymers) and 2 directions (surfaces). An overall good agreement between computed and reference solvation energies of analytically solvable model systems and a large test set of 501 neutral molecules with available experimental solvation energies has been obtained. In addition, data computed on infinite periodic systems are in line with previous works. We stress that the proposed scheme is general and can be applied both at the HF and DFT levels, using either pure or hybrid functionals. These encouraging results constitute a first step to take into account complex environ28


Page 29 of 48 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


ments in the CRYSTAL code, potentially allowing for a more accurate modeling of complex processes for both periodic and non-periodic systems. Implementation of atomic and cell gradients for the present implicit solvation model is in progress and will be discussed in a forthcoming paper.

A

Electrostatic potential at a grid point

A.1

Cubic grid

Integrating the PE (Eq. (2)) over volume yields: ZZZ

3

V

∇ · [(r) · ∇φ(r)]d r +

ZZZ

(25)

4πρint (r)d3 r = 0 V

For a cubic grid with grid spacing h, the second term can be approximated by 4πq0 . By using the Gauss’s theorem, the first volume integral can be transformed to a surface integral, which is in turn expanded as a sum of six surface integrals corresponding to the six facets of the cube (S1 → S6 ):

ZZZ

3

V

∇·[(r)·∇φ(r)]d r =

ZZ

(r)·∇φ(r)dS = S

ZZ

S1

(r)·∇φ(r)dS1 +· · ·+

ZZ

(r)·∇φ(r)dS6 S6

(26)

Evaluation of the first surface integral can be performed by considering surface S1 and a function g : (y, z) → x: ZZ

S1

(r) · ∇φ(r)dS1 = =

ZZ

S1

ZZ

D1

(x, y, z) · ∇φ(x, y, z)dS1

s

(g(y, z), y, z) · ∇φ(g(y, z), y, z) 1 +

∂g ∂y

2

+

∂g ∂z

2

dA (27)

29



Page 30 of 48

From Fig. 11: g(y, z) = x = h since the plane associated to S1 has equation x = h, hence ∂g ∂g = = 0. ∂y ∂z h

z

φ0 , q0 y

φ1 1

S1 x

Figure 11: Surface S1 for a cubic grid of spacing h. Only one neighboring cube is shown as thin gray lines. By using a centered finite-difference form for ∇ at (i, j, k), Eq. (27) can then be rewritten as: ZZ

S1

(r) · ∇φ(r)dS1 =

ZZ

D1

ZZ

φ(h + h/2) − φ(h − h/2) dA 2 · h/2 D1 ZZ φ1 − φ0 = (h, y, z) · dA h D1 ZZ φ1 − φ0 = (h, y, z)dA h D1 ZZ φ1 − φ0 = 1 dA h D1 Z hZ h φ1 − φ0 1 dydz = h 0 0 φ 1 − φ0 2 = 1 h = h1 (φ1 − φ0 ) (28) h

(h, y, z) · ∇φ(h, y, z)dA =

(h, y, z) ·

Doing similarly for the other five facets of the cube and summing: ZZ

S

(r) · ∇φ(r)dS = h1 (φ1 − φ0 ) + · · · + h6 (φ6 − φ0 ) = h

30


6 X n=1

n (φn − φ0 )

(29)

Page 31 of 48 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


By using Eq. (25), solving for φ0 and reorganizing, Eq. (3) is obtained:

h

6 X n=1

6 X n=1

hence

φ0 =

n φn + 4πq0 /h − φ0

X

!

n φn + 4πq0 /h /

n

A.2

(φn − φ0 ) + 4πq0 = 0 6 X

n = 0

n=1

X

n

n

Extension to any uniform grid shape

For a uniform grid with spacings (hx , hy , hz ) and inclination angles (α, β, γ) of the three axes along a, b and c, the elementary volume is: p V = (hx hy hz ) (1 − cos2 α − cos2 β − cos2 γ) + 2 cos α cos β cos γ

(30)

Eq. (27) now becomes: ZZ

S1

(r) · ∇φ(r)dS1 =

ZZ

ZZ

D1

(hx , y, z) · ∇φ(hx , y, z)dA

φ(hx + hx /2) − φ(hx − hx /2) dA 2 · hx /2 D1 ZZ φ1 − φ0 = dA (hx , y, z) · hx D1 ZZ φ1 − φ 0 = (hx , y, z)dA hx D1 ZZ φ1 − φ0 = 1 dA hx D1 φ 1 − φ0 = 1 · · sin α hy hz hx =

(hx , y, z) ·

(31)

since we are working in the (yz) plane. Similarly, S2 is also obtained in the (yz) plane, S3 and S4 are in the (xz) one, while S5 and S6 are in the (xy) one. Hence, Eq. (25) can be rewritten as:

31



Page 32 of 48

hx hz hx hy hy hz + (3 φ3 + 4 φ4 ) sin β + (5 φ5 + 6 φ6 ) sin γ 0 = (1 φ1 + 2 φ2 ) sin α hx hy hz hx hz hx hy hy hz + (3 + 4 ) sin β + (5 + 6 ) sin γ − φ0 (1 + 2 ) sin α hx hy hz (32)

+ 4πq0 = 0 Solving for φ0 : "

2 4 6 hy hz X hx hz X hx hy X φ0 = sin α i φi + sin β i φi + sin γ i φi + 4πq0 hx i=1 hy i=3 hz i=5 # " 2 4 6 hx hz X hx hy X hy hz X i + sin β i + sin γ i / sin α hx i=1 hy i=3 hz i=5

# (33)

For a cubic grid (hx = hy = hz and α = β = γ = 90◦ ), this reduces to Eq. (2).

B B.1

Apparent surface charge derivation Cubic grid

From Gauss’s law: ∇E = 4πσ tot

(34)

where E is the electric field (= −∇φ) and σ tot is the total charge density. Developing: E = −∇φ → ∇E = −∇2 φ = −

32

∂ 2φ ∂ 2φ ∂ 2φ − 2 − 2 = 4πσ tot ∂x2 ∂y ∂z


(35)

Page 33 of 48 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


Using a centered finite-difference form for ∇2 at (i, j, k) and Fig. 11: ∂ 2φ 2φ0 − φ1 − φ2 = ∂x2 h2 2 ∂ φ 2φ0 − φ3 − φ4 − 2 = ∂y h2 ∂ 2φ 2φ0 − φ5 − φ6 − 2 = ∂z h2 −

(36)

Hence, for any boundary grid point (i, j, k):

tot

2

4πσ (i, j, k)h = 6φ0 − 4πq tot (i, j, k)h2 /h3 = 6φ0 − h q tot (i, j, k) = 4π 3h q tot (i, j, k) = 4π q tot (i, j, k) =

3h 2π

6 X

n=1 6 X

φn φn

n=1

6φ0 −

6 X n=1

φn

!

! 6 1X 2φ0 − φn 3 n=1 ! 6 1X φ0 − φn 6 n=1 (37)

Substracting the atomic charges component (q0 (i, j, k)) obtained by the grid charging step, the apparent component of the charge at boundary grid point (i, j, k) is obtained: 3h q asc (i, j, k) = q tot (i, j, k) − q0 (i, j, k) = −q0 (i, j, k) + 2π

33


6

1X φ0 − φn 6 n=1

!

(38)


B.2

Page 34 of 48

Extension to any uniform grid shape

Since Eq. (36) is now: ∂ 2φ 2φ0 − φ1 − φ2 = ∂x2 h2x 2 ∂ φ 2φ0 − φ3 − φ4 − 2 = ∂y h2y −

−

∂ 2φ 2φ0 − φ5 − φ6 = , 2 ∂z h2z

(39)

for any boundary grid point (i, j, k), we have: 4π · σ tot (i, j, k) · h2x h2y h2z = h2y h2z (2φ0 − φ1 − φ2 ) + h2x h2z (2φ0 − φ3 − φ4 ) + h2x h2y (2φ0 − φ5 − φ6 ) 4π · h2x h2y h2z · q tot (i, j, k)/V = h2y h2z (2φ0 − φ1 − φ2 ) + h2x h2z (2φ0 − φ3 − φ4 ) + h2x h2y (2φ0 − φ5 − φ6 ) q tot (i, j, k) =

V h2y h2z (2φ0 − φ1 − φ2 ) + h2x h2z (2φ0 − φ3 − φ4 ) + h2x h2y (2φ0 − φ5 − φ6 ) · 4π h2x h2y h2z (40)

The apparent component of the charge at boundary grid point (i, j, k) is therefore: q asc (i, j, k) =

V 2 2 2 2 2 2 · h h (2φ − φ − φ ) + h h (2φ − φ − φ ) + h h (2φ − φ − φ ) −q0 (i, j, k). 0 1 2 0 3 4 0 5 6 x z x y 4πh2x h2y h2z y z (41)

For a cubic grid (hx = hy = hz and α = β = γ = 90◦ ), this reduces to Eq. (12).

Supporting Information Available The following file is available free of charge: • paper_fdpb.R1.SI.pdf: additional calculations performed with different grid sizes on the five Kirkwood model systems and the CH3 CONH2 molecule. This material is available free of charge via the Internet at http://pubs.acs.org/.

34


Page 35 of 48 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) Adamo, C.; Jacquemin, D. The calculations of excited-state properties with TimeDependent Density Functional Theory. Chem. Soc. Rev. 2013, 42, 845–856. (2) Traiphol, R.; Sanguansat, P.; Srikhirin, T.; Kerdcharoen, T.; Osotchan, T. Spectroscopic Study of Photophysical Change in Collapsed Coils of Conjugated Polymers: Effects of Solvent and Temperature. Macromolecules 2006, 39, 1165–1172. (3) Dyson, P. J.; Jessop, P. G. Solvent effects in catalysis: rational improvements of catalysts via manipulation of solvent interactions. Catal. Sci. Technol. 2016, 6, 3302– 3316. (4) D’Agostino, C.; Feaviour, M. R.; Brett, G. L.; Mitchell, J.; York, A. P. E.; Hutchings, G. J.; Mantle, M. D.; Gladden, L. F. Solvent inhibition in the liquid-phase catalytic oxidation of 1,4-butanediol: understanding the catalyst behaviour from NMR relaxation time measurements. Catal. Sci. Technol. 2016, 6, 7896–7901. (5) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H.-B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678–5796. (6) Curutchet, C.; Muñoz Losa, A.; Monti, S.; Kongsted, J.; Scholes, G. D.; Mennucci, B. Electronic Energy Transfer in Condensed Phase Studied by a Polarizable QM/MM Model. J. Chem. Theory Comput. 2009, 5, 1838–1848. (7) Olsen, J. M.; Aidas, K.; Kongsted, J. Excited States in Solution through Polarizable Embedding. J. Chem. Theory Comput. 2010, 6, 3721–3734. (8) Pisani, C.; Corà, F.; Nada, R.; Orlando, R. Hartree-Fock perturbed-cluster treatment of local defects in crystals: I. The EMBED program: general features. Comput. Phys. Comm. 1994, 82, 139 – 156. 35



(9) Wilbraham, L.; Adamo, C.; Labat, F.; Ciofini, I. Electrostatic Embedding To Model the Impact of Environment on Photophysical Properties of Molecular Crystals: A Self-Consistent Charge Adjustment Procedure. J. Chem. Theory Comput. 2016, 12, 3316–3324. (10) Wilbraham, L.; Louis, M.; Brosseau, A.; Guillot, R.; Ito, F.; Labat, F.; Métivier, R.; Allain, C.; Ciofini, I. Revealing the Origins of Mechanically Induced Fluorescence Changes in Organic Molecular Crystals. Adv. Mater. 2018, 0, 1800817. (11) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999–3094. (12) Mennucci, B., Cammi, R., Eds. Continuum Solvation Models in Chemical Physics: From Theory to Applications; John Wiley & Sons, Ltd, 2007. (13) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Acc. Chem. Res. 2008, 41, 760–768. (14) Klamt, A.; Mennucci, B.; Tomasi, J.; Barone, V.; Curutchet, C.; Orozco, M.; Luque, F. J. On the Performance of Continuum Solvation Methods. A Comment on "Universal Approaches to Solvation Modeling". Acc. Chem. Res. 2009, 42, 489–492. (15) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Performance of SM6, SM8, and SMD on the SAMPL1 Test Set for the Prediction of Small-Molecule Solvation Free Energies. J. Phys. Chem. B 2009, 113, 4538–4543. (16) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Focusing of electric fields in the active site of Cu-Zn superoxide dismutase: effects of ionic strength and amino-acid modification. Proteins: Struct., Funct., Genet. 1986, 1, 47–59. (17) Gilson, M. K.; Sharp, K. A.; Honig, B. H. Calculating the electrostatic potential of

36


Page 36 of 48

Page 37 of 48 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


molecules in solution: Method and error assessment. J. Comput. Chem. 1988, 9, 327–335. (18) Gilson, M. K.; Honig, B. Calculation of the total electrostatic energy of a macromolecular system: Solvation energies, binding energies, and conformational analysis. Proteins: Struct., Funct., Bioinf. 1988, 4, 7–18. (19) Nicholls, A.; Honig, B. A rapid finite difference algorithm, utilizing successive overrelaxation to solve the Poisson-Boltzmann equation. J. Comput. Chem. 1991, 12, 435–445. (20) Luty, B. A.; Davis, M. E.; McCammon, J. A. Electrostatic energy calculations by a Finite-difference method: Rapid calculation of charge-solvent interaction energies. J. Comput. Chem. 1992, 13, 768–771. (21) Sitkoff, D.; Sharp, K. A.; Honig, B. Accurate Calculation of Hydration Free Energies Using Macroscopic Solvent Models. J. Phys. Chem. 1994, 98, 1978–1988. (22) 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. (23) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: Applications to the molecular systems and geometric objects. J. Comput. Chem. 2002, 23, 128–137. (24) Cortis, C. M.; Friesner, R. A. An automatic three-dimensional finite element mesh generation system for the Poisson-Boltzmann equation. J. Comput. Chem. 1997, 18, 1570–1590. 37



(25) Cortis, C. M.; Friesner, R. A. Numerical solution of the Poisson-Boltzmann equation using tetrahedral finite-element meshes. J. Comput. Chem. 1997, 18, 1591–1608. (26) Baker, N.; Holst, M.; Wang, F. Adaptive multilevel finite element solution of the Poisson-Boltzmann equation II. Refinement at solvent-accessible surfaces in biomolecular systems. J. Comput. Chem. 2000, 21, 1343–1352. (27) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. USA 2001, 98, 10037–10041. (28) Holst, M.; Baker, N.; Wang, F. Adaptive multilevel finite element solution of the Poisson-Boltzmann equation I. Algorithms and examples. J. Comput. Chem. 2000, 21, 1319–1342. (29) Dyshlovenko, P. Adaptive numerical method for Poisson-Boltzmann equation and its application. Comput. Phys. Comm. 2002, 147, 335 – 338, Proceedings of the Europhysics Conference on Computational Physics, Computational Modeling and Simulation of Complex Systems. (30) Miertuš, 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. (31) Zauhar, R. J.; Morgan, R. S. The rigorous computation of the molecular electric potential. J. Comput. Chem. 1988, 9, 171–187. (32) Juffer, A.; Botta, E. F.; van Keulen, B. A.; van der Ploeg, A.; Berendsen, H. J. The electric potential of a macromolecule in a solvent: A fundamental approach. J. Comput. Phys. 1991, 97, 144 – 171.

38


Page 38 of 48

Page 39 of 48 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


(33) Allison, S.; Tran, V. Modeling the electrophoresis of rigid polyions: application to lysozyme. Biophys. J. 1995, 68, 2261–2270. (34) Bharadwaj, R.; Windemuth, A.; Sridharan, S.; Honig, B.; Nicholls, A. The fast multipole boundary-element method for molecular electrostatics: An optimal approach for large systems. J. Comput. Chem. 1995, 16, 898–913. (35) Cammi, R. Molecular response functions for the polarizable continuum model: physical basis and quantum mechanical formalism; SpringerBriefs in Molecular Science; Springer: Dordrecht, 2013. (36) Scalmani, G.; Frisch, M. J. Continuous surface charge polarizable continuum models of solvation. I. General formalism. J. Chem. Phys. 2010, 132, 114110. (37) Lange, A. W.; Herbert, J. M. A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach. The Journal of Chemical Physics 2010, 133, 244111. (38) Scalmani, G.; Frisch, M. J. Comment on "A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach" [J. Chem. Phys. 133, 244111 (2010)]. The Journal of Chemical Physics 2011, 134, 117101. (39) Lange, A. W.; Herbert, J. M. Polarizable Continuum Reaction-Field Solvation Models Affording Smooth Potential Energy Surfaces. The Journal of Physical Chemistry Letters 2010, 1, 556–561. (40) Cossi, M.; Barone, V. Solvent effect on vertical electronic transitions by the polarizable continuum model. The Journal of Chemical Physics 2000, 112, 2427–2435. (41) You, Z.-Q.; Mewes, J.-M.; Dreuw, A.; Herbert, J. M. Comparison of the Marcus and

39



Pekar partitions in the context of non-equilibrium, polarizable-continuum solvation models. The Journal of Chemical Physics 2015, 143, 204104. (42) Improta, R.; Rega, N.; Aleman, C.; Barone, V. Conformational Behavior of Macromolecules in Solution. Homopolypeptides of α-Aminoisobutyric Acid as Test Cases. Macromolecules 2001, 34, 7550–7557. (43) Cancès, E.; Maday, Y.; Stamm, B. Domain decomposition for implicit solvation models. J. Chem. Phys. 2013, 139, 054111. (44) Lipparini, F.; Stamm, B.; Cancès, E.; Maday, Y.; Mennucci, B. Fast Domain Decomposition Algorithm for Continuum Solvation Models: Energy and First Derivatives. J. Chem. Theory Comput. 2013, 9, 3637–3648. (45) Gatto, P.; Lipparini, F.; Stamm, B. Computation of forces arising from the polarizable continuum model within the domain-decomposition paradigm. J. Chem. Phys. 2017, 147, 224108. (46) Klamt, A.; Schüürmann, 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. (47) Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995– 2001. (48) Cossi, M. Continuum solvation model for infinite periodic systems. Chem. Phys. Lett. 2004, 384, 179 – 184. (49) Delley, B. The conductor-like screening model for polymers and surfaces. Mol. Simul. 2006, 32, 117–123.

40


Page 40 of 48

Page 41 of 48 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


(50) Gale, J. D.; Rohl, A. L. An efficient technique for the prediction of solvent-dependent morphology: the COSMIC method. Mol. Simul. 2007, 33, 1237–1246. (51) Fattebert, J.-L.; Gygi, F. Density functional theory for efficient ab initio molecular dynamics simulations in solution. J. Comput. Chem. 2002, 23, 662–666. (52) Fattebert, J.-L.; Gygi, F. First-principles molecular dynamics simulations in a continuum solvent. Int. J. Quantum Chem. 2003, 93, 139–147. (53) 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, 014103. (54) Fisicaro, G.; Genovese, L.; Andreussi, O.; Mandal, S.; Nair, N. N.; Marzari, N.; Goedecker, S. Soft-Sphere Continuum Solvation in Electronic-Structure Calculations. J. Chem. Theory Comput. 2017, 13, 3829–3845. (55) Yin, W.-J.; Krack, M.; Li, X.; Chen, L.-Z.; Liu, L.-M. Periodic continuum solvation model integrated with first-principles calculations for solid surfaces. Prog. Nat. Sci.: Mater. Int. 2017, 27, 283 – 288. (56) Sinstein, M.; Scheurer, C.; Matera, S.; Blum, V.; Reuter, K.; Oberhofer, H. Efficient Implicit Solvation Method for Full Potential DFT. J. Chem. Theory Comput. 2017, 13, 5582–5603. (57) 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. Europhys. Lett. 2011, 95, 43001. (58) Dziedzic, J.; Fox, S. J.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Large-scale DFT calculations in implicit solvent – A case study on the T4 lysozyme L99A/M102Q protein. Int. J. Quantum Chem. 2013, 113, 771–785. 41



(59) 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, 074103. (60) Scherlis, D. A.; Fattebert, J.-L.; Marzari, N. Stacking of oligo- and polythiophene cations in solution: Surface tension and dielectric saturation. J. Chem. Phys. 2006, 124, 194902. (61) 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. (62) Wang, H.-F.; Liu, Z.-P. Formic Acid Oxidation at Pt/H2 O Interface from Periodic DFT Calculations Integrated with a Continuum Solvation Model. J. Phys. Chem. C 2009, 113, 17502–17508. (63) Fishman, M.; Zhuang, H. L.; Mathew, K.; Dirschka, W.; Hennig, R. G. Accuracy of exchange-correlation functionals and effect of solvation on the surface energy of copper. Phys. Rev. B 2013, 87, 245402. (64) Mathew, K.; Sundararaman, R.; Letchworth-Weaver, K.; Arias, T. A.; Hennig, R. G. Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways. J. Chem. Phys. 2014, 140, 084106. (65) Cai, Q.; Hsieh, M.-J.; Wang, J.; Luo, R. Performance of Nonlinear Finite-Difference Poisson-Boltzmann Solvers. J. Chem. Theory Comput. 2010, 6, 203–211. (66) Boschitsch, A. H.; Fenley, M. O. Hybrid boundary element and finite difference method for solving the nonlinear Poisson-Boltzmann equation. J. Comput. Chem. 2004, 25, 935–955.

42


Page 42 of 48

Page 43 of 48 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


(67) Grant, J. A.; Pickup, B. T.; Nicholls, A. A smooth permittivity function for PoissonBoltzmann solvation methods. J. Comput. Chem. 2001, 22, 608–640. (68) Li, L.; Li, C.; Sarkar, S.; Zhang, J.; Witham, S.; Zhang, Z.; Wang, L.; Smith, N.; Petukh, M.; Alexov, E. DelPhi: a comprehensive suite for DelPhi software and associated resources. BMC Biophys. 2012, 5, 9. (69) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; Andrew, M. J. Electrostatics and diffusion of molecules in solution: simulations with the University of Houston Brownian Dynamics program. Comput. Phys. Commun. 1995, 91, 57 – 95. (70) Bashford, D. An object-oriented programming suite for electrostatic effects in biological molecules An experience report on the MEAD project. Scientific Computing in Object-Oriented Parallel Environments. Berlin, Heidelberg, 1997; pp 233–240. (71) Young, D. Iterative Solution of Large Linear Systems; Dover Books on Mathematics Series; Dover Publications, 2003. (72) Hadjidimos, A. Successive overrelaxation (SOR) and related methods. J. Comput. Appl. Math. 2000, 123, 177 – 199, Numerical Analysis 2000. Vol. III: Linear Algebra. (73) Tornberg, A.-K. The Ewald sums for singly, doubly and triply periodic electrostatic systems. Advances in Computational Mathematics 2016, 42, 227–248. (74) Parry, D. E. The electrostatic potential in the surface region of an ionic crystal. Surf. Sci. 1975, 49, 433–440. (75) Zhan, C.-G.; Bentley, J.; Chipman, D. M. Volume polarization in reaction field theory. J. Chem. Phys. 1998, 108, 177–192. (76) Cossi, M.; Menucci, B.; Pitarch, J.; Tomasi, J. Correction of Cavity-Induced Errors

43



in Polarization Charges of Continuum Solvation Models. J. Comput. Chem. 1998, 19, 833–846. (77) Decherchi, S.; Colmenares, J.; Catalano, C. E.; Spagnuolo, M.; Alexov, E.; Rocchia, W. Between Algorithm and Model: Different Molecular Surface Definitions for the Poisson-Boltzmann Based Electrostatic Characterization of Biomolecules in Solution. Comm. Comput. Phys. 2013, 13, 61–89. (78) Pang, X.; Zhou, H.-X. Poisson-Boltzmann Calculations: van der Waals or Molecular Surface ? Comm. Comput. Phys. 2013, 13, 1–12. (79) Connolly, M. L. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983, 221, 709–713. (80) Connolly, M. L. Analytical molecular surface calculation. J. Appl. Crystallogr. 1983, 16, 548–558. (81) Shrake, A.; Rupley, J. A. Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J. Mol. Biol. 1973, 79, 351–371. (82) Cramer, C. J.; Truhlar, D. G. General parameterized SCF model for free energies of solvation in aqueous solution. J. Am. Chem. Soc. 1991, 113, 8305–8311. (83) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378–6396. (84) Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G. Charge Model 5: An Extension of Hirshfeld Population Analysis for the Accurate Description of Molecular Interactions in Gaseous and Condensed Phases. J. Chem. Theory Comput. 2012, 8, 527–541. 44


Page 44 of 48

Page 45 of 48 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


(85) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; Pierre, M. D. L.; D’Arco, P.; Noël, Y.; CausàÂă, M.; Rérat, M.; Kirtman, B. CRYSTAL14: A program for the ab initio investigation of crystalline solids. Int. J. Quantum Chem. 2014, 114, 1287–1317. (86) Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rérat, M.; Casassa, S.; Baima, J.; Salustro, S.; Kirtman, B. Quantum-mechanical condensed matter simulations with CRYSTAL. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8, e1360. (87) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623. (88) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (89) Bondi, A. van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441–451. (90) Wang, B.; Wei, G. Accurate, robust and reliable calculations of Poisson-Boltzmann solvation energies. ArXiv e-prints 2016, 1607.04594. (91) Li, A. Performance comparison of Poisson-Boltzmann equation solvers DelPhi and PBSA in calculation of electrostatic solvation energies. J. Theor. Comput. Chem. 2014, 13, 1450040. (92) Born, M. Volumen und Hydratationswärme der Ionen. Zeitschrift für Physik 1920, 1, 45–48. (93) Lu, B. Z.; Zhou, Y. C.; Holst, M.; McCammon, J. A. Recent Progress in Numerical Methods for the Poisson-Boltzmann Equation in Biophysical Applications. Comm. Comput. Phys. 2008, 3, 973–1009. 45



Page 46 of 48

(94) Zhou, Y. C.; Feig, M.; Wei, G. W. Highly accurate biomolecular electrostatics in continuum dielectric environments. J. Comput. Chem. 2007, 29, 87–97. (95) Tomasi, J.;

Mennucci, B. In Self-consistent Reaction Field Methods;

von

Ragué Schleyer, P., Allinger, N., Clark, T., Gasteiger, J., Kollman, P., Schaefer, H., Schreiner, P., Eds.; Encyclopedia of Computational Chemistry; 2002. (96) Mobley, D. L.; Dill, K. A.; Chodera, J. D. Treating Entropy and Conformational Changes in Implicit Solvent Simulations of Small Molecules. J. Phys. Chem. B 2008, 112, 938–946. (97) Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J. Comput.-Aided Mol. Des. 2014, 28, 711– 720. (98) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. Correction to Small Molecule Hydration Free Energies in Explicit Solvent: An Extensive Test of Fixed-Charge Atomistic Simulations. J. Chem. Theory Comput. 2015, 11, 1347–1347. (99) 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. (100) Buša, J.; Hayryan, S.; Hu, C.-K.; Skřivánek, J.; Wu, M.-C. Enveloping triangulation method for detecting internal cavities in proteins and algorithm for computing their surface areas and volumes. J. Comput. Chem. 2009, 30, 346–357. (101) 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.

46


Page 47 of 48 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


(102) Durian, D. J.; Franck, C. Wetting phenomena of binary liquid mixtures on chemically altered substrates. Phys. Rev. Lett. 1987, 59, 555–558. (103) Bonn, D.; Eggers, J.; Indekeu, J.; Meunier, J.; Rolley, E. Wetting and spreading. Rev. Mod. Phys. 2009, 81, 739–805. (104) Rimola, A.; Costa, D.; Sodupe, M.; Lambert, J.-F.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 2013, 113, 4216–4313. (105) 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. (106) Rignanese, G.-M.; De Vita, A.; Charlier, J.-C.; Gonze, X.; Car, R. First-principles molecular-dynamics study of the (0001) α-quartz surface. Phys. Rev. B 2000, 61, 13250–13255. (107) Tosoni, S.; Civalleri, B.; Ugliengo, P. Hydrophobic Behavior of Dehydroxylated Silica Surfaces: A B3LYP Periodic Study. J. Phys. Chem. C 2010, 114, 19984–19992.

47



Graphical TOC Entry Ωout

Ωout

Ωin

Ωin

∇ · [(r) · ∇φ(r)] + 4π · ρ(r) = 0 1

48


Page 48 of 48

Implicit solvation using a generalized finite-difference approach in

Recommend Documents