ARGOS: An adaptive refinement goal‐oriented solver for the linearized Poisson–Boltzmann equation

An adaptive finite element solver for the numerical calculation of the electrostatic coupling between molecules in a solvent environment is developed and tested. At the heart of the solver is a goal‐oriented a posteriori error estimate for the electrostatic coupling, derived and implemented in the present work, that gives rise to an orders of magnitude improved precision and a shorter computational time as compared to standard finite difference solvers. The accuracy of the new solver ARGOS is evaluated by numerical experiments on a series of problems with analytically known solutions. In addition, the solver is used to calculate electrostatic couplings between two chromophores, linked to polyproline helices of different lengths and between the spike protein of SARS‐CoV‐2 and the ACE2 receptor. All the calculations are repeated by using the well‐known finite difference solvers MEAD and APBS, revealing the advantages of the present finite element solver.


| INTRODUCTION
Electrostatics is the key that allows proteins and other biological macromolecules to perform their multiple functions [1][2][3][4]. At the heart of many electrostatic theories lies the Poisson-Boltzmann equation (PBE) providing the electrostatic potential (ESP) of a macromolecule in an ionic solution [2,[5][6][7]. The PBE is widely applied in different areas of research, such as biomolecular recognition, kinetics of molecular association and enzyme catalysis (see [7] and the references therein).
It treats the environment of the macromolecule as a polarizable continuum containing a Boltzmann-distributed continuous charge density of the ions. The macromolecule itself is described by atomic partial charges that are embedded in a molecule-shaped cavity with a continuous polarization, characterized by a dielectric coefficient.
Analytical solutions of the PBE are available only for specific simple systems. For complex biomolecular surfaces and arbitrary distribution of charges, numerical methods have to be employed to solve the PBE. From a mathematical point of view the PBE represents a second-order semilinear elliptic partial differential equation with a measure-valued right hand side. The PBE imposes several challenges impacting the accuracy of numerical algorithms, including discontinuous dielectric coefficients representing interfaces (e.g., solventbiomolecule), singularities as a result of the representation of the partial charges by delta functions, and the complicated geometry of molecular domains in three spatial dimensions. To effectively handle the above challenges, a PBE solver must employ some kind of adaptivity. While the error indicators controlling the mesh refinement steps in the existing adaptive PBE solvers attempt to detect the regions where the error in the electrostatic potential (ESP) is the highest, this is far from the best strategy when it comes to the calculation of the electrostatic interaction between molecules. Since the interaction energy E is a number that depends on the ESP, the purpose of the error indicator is to detect the regions for which the error in the ESP contributes the most to the error in the goal quantity E. Those are the regions that should be refined and are not necessarily the regions, where the error in the ESP itself is the highest. The use of standard error indicators results in much denser adaptive meshes with a huge number of nodes and inherently unreliable evaluation of the goal quantity E.
Most of the common approaches for the numerical solution of the PBE fall in one of the following three categories: finite difference (FD), boundary element (BE) and finite element (FE) methods. The majority of the PBE solvers like DelPhi [8], GRASP [9], MEAD [10], UHBD [11], the PBEQ module of CHARMM [12], and the PMG solver in the APBS software suite [13] are based on FD methods, utilizing regular 3D lattices, due to the ease of numerical implementation.
Finite difference schemes rely on solving the PBE in a lattice representation. The charges and dielectric permittivity are mapped onto a rectangular grid and a discrete approximation employing difference quotients is applied to the partial differential equation. As a result, a system of linear/nonlinear algebraic equations is build up, which is often solved by iterative methods, for example, the Conjugate Gradient, multigrid, or the Newton-Raphson method.
To speed up the convergence and reduce the memory requirements, a focusing technique was developed in Refs. [8,14]. This procedure involves solving the equation on a coarse grid, covering a large region, followed by a solution on a finer grid on a smaller region, with a boundary condition, interpolated from the values of the coarse grid solution. This approach was further improved in Ref. [15] by implementing it in parallel and then was applied to the electrostatics of microtubules and ribosoms. Due to the highly improved performance, focusing schemes are utilized in most of the FD-based solvers for the PBE which use regular cubic lattices.
As the underlying grid in regular 3D lattices is not adapted to the molecular surface, the accuracy of the solution is decreased. The correct assignment of the dielectric coefficient in points close to the solute-solvent interface is crucial for an optimal convergence of FD methods. One way to achieve it is the use of certain averaging procedures. More precisely, for each integer or half-integer grid point the values of the dielectric coefficient ϵ, mapped on the grid, are defined by averaging its actual values at the neighboring points, which depend on the modeling assumptions concerning the definition of the solutesolvent interface and the dielectric permittivity in the molecular and solvent regions [14,16]. The inverse Debye length k is defined in a similar fashion. Due to the averaging, there is a region around the true interface in which these model parameters vary between their respective values on both sides of the interface. The thickness of this layer depends on the grid spacing h and this is where the discretization error of the FD method is typically much larger than that produced by the FE method, discussed in the end of this section.
Finite difference methods also exhibit difficulties caused by nonsmooth source terms, such as the delta distributions introduced as a result of the mathematical modeling of point charges. A redistribution of the charges around their closest grid points is necessary in order to minimize the errors occurring in the electrostatic potential at small distances from the charges (see Refs. [14,16] for different approaches to such a redistribution). Another problem, even though not unique to the FD method and also present in finite element implementations, that arises in the computation of free energy differences (e.g., the calculation of solvation energies) is the so called grid artifact [17][18][19].
The reason behind it is the singular charge distribution in the molecular system resulting from the modeling of the atomic charges as discrete point charges instead of using alternative charging models, where the atomic charge is represented by a more regular charge density function (see [20] for an example of such alternative charging model, where each atom is represented as a uniformly charged ball).
In the case of such a singular charge distribution, the Coulombic selfinteraction energy of the point charges becomes infinite and under certain assumptions (see assumptions 1-3 in section 5 from [18]) it cancels out in analytical expressions for the free energy difference.
However, in the discrete FD analog, the Coulombic self-interaction energy participates in the free energy as a finite, but large contribution that is inversely proportional to the grid spacing. Therefore, the discrete approximation of the (infinite) Coulombic self-interaction energy is not guaranteed to cancel out in the difference of the numerically computed electrostatic energies unless one uses exactly the same grids and redistribution of the point charges. Another way to avoid the spurious grid artifact in the calculation of the energy difference is to apply proper splitting techniques for the full potential in the PBE. By applying such decompositions of the full potential, the source term or, equivalently, the right hand side in the PBE, that is the delta distributions describing the discrete charge density due to the point charges in the solute, are transformed to a more regular distribution which becomes the source term in the equation defining the regular part of the potential in all these decompositions. In the so-called weak formulation 1 of the problem defining the regular part of the potential, this better behaved right hand side can be represented as a surface integral over the solute-solvent interface of a (polarization) charge density function, see Equation (21) in [21] and the interpretation after it. Thus, mathematically speaking, this right hand side can be interpreted as a certain charge density concentrated on the molecular surface (see Remarks 3.1 and 3.2 in [22] or Remarks 2.1 and 2.2 in [23]). We note that these splitting techniques can be applied in conjunction with any numerical procedure to solve the PBE, and in particular with the FD method.
However, the surface charge density on the interface between the molecular region and the solution region remains hard to handle with the FD method since no explicit molecular surface is constructed (see, e.g., Refs. [21,24]).
Another disadvantage of FD methods is that the errors, caused by the specific way of handling interfaces and the partial charge distribution in the solute, strongly depend on the particular position and orientation of the molecular system relative to the FD grid. One approach to cope with this problem is to use finer grids, which naturally leads to a computationally more expensive problem. Another approach is to use a rotational averaging scheme. The idea is to make several calculations with slightly different mutual orientation between the molecular system and the FD grid, and finally take an average (see Refs. [14,21]).
A third approach, which could cope better with curved interfaces and regions where the solution is sharply varying (mainly around the partial charges and across interfaces), is to involve some kind of adaptivity in the FD method in conjunction with splitting techniques for the full potential. Finite difference methods based on adaptive Cartesian grids (ACG) have been developed for solving nonlinear equations in domains with curved boundaries [25][26][27]. Here, the grid representation is based on quadtrees in two dimensions and octrees in three dimensions. The ACG is recursively adapted by identifying mesh cells, which intersect the molecular surface, and test whether they satisfy a certain set of conditions. Using a coarser mesh away from the surface reduces the total number of grid points, which for ACG often is several orders of magnitudes less than that required in a conventional lattice grid. However, the underlying mesh still does not conform to the solute surface and thus an accurate evaluation of the ESP near the interface remains problematic. In an attempt to avoid this shortcoming, the technique proposed in Ref. [28] introduces additional mesh points at the Cartesian grid-surface intersection and uses a reconstruction of the ESP from the solution obtained in the exterior and interior domains. Note that, in this approach the continuity of the normal component of the displacement field is not enforced in a strong (pointwise) sense, which is a potential source of error. Moreover, the accuracy typically deteriorates when a charge approaches the surface and the mesh does not adequately resolve the rapid variations of the potential induced by such a charge. Recently, an FD method based on ACG was implemented in the finite difference solver CPB, particularly aiming at the approximation of the PBE for complex biomolecular structures [29]. This solver also utilizes a splitting technique for the full potential to eliminate the singular behavior at the charge sites.
The boundary element method (BEM) is another alternative discretization technique for partial differential equations. It was first adapted to the computation of the ESP of macromolecules in Ref. [30]. Here, the PBE is reformulated into boundary integral equations for the boundary of the solute domain by employing Green's theorem. Therefore, only the two-dimensional molecular surface needs to be discretized, and thus the number of degrees of freedom (DOFs), that is, the number of unknowns is greatly reduced. Moreover, the interior point charges can be placed at their exact locations and no redistribution is needed as it is the case with the FD method. Consequently, to obtain reliable results, BEM does not need focusing or rotational averaging, unlike the FD method. The grid artifact, as a result of the infinite Coulombic self-interaction energy in the case where the atomic charge density is modeled by discrete point charges, also does not exist in BEM [31].
The resulting boundary integral equations are solved by using collocation, Galerkin, or least squares methods. The first of these approaches is the most commonly used by the PBE community since it is the simplest to apply in practice. Further developments of the BEM with the same focus of application can be found in Refs. [32,33] and deal mainly with algorithmic improvements. The BEM has some advantages over other numerical methods like the FE and FD methods, which involve volume-domain discretization. Since only the molecular surface is discretized, the number of unknowns is greatly reduced. Moreover, boundary conditions at infinity are treated exactly as opposed to the FE and FD methods. In contrast to the FD method, the jump in the normal component of the electrostatic field across the molecular surface (interface) is explicitly accounted for, which results in a more accurate solution near the solute surface.
The BEM also possesses certain disadvantages, mainly related to the requirement of the efficient evaluation of singular and hypersingular integrals and numerous other integral operations, a dense nonsymmetric coefficient matrix for the linear system of equations [31,34]. Moreover, in the presence of an ionic charge density of the solvent, BE methods usually neglect the ion exclusion (Stern) layer and thus only take into account a single surface separating the solute and the solvent.
One of the promising strategies to improve the efficiency of the BEM is to combine it with a fast multipole method (FMM) [35,36]. Further improvements can be achieved by using an adaptive FMM for accelerating the convolution type matrix-vector multiplications [37]. Solvers based on these concepts benefit from a more efficient use of the memory and improved load-balance between the local-and far-field calculations, resulting in a significant speed up of the computations.
The third approach to the numerical solution of the PBE is the finite-element (FE) method, which approximates the potential by a superposition of a set of basis functions. The FE method, which is employed in the present work, is the most successful for elliptic PDEs, since it combines geometrical flexibility and satisfactory convergence analysis with the ability to handle nonlinear problems involving interface jump conditions and nonsmooth source terms. Moreover, it enjoys a wide range of efficient iterative solvers for the resulting sparse linear systems, where the PBE is typically solved on unstructured finite element meshes. The finite element discretization results in a linear or nonlinear algebraic system of equations for the expansion coefficients in the linear combination of basis functions that represents the approximate potential. This system is then solved usually by means of iterative methods such as the Conjugate Gradient or the Newton-Raphson method.
One of the first applications of the FE method to continuum electrostatics of solubilized biomolecules appeared in Ref. [38]. Since then, the popularity of this approach has increased synchronically with the increase of computational power and advancement of FE mesh generation techniques, see, for example, Refs. [13,[39][40][41][42]. As representatives of the solvers, implementing the FE method for the numerical approximation of the PBE and the linearized PBE (LPBE), we mention APBS [13] (utilizing adaptive mesh refinement based on a residual type error indicator) and mFES [17], respectively.
However, application of FEM solvers for the PBE and the LPBE to a reliable calculation of the electrostatic interaction between molecules becomes extremely difficult in the absence of appropriate error estimators and adaptive mesh refinement techniques. The two FEM solvers mentioned above face the same issues. More precisely, mFES does not implement error estimation and adaptive mesh refinement techniques, besides the fact that the initial mesh is somehow adapted to the shape and position of the solute. The tetrahedral volume mesh inside the solute has approximately constant edge size related to the mesh size on the solute-solvent interface, whereas the mesh size outside the solute is gradually increasing up to the boundary of the computational domain. On the other hand, APBS utilizes adaptive mesh refinement based on a residual type error indicator. However, in APBS, the error is controlled in the so called energy norm, which is a global (integral) quantity, whereas for the evaluation of the electrostatic coupling between molecules the ESP needs to be evaluated only at particular points of interest. Therefore, adaptive mesh refinement based on energy norm error indicators attempts to enhance the global quality of the solution. This approach typically yields optimal rates of convergence with respect to the global energy norm, while the improvement in the potential at the points of interest happens at a much slower rate. Thus, an adaptive FE solver which is based on energy norm error indicators is inherently unreliable and ineffective for the computation of the electrostatic interaction between molecules.
In the present paper, we report an Adaptive Refinement Goal-Oriented Solver [43] (ARGOS) for the calculation of the electrostatic interaction between molecules. In order to have a consistent description of the electrostatic interaction, which is symmetric in nature, one needs to consider the linearized PBE with a homogeneous Dirichlet boundary condition as the primal problem governing the ESP. Besides the smallness of the ESP, this is why ARGOS, in its present implementation, is based on the linearized PBE. The most notable feature of this solver is the goal-oriented a posteriori error control, which allows adaptive refinement of the mesh specifically aiming at the reduction of the error in the electrostatic coupling, herewith significantly improving the accuracy of the numerical solution and the computational efficiency. We compare the accuracy of computed interaction energies obtained via this new approach with those from FD solvers (MEAD and APBS) both for practically relevant problems and for model problems with known analytical solutions.
The structure of the paper is as follows. In Section 2, the theory of a goal-oriented error estimator for the electrostatic coupling between molecules and its algorithmic realization are discussed. In Section 3, we start with an application of our new ARGOS to three test cases with analytically known solutions and evaluate the accuracy of ARGOS in comparison to the FD solvers MEAD and APBS. We proceed with the application of the different solvers to the computation of the interaction energy between two chromophores, that are attached to a polypeptide. Besides reevaluating the coupling between the transition densities of the chromophores [44], we study the interaction between their ground state charge densities for this system as well as the SARS-CoV-2 and the ACE2 cell receptor. The superior performance of ARGOS is used to evaluate the accuracy and numerical efficiency of the FD solvers MEAD and APBS. φ of a discrete charge distribution that is placed in a dielectric continuum containing a continuous ionic charge distribution. In the special case of only two univalent ion species in the solvent with the same concentration but opposite charges, the PBE is given by where k B is the Boltzmann constant, T the temperature, e 0 the elementary charge, N the number of charges in the molecule, z i is the numerical value of the i-th partial charge (in units of the elementary charge), and x i its coordinate in the molecule. The dielectric coefficient ϵ is piecewise constant and defined by where Ω m is the molecular domain and Ω s is the solvent domain composed of Ω ions and the ionic exclusion layer (IEL) Ω IEL . The latter is a layer around the molecules which no ions can penetrate. It is defined as the difference between the union of the van der Waals balls of the atoms, where the atomic radii have been increased by a counterion radius R ion , and the molecular region defined by the solvent excluded surface (SES) formed by the contact points of the van der Waals surface and a (solvent) probe sphere, rolling over it (see Refs. [45][46][47][48]). In 1983 Connolly [49] gave an analytical description of the SES, which is, therefore, also known as the Connolly surface. We set the probe radius in the definition of the SES to 1 Å in the calculations on our first test system (Alexa 488 and Alexa 594 dyes in Section 2.4) and to 1.4 Å in our second test system (SARS-CoV-2 and ACE2 cell receptor in Section 2.4), while the counterion radius R ion in the definition of the IEL is set to 2 Å in all calculations. It is known that the electrostatic interaction energy depends on the definition of the solute surface (see, e.g., Ref. [50]), but the study of this problem is beyond the scope of this paper. For a system of N mol molecules, Ω m consists of the domains Ω m1 ,…, Ω mN mol , where Ω mi corresponds to the domain of the ith molecule (see Figure 1). The coefficient k 2 is also piecewise constant and given by where the ionic strength I s is defined as with the molar concentration of each ion type c 1 = c 2 and N A denoting Avogadro's number. Assuming the potential φ to be small (φ ( kBT e0 ), the PBE (1) can be linearized. The resulting equation is referred to as the linearized PBE (LPBE) and will be used in this work. For the special case I s = 0, Equation (4) becomes the Poisson equation: Equations (1), (4), and (5) include a natural interface Γ between the molecular and the solvent domains due to different dielectric coefficients ϵ = ϵ(x) (see (2)). The potential φ and the normal component of the displacement field ϵrφ must be continuous across the interface, that is, where Á ½ Γ denotes the jump of the enclosed quantity across Γ. In order to find a unique solution of (1) and (4), a Dirichlet boundary condition for φ must be imposed. The usual approach in physics prescribes a vanishing potential at infinite distance from the boundary of Ω m . In practical applications, one works in a bounded computational domain Ω and imposes the boundary condition where g is usually close to zero and calculated accurately enough by solving a simpler problem, possibly with a known analytical solution.

| Integral formulation of the problem
The right-hand side of (4) is a functional, consisting of a linear combination of delta functions, while the left hand side defines a proper function everywhere except at the interface Γ, where ϵ is discontinuous and rÁ (ϵrφ) is not well defined. A natural way to interpret (4) is to rewrite it in its integral (weak) formulation, where the left-hand side and the right-hand side are comparable objects with well-defined meaning. To obtain the weak formulation, we multiply both sides of (4) by a function v equal to zero on ∂Ω, also called a test function, and then integrate the equation over the domain Ω: On the left-hand side we apply the divergence theorem for integration by parts together with the condition (7) and obtain (for details see Ref. [51]) It can be rigorously shown that under certain conditions on the regularity of the interface Γ there exists a unique potential φ that satisfies the integral formulation (10) for any test function v in a certain function space [22,23]. By directly discretizing (10), one obtains an approximation φ h of the exact potential φ.
Another widespread approach, from which both the numerical treatment and analysis of the LPBE (also the PBE) benefits, involves splitting the potential φ into two parts G and u, that is, φ = G + u (see, e.g., Refs. [21,24,52,53] is the singular Coulomb potential of the point charges in a uniform dielectric medium with ϵ = ϵ m without any molecular cavities and u is the reaction field potential that takes into account the polarization of the solvent.
Thus, one completely avoids the above mentioned grid artifact by computing only the reaction field potential u h instead of the potential φ that has singularities at the positions of the charges. This approach is particularly beneficial when one needs the reaction field potential only or when calculating the solvation energy, that is, the free energy difference between the solute in vacuum and in a solvent environment where q i = z i e 0 , u vac and u solv are the reaction field potentials for the molecule in vacuum and the solvent, respectively.

| Goal-oriented a posteriori error estimates
The purpose of adaptive finite element methods (AFEM) is to achieve high accuracy at reduced computational costs. The general idea In many applications one is interested not in the solution itself, but rather in a specific quantity that depends linearly on the solution.
Thus, the target quantity is expressed in terms of a linear functional, called a goal functional. In this case, the error indicators are provided only for the goal quantity. For example, if one is interested in the average value of a solution w, the goal functional acting on w can be represented by To outline the method, let us consider a general case of a linear differential equation with a solution w (see, e.g., Ref. [54]): where A is a linear differential operator and f is a generalized force, exerted by the environment. Problem (14) is called primal. Let ℒ w ð Þ be a quantity of interest obtained from the solution w by applying a goal functional ℒ and w h an approximation of w that depends on a discretization mesh parameter h. The strategy involves controlling the error in terms of the local residuals ρ K (w h ) computable on each mesh cell K (e.g., tetrahedrons or patches of tetrahedrons): The idea is to consider the adjoint (or dual) problem where A * is the adjoint operator defined via hAu, vi = hA * v, ui for all functions u and v with the appropriate boundary conditions, where hÁ, Ái is the standard L 2 integral inner product, and j is the density function associated with ℒ. Taking into account the definition of the adjoint operator, (17), and (15), the error in the goal quantity is obtained as where z is unknown and has to be substituted with an approximate solution. In this case, we have where s h (z) is a sensitivity factor obtained from the approximate solution of the adjoint problem (17). Thus, the adjoint solution demon- where q i,2 = z i,2 e 0 and N 2 is the number of charges in M 2 . The energy can also be expressed by the potential φ 2 of M 2 at the position of the charges q i,1 of M 1 , that is, where q i,1 = z i,1 e 0 and N 1 is the number of charges in M 1 . Here, we assume φ 1 and φ 2 to be the solutions of the LPBE (4) together with the interface conditions (6) and (7), and a homogeneous Dirichlet boundary condition g = 0 in (8). Note that both, the linearity of the LPBE and the homogeneous Dirichlet boundary condition on ∂Ω, are necessary in order to reflect the symmetry of the interaction, that is, In what follows, without loss of generality, we define the primal problem to be It describes the potential created by the charges of molecule M 2 . For convenience, we denote the right-hand side of (22) as the functional Applying (10), we obtain the integral formulation of the primal problem (22), that is, It is easy to see that the goal functional related to our target quantity (21) is given by Therefore, the interaction energy (21) can be written as For further analysis, it is necessary to express the adjoint problem in an integral formulation similar to (24). Recalling the definition (17) and taking into account that in our case A * = A we obtain where z is the solution of (27) and is defined to be equal to zero on the boundary ∂Ω. Here, we would like to emphasize the difference between the primal and adjoint problems: the solution φ 2 of the primal problem is generated by the charges of the molecule M 2 , whereas the solution z of the adjoint problem is determined by the charges of M 1 (see (23) and (25)). Now, by taking into account the definitions (23) and (25) of ℱ and ℒ, respectively, and of the potential z, it is clear that the interaction energy E 1-2 (=E 2-1 ) can also be calculated as The optimal adaptive algorithm balances the error distribution of the calculated potential over the domain Ω in such a way that for a given number of DOFs the error in the goal quantity with e: = φ 2 À φ 2,h and e i : = φ 2 (x i,1 ) À φ 2,h (x i,1 ) is minimal. In general, this minimization problem cannot be solved exactly, but adaptive algorithms aim at approximating its solution well.
A difficulty in the problem under consideration is that the goal functional defining the quantity of interest, that is, the electrostatic interaction, requires point evaluations. Hence, it is not representable in the form (13) through a regular density function j, but defined in terms of the delta functions δ(x À x i,1 ), cf. (25). A standard approach in such a situation employs regularization of the goal functional by replacing the point evaluations with averages over small balls B(x i,1 , ρ), where ρ is the radius around the points of interest x i,1 . The application of this approach, however, changes the goal functional when the averaged functions are not harmonic and requires the evaluation of integrals of discontinuous functions over balls. For achieving high enough accuracy, the mesh has to be a priori adapted around the points of interest and special integration rules have to be used.
In order to avoid these conceptual shortcomings, here, we derive a nontrivial representation of the error in the goal quantity which does not require averaging and directly exploits the original goal functional.
The purpose of our goal oriented error indicator is to detect the regions for which the error in the electrostatic potential contributes the most to the error in the goal quantity, that is the interaction energy. Those regions are then being refined and are not necessarily the regions where the pointwise error or the so-called energy norm of the error in the potential are the highest. This is in contrast to the functional type a posteriori error estimates developed in Ref. [23] or the residual-based error indicator used, for example, in the fe-manual routine of APBS, which aims at the reduction of the global (integral) energy norm of the error in the potential.

Error indicator
The following approach does not involve any splitting of the full potential and is based on a direct discretization of the problem defining φ 2 . Hence, the primal and adjoint problems are given by the expressions (24) and (27), respectively. The goal quantity is defined in (26), and the approximate interaction energy is simply where φ 2,h is the solution of the discrete primal problem Here, φ 2,h is a continuous piecewise linear function from the finite element space V 1 h defined on the mesh T h . By taking into account (26) and (30), and the fact that the goal functional ℒ is linear, the error in the interaction energy is obtained as Furthermore, to obtain an error estimator in terms of the solution to the adjoint problem, we test (27) with the function φ 2 À φ 2,h and (24) with the adjoint solution z (for a proof that one can do this, see Proposition 5.12 in Ref. [51]), and similarly to the general approach presented by (18) and (19), we obtain where in our particular case we have In (34), h is a FE approximation of the adjoint solution z found in a richer space than the one which φ 2,h belongs to. In particular, we find it in the space V 2 h consisting of continuous piecewise quadratic functions defined on the same mesh T h (see Eq. (S1) in the Supporting Information). Furthermore, I h is the nodal interpolation operator in V 1 h , N V is the number of nodes (vertices) in the mesh T h , and ψ i f g NV i¼1 are the nodal P 1 basis functions in V 1 h equal to 1 at the i-th node and zero elsewhere (see Figure 3). The nodal interpolant I h f of a continuous h can be taken as a test function in (31)). The functions ψ i f g NV i¼1 form the so-called partition of unity (the "PU" in η PU i ), that is, their sum equals one in Ω.
cators that are used to drive the adaptive mesh refinement. Note, that the first term in the definition (34) of the error indicator on the ith node patch can be nonzero only if the support of ψ i contains a partial charge from molecule M 2 (recall that the charges of M 2 define the primal problem). In this case, the mesh cells K in the notation of (16) and (19) are just nodewise patches of tetrahedrons corresponding to the supports of the P 1 basis functions ψ i f g NV i¼1 . So far, we have constructed the error indicators η PU i φ 2,h , z 2 ð Þ h that can guide our adaptive mesh refinement. To obtain an approximation of the goal quantity, that is, the electrostatic interaction between M 1 and M 2 , we could directly use the solution φ 2,h of the primal problem and the expression (30). However, it is better to make use of the approximate solution z 2 ð Þ h of the adjoint problem and the expression (28). Indeed, one can show the following equality: This expression shows that the interaction E Adj 1À2 computed by means of the more accurate z 2 ð Þ h is equal to the interaction E Prim 2À1 corrected with the approximately evaluated error estimator hρ(φ 2,h ), s h (z)i. The quantity E Adj 1À2 is a more accurate approximation of the electrostatic interaction E 1-2 (see (28)), which will be used in our applications. It can be shown that, the theoretically optimal convergence rates for the error in the goal quantity using the primal and adjoint solutions are (see eqs. (5.127) and (5.128) in [51]) where #DOFs denotes the number of degrees of freedom, that is, the number of unknowns for the primal or adjoint problems. From this expression, it is clear that the interaction energy obtained with the adjoint solution z 2 ð Þ h converges faster to the exact interaction energy than the one with the primal solution. We will demonstrate these convergence rates in Section 3.
We will refer to the corresponding solver as ARGOS and the solution obtained by ARGOS will be denoted as E ARGOS ≔ E Adj 1À2 . If, for the definition of the primal problem (22), we use the charges of molecule M 1 instead, then the definition of the adjoint problem (27) will involve the charges of molecule M 2 . In this case, the numerical approximations of the interaction energies E 1-2 and E 2-1 obtained by ARGOS are denoted as E Prim 1À2 and E Adj 2À1 , respectively. The reasoning here is that this time the energy E 1-2 is computed by using the solution of the primal problem, whereas E 2-1 is obtained from the more accurate solution of the adjoint problem.
We note that, a similar approach to the estimation of the error in the goal quantity E can be also executed in conjunction with the splitting of the full potential φ 2 = G 2 + u 2 , where G 2 is the Coulomb potential of the charges of M 2 in a uniform dielectric medium with ϵ = ϵ m , having the form (11), and u 2 is the corresponding reaction field potential. An adaptive algorithm based on a goal oriented error estimation together with a splitting of the potential which eliminates the grid artifact due to the infinite Coulomb self-interaction can be very reliable for the calculation of the solvation energy (12). This feature will be included in future releases of ARGOS.

| Adaptive algorithm
In general, adaptive algorithms can be represented by the following flow chart ( Figure 2). Instead of the mesh parameter h, we will use a subindex ℓ = 1, 2,. … to denote the refinement level and we recall that The adaptive finite element solver based on the error estimator described above can be summarized as follows.

Solve
The solution φ 2,ℓ from the space V ℓ 1 , defined on the current mesh T ℓ , is found by solving (31).

Estimate
Next, we find the solution z ℓ (2) of the adjoint problem by solving the finite dimensional problem (S1). From z ℓ , we compute the nodewise

Mark
Based on the values of η PU i , a set of node patches N ℓ is marked for refinement. This set can be chosen through different marking strategies (see the Supporting Information).

Refine
Having performed a refinement based on the marked patches N ℓ , a new mesh T ℓþ1 is obtained.
The detailed algorithm of ARGOS reads: : until convergence or maximum number of refinement levels reached

| Finite element method
The finite element (FE) method aims at approximating the exact solution of a differential equation by a piecewise polynomial function. For this purpose, the computational domain is split into individual small (finite) elements, such as triangles in 2D or tetrahedrons in 3D. The construction is such that the FE space used for the approximation is finite dimensional with a basis consisting of locally (on patches of elements) supported piecewise polynomials. The latter have to satisfy certain interpolation conditions, for example, continuity across element interfaces, as used in the adaptive algorithms presented in this paper.
In case of the lowest approximation order, that is, piecewise linear polynomials, each basis function is associated with a node in the mesh, that is, a vertex of an element. The approximate solution is represented in the chosen basis with unknown expansion coefficients, for example, the nodal values of the approximate solution. Based on a variational formulation of the differential equation, taking into account boundary conditions, a system of algebraic equations is derived for these expansion coefficients.
In our case, the FE space V 1 h for the primal problem is based on the Lagrange P 1 element and consists of continuous piecewise linear functions. The FE space V 2 h for the adjoint problem is based on the Lagrange P 2 element and consists of continuous piecewise quadratic functions. In the case of V 1 h , there is one basis function ψ i for each vertex (node) V i of the mesh T h , which is equal to 1 at V i and 0 at all other vertices, see Figure 3. Here, we will outline the construction of the system of linear algebraic equations arising from (31). The system for the adjoint problem can be constructed in a similar way. Recall, that Equation (31) has to be satisfied for all test functions v V 1 h . Equivalently, this means that the equation must hold for all basis functions ψ i , i = 1, …, N Ω that are zero on the boundary ∂Ω. In this notation, the number of nodes N V in the mesh T h equals the number of degrees of freedom N Prim DOFs (=N Ω ) plus the number of nodes N ∂Ω on the boundary ∂Ω. Hence, in the primal problem, we have N Prim DOFs ¼ N Ω and we obtain N Prim DOFs equations for φ 2,h : Moreover, φ 2,h belongs to V 1 h and φ 2,h = 0 on ∂Ω. Therefore, φ 2,h can be represented as a linear combination of the basis functions ψ j which are equal to zero on the boundary ∂Ω with unknown expansion coefficients φ j 2 , that is, The goal is to find the vector of unknowns corresponding to the interior vertices. For this purpose, we substitute expression (39) for φ 2,h in (38) and use the linearity of the integral, gradient operator, and right-hand side to obtain the following system of linear equations F I G U R E 2 Schematic representation of adaptive refinement as explained in the text F I G U R E 3 Interior P 1 basis function in 2D In matrix notation, the above equations read where A is the symmetric positive definite matrix with entries The system (40) is then solved either by a direct or by an iterative method. In our numerical tests we have used the conjugate gradient method with a sufficiently high accuracy in order to make the error due to the iterative solution of the linear system subordinate to the overall error. In this way, it is possible to study the effect of the adaptive mesh refinement in a reliable way.

| Mesh construction
A very important step in the modeling of the electrostatic potential in biomoleclar systems through the PBE is the definition of the solutesolvent interface. Here, two approaches are possible. The first one, and much more commonly used [8, 9, 11-13, 17, 18], is where the transition between the low dielectric region inside the solute cavity and the high dielectric region of the surrounding continuum is sharp and discontinuous. In the second approach, this transition is smooth and depends on the (smoothly varying) electronic charge density of the solute [50,[55][56][57]. As mentioned in Section 2.1.1, in the present work, we utilize the first approach in combination with a solute charge density described by discrete point charges at atomic positions. In this approach, one needs to define the solute-solvent interface, that is, the molecular surface separating the solute interior from the solvent.
There are many possible ways to define this surface. In all our calculations, as a definition of the molecular surface we take the SES with a To generate a triangulated surface mesh which approximates the SES and/or the van der Waals surface, we use the meshing software NanoShaper [58]. It also has the functionality to detect pockets and cavities in proteins which are defined according to their volumes, compared to the volume of a single water molecule. If the cavity has a volume less than that of a water molecule modeled as a sphere with a radius of 1.4 Å, then it is assigned the dielectric coefficient of the protein. Similarly, the pockets are regions in space that can be accessed by a water molecule, but cannot be accessed by a virtual probe of bigger size, for example, 3 Å. The dielectric constant of water is assigned to these pockets. The

Volume mesh generation and mesh refining
Once the surface meshes are generated, we use TetGen [61] to create the initial tetrahedral meshes, where one can optionally prescribe a uniform mesh size in the molecular domain, related to the average edge length on the molecular surface. Since, our FE algorithm is adaptive, prescribing a uniform mesh size inside the protein is not necessary in order to properly resolve the singularities associated with the fixed point charges. Instead, the mesh will be appropriately refined after a few mesh refinement steps. For the adaptive mesh refinement (AMR) we use mmg3d [59,60]. The implementation of the prototype version of the adaptive solver is realized in FreeFem [62]. Here, we want to point out that the particular choice of a molecular surface definition and of the IEL, as well as of the utilized meshing software is not essential and other alternatives can also be used.

| MEAD solver
Here, we will give a brief overview on the Macroscopic Electrostatics In our calculations, we used MEAD version 2.2.8a and its programs potential and multiflex. The former is applied to calculate the ESP at the coordinates specified by the user. We used multiflex to calculate the interaction energy. An important issue in the MEAD procedure involves calculating both energies E 1-2 and E 2-1 , expressed through potentials φ 1 (see (20)) and φ 2 (see (21)). The final result is obtained as the average of the two, that is, The implementation of MEAD allows the user to set model parameters including the dielectric coefficients inside and outside the molecular cavities, the ionic strength of the solvent, the grids dimensions and spacings, as well as to revert to the old way of testing for convergence of the iterative solver Successive Over Relaxation (SOR).
The latter is a variant of the Gauss-Seidel method for solving a linear system of equations. The algorithm is controlled by the following parameters, hidden from the user: minits, maxrmsdiff and maxiter. The parameter minits specifies the minimum number of iterations and is defined as 30(N d À 1)/2, where N d is the number of grid points in each coordinate direction for both the coarse-grained and the fine mesh. The parameter maxrmsdiff is the tolerance for the maximum root mean square of the difference (rmsdiff) between two consecutive iterates φ k and φ k + 1 generated by the SOR method.
By default, maxrmsdiff is set to 2 Â 10 À5 /(N d À 1). The parameter maxiter denotes the maximum number of iterations and is set to 10Âminits. The SOR iterative method is terminated when either the maximum number of iterations is reached, or the current iteration number is greater than or equal to minits, while rmsdiff is less than or equal to maxrmsdiff. It turns out that the default values of these parameters are inappropriate in the majority of our tests. Therefore, in this paper, we also present results of MEAD for a set of manually improved values of SOR parameters, which requires a recompilation of the source code. The selected values of these and other parameters will be given for every problem separately. It is important to note that MEAD performs the calculations in singleprecision arithmetics.

| APBS solver
Similar to MEAD, the Adaptive Poisson-Boltzmann Solver (APBS) solves the linearized PBE 2 , applying an FD method [13]. However, unlike MEAD, it employs double precision arithmetics in combination with a multigrid iterative solver and a better stopping criterion.
In analogy with the MEAD calculations, for each frame we perform two computations with APBS: one to find the potential φ APBS  (20) and (21), respectively. The final result is obtained by the average of the two, in analogy with (41): To present the best case scenario for APBS, we use the finite difference solver mg-auto with parameters which provide a relatively good computational efficiency and a very high accuracy.
Below, we list the keywords, most relevant to our computations, from the input ".in" file of APBS, where "#" denotes the beginning of a comment in the input file. The full list of keywords together with particular values of the parameters they define is given in the Supporting Information.
mg-auto # specifies the type of calculations that the finite difference multigrid PMG code performs.
dime nx ny nz # are the number of grid points in the x-, y-, and zdirection, respectively, for all focused FD grids.
cgcent mol 1# indicates that the coarse grid is centered at the center of the minimal box, with sides parallel to the coordinate axes, such that it contains all the atoms' centers of the system under study (mol 1 is its internal number).
cglen xlen ylen zlen # defines the computational domain which is a rectangular parallelepiped with side lengths xlen, ylen, and zlen in Å, in the x-, y-, and z-directions, respectively.
fgcent mol 1# indicates that all focused grids, including the finest one, are centered at the center of the minimal box, with sides parallel to the coordinate axes, such that it contains all the atoms' centers of the system under study (mol 1 is its internal number).
fglen xlen ylen zlen # specifies the fine mesh domain, a rectangular parallelepiped with side lengths xlen, ylen, and zlen in Å, in the x-, y-, and z-directions, respectively.  The standard CHARMM force field (version v35b3) [64] was used for polyproline. The force field parameters of the Alexa chromophores were created by an analogy approach from that of similar chemical groups in the CHARMM force field (for details on the MD simulations including force field parameters see Ref. [44]). The water molecules in the MD simulations were described, using the TIP3P model [65].

| Systems under study
The solver will also be applied to compute the interaction The preparation of the system included the following steps. First, CHARMM GUI was used to add hydrogen atoms to the heavy atoms resolved in the crystal structure 6m0j. In addition, the system was placed in a box of water molecules and explicit ions were added to obtain an overall charge neutral system. The geometry of the system was optimized by using an energy minimization constraining the heavy atom coordinates of the glycans and the protein backbone as in the crystal structure. The calculation of the interaction between the Alexa chromophores and the RBD and ACE2 proteins by MEAD, APBS and ARGOS follows the same procedure, starting with computing the ESP from the Poisson or LPBE (see (5) and (4)). This potential is then used to calculate the coupling (see (20) and (21)).

| Evaluation of the accuracy of the numerical methods
We start with three models, for which analytical solutions for the ESP exist. By comparing the numerical with the analytical solutions, the sources of numerical errors and their magnitude will be determined.
The numerical accuracy of the solvers (MEAD, APBS and ARGOS) is evaluated by the relative errors (in %) defined as where E is the exact (analytical) solution, Err Adj = Err ARGOS and Err Prim correspond to the relative errors of the numerical solutions  Table 1 Table 2.
The computational domain for ARGOS is a cube with a side length of around 160,000 Å for all frames with POLY6 and around 400,000 Å for all frames with POLY20. To demonstrate the convergence of the ARGOS results with respect to the number of degrees of freedom (DOFs), we compute the relative errors for the interaction energies obtained from the solutions of the adjoint and primal problems, Err ARGOS and Err Prim defined in (43), and average them over all frames ( Figure 5). Both Err Adj and Err Prim exhibit the order of convergence, predicted by (36) and (37). The solution of the adjoint problem exhibits a higher convergence order O(#DOFs À1 ) (see (37)), and thus, much smaller relative errors, ranging from 2 Â 10 À1 % down to 2 Â 10 À4 % for the largest number of DOFs, that is, 10 7 . For a nonhomogeneous Dirichlet boundary condition, the relative errors   (43)) computed by ARGOS as a function of the number of DOFs (averaged over the same 19 frames). The dashed lines represent the optimal convergence rates for the finite element spaces used, that is, P 1 for the primal and P 2 for the adjoint problem (see (36) and (37) Tables 3 and 4, respectively. All grids (coarse-grained and focused) for both MEAD and APBS are centered at (0, 0, 0), that is, slightly offset from the center of the ion. In the case of MEAD, the calculations were performed for two different sets of SOR parameters (Table 3). APBS yields very similar results for calculations with etol 1e À 06 and etol 1e À 09.
Therefore, we present only the latter results.
In the ARGOS calculations, we used GridScale An illustration of the mesh for GridScale = 2 before and after refinement is presented in Figure 7. Comparison of the initial and final mesh topologies shows significant refinement in the critical areas: around the ion surface, the points, where the ion charge and probe charge are located, and in between the two charges. The convergence of the relative errors with increasing refinement levels of ARGOS is shown in Figures S12 and S15.
The analytical solution for the interaction energy between the ion and the probe charge follows from the well-known expression for the electrostatic potential φ: The relative errors of the calculations obtained with MEAD, APBS and ARGOS using different solver configurations are given in For an ion radius R = 1 Å, the relative errors of MEAD and APBS are approximately 60 % at the interface between the two dielectrics (see Figure S17). This result shows that the presence of charges near sharp interfaces has a strong impact on the accuracy of the FD approximation of the electrostatic interaction. The relative error for ARGOS is also maximal at the interface, but it does not exceed 1 % F I G U R E 7 ARGOS (GridScale = 2) mesh refinement (from left to right) for the solution of the primal and adjoint problems for Model-2. Here, the ion sphere radius is 2 Å and the probe charge is at a distance r = 15 Å from the center of the spherical ion. The color code corresponds to the solution of the adjoint problem (27) for the ESP in units of 10 À3 e 0 /Å for, both, GridScale = 2 and 4. The error for r > 1 Å is around 10 À2 % for, both, GridScale = 2 and 4. The convergence of the relative error for ARGOS with respect to the mesh refinement levels is shown in Figures S12, S15, and S16.
We also recall that the results obtained by MEAD and APBS for a   (43)) for MEAD and APBS averaged over 10 directions d ! i between ion and probe charge (Eq. (S4)) together with the relative errors for ARGOS as a function of the distance r between the ion and the probe charge. The ARGOS calculations were performed with MRL = 5 and MRL = 3 for GridScale = 2 and 4, respectively F I G U R E 9 Same as in Figure 8, but for Model-3 and performing the ARGOS calculations with MRL = 3 and 4 for GridScale = 2 and 4, respectively by MEAD and APBS. Similarly to the results presented in Figures 8   and 9, the relative errors for ARGOS are more than two orders of magnitude smaller around the interface than those of MEAD and APBS ( Figure S18).

| Model-3: Presence of interfaces and ions
We supplemented Model-2 with a continuous charge density of ions outside the cavity, described by an ionic strength of 0.3 M and taking into account the ion exclusion layer (IEL). Here, the IEL is enclosed between the ion sphere and another concentric sphere with a radius enlarged by that of a probe ion's radius R ion = 2 Å. For MEAD and APBS we used the configurations, defined in Tables 3 and 4, respectively. The APBS results for etol 1e À 06 and etol 1e À 09 again do not differ, and thus we report only the latter results. The computational domain for ARGOS is the same as for Model-2. We compare the solutions of ARGOS obtained at MRL = 3 and 4 for Grid-Scale = 2 and 4, respectively, with those obtained by MEAD and APBS. The convergence of the ARGOS relative errors with increasing mesh refinement levels is shown in Figure S22. The analytical solution for the ESP φ in the present model is given as: where r is the distance of the probe charge from the center of the Born ion and a = R + R ion and k ¼ In the next sections we will use ARGOS as a reference to explore the accuracy of the MEAD and APBS solvers in computations of interaction energies in real systems.

| Interactions between chromophores in a dielectric environment
The interaction energy between two chromophores, described by molecule-shaped cavities with atomic partial charges, in a dielectric environment cannot be computed analytically. Therefore, we exploit the high accuracy of the ARGOS results and determine the deviations, D AM and D AA , of the MEAD and APBS calculations, respectively, relative to the solution E Adj 1À2 of the adjoint problem of ARGOS as: Refinement level ℓ = 4 and GridScale = 2 are employed in the ARGOS calculations. The computational domain is a cube with a side length between 160,000 and 400,000 Å depending on the size of the smallest box containing the molecular system. For MEAD and APBS we used configurations M1-2 and A1 given in Tables 1 and 2, respectively.
We start with the coupling between the charge densities of the electronic ground states of the two chromophores (see Tables S3 and   S4 for the respective partial charges) in equilibrium with the solvent. To study the impact of the description of interfaces, we investigate different MD frames for POLY6 and POLY20, as in Model-1 studied above, but using the following models with different types of interfaces: •  Figure 10C).
The Alexa chromophores in their electronic ground state are negatively charged by two elementary units (À2e 0 ) [44]. In order to investigate, how the interchromophore charge density coupling changes, if no monopole moments are present, we artificially neutralize the chromophores by adjusting the charges within the SO 3 groups such that the overall charge of the chromophores is zero, as described in Tables S5 and S6. Please note that, this neutralization is not realistic for the present chromophores but represents a common situation for other chromophores, which would, of course, differ in the details of the charge density. Note also that, these charges have been employed in Model-1 above.

| Charged dyes
The deviations D AM and D AA of the electrostatic interaction energies between the charged chromophores obtained with MEAD and APBS and the ARGOS reference solutions are presented in Figure 11 for three different models described above.
Examination of the MEAD results shows, that including the polyproline domain can notably change the error (e.g., frames 5, 9 and 10 in Figure 11), but it still stays small. The maximum value of D AM does not exceed 3 %. Furthermore, for all frames, except frame 16, including the ionic strength of the solvent results in an increase of D AM . For example, the relative deviations for frames 5 and 9 have increased 100 times. To understand these effects we investigate the reference interaction energies obtained with ARGOS (see Figure 12, and for logarithmic scale of the y-axis see Figure 13).
As evident from Figure 12, the difference between interaction energies, computed for different models, is more pronounced for POLY6 (frames 1 to 10). Introducing molecular cavities for the dyes and polyproline, increases the interaction energy in comparison to the simple model of a homogeneous dielectric. This effect is due to the lower dielectric coefficient inside the cavities. For larger interchromophore distances this effect gets smaller.
F I G U R E 1 0 Interfaces in different models explained in Section 3.2 in the text. The conformation of the two dyes and the shorter polyroline helix POLY6 corresponds to frame 10 from the MD simulation mentioned in the text. Figures generated using VisIt [67] F I G U R E 1 1 Relative deviations of the interaction energies of charged chromophores, obtained with MEAD (green symbols) and APBS (blue symbols) with respect to the reference ARGOS solutions. The electrostatic energies are calculated for MEAD configuration M1-2 from Table 1 and APBS configuration A1 from Table 2. The circles, squares, triangles correspond to the models DYES, DYES-POLYPRO, DYES-POLYPRO-IONS, respectively, that are explained in the text. The different frames correspond to different molecular geometries taken from MD simulations The results, shown in Figure 13, demonstrate a significant screening of the interaction by the ions, resulting in a decrease of the interaction energy for all frames. The magnitude of the screening effect for POLY6 can vary significantly, reducing the interaction energy by factors between 1.5 and 100, depending on the conformation of the system (see Figure 13). For the larger distances between the chromophores in POLY20, the magnitude of the screening effect varies stronger: between 10 2 to 10 6 times (see the results for frames 11 and 12 in Figure 13). The large increase in screening for POLY20 with respect to POLY6 can be explained in the following way. When the charged chromophores are far apart, they can be approximated by Born ions, immersed in the dielectric continuum with ionic strength. In this case, the interaction energy can be obtained analytically from the potential in (45). It includes an exponential factor, which imposes a highly nonlinear decline with increasing distance. To estimate its influence, we compare the screening effect for frames 12 and 14, which correspond to interchromophore distances of 78 and 42 Å, respectively. The ratio between the interaction energies, calculated by ARGOS, is about 10 À4 (see Figure 13 for the model DYES-POLYPRO-IONS). The estimation of this ratio, using (45) with k 2 ¼ 2:529 Å À2 , gives 9 Â 10 À4 , providing a qualitative explanation of the effect of the exponential screening factor. Without exponential factor we obtain a ratio of 0.5. An examination of the results in Figure 13 shows a significant drop in the interaction energies due to ion screening for frames 12, 18 and 19. At the same time these frames exhibit the largest relative deviations for MEAD ( Figure 11).  Figure 12, but with logarithmic scale of vertical axis. Note that, the interaction energies of the first three models cannot be distinguished in the present representation the single values (Table 5). Also, the average over different grid positions improves the accuracy, as seen from Table S9 and from the standard deviations in Table 5. Decreasing the value of the parameter etol improves the accuracy as well. All these effects reduce the relative deviation of the APBS result for frame 12 from the reference (ARGOS) value À5.5959e -10 e 2 0 =Å in Figures 11 to 3 % in Table 5. Reduction of the grid spacing of APBS or equivalently an increase of the number of grid points is expected to further decrease the relative deviation D AA . However, a limiting factor for this strategy would be primarily the memory capacity required to store the coefficients and the potential values at all grid points.

| Neutral dyes
The relative deviations for the interaction energies between ground state charge densities of neutral chromophores, presented in Figure 14, reveal a drastic increase compared to the previous case of charged chromophores in Figure 11. For MEAD computations, the relative deviations exceed 10 % for most of the models and frames.
Moreover, for some frames MEAD produces a different sign of the interaction energies compared to ARGOS (see frames 15 and 10 for M1-1 in Table 6), and even a different sign of E MEAD 1À2 and E MEAD 2À1 (see frames 15 and 16 for M1-1 in Table 6). Obviously, a strong asymmetry is an indicator of a poor accuracy. In contrast to the MEAD calculations, the energies E Adj 1À2 computed by ARGOS match the energies E Adj 2À1 very well. The respective values are even computed on two different sequences of adapted meshes, both starting from one and the same initial mesh, demonstrating again a very good symmetry, and thus the superior performance of the FE-solver. The Table 6 were performed with GridScale = 2 at MRL = 4. Increasing the GridScale to 4 in the calculations for frames 15 and 10 leads to very similar energies E Adj 1À2 and E Adj 2À1 of 1.9649e À 07 and 1.9648e À 07, and 2.7915e À 07 and 2.7911e À 07, respectively.

ARGOS calculations in Figures 14 and 15 and in
Similarly to the previous case, the APBS calculations deviate The interaction energies for frames 6 and 10 are analyzed in Tables 7 and 8, respectively (for more details see Tables S10 and S11).

Like in the case of charged dyes, the values of E APBS
1À2 and E APBS 2À1 for all the tested grid spacings do not exhibit the required symmetry The effect is even more pronounced for particular relative positions of the molecular system and the grid, when E APBS 1À2 and E APBS

2À1
can have different sign (see the Supporting Information, grid centering perturbation i = 10 in Table S10 for h = 0.2 Å of frame 6). Moreover, the asymmetry of the solutions, as well as the standard deviations Even at the smallest grid spacing the variation of computed energies with respect to small perturbation of the grid position remains substantial. On the other hand, the convergence also indicates that the ARGOS solution is indeed accurate. Similar to the case of charged dyes, our analysis demonstrates that one calculation per frame is in general not sufficient to obtain a reliable result. Instead, an averaged T A B L E 5 Sensitivity of the APBS results to the etol parameter and the relative position of the grid and the molecular system for frame 12 for model DYES-POLYPRO-IONS with charged chromophores 0.286 (etol 1e-06) À4.82e À 10 1.80e À 11 À6.94e À 10 1.23e À 11 À5.88e À 10 5.08 0.286 (etol 1e À 09) À5.79e À 10 1.79e À 11 À5.79e À 10 1.16e À 11 À5.79e À 10 3.45 Note: σ 1À2 and σ 2À1 are standard deviations of interaction energies obtained for different grid positions. All energies are in units of e 2 0 =Å.

F I G U R E 1 4
Relative deviations of interaction energies of chromophores in the charge-neutral ground states, obtained with MEAD (green) and APBS (blue) with respect to the ARGOS results. MEAD configuration M1-2 from Table 1 and APBS configuration A1 from Table 2 were used. The circles, squares, triangles correspond to the models DYES, DYES-POLYPRO, DYES-POLYPRO-IONS, respectively. The different frames correspond to different molecular geometries taken from MD simulations value, obtained from a series of calculations with variations of the grid centering and grid spacing h is required. This example also shows that if one wishes to use a finer final focused grid, it is required to considerably increase the values of the dime parameter of APBS (see Table S10) which, on the other hand, increases the computational costs and even more so the memory requirements.
Both FD solvers do not ensure the natural symmetry condition.
Reducing the grid spacing helps to improve the accuracy, but not to the level obtained for the charged chromophores. In contrast to the latter, in the present case, the leading term in the interactions is due is, ϵ s À ϵ m j j .
In Table 9, we summarize five setups for the DYES model with different dielectric coefficients ϵ m and ϵ s where the Alexa chromophores are prescribed neutral ground state or transition charges (see Tables S1 and S2 for Table S16) cut down the relative deviation D AA between APBS and ARGOS to 0.5 %, which is again an evidence that the ARGOS reference solutions are indeed accurate.
The smaller relative deviations obtained for the coupling between transition charges can be explained in the following way: in the case of transition charges the potential is created by a smaller number of nonzero charges (see Tables S3 and S4, and S1 and S2 for the groundstate and transition charges, respectively, for the two dyes) compared to the case of ground-state charges, and then multiplied again by a smaller number of nonzero charges 4 to obtain the interaction energy.
In order to describe a situation relevant for FRET, the dielectric coefficients ϵ m and ϵ s of the dyes with transition charges and the solvent, respectively, have to be decreased from ϵ m = 4 and ϵ s = 80 in Setup 2, considered so far, to ϵ m = 1 and ϵ s = 2 (Setup 5 in Table 9).
The electrostatic interaction energies in this case have been assessed with MEAD in our previous publication [44]. Here, we reinvestigate these results for a large number of frames of POLY6 ( Figure S27). For most frames of POLY6, the relative deviation of MEAD from ARGOS stays below 3 %, while for some the error can reach up to 20 %. The large deviations, however, occur only for small couplings and, therefore, do not affect the FRET efficiencies calculated previously [44].
Hence, it is not the calculation of the coupling between transition densities but that between the ground state charge densities of the chromophores that can be critically improved by ARGOS. The latter aspect will be crucial in an implicit description of the solvent in MD simulations that will allow one to simulate much larger time scales than by using an explicit description. In this way the sampling of chromophore configurations can be improved which is critical for the interpretation of FRET experiments. An important aspect of an implicit solvent description concerns the screening of the interchromophore coupling by the ions, which is investigated in the following.

Screening by ions
In order to investigate the ionic screening effects in more detail, we applied ARGOS to calculate the reference interaction energies for all interface models. Two additional models without dielectric boundaries (ϵ m = ϵ s = 80), but including the ion exclusion layer (IEL) around different parts of the molecular system, were introduced: • DYES-IONS includes ionic charge density in the aqueous solvent,

| Interaction between SARS-CoV-2 protein and cell receptor ACE2
Now, we present the electrostatic interaction between the RBD of the spike protein of SARS-CoV-2 and the cell receptor ACE2 (see Section 2.4) calculated with APBS and ARGOS. In order to investigate how the electrostatic coupling between RBD and ACE2 changes if LYS417-ASP30 or ASN501-LYS353 interactions are absent (see Figure 4B), we artificially set to zero the charges of the corresponding residues. Thus, in our work we have used three different setups for the charge distribution: • SARS2-WT corresponding to the crystal structure; • SARS2-NK with neutralized charges on ASN501 and LYS353; • SARS2-KD with neutralized charges on LYS417 and ASP30.
Without loss of generality, let M 1 denote the ACE2 receptor and M 2 denote the RBD of SARS-CoV-2. The total charge in M 1 and M 2 for the three setups SARS2-WT, SARS2-NK, and SARS2-KD is À27 e 0 and 2 e 0 , À28 e 0 and 2 e 0 , and À26 e 0 and 1 e 0 , respectively. We have calculated the electrostatic couplings E 1-2 and E 2-1 (Eqs. (20) and (21) Table S17.

Surface and volume meshes in ARGOS
Concerning the surface mesh in ARGOS, now we make use of the external routine mmgs_O3 to simplify the molecular model, unlike in the case of the Alexa dyes where we directly used the meshes generated by NanoShaper using its built in Laplacian smoothing. It would be best to apply mmgs_O3 to the nonsmoothed by NanoShaper raw surface mesh ( Figure 20A) since this raw mesh has its vertices analytically sampled from the SES. However, with the current version 5.5.2 of mmgs_O3 this was not possible due to some really bad quality, needle-like triangles ( Figure S28) that prevented the completion of the procedure. Therefore, we simplify the surface already smoothed by NanoShaper. Nevertheless, the final mesh ( Figure 20C) after the application of mmgs_O3 is still very close to the raw surface mesh on Figure 20A. This is seen by computing the (symmetric) Hausdorff (H-) distance (0.16 Å) and the mean H-distance (0.018 Å) between these two meshes (see Eqs. (S12) and (S13)). For the case of nonzero ionic strength I s , the outer surface of the IEL is constructed following a similar procedure as for the SES, that is, we first construct the raw van der Waals surface with R ion = 2 Å in NanoShaper to which we apply its built-in Laplacian smoothing and then we simplify the result with mmgs_O3 ( Figure S31).
In this example, the molecular system contains a total of 7 cavities, detected by NanoShaper in the surface mesh construction phase (see Figure S29). Since the formed cavities are relatively small, the computed by ARGOS interaction energies are practically the same for an assignment of their dielectric coefficient as the dielectric coefficient ϵ s of the solvent or ϵ m of the protein (see Tables S19 and S20).
The tetrahedralization of the inside and outside of the molecular surface mesh is performed using TetGen without an edge size constraint inside or outside the molecular model. This means that our volume mesh only conforms to the surface mesh on Figure 19A, but is irregular inside and outside, becoming coarser and coarser away from the molecular surface ( Figure 19B). We point out that the coarse initial mesh inside and outside the molecular region does not affect the accuracy of ARGOS, because it is being refined by the adaptive algorithm precisely where it is needed to resolve the problematic regions contributing the most to the error in the electrostatic coupling ( Figure 21). The CPU times to construct and simplify the surface and volume meshes are given in Table 13 in Section 3.4.

Results of APBS and ARGOS
The results of APBS and ARGOS for the three setups SARS2-WT, SARS2-NK, SARS2-KD with I s = 0 M and I s = 0.05 M are given only for ARGOS in Table 10 (for the detailed results, see Tables S18,   S20-S22). The comparison of the data demonstrates a high accuracy of the averaged APBS values that is likely due to the fact, that both F I G U R E 1 9 (A) Surface mesh S 3 of the molecular system RBD of SARS-CoV-2 and the ACE2 receptor. The red rectangle represents the zoomed-in part of the system shown on Figure

| CPU times
Although ARGOS [43] is currently in a prototype version (implemented in FreeFem [62]) and is currently undergoing a substantial performance improvement, we present first CPU times in comparison to APBS obtained on a machine with Intel(R) Core(TM) i5-7500 CPU @ 3.40 GHz and 16 GB of RAM. Since the primal and adjoint problems are discretized using the linear and quadratic Lagrange finite elements, respectively, the number of DOFs, N Adj DOFs , in 3D for the adjoint problem is approximately 8 times larger than that for the primal. Here, we note that for a Dirichlet BC the number of DOFs for the FE space V 1 h , N Prim DOFs , is exactly equal to N V À N ∂Ω , where N ∂Ω is the number of mesh vertices on the boundary. Consequently, the major part of the CPU time per MRL for ARGOS is spent for the solution of the discrete adjoint problem and the refinement of the mesh for the next MRL. In all the presented examples related to the Alexa 488 and Alexa 594 chromophores, we use a basic Jacobipreconditioned conjugate gradient (CG) method (with a stopping criterion ensuring a very high precision 6 ) as an iterative solver for the resulting linear systems for, both, the primal and adjoint problem. For the system consisting of the RBD of the spike protein of SARS-CoV-2 and the ACE2 receptor, we also use the CG method, but this time with a block-Jacobi preconditioner for, both, the primal and adjoint problem (see Footnote 6). As it was mentioned above, the mesh refining is done using mmg3d v4.0 [60], but a newer version of it can also be used. Any other remeshing software can be utilized as well, which could drastically reduce the CPU time for the mesh refining steps. In Tables 11 and 12  As a first example, in Table 11, we consider Model-2 where the probe charge is placed at a distance r = 10 Å from the center of the Born ion. The total execution time for ARGOS with GridScale = 2 and a CG tolerance of 1e À 09 is 46.9 sec (the corresponding time for a CG tolerance of 1e À 12 is 51.14 sec, see Footnote 6). In comparison, APBS with configuration A2 (see Table 4) for etol 1e-06 takes 658 sec on average.
Note that from this one computation with APBS one can extract the interaction energies for all distances r ≤ 48 Å. However, the size of the final focused grid does not allow to place the test charge at a larger distance, whereas with ARGOS the size of the molecular system is not an issue. A further comparison of the CPU times for GridScale = 2 and 4 in ARGOS is given in Tables S23 and S24.
In Table 12 we give the CPU times of ARGOS for frame 10 of model DYES ( Figure 10A) with neutral ground-state charges (see Figure 17). The total execution time for GridScale = 2 and with a CG tolerance of 1e À 09 is 99 sec compared to 109 and 143 sec of APBS for configuration A1 (see Table 2) with etol 1e À 06 and etol 1e À 09, respectively. However, those CPU times of APBS are only for the computation of E APBS 1À2 , which alone is not enough to obtain a somehow reliable interaction energy and do not include the time of 17 sec needed by the plugin multivalue to read each of the potential ".
dx" files and return the potential at the atom's centers. Taking only

E APBS
1À2 or E APBS 2À1 can result in a very large error (see Table S11, together with the other tables concerning the inspection of problematic APBS frames). Therefore, to compute E APBS (see (42)), it takes in total 218 and 286 sec for etol 1e À 06 and etol 1e À 09, respectively.
For comparison, we report the CPU times for a calculation with nonzero ionic strength of 0.3 M involving frame 12 of model DYES-POLY-PRO-IONS in the case of charged dyes (see Tables 5 and S9). ARGOS with a CG tolerance of 1e À 09 and three MRLs takes 616 sec compared T A B L E 1 0 Electrostatic coupling, obtained with ARGOS for MRL = 2, error tolerance in the CG of 1e À 14, and ϵ = ϵ m inside the cavities, between the RBD of the spike protein of SARS-CoV-2 and the cell receptor ACE2 for the wild type (WT) and two modified systems (NK, KD) whose charges of certain amino acids were set to zero, as explained in the text Note: The results obtained for ϵ = ϵ s inside the cavities are practically the same and are given in Table S19. All energies are calculated in units of e 2 0 =Å. The values of E Adj 1À2 and E Adj 2À1 are obtained on two different sequences of refined meshes, both starting from the same initial mesh in Figure 19B.
to 837 sec necessary to compute E APBS 1À2 with configuration A1 for etol 1e À 09 (recall that 1e À 06 is not enough for this frame), and 747 sec for etol 1e À 06. Again, those times do not include the time of around 32 sec spent by multivalue to read the ".dx" file for each run and return the potential at the specified atom's centers. If the CG tolerance in ARGOS is decreased to 1e À 12, then the time needed increases to 1114 sec, even though the resulting interaction energy is practically the same.
Finally, in Table 13, we report the CPU times for the calculation of the electrostatic coupling for the setup SARS2-KD with vanishing ionic strength. The corresponding times for I s = 0.05 M are given in Table S25. From our experience with the calculations presented in this work, an error tolerance of 1e À 9 in the CG method and 2 refinement steps are always enough to achieve an accuracy 7 within 0.5 %. The solution of the adjoint problem (using a block-Jacobi preconditioned CG with a very low error tolerance of 1e À 14) takes the largest part of the CPU time per MRL, followed by the refinement of the mesh and the assembling of the corresponding stiffness matrix. It is important to note that the mesh quality can be crucial for the convergence and stability of the direct or iterative method used to solve the primal and adjoint linear systems. Therefore, the mesh refining routine should be chosen with care. We also note that using a more reasonable error tolerance of 1e À 9 in the CG, exploiting parallelism and exchanging the Jacobi-preconditioned CG algorithm for a direct solver such as MUMPS [72,73] or an advanced iterative method, for example, a CG with an appropriate preconditioner or a multigrid method, will dramatically improve the performance of ARGOS, which is already comparable or higher now than that of APBS. The CPU time for APBS with the configuration A4 in Table S17 for this same problem amounts to two times (770 + 32) = 1604 sec which is the time for the computation with the charges placed in each of M 1 and M 2 plus the time for two calls of the subroutine multivalue). For comparison, the CPU time for ARGOS with a CG tolerance of 1e À 9 and for MRL = 0, 1, and 2 drops down to 1448 sec.
T A B L E 1 1 CPU times of ARGOS (including the generation of the initial mesh with NanoShaper and TetGen) for Model-2 obtained for different MRLs Note: The probe charge is placed at a distance r = 10 Å from the center of the Born ion. GridScale = 2, a CG tolerance of 1e À 09, and a homogeneous Dirichlet boundary condition were used. Err ARGOS denotes the relative error (Eq. (43)) and t MRL , t Adj mat , t Adj CG , and t ref are the total CPU time per MRL, the CPU time to assemble the matrix for the adjoint problem, the CPU time to solve the adjoint problem, and the CPU time to refine the mesh, respectively, all times given in units of seconds.
T A B L E 1 2 CPU times of ARGOS (including the generation of the initial mesh with NanoShaper and TetGen) in seconds (for definition of the different times see Table 11) for frame 10 of model DYES ( Figure 10A) with neutral ground-state charges, GridScale = 2 and a CG tolerance of 1e À 09  Table S11, obtained from a better quality initial volume mesh.
T A B L E 1 3 CPU times of ARGOS (including the generation of the initial surface mesh with NanoShaper (10 sec), its simplification with mmgs_O3 (17 sec), the generation of the initial volume mesh with TetGen (26 sec), and an additional optimization with mmg3d (105 sec)) in seconds (for definition of the different times see Table 11) for the setup SARS2-KD with an initial volume mesh shown on Figure 19B and with a CG tolerance of 1e À 14  Due to the use of an appropriate (goal-oriented) error estimation in combination with an adaptive mesh refinement the present FE solver requires much less degrees of freedom than MEAD or APBS to obtain a comparable or better accuracy. This advantage, if complemented with parallelization techniques, will make ARGOS a very fast and reliable PBE solver. First benchmark studies show that ARGOS already in its present nonoptimized version is faster and more accurate than the standard FD solvers MEAD and APBS at least for medium sized systems of around 10,000 atoms.

ACKNOWLEDGMENTS
This work was supported by the Doctorate College program "Nano-

ENDNOTES
1 The finite element method can be thought of as a discretization of the weak formulation of a problem. 2 APBS can solve both the PBE and the linearized PBE with the FD method. The application of the finite element method, however, is still at a very experimental stage and in our calculations we could not use it. 3 In the case of exactly prescribed Dirichlet boundary data, the numerical solution converges to the exact solution of the LPBE for the particular triangulation of the spherical cavity. This exact solution is different from the exact solution given by (44) for a cavity which is a perfect sphere, see Figures S12 and S15 for the convergence of the relative error with respect to the MRLs in the case of nonhomogeneous Dirichlet BC. 4 Here by a smaller number of nonzero charges we also refer to the fact that many transition charges in both dyes are either zero or negligible in magnitude relative to a few nonzero ones. 5 For the case with cavities, the mesh inside the molecular region is slightly denser in order for the volume mesh to conform to the surface of those cavities, and the total number of grid points and elements in the mesh becomes 278,128 and 1,520,738, respectively. 6 The reference solutions of ARGOS presented in the paper are obtained with an error tolerance between 1e À 14 and 1e À 12 for the examples related to the Alexa chromophores and an error tolerance of 1e À 14 for the system consisting of the RBD of the spike protein of SARS-CoV-2 and the ACE2 receptor. This means that on the last CG iteration the residual in the linear system, measured in the norm generated by the inverse of the preconditioner matrix C, is smaller than that error tolerance. In particular, C is a diagonal matrix whose diagonal is equal to the diagonal of the system's matrix in the case of the Alexa chromophores and a block-diagonal matrix in the case of the SARS-CoV-2 example. 7 Here, by 0.5 % we mean the relative distance from the reference solution after 5 or 6 refinement steps when no exact solution is available.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.