{'fraction_non_alphanumeric': 0.03463012425966787, 'fraction_numerical': 0.011682731390082452, 'mean_word_length': 4.680211081794195, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Goal-Oriented Adaptivity and Multilevel Preconditioning for the Poisson-Boltzmann Equation

In this article, we develop goal-oriented error indicators to drive adaptive refinement algorithms for the Poisson-Boltzmann equation. Empirical results for the solvation free energy linear functional demonstrate that goal-oriented indicators are not sufficient on their own to lead to a superior refinement algorithm. To remedy this, we propose a problem-specific marking strategy using the solvation free energy computed from the solution of the linear regularized Poisson-Boltzmann equation. The convergence of the solvation free energy using this marking strategy, combined with goal-oriented refinement, compares favorably to adaptive methods using an energy-based error indicator. Due to the use of adaptive mesh refinement, it is critical to use multilevel preconditioning in order to maintain optimal computational complexity. We use variants of the classical multigrid method, which can be viewed as generalizations of the hierarchical basis multigrid and Bramble-Pasciak-Xu (BPX) preconditioners.


INTRODUCTION
The Poisson-Boltzmann equation (PBE) is a widely used model for electrostatic interactions of charged bodies in dielectric media, such as molecules, ions, and colloids, and thus is of importance in many areas of science and engineering, including biochemistry, biophysics, and medicine. (See the classical texts [65,74] for a derivation of the PBE.) The importance of the PBE model is reflected by the popularity of software packages such as APBS [9], CHARMM [34], DelPhi [70], and UHBD [64], within the molecular modeling communities. It provides a high fidelity mean-field description of electrostatic interactions and ionic distributions of solvated biomolecular systems in equilibrium. The partial differential equation itself is challenging to solve numerically due to singularities of different orders at the positions of permanent point charges and the presence of a dielectric interface.
In this article, we develop an adaptive multilevel finite element method for the PBE using goal-oriented a posteriori error indicators. This adaptive algorithm, which is a variant of that studied for the PBE in [55,8,9], deviates substantially from previous work in that the error indicator is based on a user defined quantity of interest or goal. This is in contrast to traditional residual-based adaptive refinement algorithms (like those developed for the PBE in [39]) that drive-refinement to minimize the global error measured in an energy-norm. The goal-oriented refinement methodology has been successfully employed in a wide range of application areas, including fluids, elasticity, and fluid structure interaction [11]. Despite these successes, we show that this methodology applied directly to the PBE does not necessarily lead to a successful adaptive algorithm. To remedy this issue we propose a novel marking strategy which recovers the performance commonly seen in other applications. This is the first time that this particular goal-oriented refinement strategy has been applied to the PBE specifically, and molecular biophysics in general.
At the core of any adaptive finite element approach are the iterative methods used to solve the discretized equation. However, due to the ill-conditioning of the linear systems arising from the discretization of the PBE, the convergence rate of traditional iterative solvers is significantly deteriorated. To remedy this, we combine modern Bramble-Pasciak-Xu (BPX)-type multilevel preconditioners with the goal-oriented adaptive algorithm mentioned above. When applied to the PBE, our results demonstrate that the overall algorithm is accurate, highly efficient and scalable with respect to the number of levels in the adaptive hierarchy.
An outline of the article is as follows. In section 2, we give a brief overview of the Poisson-Boltzmann equation, and describe the most useful formulations for modeling and numerical simulation, such as the regularized formulations described in [39,56,37]. We also discuss the solvation free energy functional corresponding to a given reaction potential, which will form the basis of our goal-oriented error indicators developed later in the article. We describe adaptive finite element methods in section 3, including weak formulation of the regularized PBE, discretization by finite element methods, and adaptive algorithms driven by a posteriori error indicators. In section 4, we describe a particular class of error indicators known as goal-oriented indicators, and describe several indicators designed for the PBE. In section 5, we discuss a local multigrid algorithm used to precondition an iterative Krylov method for solving the linear systems arising from adaptive mesh refinement. With some care, these methods enable an algorithm whose complexity is close to optimal. The results from a sequence of numerical experiments

THE POISSON-BOLTZMANN EQUATION
The Poisson-Boltzmann equation (PBE) is a second-order nonlinear partial differential equation whose solution gives the electrostatic potential, φ(x), for a solute molecule immersed in an implicitly defined solvent. Using a mean-field approximation, the solvent is treated as a bulk medium where ions are distributed according to the Boltzmann distribution. Figure 1 is a schematic representation of the domain of the PBE, denoted by Ω. The innermost region, Ω m , contains the explicitly represented solute molecule. The outer region, Ω s , is the bulk solvent and contains the implicit solvent ions. Between Ω s and Ω m is the ion exclusion layer, which separates the solute from the solvent ions, and has a width dependent on the size of the solvent ions. For simplicity, we will assume the solvent ions are small and the ion exclusion layer can be neglected. Hence, the interface between the solute and solvent is a surface, denoted by Γ =Ω m ∩Ω s . The shape of the surface is governed by the short-range repulsive van der Waals interactions, which prevent the solvent from penetrating the solute. The precise definition of the surface varies depending on the model [10].
The PBE for a 1:1 electrolyte (e.g., sodium chloride) is where u(x) = e c φ(x)/k B T is the dimensionless potential, e c is the charge of an electron, k B is Boltzmann's constant, and T is the temperature. Here, · denotes the jump across the interface and n is the outward pointing normal of ∂Ω m . The dielectric function (x) jumps one or two orders of magnitude at the interface Γ. For example, commonly used values are (Ω m ) = m = 2 and (Ω s ) = s = 80. The modified Debye-Hückel parameter,κ, has a similar discontinuity, withκ(Ω m ) = 0 andκ(Ω s ) =κ s > 0. The fixed ions within the solute are represented by a sum of Dirac delta distributions, with fixed charge centers, x i , and charges e c q i . This charge distribution induces singularities in the electrostatic potential and has, until recently, proved to be difficult to treat numerically. To address this issue the PBE is reformulated so that the singularities are explicitly removed [48,86]. Following [39,56,37], this is accomplished by writing the potential as a sum of a singular term u c and a nonsingular remainder u r . The singular term is the Coulomb potential which satisfies the Poisson equation There are numerous fast algorithms (with linear or near linear complexity) for evaluating u c on a set of quadrature points, such as fast multipole [49], multilevel summation [72,51] and particle mesh Ewald [46]. Substituting u = u r + u c into the PBE (Eq. 2.1) gives a modified form of the PBE, which was termed in [39] as the Regularized PBE (RPBE): (2.5) Note that, because bothκ(x) and ( (x) − m ) are zero for x ∈ Ω m and the centers of the atoms in the solute are well separated from Γ, the singular function, u c , is never evaluated near the singularities. This formulation was used in [39] to develop continuous a priori L ∞ estimates of solutions, and subsequently to show existence and uniqueness of solutions of the PBE. Furthermore, the authors established discrete a priori L ∞ estimates for Galerkin solutions, making possible quasi-optimal a priori error estimates, as well as a provably convergent adaptive finite element method for the RPBE.
Recently, an alternative 3-term splitting of the PBE has been proposed in [56] which addresses the inherit subtractive cancellation in the reconstruction of the electrostatic potential used by the RPBE. In this article, the authors establish mathematical results for the alternative splitting, including continuous and discrete a priori L ∞ estimates, existence and uniqueness of solutions, quasi-optimal a priori error estimates, and a convergent adaptive finite element method (AFEM). (Whereas in [39] only AFEM convergence was shown, it was shown in [56] that AFEM is a contraction for the RPBE, using a new AFEM convergence framework for nonlinear problems developed in [57].) The 3-term splitting decomposes the electrostatic potential into where u c is the Coulomb potential, but here it is restricted to the subdomain Ω m . The harmonic term, u h , is defined as the solution to (2.8) Applying the definitions of u c , u h and substituting into the PBE (Eq. 2.1), one obtains an equation for u 3 , In contrast to the RPBE, this formulation avoids the subtractive cancellation since u 3 = u in Ω s , and the Coulomb term is not used to reconstruct the potential in the solvent subdomain.
For systems which are not highly charged, the variation in the potential is relatively small, and the hyperbolic sine term is well approximated by its linearization. This approximation, which replaces sinh(u) with u, results in what is known as the linear Poisson-Boltzmann equation (LPBE), or the linear regularized Poisson-Boltzmann equation (LRPBE) in case of the RPBE. Although this approximation reduces the ionic response of the solvent [10,47], it can significantly reduce the complexity of many numerical algorithms (e.g., boundary element and boundary integral methods) [63,25,62].
One important use for the solution to the PBE is in the calculation of solvation free energies. This quantity measures the thermodynamic work of moving the solute molecule from a vacuum to a solvent environment. The solvation free energy can be written as a sum of nonpolar and polar contributions. The nonpolar term depends on the solvent accessible surface area, excluded volume, and nonpolar forces which are typically assumed to be independent of the electrostatic potential [60,78]. The polar term, S, is a linear functional of the solution to the RPBE, u r , (also known as the reaction potential) [60] and can be expressed as For the 3-term splitting, the reaction potential is the sum of u 3 and u h .

ADAPTIVE FINITE ELEMENT METHODS
The finite element approach provides a natural framework for dealing with the complex molecular surfaces which arise in the PBE. Although there are modified finite difference methods which address this difficulty [83], finite element methods provide an attractive alternative when paired with an adaptive unstructured mesh designed to conform to the shape of the solute molecule [62]. In this section, we present a general adaptive finite element method for the regularized PBE, including the weak formulation, discretization, solution using an inexact global Newton iteration, and adaptive refinement procedure. For more details on the finite element method, see [73,32,26,45].
3.1. Weak Forms. To give a well-defined weak formulation, the nonlinearity involving exponentials must be controlled; in [39,56], a priori L ∞ estimates are obtained for any solution to the RPBE, giving almost everywhere pointwise bounds of the form: α u r β. This leads to working with a well-defined solution space that consists of a non-empty, topologically closed, convex subset of H 1 (Ω): It is shown in [39,56] that there exists a unique solution to either regularized form of the RPBE in M e ⊂ H 1 (Ω). The weak formulation is: Find u r ∈ M e such that (3.5) The linear functional L(·) is defined by integrating the right hand side of Eq. 2.5 by parts and applying the jump condition to eliminate the interface terms. The weak form for the 3-term split regularized PBE requires solving two problems: first for the harmonic term on Ω m and second for the split potential on the whole domain Ω. Define the solution space to the harmonic problem as where a m (·, ·) is the restriction of the bilinear form to the Ω m subdomain and 3.2. Solving. Due to the hyperbolic sine, the RPBE has a strong nonlinearity. The discretized nonlinear problems defined in Eqns. 3.2 and 3.6 can be solved using an inexact-Newton method [44]. For brevity, we give details for Eq. 3.2. Define the weak residual functional to be Here u h r is the discrete solution satisfying the system of nonlinear equations where V h is the space of piecewise linear functions defined by the tetrahedral mesh. Linearizing Eq. 3.9 around u h r results in In the linear RPBE, sinh(u) is replaced by u, and b (u, v) = (κ 2 u, v). Newton's method defines the nonlinear update vector s h as the solution to Given an initial guess u h r ≈ u 0 , the updated solution is defined as u 1 = u 0 + s h . This process can be repeated until a desired level of convergence is achieved. An inexact-Newton method uses an iterative solve to find an approximate solution to Eq. 3.11, with a relatively large tolerance for the linear solve when far from the nonlinear solution. However, as the exact solution to Eq. 3.9 is approached, the linear solver tolerance is tightened so that quadratic convergence is achieved.
The computational complexity of the Newton solver is dominated by the method used to solve the N linear algebraic equations [15,50] within each iteration. Multilevel methods provide an advantage in that they are provably optimal or nearly optimal methods for solving these systems [13,50,80]. The presence of geometrically complex discontinuities in the dielectric and in the Debye-Hückel parameterκ in the PBE destroy classical multilevel method efficiency, and can even cause divergence. This is analyzed at length for the PBE in [53,58,59], and various techniques based on coefficient averaging and algebraic enforcement of variational (Galerkin) conditions are examined. Algebraic multilevel methods have been used successfully for many similar problems; cf. [16,17,30,36,35,31,71,76,77]. A fully unstructured algebraic multilevel approach is taken in FETK, more details are provided in Section 5.
Starting from an initial mesh T 0 , the adaptive mesh refinement procedure builds a sequence of conforming meshes T 0 , T 1 , . . . , T l [39,56,43]. This procedure is divided into four steps: SOLVE, ESTIMATE, MARK, and REFINE. In the SOLVE step, a solution is computed on the current mesh. Using this result, the ESTIMATE step computes elementwise error indicators and an estimate of the global error. In a production environment, the procedure terminates if the global error estimate is below some prescribed tolerance. Otherwise, the MARK step selects elements for refinement. This step is crucial to the convergence of the method. Finally, REFINE subdivides the marked elements possibly subdividing additional unmarked elements in order to produce a conforming mesh. The refinement technique used in this article is longest edge bisection [69].

ERROR INDICATORS
In this section we present a posteriori error indicators for use in adaptive refinement. These estimators are typically developed by considering the residual of the weak form. For example, given a finite element solution u h r ∈ V h , the weak residual for the linear RPBE is (compare to Eq. 3.8) for a given v ∈ V . For the remainder of this article, we restrict our attention to two classes of error indicators: energy-based and goal-oriented. The first class estimates the error in the energy norm, although this idea can be generalized to other norms, (e.g., the H 1 norm). The second class, called goal-oriented, focuses on estimating the error in a user specified quantity of interest or goal functional. In the following sections, we derive error indicators from both classes for the linear RPBE. For the goal-oriented indicator, the solvation free energy is used as the target functional.

Energy Norm Indicators.
A standard a posteriori error indicator is based on bounding the error in the energy norm. It is easily derived by breaking the weak residual into its elementwise components and integrating by parts over each element [1]. This technique was used in [40] to derive the following estimator for the linear RPBE and This indicator gives a bound on the error measured in the energy norm ( Bounding the error in other norms is possible. For example, in [39,56] a similar a posteriori error indicator for the RPBE was shown to bound the error measured in the H 1 norm. Other efforts have focused on formulating the RPBE as a first-order system least squares (FOSLS) problem, which has a natural error estimate [22,38].

Goal-Oriented Indicators.
Key to the development of goal-oriented error indicators is relating the weak residual to the error in the goal functional. For symmetric linear problems, a direct application of the Riesz representation theorem shows that there exists a dual function that when paired with the weak residual gives the error in the goal [21]. The challenge is to approximate this function and utilize that approximation to develop error indicators. However, for nonlinear problems, like Eq. 2.5, the definition of the dual function is not so clear. The first part of this section discusses a strategy for defining the dual function for both the nonlinear RPBE and the three term splitting. Using this definition of the dual, two goal-oriented error indicators are proposed for the linear RPBE utilizing the solvation free energy as the quantity of interest. Following [21], a dual function for the nonlinear RPBE can be defined by considering the constrained minimization problem where a(·, ·), b(·, ·) and L(·) are specified in Eq. 3.5, and S is the goal functional. Note that because the RPBE constraint determines the solution uniquely, the minimization problem has the same solution. However, specifying the minimization provides an additional mathematical framework to define the dual function. To see this, consider the Lagrangian associated with the minimization problem where the Lagrange multiplier, w ∈ H 1 0 (Ω), is also the dual function. Taking the first variation of Θ with respect to u gives the dual problem: where DR(·)·, · was defined in Eq. 3.10. As discussed above, if b(·, ·) is linear in the first argument then Eq. 4.8 simplifies to a(v, w) + (κ 2 v, w) = S(v) ∀v ∈ H 1 0 (Ω). For the linear problem the error in a goal functional S(·), like the solvation free energy, is simply expressed in terms of the weak residual Thus, if w is known, the error in S(·) is easily calculated. In the nonlinear case the error in the goal satisfies where E is quadratic in the error in u h r [21]. For the three term splitting of the PBE, again we setup a constrained minimization problem to define the dual. Using the notation from the weak form in Eq. 3.6, the corresponding Lagrangian is where w 3 ∈ H 1 0 (Ω) and w h ∈ H 1 0 (Ω m ) are dual functions. The functional S m (·) is a restriction of the original goal functional to the Ω m domain. Taking the first variation of Θ 3-term with respect to u 3 and u h gives the dual problem To make the application of the dual functions in error indicators more concrete, we present two goal-oriented indicators for the linearized RPBE. As an example we will focus on accurate computation of the solvation free energy (see Eq. 2.12). Unfortunately, S(·) is not bounded on H 1 0 (Ω) due to the inclusion of delta distributions. A common approach to circumventing this issue is to use a mollified version of the functional [1,11]. In this case, the mollified solvation free energy is where θ is a locally supported function defined such that One possible choice for θ is the step function where B σ is the volume of a ball of radius σ.
A simple error indicator suggested by Eq. 4.9 is to first solve the dual problem using the same approximation space as the primal. This approximate dual could then be substituted for w in Eq. 4.9 to compute the value of the indicator. However, if the same finite dimensional space is used for solving both the dual problem and the primal problem, then, because of Galerkin orthogonality, Eq. 4.9 will be zero (see [11]). A remedy is to instead solve the dual problem using a finer approximation space, U h ⊂ H 1 0 (Ω). One convenient choice is to maintain the same mesh and use higher order polynomials for U h . In the examples below, V h is the space of piecewise linear polynomials and U h is the space of piecewise quadratics. Let the finer resolution solution of the dual problem be denoted w h,2 . Substituting w ≈ w h,2 into Eq. 4.9 yields h,2 is a convenient projection (e.g., nodal injection) of the fine space U h onto V h . The choice of the projection operator will affect the quality of the indicator. Decomposing the error into its elementwise contributions gives where the subscript K indicates the restrictions of the linear functional, bilinear functional or inner product to element K and Here, the absolute value of the integrand in Eq. 4.18 has been taken over each element. In the numerical experiments in section 6, Eq. 4.19 is referred to as the the "goal-quadratic" error estimator used by the various adaptive refinement marking strategies. The error indicator in Eq. 4.19 requires solving a dual problem which is substantially larger than the primal problem. To alleviate this issue, we develop a second goal-oriented error estimator that finds an approximation to the dual in V h (which is the same space as the primal problem). The error in the goal is estimated by solving many local elementwise boundary value problems. The technique proposed here is similar to the development of the equilibrated residual method for computing goal-oriented estimators [66,67]. However, the less accurate but simpler element residual method (ERM), as discussed in [1], is used. Using the parallelogram law, the error in the linear functional can be rewritten [1,66,67] as where again |||v||| 2 = a(v, v) + (κ 2 v, v) is the square of the energy norm. Define elementwise error functions φ K = u r −u h r | K and ψ K = w −w h | K to be computed by the element residual method. The error in the primal problem on element K is approximated by the solution to and (4.23) Similarly, the error in the dual problem on element K is approximated by the solution to (4.25) For a derivation of these equations see [40]. Explicitly stated, the "goal-linear" error estimator used in the numerical experiments in section 6 is given by (4.26)

MULTILEVEL PRECONDITIONING
As the mesh is refined, the conditioning of the linear system deteriorates, and preconditioning is necessary to accelerate the convergence of iterative solvers (e.g., the conjugate gradient method). The challenge in designing an efficient preconditioner is balancing the cost of applying the preconditioner with its effectiveness in improving the conditioning of the underlying system.
In a given finite element mesh, at level j, we denote the set of nodes and its cardinality by N j and N j , respectively. We call the set of nodes introduced precisely at level j the fine nodes, and denote them N f j . As the mesh is refined, N f j is appended to N j−1 , leading to the following hierarchy of nodes: In the local mesh refinement setting, the way the coarse and fine nodes are processed plays a central role in determining the overall efficiency of a preconditioner. If the computational cost per level can be maintained proportional to N j − N j−1 , or slightly larger, then total cost will be order N , and the resulting preconditioner is said to have optimal computational complexity per iteration. If the resulting preconditioned system has a bounded condition number (independent of problem size), a solution can be obtained using an iterative method with a bounded number of iterations. Hence, the combination of optimal per iteration complexity with a bounded condition number leads to a solver with optimal overall complexity.
In this article, we restrict the presentation of local multilevel preconditioning to a purely geometric (node based) perspective because the computational complexity is exactly governed by the number of nodes processed by the preconditioner at each level. The local multilevel preconditioners of interest can be classified into two groups: multiplicative local multigrid (MG) [18,20,28,68] and additive local MG [28,29]. The additive local MG preconditioner is often called the Bramble-Pasciak-Xu (BPX) preconditioner in the literature. In this article, we report on only the multiplicative variants. We use the term "classical" to refer to the application of a preconditioner in the uniform refinement setting. There is abundant literature on MG preconditioners [33,75,80,81] and local MG (e.g., see [19,79] for a review).
Proofs demonstrating the optimality of classical MG and classical additive MG preconditioners rely on a geometric increase in the number of nodes per level. This is because the cost per iteration of these classical preconditioners is proportional to N j (not N j −N j−1 ) per level, which results in suboptimal complexity if J j=1 c j N j is not O(N ). This occurs frequently in the local refinement setting due to the slow increase in the number of nodes between levels.
The family of hierarchical basis (HB) preconditioners, developed by Bank, Dupont, and Yserentant [12,14,82], maintain a per-level cost proportional to N j − N j−1 , by only processing (smoothing) the fine nodes at each level. Although the cost per iteration is optimal, HB preconditioners do not achieve a uniformly bounded condition number, and suffer from O(J 2 ) and O(2 J ) iterations in two-and three-dimensions, respectively. To address this deficiency, we investigate local MG preconditioners which process a larger set of nodes, X j , but still maintain a cost which is proportional to N j − N j−1 at each level. Hence, we seek a set X j such that N f j ⊂ X j ⊂ N j , with cardinality, X j , which is proportional to N j − N j−1 . At the same time, X j , should be large enough that the resulting system has a bounded condition number, leading to a solver with optimal overall complexity. Aksoylu and Holst [6] showed that this is possible even for three-dimensional local refinement routines.
In the local mesh refinement setting, Aksoylu, Bond, and Holst [2] studied the implementation and algebraic aspects (e.g., matrix representations) of multilevel preconditioners. Subsequent articles provide a comprehensive overview of local MG preconditioners with various emphases: for a theoretical treatment, see [3,4,5]; for optimality analysis in three-dimensional local refinement routines, see [6]; for surface mesh applications in computer graphics, see [7]. 5.1. Local multigrid preconditioners. As mentioned in the previous section, the fundamental difference between classical and local MG preconditioners is the smoothing operation. In classical MG, the smoother acts on all degrees of freedom on every level. In contrast, local MG only smooths a small subset, typically a neighborhood, of the fine degrees of freedom. Pseudo-code for a local MG V-cycle is provided in Algorithm 5.1.
The set of nodes used by the underlying method's smoothing operation defines the type of preconditioner. In this article, we introduce three different sets of nodes to be used in the local MG preconditioner. On a given refinement level, the marked region is the set of elements which have been marked for refinement. The refinement region is union of the newly introduced elements as a result of the refinement and closure procedure. Figure 2 depicts a two-dimensional refinement region formed by a local refinement routine which consists of quadrasection closed by bisection (the so-called red-green refinement). The sets of nodes used to define the preconditioners are as follows: HB: •. BPX: •, . BEK: •, , +. ONERING: •, , +, . FIGURE 2. An example of red-green refinement (quadrasectionbisection) in two dimensions. On the left, the mesh before refinement, with the marked region shaded in gray. On the right, the mesh after refinement (red edges) and closure (green edges), with the refinement region shaded in gray. The labels indicate which nodes are assigned to each of the preconditioner smoothing sets, X j .
• X j -HB: The set of fine nodes on level j, i.e., N f j [12,14,82]. • X j -BPX: The set of nodes whose corresponding basis functions have support entirely contained in the refinement region [28]. • X j -BEK: The set of nodes whose corresponding basis functions have non-empty intersection with the marked region. This set is named after the Bornemann-Erdmann-Kornhuber type refinement [23] routine. It consists of fine nodes and their immediate neighboring coarse nodes in the marked region. For P1 elements, this set can be inferred from the nonzero pattern of the prolongation operator. • X j -ONERING: The set of nodes whose corresponding basis functions have nonempty intersection with the refinement region [2,24,27,41]. This set consists of fine nodes and their immediate neighboring coarse nodes in the refinement region. For P1 elements, this set can be inferred from the nonzero pattern of the coarse-fine subblock of the stiffness matrix. As an example, we have labeled the nodes in Figure 2 to show which nodes are in each of the sets described above. The set corresponding to the classical multigrid preconditioner contains all of the nodes on each level.
We should note that the practical implementation of the various local MG preconditioners varies significantly depending on the particular preconditioner. Special care must be taken in order to achieve optimal computational as well as storage complexities. The implementation aspects of how to construct optimal complexity preconditioners are studied in more detail in [2]. 6. NUMERICAL EXPERIMENTS 6.1. Adaptive Refinement. In this section, the effectiveness of goal-oriented mesh refinement is compared to refinement using the energy-based error indicator. Of interest is  computing the solvation free energy of the 921-atom Fasciculin-1 protein [61], using the solution of the linear RPBE. All tests were performed using FETK [54].
In order to solve the RPBE, a definition of the molecular surface and a mesh conforming to that surface is needed. Various definitions of the molecular surface have been proposed in the literature, and the particular value obtained for the solvation free energy will depend strongly on the surface geometry [10]. However, the performance of the algorithms proposed here are insensitive to the choice of molecular surface, as long as it is sufficiently smooth, and the underlying mesh is conforming. Historically, generation of the mesh conforming to the surface was a great impediment to using finite elements for solving the PBE, and only recently, with the development of tools like PDB2PQR and GAMer, has molecular meshing become a routine task. The first step is to prepare the structure using PDB2PQR [42], which adds missing hydrogens, assigns charges, and specifies a radius for each atom in the protein. Next, the resulting PQR file is passed to GAMer [56,52,84,85], which produces a tetrahedral mesh conforming to the shape of the protein. Finally, this mesh is used by FETK to solve the linear RPBE, and compute the solvation free energy. 6.1.1. Global Marking Strategy. The adaptive refinement algorithm creates a sequence of grids T 0 , T 1 , . . . T l . . . based on an error indicator. Critical to this algorithm is the third step MARK. The goal of this step is to select elements for refinement. There are several different choices for marking strategies [1,11]. The strategy used here is to mark all elements in the l th refinement level that satisfy where γ ∈ (0, 1). This criteria yields no refinement for γ = 1 and uniform refinement for γ = 0. The convergence of the solvation free energy for Fasciculin-1 using the global marking strategy with both energy-based (Eq. (4.2)) and goal-oriented indicators (Eq. (4.19)) can be seen in Fig. 3. The figure shows the relative error in the solvation free energy as a function of the number of unknowns in the primal problem. The error in the solvation free energy is estimated by computing a high resolution solution to the PBE using uniform mesh refinement. The vertical line near 1.5 × 10 6 unknowns marks the size of this reference solution. Notice that this figure differs from traditional finite element convergence plots which show error as a function of element radius. Since the resolution of Energy-Based Goal-Oriented an adaptively refined mesh can greatly vary over the spatial domain, the element radius is not an appropriate measure of the resolution of the mesh. Furthermore, the order of accuracy of the approximation is based on the type of basis functions used, and thus is the same for both uniform and adaptive refinement. The benefit of adaptive refinement is that the mesh can be refined in regions that heavily contribute to the error, resulting in higher accuracy with fewer total degrees of freedom. In Fig. 3, the lower line shows the convergence of the solvation free energy for several levels of adaptive refinement using the energy-based indicator with the global marking strategy. This scheme makes steady progress to the correct solvation free energy. On the other hand, the upper line shows the results using goal-oriented indicators with the global marking strategy. This scheme performs very poorly.
The reason for the poor performance can be explained by looking at Fig. 4. The images are cut-aways of the 3D mesh, colored to indicate the distribution of marked elements. The left image shows elements selected by the global marking strategy using the energybased indicator. While the image on the right uses the goal-oriented indicator. The white elements are unmarked elements in the solvent subdomain and the gray elements are unmarked elements in the solute subdomain. Elements colored red are marked solvent elements, while blue elements are marked solute elements. Notice for the energy-based indicator only elements in the solvent subdomain are marked. This indicates that the reaction potential in the solute domain is relatively well approximated compared to the solution in the solvent subdomain. The distribution of the elements marked using a goaloriented indicator is focused on a few locations in the inner subdomain around the solute atoms. This explains the poor convergence for the goal-oriented indicator in Fig. 3. Since the strategy does not indicate there is error in the solvent domain, no refinement takes place.
To gain further insight into why the global goal-oriented strategy only marks elements in the solute domain, we refer the reader to Fig. 5. This image shows a cut-away with elements colored by their approximate signed contribution to the error in the goal. This where γ ∈ (0, 1). The split-marking strategy marks a significantly different group of elements, especially for the goal-oriented indicators. The color coding of Fig. 6 is the same as in Fig. 4. Again the energy-based refinement selects elements primarily in the solvent subdomain. However, because split-marking forces refinement in both subdomains, a few elements along the interface in the solute domain are also selected. In contrast, split-marking using a goal-oriented indicator marks a few elements in the solvent subdomain, while also marking elements surrounding the solute atoms.

Energy-Based
Goal-Oriented   . Convergence of the solvation free energy for Fasciculin-1 with both goal-oriented and energy-based indicators using the split marking strategy. The convergence measured against the size of the primal/dual problem is on the left/right. Figure 7 shows the relative error of the solvation free energy in Fasciculin-1 as a function of problem size for several indicators using the split-marking strategy. In the figure, there are two plots. The plot on the left shows the relative error in the solvation free energy as a function of the number of unknowns in the primal problem. The second plot, with the exception of the energy-based refinement strategy, shows the relative error as a function of the size of the dual problem. For energy-based refinement, since no dual is needed, the horizontal axis is the size of the primal problem.
In the figure, the split marking strategy using the energy-based indicator converges steadily. This is similar to the convergence of energy-based refinement using global marking. On the other hand, there is a dramatic improvement in the convergence of both the linear and quadratic goal-oriented refinement techniques, with the rate even increasing slightly as the size of the problem increases. Compare this with the poor results for goal-oriented refinement with the global marking strategy from Fig. 3. The Number of levels (r) improvement comes from the split marking strategy explicitly taking into account the error in the solvent and trying to control where it is large. The second plot in Fig. 7 shows that solving the dual problem using piecewise quadratic elements has substantial additional cost. Although, it is likely that this cost would be mitigated by the additional work needed when solving the nonlinear RPBE (see Eq. 2.5). In contrast, the goal-oriented strategy using an indicator constructed from a linear dual problem (and split marking) does not suffer from the same problem, and is the most efficient method when the total cost is taken into account. 6.2. Performance of Preconditioners. As was discussed in section 5, classical multigrid (MG) preconditioners perform best in the uniform refinement setting, where there is a rapid geometric growth in the number of unknowns as the mesh is refined. In the local refinement setting, this growth is much slower, and is frequently subgeometric when the refinement is concentrated in the neighborhood of a low dimensional feature (e.g., a point or a line). As a result, the per-iteration complexity of classical MG may fail to scale linearly (or suffer from a large scaling constant) as the number of unknowns increases. In contrast, many local MG preconditioners do not have the same restriction, and maintain optimal per-iteration complexity in both the local and uniform refinement settings.
In Figure 8, the number of primal unknowns is shown as a function of the refinement level for three different refinement strategies. For the energy-based marking/refinement strategy, the growth in the number of unknowns is geometric, but subuniform. For goaloriented refinement the growth is even slower, but still geometric. The dashed reference line shows a geometric growth rate proportional to 2 r .
As a measure of the relative locality of each preconditioner, we compute the ratio of the number of unknowns processed by the smoother, shown in Figure 9. Classical MG is a global method, so the ratio of smoothed to total unknowns is 1 in both the energy based and goal oriented refinement cases. As expected, the BEK preconditioner consistently FIGURE 9. Ratio of the number smoothed unknowns, X j , to total unknowns, N j , on each level for different multilevel preconditioners FIGURE 10. Conjugate gradient method iteration counts for the preconditioners used. smooths more unknowns than HB, but fewer than MG. This is because the set processed by the BEK is a superset of the set processed by HB.
By design, the local multilevel preconditioners considered here maintain an optimal per-iteration complexity in both the local and global refinement settings. The primary challenge for these preconditioners is achieving a bounded condition number, independent of problem size. It can be shown that for a given tolerance, the number of conjugate gradient (CG) iterations can be bounded by a function of the condition number. Hence, if the condition number is bounded, so is the number of CG iterations. In Figure 10, we report the number of CG iterations as a function of refinement level for MG, BEK, and HB. As predicted by theory, application of the HB preconditioner leads to a slow growth in the number of CG iterations as the mesh is refined regardless of indicator type. In contrast, for both goal-oriented and energy-based refinement the iteration count is bounded for both the classical MG and BEK preconditioners. The iteration counts for BEK are modestly higher than classical MG, but the work per level is reduced, since BEK smooths only a fraction of the unknowns smoothed by classical MG. For this reason, BEK is a compelling alternative to classical MG and HB.

CONCLUSION
In this article, we developed goal-oriented error indicators for accurate computation of the solvation free energy from solutions of the regularized Poisson-Boltzmann equation. We found that due to oscillations and imbalanced cancellation in the error contributions, global marking strategies based on goal-oriented error indicators were not viable for driving adaptive mesh refinement. To address this problem, we developed a split marking strategy based on considering each subdomain individually. In numerical experiments, we calculated the solvation free energy for a 921-atom Fasciculin-1 protein.
Through these experiments, we showed that the new marking strategy, combined with goal-oriented refinement, is more efficient than energy-based refinement in the context of solvation free energy calculations.
The use of adaptive mesh refinement puts a greater burden on the preconditioner to maintain optimal runtime efficiency. To address this issue, we investigated the use of local multigrid methods, which have a lower per-iteration complexity compared to classical global multigrid. In particular, the BEK variant proved to be a compelling alternative to classical multigrid since it has optimal per-iteration complexity, while still maintaining a bounded iteration count as the mesh is refined. The result is an iterative solver with an optimal overall complexity, scaling linearly with problem size.

ACKNOWLEDGMENTS
BA was supported in part by NSF Award 1016190. SB was supported in part by NSF Award 0830578. EC was supported in part by a University of Illinois CSE Fellowship and by a DOE Office of Science ASCR-UQ effort at Sandia National Laboratory under contract DE-AC04-94AL85000. MH was supported in part by NSF Awards 0715146 and 0915220, and by DOD/DTRA Award HDTRA-09-1-0036. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.