A high order compact multigrid solver for implicit solvation models

Subscriber access provided by Iowa State University | Library

Condensed Matter, Interfaces, and Materials

A high order compact multigrid solver for implicit solvation models Arcesio Castaneda Medina, and Rochus Schmid J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00774 • Publication Date (Web): 16 Jan 2019 Downloaded from http://pubs.acs.org on January 17, 2019

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

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

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

A high order compact multigrid solver for implicit solvation models Arcesio Castañeda Medina∗ and Rochus Schmid∗ Ruhr-University, Bochum, Germany E-mail: [email protected]; [email protected] Abstract The electrostatic problem defined by the continuum solvation models used in molecular mechanics and ab initio molecular dynamics is solved in real space through multiscale methods. First, the Poisson equation is rewritten as a stationary convection-diffusion equation and discretized by a general mesh size fourth-order compact difference scheme. Then, the linear system associated to such discrete version of the elliptic partial differential equation is solved by a parallel (geometric) multigrid solver whose convergence rates and robustness are improved by an iterant recombination technique in which the multigrid acts as a preconditioner of a Krylov subspace method. The numerical tests performed on ideal and physical systems described by linear Poisson equations under different boundary conditions show the good performance of this accelerated multigrid solver. Furthermore, non-linear Poisson equations, like the regular modified PoissonBoltzmann equation, can also be solved by using in addition simple iterative schemes.

1

Introduction

The first principles investigation of solvation processes is of fundamental relevance for different fields of chemistry, physics and biology, not only because their presence in many daily 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

life phenomena, but also due to the technological implications. Despite of the enormous increase in the computational resources and the theoretical and numerical approaches of the last decades, the atomistic study of solvation poses several challenges. In particular, it is well known 1 that thermodynamic properties are physically meaningful only if a rather large number of solvent molecules and sufficiently long molecular dynamics simulations are considered. This imposes serious limitations on ab initio molecular dynamics (AIMD) investigations of solvation processes, and to a lesser extent, also on its classical force-field based counterpart. It has been realized, however, that a mean field approach can be a reasonable and much cheaper alternative to explicit solvation models. In these so-called implicit models, the molecular electrolyte is replaced by a continuum, which is defined by the dielectric properties of the embedding solvent and by the concentrations of the ionic species. The solvation process is visualized as the interplay between the polarization charge induced in the solute-electrolyte interface and the free charge density inside the solute cavity. In the last years, such continuum solvation models (CSMs) have became popular in the field of AIMD, since accurate thermodynamic properties can be predicted with a moderate numerical effort and a relatively small number of adjustable parameters. 2,3 The main task in CSMs is the solution of the electrostatic problem, defined by the charge densities and the dielectric permittivities of the solute-electrolyte system. Different levels of description of the electrostatic effects are formally accounted for by a family of linear and non-linear elliptic partial differential equations (EPDE) similar to the stationary convectiondiffusion equation. Thus, for a pure solvent, the electrostatic potential is defined by a linear generalized Poisson equation (GPE), whereas for an electrolyte, a non-linear modified Poisson-Boltzmann equation (MPBE), or its linearized version (LMPBE), are commonly employed. 4 The dielectric media, on the other hand, are assumed to be static, linear and isotropic. Typically, the permittivity is an ad hoc function, changing smoothly or discontinuously from one (inside the solute) to the bulk’s permittivity of the solvent (deep inside the electrolyte). In the AIMD context, such transition is modeled by a functional dependence

2


Page 2 of 26

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


with respect to the electronic charge density and/or the ionic positions of the solute. The proper physical description of the dielectric function, however, is a problem on its own, as the permittivity can depend on the potential, the dipole distribution of the solvent in the interface region, and the concentrations of the ionic species in the electrolyte. 5 In the field of numerical computation, the traditional approaches to solve the convectiondiffusion equation are based on finite differences (FD), finite elements or finite volumes methods. 6 In the area of first principle molecular dynamics, on the other hand, several methods of solution have been developed, commonly, depending on the approach followed to perform the ab initio calculations. Different variants of the finite difference approach have been proposed, 2,7–9 as well as, functional-minimization, 4,10 multiwavelet based 11 and semi-analytical methods. 12,13 In a FD approach, the approximations employed to represent the EPDE in a grid determine the accuracy of any numerical solver. For a long time, second-order FD schemes have been used to solve the convection-diffusion equations, although, it is well known special treatments can be required, in particular, due to the convection term. 6 It has been shown that higher order discretization methods avoid such kind of caveats, besides being more accurate and stable. 14 Two basic methods can be used to increase the accuracy of the discretization. 6 In the so-called defect correction, this is achieved by an improvement of the low order solution through a high order representation of the differential operators. A closely related method, the τ −extrapolation, has been used to increase the accuracy of a multigrid solver for the standard Poisson equation. 15 Very recently, a high-performance multigrid solver coupled to a defect correction scheme, aiming to solve the more general Poisson equations of CSMs, has been presented. 9 On the other hand, accuracy can also be increased by higher order discretizations of the EPDE. In principle, this just requires the extension of the low order finite difference expressions for the partial derivatives appearing in the EPDE. However, such an approach has the disadvantages of a more complicated treatment of the boundary conditions and an increase in the communication cost for parallel

3



implementations. The so-called high order compact difference schemes avoid such problems at the expense of more elaborated expressions in the discretization of the EPDE. 16–18 The use of a fourth-order compact scheme to solve the GPE in the AIMD context was already proposed in reference. 2 There, the original GPE was transformed, by a change of variables, into a Helmholtz equation with variable wavenumber. It is known, however, that the convergence of iterative methods for such an equation can only be guaranteed under very specific conditions. 6 As any finite difference approach, the high order compact discretization scheme leads to a system of linear equations, whose direct solution (e.g., by Gaussian elimination) is hardly attainable for real applications. The classical iterative schemes (e.g., by Gauss-Seidel) are not optimal options, either. Nevertheless, when they are combined with multiscale methods, as in the (geometric) multigrid technique, they lead to quite general and highly efficient solvers. 6 On the other hand, as the complexity of a problem increases, as it is the case for the GPE and the MPBE, the issue of selecting the multigrid parameters properly is not trivial. Complications as anisotropies, convection-dominance and nonlinearities can appear at the same time, affecting the central idea of multigrid. The combination of multigrid with Krylov subspace approaches can lead to robust and efficient numerical solvers in such difficult cases. 19,20 This procedure, also known as iterant recombination, 6 leads to a reduction of the number of steps required to achieve convergence and allows to solve also realistic convectiondiffusion problems. Although it has been already applied to solve 2D linear and nonlinear elliptic problems with the so-called upwind discretization schemes, 19,20 to our knowledge, its potential to solve the GPE and the MPBE in conjunction with high order compact discretizations has not been explored, yet. The paper is organized as follows. In Section 2 the continuum solvation models in density functional theory (DFT) calculations are briefly described. The electrostatic problem to be solved in the DFT-CSM and the main features of the actual implementation are discussed in Section 3, leaving out some technical and numerical details for the Appendix 6 and the

4


Page 4 of 26

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


Supplemental Information. The validation and computational efficiency of the accelerated multigrid technique are presented in Section 4. In Section 5 the performance of the solver in a Car-Parrinello approach to AIMD is analyzed. Finally, some conclusions are drawn in 6.

2

DFT with a continuum solvation model

Many of the modern continuum solvation models in density functional theory are based on the Fattebert-Gygi (FG) implicit solvation model. 2 In the FG-CSM, the solute cavity is defined by the electronic density and, in contrast to the conductor-like screening (COSMO) 21 or the polarizable continuum (PCM) models, 22 the dielectric transition at the solute-solvent interface is smooth. As a result, the ionic forces can be accurately computed and (microcanonical) molecular dynamics simulations are energy conserving. Alternative functional forms for the definition of the solute cavity, or equivalently a smooth dielectric permittivity, have been proposed, after the FG model 3,7,23–25 Some of these implicit models also take into account electrolyte solutions 24,26 and so-called cavitation energy terms. 23,27 The energy density functional for a system of electrons and nuclei, with densities ρe and ρI , respectively, embedded in a continuum electrolyte, with charge density ρc , is given by: 7

E[ρe ] = Ts [ρe ] + Ene [ρe ] + Exc [ρe ] + Ees [ρ] ,

(1)

where Ts , Ene and Exc are the in vacuo electronic kinetic, electron-nuclei and exchangecorrelation energies, respectively. The CSM is accounted for by the electrostatic energy term 1 Ees [ρ] = 2

Z drρ(r)φ[ρ(r), (r)] ,

(2)

defined by the total charge density ρ = ρe + ρI + ρc , the dielectric permittivity and the

5



Page 6 of 26

potential φ, which is the solution of the Poisson equation:

∇((r)∇φ(r)) = −4πρ(r) ,

(3)

under suitable boundary conditions. In a non-electrolytic embedding media (ρ = ρe + ρI ) only the dielectric effects need to be taken into account and the equation (3) is then referred as the generalized Poisson equation (GPE). 4 In the general case, for an electrolyte with M ionic species, however, the source term of this equation has the additional contribution:

ρc (r) = γ(r)

M X

zj cj (r) ,

(4)

j=1

where cj is the concentration of the species j with charge zj (in units of the elementary electronic charge e), and γ is a smooth arbitrary function introduced to avoid the existence of non-null cj (r) inside the solute’s cavity. For a Maxwell-Boltzmann distribution with bulk concentrations c∞ j and temperature T :

cj (r) =

c∞ j

zj φ(r) , exp − kB T

(5)

the expression in the eq. (3) leads to the Poisson-Boltzmann equation (PBE) of the GouyChapman model. 28 At low potentials, i.e. |zj φ(r)| kB T , a linearized form of the PoissonBoltzmann equation (LPBE), also known as the Debye-Hückel equation, is obtained:

∇((r)∇φ(r)) + (γ(r)∞ λ−2 D )φ(r) = −4πρ(r) + γ(r)Q ,

1 ∞ kB T

(6)

PM

2 ∞ j=1 zj cj

−1/2

where ∞ is the permittivity in the bulk electrolyte, λD = is the PM Debye length and Q = j=1 zj c∞ j is the total charge in the electrolyte. The PBE describes important aspects of the solute-electrolyte system, however, it leads to an overestimation of

6


Page 7 of 26 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 cj ’s in the interface region as the ions are considered as point particles. Improved models can be constructed by accounting for size correlation effects. For example, from a lattice-gas model of an electrolyte with equal sized ionic species one can derive: 29 zj φ(r) exp − kB T , c∞ zi φ(r) i exp − kB T + 1 cmax

c∞ j

cj (r) = 1+

PM

i=1

(7)

i

where cmax is the maximum allowed local concentration for the j-th species. In this case, j eq. (3) is referred to as modified Poisson-Boltzmann equation (MPBE). Both PBE and MPBE lead to a non-linear self-consistent problem. Once the electrostatic potential is known, the total energy of the system (eq. (1)) can be calculated straight forward. On the other hand, during an electronic ground state minimization or a molecular dynamics simulation, the forces acting on the Kohn-Sham (KS) single particle states and the (solute) ions need also to be computed. In general, the contributions of the CSM depend on the assumed functional form of the dielectric permittivity. Thus, in a soft sphere PCM there is no contribution to the effective KS potential, but just to the ionic forces. For a permittivity which depends only on the electronic density, as in the FG-CSM, there are only contributions to the KS potential. In this case, the variation f eq. (1), under the constraint imposed by the expression (3), leads to the effective KS potential:

veff (r) = v(r) + vxc (r) + φ(r) + v (r) ,

(8)

where v and vxc are the in vacuo external and exchange-correlation potentials and

v (r) = −

1 d |∇φ(r)|2 , 8π dρe

is the additional potential due to the CSM.

7


(9)


3

Page 8 of 26

High order compact discretization and accelerated parallel multigrid solver

The main computational task in a CSM is the solution of the generic Poisson equation (3), or equivalently, the (stationary) convection-diffusion equation:

(r)∇2 φ(r) + ∇(r) · ∇φ(r) = −4πρ(r) .

(10)

The compact discretization of an EPDE with constant coefficients (e.g., constant and ∇ in the above equation) is not much more complex than the well known non-compact FD schemes. Except for a slight modification (“smoothing”) of the source term, the finite difference representation of the EPDE (i.e., the so-called stencil) involve just the mesh spacings, h ≡ (hx , hy , hz ), and some predefined rational constants. The situation changes dramatically, however, if the coefficients (functions) appearing in the EPDE are not constant. In this case, the coefficients of the high order compact stencils depend also on the grid values of these functions. This is due to the principle behind the Mehrstellen discretization: 16 high order accuracy is gained at the expenses of approximating the partial derivatives by using the EPDE itself. In fact, the derivation of a fourth order compact stencil for the general convection-diffusion equation in grids with unequal meshsizes is much more elaborated than that of a non-compact scheme. 30 Thus, any possible reduction in the number of these functions results in a simpler high order compact discretization. As it is proposed in reference, 2 a substantially simplified version of eq. (10) can be obtained by performing the change of √ variable u = φ, 31 leading to

−∇2 u(r) +

∇2 1/2 (r) ρ(r) u(r) = 4π 1/2 . 1/2 (r) (r)

(11)

Although easier to discretize, this type of Helmholtz equation is numerically hard to solve by the standard iterative methods. To be more precise, its convergence can be guaranteed 8


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


only for the so-called positive definite case, 32 i.e. only for ∇2 1/2 (r) ≥ 0. As a consequence, the stability of a FD based algorithm depends strongly on the functional form of and the spacing h. 31 On the other hand, since (r) ≥ 1, a slight modification of eq. (10), is possible before its discretization:

Lφ(r) ≡ [−∇2 + ∇ ln −1 (r) · ∇]φ(r) = f (r) ≡ 4π

ρ(r) . (r)

(12)

This rather simple change reduces the complexity of the high order compact stencil (HOCS) considerably, which can be further manipulated algebraically in order to diminish the computational effort. This is of particular importance for the non-linear (self-consistent) problem, where the HOCS must be computed more than once. Detailed definitions and formulae of the fourth order compact stencil for eq. (12) are presented in the Supplemental Information (SI). Although slightly more involved, the compact scheme has important formal and practical benefits in comparison with non-compact high order discretizations. In particular, it leads to linear systems with lower bandwidth matrices, as well as, to a reduction of the communication load in parallel implementations and an easier treatment of the boundary conditions. 17 Furthermore, in the context of general multigrid solvers it allows to cover a wider range of convection dominated problems, which are difficult to deal with, when using non-compact low order discretization schemes. 14 In our case, the linear system obtained by the fourth order compact difference scheme is solved by a geometric multigrid solver (MG), 15 which is further accelerated by the iterant recombination technique (IR). 19,20 A brief description of the implemented MG and IR is presented in the Appendix 6. In a nutshell, the iterant recombination is a Krylov subspace technique in which the MG acts as a preconditioner. Its main purpose is to obtain an optimal input for the multigrid’s V-cycle using a linear combination of the actual solution and the most recent previous solutions. The IR can also be applied at a given stage of the V-cycle, i.e., at the coarser MG levels, diminishing the

9



Page 10 of 26

storage requirements and the computational load. 20 Here, however, the IR is performed only at the finest MG level. For a self-consistent optimization of the electronic states or an AIMD simulation, on the other hand, good initial inputs for the accelerated multigrid are computed by a second order extrapolation of two previous solutions. 33 On the practical side, the the driving routines of the solver are implemented in Python, whereas the heavy MG operations (smoothing, interpolation, restriction and overlapping) are written using the compiled language C using the MPI standard for communication. The parallel implementation follows a grid partitioning strategy, splitting the computational domain along the Cartesian directions and using one ghost layer to do the exchange of values between neighborhood processes.

4

Validation on test cases

In this section the accelerated multigrid solver is tested against some analytical solvable cases of the GPE, LMPBE and MPBE. The dielectric media is modelled by the permittivity: 4 ∞ − 1 (r) = 1 + 2

1 + erf

|r − r0 | − d0 ∆

,

(13)

where r0 is the origin, ∞ = 78.36, d0 = 1.7 and ∆ = 0.3. The source term in equation (12) is computed from the action of the differential operator L onto a normalized Gaussian potential (σ = 0.5): φexact (r) =

1 2πσ 2

3/2

(r − r0 )2 exp − . 2σ 2

(14)

The GPE under homogeneous (null) Dirichlet conditions is discretized in a cubic domain ((0, L)3 , L = 10) by a HOCS with h ≡ hx = hy = hz . Three grid-levels (ng = 3) are used in V(1,1) multigrid cycles whose coarsest grid defect equation is, approximately, solved by Nb = 100 Gauss-Seidel relaxation steps (see Appendix 6). The IR optimal inputs for each

10


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


V-cycle are obtained from five iterants, l = 5, in the equation (18). A numerical solution at a given spacing, φh , is accepted as soon as the `2 -norm of the residual is below 10−12 . The corresponding error is defined by δh = ||φh − φexact ||∞ . Table 1 shows the results obtained for increasingly finer grid meshes. The errors obtained for both the accelerated (MGIR) and non-accelerated (MG) solver are equal within numerical precision. The accuracy order of the discretization, estimated from the errors at two consecutive mesh sizes ((hi /hi+1 )p ≈ δhi /δhi+1 ), is always close to the theoretical value of four. The time spent in the calculation of the stencil, tstencil , and the smoothing of the right-hand side of the GPE trhs , are independent of the IR and just the averages are shown in the table. The number of iterations in the MGIR, on the other hand, is reduced by a factor of two, approximately. More importantly, with the IR, the time to solve the GPE, tsolver , is also decreased, most noticeably for the finest meshes. Similar results are obtained when the grid is decomposed into more than one MPI task. In particular, apart from the solver, the most expensive part of the algorithm is the calculation of the high order stencil. For the MGIR, it amounts around 10% of the total computation time, ttotal = trhs + tstencil + tsolver . As it was mentioned before, this stems from the algebraic complexity of the HOCS. Note, however, this calculation scales such that the ratio tstencil /ttotal is more or less independent of the number of MPI processes. The LPBE (eqs. (6) and (12)) is tested for a neutral monovalent 1:1 (Q = 0) aqueous electrolyte with c∞ i = 0.1 mol/L and T = 300 K. The source term is computed from the LPBE, the Gaussian potential in eq. (14), and the exclusion function appearing in the extra linear term: 4

γ(r) =

(r) − 1 , ∞ − 1

(15)

The same general trends observed for the GPE solver are also found from the numerical solution of the LPBE (table 2). In particular, the IR reduces the number of iterations to

11



Page 12 of 26

Table 1: Error (δh ), accuracy order (p), number of required V-cycles (Ih ) and computational times of the HOCS-MG algorithm for the GPE (eqs. (13) and (14)) using nc MPI processes. Both, the standard (MG) and the multigrid preconditioned Krylov method (MGIR) are used. The discretization is done on grids with equal mesh size (N points along each Cartesian direction). nc 1

N 101 161 221 281 341

δh 3.72e-05 5.74e-06 1.64e-06 6.13e-07 2.64e-07

p – 4.0 4.0 4.1 4.4

trhs (s) 0.04 0.12 0.30 0.60 1.04

tstencil (s) 0.33 1.34 3.50 7.50 14.27

Ih 23 18 19 15 22

4

101 161 221 281 341

3.72e-05 5.74e-06 1.64e-06 6.09e-07 2.64e-07

– 4.0 4.0 4.1 4.3

0.017 0.042 0.088 0.17 0.29

0.086 0.34 0.89 1.88 3.56

23 18 19 15 22

non-IR tsolver (s) 3.71 9.98 31.86 52.67 132.41 1.39 3.42 8.50 13.50 37.12

Ih 9 9 9 9 10

IR tsolver (s) 1.95 8.12 20.30 40.66 81.84

9 9 9 9 9

0.83 2.60 6.25 12.53 22.16

one half and the corresponding computational times for the solver are below those of the non-accelerated MG. Actually, the solution of the LPBE is computationally not much more expensive than solving the GPE. The non-linear MPBE equation, with an ionic distribution Table 2: Error (δh ), accuracy order (p), number of required V-cycles (Ih ) and computational times of the HOCS-MG algorithm for the LPBE (eqs. (13), (14) and (15)) using four MPI processes (nc = 4). Both, the standard multigrid (MG) and the multigrid preconditioned Krylov method (MGIR) are used. The discretization is done on grids with equal mesh size (N points along each Cartesian direction). N 101 161 221 281 341

δh 3.71e-05 5.73e-06 1.63e-06 6.24e-07 2.86e-07

p – 4.0 4.0 4.0 4.0

trhs (s) 0.018 0.046 0.099 0.187 0.322

tstencil (s) 0.099 0.396 1.024 2.101 3.804

Ih 23 18 19 15 20

MG tsolver (s) 1.73 4.11 10.92 21.57 36.79

MGIR Ih tsolver (s) 11 0.92 10 2.80 10 6.67 10 13.24 10 24.89

given by eq. (7), is tested following the same approach as for the GPE and the LPBE. The maximum local concentration is computed from cmax = p( 4π Ri3 NA )−1 , where p = 0.74 is the i 3 close packing factor for particles of effective radius Ri = 3Å, and NA is Avogradro’s number. 34

12


Page 13 of 26 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 methods to solve the MPBE are by themselves a field of investigation, nevertheless, two main approaches can be distinguished: (i) those based on the multiscale approach (like the Full Approximation Scheme (FAS) 6 ), and (ii) those which rely on the solution of some kind of linearization of the full problem and a Newton-iterative scheme. In the last case, the accuracy of the method is delegated to a large extent to the solver for the linear problem. In this work, a rather simple iterative mixing scheme is used (cf. algorithm 3 in the reference 4 ). Given and ρm , at the m-th step, a solution φm+1 of eq. (12) is obtained through the GPE solver. Then, new particle distributions in the electrolyte are computed, cm+1 ≡ cj (φm+1 ), j and a linear mixing of the new and old ionic charge distributions is used to calculate the new source term ρm+1 . The cycle is finalized when the difference between consecutive values of the electrolyte concentrations, ||cm+1 − cm ||∞ , are bellow 10−10 . The results are shown in table 3. For a given mesh size (N ), the error in the potential is consistent with the HOCS order. On the other hand, for both the accelerated and the non-accelerated MG solvers, the total number of linear-mixing steps (I) is rather low. At each of these steps, the iterations required by the GPE solver decreases rapidly, and after the fourth mixing step just one iteration is needed. Note, however, at the first steps of the linear-mixing scheme, the MGIR solver needs less iterations. As a result, the total computational time for the scheme with a MGIR solver is lower than that with the standard MG solver, in particular, for large N . As a final validation test, the parallel performance of the GPE solver has been checked by using very fine discretization grids (table 4). The speedup sc is estimated from the ratio T1 /Tnc , where T1 and Tnc are the times required to solve with 1 and nc processes, respectively. Although, the speedup does not scale perfectly linear, which is expected due the communication overhead for the coarse grid correction steps, the performance of the algorithm increases significantly with the number of parallel tasks. As it is expected, the parallel scaling is better for finer grids, since the ratio of computation to communication is increasing. Also, the speedup of the multigrid preconditioned Krylov algorithm is higher in comparison with that of standard multigrid.

13



Page 14 of 26

Table 3: Error (δh ), accuracy order (p), number of iterations required by the linear-mixing scheme (I), GPE iterations (IGP E ) and total computational time (tT ) of the HOCS-MG algorithm for the MPBE (eqs. (3) and (4)) using four MPI processes. The dots in IGP E denote that just one iteration of the GPE is needed. Both, the standard (MG) and the multigrid preconditioned Krylov method (MGIR) are used. The discretization is done on grids with equal mesh size (N points along each Cartesian direction). MG N 101 161 221 281 341

δh 3.72e-05 5.74e-06 1.64e-06 6.26e-07 2.88e-07

p – 4.0 4.0 4.0 4.0

I 8 8 8 10 14

IGP E [23,12,7,4,2,. . . ] [18,6,5,3,. . . ] [19,7,4,3,. . . ] [15,8,5,2,. . . ] [22,11,5,2,. . . ]

tT (s) 3.99 10.59 26.49 61.96 154.65

I 8 8 8 10 12

MGIR IGP E [9,5,4,3,2,. . . ] [9,4,3,2,2,. . . ] [9,5,3,2,. . . ] [9,5,4,2,. . . ] [9,5,4,2,. . . ]

tT (s) 2.73 9.20 23.01 56.48 125.35

Table 4: Computation time for the HOCS-MG algorithm for the GPE and the corresponding speedup sc for two different grids N and different number of parallel tasks nc . Both, the standard (MG) and the multigrid preconditioned Krylov method (MGIR) are used. The calculations are done on a HPC cluster with nodes (16-core Intel-Xeon CPU, 64GB RAM) connected by an Intel OPA interface (100 Gbps-bandwidth). N

nc

341

16 32 48 64 80

401

16 32 48 64 80

tsolver (s) MG 27.05 12.90 15.50 12.11 8.77 51.58 27.05 26.05 25.84 20.00

sc 5.2 10.8 9.0 11.5 15.9 5.8 11.1 11.5 11.6 15.0

14

tsolver (s) MGIR 14.94 9.98 7.22 6.56 5.00 25.24 13.33 14.24 17.01 8.79


sc 5.9 8.9 12.3 13.5 17.7 7.1 13.4 12.5 10.5 20.3

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

AIMD with a continuum solvation model

In this section the performance of the accelerated multigrid solver in actual AIMD simulations is studied. Only pure solvents are considered and the dielectric effects are modelled by the permittivity: 2

FG [ρe (r)] = 1 +

∞ − 1 , 1 + (ρe (r)/ρe0 )2β

(16)

where ρe0 = 4.0 × 10−4 and β = 1.3 determine the cavity size and the smooth modulation of the dielectric, respectively. As this expression does not depend on the nuclei positions, there is no contribution of the CSM to the ionic forces, and the coupling between the quantum solute and the classical solvent is entirely due to the extra potential term in eq. (9). Both the optimization of the electronic ground state and molecular dynamics simulations are performed following a Car-Parrinello molecular dynamics (CPMD) scheme. 35 The required energies and forces are computed from a real space approach to KS-DFT, where the wavefunctions and electronic densities are represented on a discrete Cartesian grid with mesh spacings h = (hx , hy , hz ). 33 A six-order finite difference approximation is used for the Laplacian in the electronic kinetic energy operator. Troullier-Martins norm-conserving pseudopotentials 36 in their separable form 37 are used to describe the electron-ion interactions and the PBE functional 38 is used to account for the exchange-correlation energy. As a first test, the ground state of the protein 1y49 4 is computed using a damped dynamics method. The KS equations and the GPE are discretized in a rectangular domain such that the mesh sizes along each Cartesian direction are as similar as possible to each other, in order to keep the convergence properties of the multigrid technique. 6 Homogeneous Dirichlet conditions are used at the boundaries of the domain. Table 5 shows the computed electrostatic energy, Ees , and the electrostatic contribution to the solvation energy, ∆Gel sol , for several grid sizes. The last quantity is defined as the difference between the total energy of the embedded system and the energy of the same system in vacuum. As in the validation

15



Page 16 of 26

tests, the GPE solver is assumed to converge when the `2 -norm of the residual is below 10−12 . The optimization of the KS states, on the other hand, is completed when the change in the total energy is below 10−6 Ha. As it is expected, for all mesh sizes, the protein is stabilized by the embedding solvent. In each optimization step, 20% of the time is spent in the solution of the GPE. Note, however, this an average quantity and at the beginning of the optimization the number of iterations required by the solver are comparatively higher than those needed when the ground state is almost reached. Table 5: Electrostatic energy, Ees , and electrostatic contribution to the solvation energy, ∆Gel sol , for the 1y49 protein after optimization of the KS states. The GPE is solved by the MGIR in a rectangular domain with grids of mesh sizes h (numbers of points N). Also the ratio between the (mean) time spent by the accelerated multigrid algorithm and the (mean) total time of each SCF cycle is shown. The calculations have been done using 48 MPI tasks. All quantities are given in atomic units. N

h

Ees

∆Gel sol

t¯GP E /t¯step

(169,177,169) (137,137,129) (113,113,113) (89,89,81)

(0.205,0.198,0.197) (0.252,0.254,0.257) (0.305,0.308,0.293) (0.387,0.391,0.409)

3.807 3.813 3.842 3.929

-0.159 -0.164 -0.240 -0.255

0.187 0.197 0.168 0.230

Previous studies on the standard Poisson equation showed the relevance of a sufficiently accurate solution of the electrostatic problem to obtain proper conserved dynamics. 39 Therefore, as a second test, the dynamics of a Germanium (100) surface in contact with water is considered. The surface is modelled by a slab with three monolayers of Ge. Both, the bottom and the top of the slab are saturated by Hydrogen, while the top is also explicitly and implicitly hydrated (fig. 1). Only the topmost layer of the hydrogen terminated slab and the explicit water is allowed to move. The GPE is discretized using periodic boundary conditions in the xy plane and homogeneous Dirichlet conditions in the z direction. It must emphasized that these conditions are also applied to the KS equations, i.e. to the electronic wavefunctions and densities. The real space approach to KS-DFT allows to deal with arbitrary boundary conditions, and for 2D periodic geometries there is no need to consider the interaction of replicas along the non-periodic direction. After minimization of the KS energy 16


Page 17 of 26

functional, microcanonical (NVE) simulations are performed. The electronic fictitious mass in the CPMD is 5000 a.u. and the equations of motion are integrated with a time step of 0.12 fs. Figure 1 shows that the MGIR solver leads to a perfectly energy conserving dynamics. In contrast to the optimization of the KS states, the number of iterations required to solve the GPE at each molecular dynamics step do not change too much during the whole simulation. For this test, in particular, the ratio t¯GP E /t¯step ≈ 10%.

0.08

Kions

0.06

U

0.04 0.02 E/Ha

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.00 0.02 0.04 0.06 0.08 0.0

0.5

1.0

1.5

t/ps

2.0

2.5

3.0

3.5

Figure 1: An explicitly/implicitly hydrated Ge(100) slab with Hydrogen terminations. On the left, the isosurface of the electronic density and a slice plot of the dielectric permittivity (eq. (15)) for an arbitrary molecular dynamics step are shown. Dirichlet (periodic) boundary conditions are used in the z (x, y) direction(s). The GPE is discretized in a grid with h = (0.315, 0.315, 0.309). On the right side, the energy components of a CPMD run are given: ionic kinetic energy (Kions ), potential (U ) and total energy (E). The last two are measured with respect to their initial values.

6

Conclusions

The common continuum solvation models used in first-principle methods lead to an electrostatic problem in a linear, but inhomogeneous dielectric medium. Formally, this problem is described by a generic Poisson equation with the same structure as that of the stationary convection-diffusion equation. Convection-dominated problems, or equivalently situa17



Page 18 of 26

tions with high polarization charges, are difficult to handle by the standard finite difference schemes and commonly require a special treatment. We have shown that the combination of a high order compact discretization and a multigrid preconditioned Krylov method lead to a general, stable and robust algorithm to solve such problems. The numerical results demonstrate that the solver is computationally efficient and sufficiently accurate to perform reliable ab initio molecular dynamics simulations. Especially in the context of a real space approach with numerically discretized wavefunctions and densities it can directly replace the multigrid solvers for the electrostatic problem and allow the investigation of solvated processes without a substantial additional numerical effort.

Appendix A: The accelerated multigrid technique The geometric multigrid technique is based on two basic ideas: 6 (i) the standard iterative methods used to solve linear system are good smoothers, in the sense they are able to reduce short-wavelength errors in the sought solution, and (ii) by employing a hierarchy of coarser grids, such that errors can systematically be reduced. At the m-th iteration of Lh φh = fh , m the error vhm is determined from the defect dm h , and then added to the actual φh in order to

obtain the next solution guess:

m+1 m m m m m φm = φm h → dh = fh − Lh φh → Lh vh = dh → φh h + vh .

(17)

Of course, the exact solution of the defect equation is not known, however, if it is restricted to a coarser grid, e.g. H1 = 2h, and the approximated smoothed error can then be interpolated back from H1 to H0 = h. Thus, it is possible to restrict (and then smooth and interpolate) n times (Hn = 2Hn−1 ) until an enough accurate solution at the coarsest grid is available. From this coarsening principle emerges the so-called V-cycle of the multigrid method (Fig. 2). The procedure starts by defining the discretized versions of the LHk operator and fH0 , as well as, an initial guess φ0H0 for the sought solution. The number of presmoothing and postsmoothing 18


Page 19 of 26 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


steps at each grid level are also predefined parameters, and at the coarsest grid, a sufficiently accurate solution is found by relaxation. Finally, the solution is accepted as soon the residual is below a given threshold, otherwise it is used as input for a new V-cycle. In the present implementation, a non-lexical sorted Gauss-Seidel procedure is employed as smoother for the finest and the coarser grids, using four and two colors, respectively. A second-order seven-point stencil representation of the partial differential equation is used h as the coarse grid operator. Full-weighting restriction IhH and (tri)linear interpolation IH

are used as transfer operators in the V(m,n) cycle with ng grid levels, and m (n) presmoothing (postsmoothing) steps. At the coarsest grid, the linear systems are solved by an adequate number of smoothing steps (Nb ). The convergence rates are enhanced by scaling the correction (error) at the finest grid. 14 When it is applied at the finest grid, the iterant recombination procedure is an outer loop of the multigrid V-cycle chain (Fig. 3). The optimal input solution for the next coarse grid correction is determined as a linear combination of the actual solution and the most recent l previous solutions: 6

φh,opt = φm h +

l X i=1

αi (φm−i − φm h ). h

(18)

The coefficients αi are determined from a defect minimization problem, which reduces to the solution of the linear system Hβ = α, defined by the (Euclidean) inner products of the residuals:

m m m−i i, βi =hdm h , dh i − hdh , dh m−k , dm−k i − hdm i. Hik =βi + hdm−i h , dh h h

(19)

In general, just a few of the previous iterants (2 ≤ l ≤ 15) are enough to obtain a well conditioned linear system and an optimal solution. 19,20 At the early stages, the recombination process can be performed with the iterants already available. The technique also fits to the 19



← φm φm−1 H0 H0 FALSE

LHk , φ0H0 , fH0

||fH0 − LH0 φm H0 ||∞ ≤ THRSH

m−1 = fH0 ) φm H0 ← PRESMOOTH(LH0 φH0

Page 20 of 26

TRUE

END

m POSTSMOOTH(LH0 φm H0 = fH0 )→ φH0

m dm H0 ← fH0 − LH0 φH0

m dm H1 ←RESTRICT(dH0 )

m m φm H0 + α×INTERPOLATE(vH1 )→ φH0

m vH =0 1 m m m vH ←PRESMOOTH(L H1 vH1 = dH1 ) 1 m m dm ← d − L v H 1 H1 H1 H1

m m POSTSMOOTH(LH1 vH = dm H1 )→ vH1 1

m m m vH +INTERPOLATE(vH )→ vH 1 2 1

m dm H2 ←RESTRICT(dH1 )

m vH 2

m vH =0 2 m ←SOLVE(LH2 vH = dm H2 ) 2

Figure 2: Geometric multigrid iterative solver algorithm with three grid levels (Hn = 2Hn−1 for 1 ≤ n ≤ 2 and H0 = h). The plots show 2D projections of the potential (after the postsmoothing steps) for the ground state of a molecule of H2 O.

20


Page 21 of 26 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


parallel features of the multigrid solver, as the required overlaps are computed following the grid decomposition scheme, and the data is accessible to each one of the MPI processes. LHk , φ0H0 , fH0 ld

ld

ld

b

b

⊕ qp

b

ld

⊕

qp

ld

⊕ qp

b

b

b

l+1

Figure 3: Schematic representation of the multigrid grid iterant recombination process. The actual and the most recent l solutions () are linearly combined (⊕) to be used as an optimal input guess (D) for the m-th V-cycle (m ≥ l) with smoothing steps ( ) and coarsest grid solution (#).

Acknowledgement This work is supported by the Cluster of Excellence RESOLV (EXC 1069) funded by the Deutsche Forschungsgemeinschaft (DFG) within the framework of the German Excellence Initiative and by the DFG priority program 1928 (COORNETs).

Supporting Information Available Some more detailed information about the high order discretization scheme can be found in the supplemental material. • hocs.pdf: Definitions used for the O(h4 ) compact discretization and stencils formulae.

References (1) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press, 2009. 21



(2) Fattebert, J.-L.; Gygi, F. Density Functional Theory for Efficientab Initio Molecular Dynamics Simulations in Solution. J. Comput. Chem. 2002, 23, 662–666. (3) Petrosyan, S. A.; Rigos, A. A.; Arias, T. A. Joint Density-Functional Theory: Ab Initio Study of Cr

2

O

3

Surface Chemistry in Solution. J. Phys. Chem. B 2005, 109,

15436–15444, 00000. (4) 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. (5) Pache, D.; Schmid, R. Molecular Dynamics Investigation of the Dielectric Decrement of Ion Solutions. ChemElectroChem 2018, 5, 1444–1450. (6) Trottenberg, U.; Oosterlee, C. W.; Schüller, A. Multigrid ; Academic Press, 2001. (7) 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. (8) Garcia-Ratés, M.; López, N. Multigrid-Based Methodology for Implicit Solvation Models in Periodic DFT. J. Chem. Theory Comput. 2016, 12, 1331–1341. (9) Womack, J. C.; Anton, L.; Dziedzic, J.; Hasnip, P. J.; Probert, M. I. J.; Skylaris, C.-K. DL_MG: A Parallel Multigrid Poisson and Poisson-Boltzmann Solver for Electronic Structure Calculations in Vacuum and Solution. J. Chem. Theory Comput. 2018, (10) Bani-Hashemian, M. H.; Brück, S.; Luisier, M.; VandeVondele, J. A Generalized Poisson Solver for First-Principles Device Simulations. The Journal of Chemical Physics 2016, 144, 044113. (11) Fosso-Tande, J.; Harrison, R. J. Implicit Solvation Models in a Multiresolution Multiwavelet Basis. Chem. Phys. Lett. 2013, 561-562, 179–184. 22


Page 22 of 26

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


(12) Otani, M.; Sugino, O. First-Principles Calculations of Charged Surfaces and Interfaces: A Plane-Wave Nonrepeated Slab Approach. Phys. Rev. B 2006, 73 . (13) Ringe, S.; Oberhofer, H.; Hille, C.; Matera, S.; Reuter, K. Function-Space-Based Solution Scheme for the Size-Modified Poisson–Boltzmann Equation in Full-Potential DFT. J. Chem. Theory Comput. 2016, 12, 4052–4066. (14) Gupta, M. M.; Zhang, J. High Accuracy Multigrid Solution of the 3D Convection–Diffusion Equation. Applied Mathematics and Computation 2000, 113, 249–274. (15) Köstler, H.; Schmid, R.; Rüde, U.; Scheit, C. A Parallel Multigrid Accelerated Poisson Solver for Ab Initio Molecular Dynamics Applications. Comput. Visual Sci. 2008, 11, 115–122. (16) Collatz, L. The Numerical Treatment of Differential Equations; Springer Science & Business Media, 2013. (17) Spotz, W. F.; Carey, G. F. High-Order Compact Scheme for the Steady StreamFunction Vorticity Equations. Int. J. Numer. Meth. Engng. 1995, 38, 3497–3512. (18) Zhang, J. An Explicit Fourth-Order Compact Finite Difference Scheme for ThreeDimensional Convection–Diffusion Equation. Commun. Numer. Meth. Engng. 1998, 14, 209–218. (19) Washio, T.; Oosterlee, C. W. Krylov Subspace Acceleration for Nonlinear Multigrid Schemes. Electron. Trans. Numer. Anal. 1997, 6, 3–1. (20) Oosterlee, C. W.; Washio, T. Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows. SIAM J. Sci. Comput. 2000, 21, 1670–1690. (21) Klamt, A.; Schuurmann, G.; AG, B. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. 7.

23



(22) 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. (23) Andreussi, O.; Dabo, I.; Marzari, N. Revised Self-Consistent Continuum Solvation in Electronic-Structure Calculations. J. Chem. Phys. 2012, 136, 064102, 00018. (24) Jinnouchi, R.; Anderson, A. Electronic Structure Calculations of Liquid-Solid Interfaces: Combination of Density Functional Theory and Modified Poisson-Boltzmann Theory. Phys. Rev. B 2008, 77 . (25) 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. (26) 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. (27) 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. (28) Newman, J.; Thomas-Alyea, K. Electrochemical Systems; Wiley, 2012. (29) Tresset, G. Generalized Poisson-Fermi Formalism for Investigating Size Correlation Effects with Multiple Ions. Phys. Rev. E 2008, 78, 061506. (30) Medina, A. C.; Schmid, R. Solution of High Order Compact Discretized 3D Elliptic Partial Differential Equations by an Accelerated Multigrid Method. Journal of Computational and Applied Mathematics 2019, 350, 343–352. (31) Concus, P.; Golub, G. Use of Fast Direct Methods for the Efficient Numerical Solution of Nonseparable Elliptic Equations. SIAM J. Numer. Anal. 1973, 10, 1103–1120. 24


Page 24 of 26

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


(32) Ernst, O. G.; Gander, M. J. Why It Is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. Numer. Anal. Multiscale Probl. 2012, 83, 325–363. (33) Schmid, R. Car-Parrinello Simulations with a Real Space Method. J. Comput. Chem. 2004, 25, 799–812, 00017. (34) López-García, J. J.; Horno, J.; Grosse, C. Poisson–Boltzmann Description of the Electrical Double Layer Including Ion Size Effects. Langmuir 2011, 27, 13970–13974. (35) Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A.; Joannopoulos, J. D. Iterative Minimization Techniques for \textit{ab Initio} Total-Energy Calculations: Molecular Dynamics and Conjugate Gradients. Rev. Mod. Phys. 1992, 64, 1045–1097. (36) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for Plane-Wave Calculations. Phys. Rev. B 1991, 43, 1993, 10591. (37) Kleinman, L.; Bylander, D. M. Efficacious Form for Model Pseudopotentials. Phys. Rev. Lett. 1982, 48, 1425–1428. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (39) Schmid, R.; Tafipolsky, M.; König, P. H.; Köstler, H. Car–Parrinello Molecular Dynamics Using Real Space Wavefunctions. Phys. Status Solidi B 2006, 243, 1001–1015, 00012.

25


1 2 3 e4h (ǫ(rh )∇ e h φ(rh )) = −4πe ∇ ρ(rh ) 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Multigrid solver + and iterant recombination Journal of Chemical Theory Computation

φ1

φ2

φi +=

l X

αj (φi−j − φi )

j=0

Grid size

High order compact discretization

Iterations ACS Paragon Plus Environment

φN

Potential Page 26 of 26

A high order compact multigrid solver for implicit solvation models

Recommend Documents