Improvements to the APBS biomolecular solvation software suite

The Adaptive Poisson-Boltzmann Solver (APBS) software was developed to solve the equations of continuum electrostatics for large biomolecular assemblages that has provided impact in the study of a broad range of chemical, biological, and biomedical applications. APBS addresses three key technology challenges for understanding solvation and electrostatics in biomedical applications: accurate and efficient models for biomolecular solvation and electrostatics, robust and scalable software for applying those theories to biomolecular systems, and mechanisms for sharing and analyzing biomolecular electrostatics data in the scientific community. To address new research applications and advancing computational capabilities, we have continually updated APBS and its suite of accompanying software since its release in 2001. In this manuscript, we discuss the models and capabilities that have recently been implemented within the APBS software package including: a Poisson-Boltzmann analytical and a semi-analytical solver, an optimized boundary element solver, a geometry-based geometric flow solvation model, a graph theory based algorithm for determining p$K_a$ values, and an improved web-based visualization tool for viewing electrostatics.


Introduction
Robust models of electrostatic interactions are important for understanding early molecular recognition events where long-ranged intermolecular interactions and the effects of solvation on biomolecular processes dominate. While explicit electrostatic models that treat the solute and solvent in atomic detail are common, these approaches generally require extensive equilibration and sampling to converge properties of interest in the statistical ensemble of interest [1]. Continuum approaches that integrate out important, but largely uninteresting degrees of freedom, sacrifice numerical precision in favor of robust but qualitative accuracy and efficiency by eliminating the need for sampling and equilibration associated with explicit solute and solvent models.
While there is a choice among several implicit solvation models [1][2][3][4][5][6], one of the most popular implicit solvent models for biomolecules is based on the Poisson-Boltzmann (PB) equation [7][8][9]. The PB equation provides a global solution for the electrostatic potential (φ) within and around a biomolecule by solving the partial differential equation (1) The solvent is described by the bulk solvent dielectric constant ε s as well as the properties of i = 1, . . . , M mobile ion species, described by their charges q i , concentrations c i , and steric ion-solute interaction potential V i . The biomolecular structure is incorporated into the equation through V i , a dielectric coefficient function , and a charge distribution function ρ.
The dielectric is often set to a constant value ε m in the interior of the molecule and varies sharply across the molecular boundary to the value ε s which describes the bulk solvent. The shape of the boundary is determined by the size and location of the solute atoms as well as model-specific parameters such as a characteristic solvent molecule size [10]. The charge distribution ρ is usually a sum of Dirac delta distributions which are located at atom centers. Finally, β = (kT ) −1 is the inverse thermal energy where k is the Boltzmann constant and T is the temperature. The potential φ can be used in a variety of applications, including visualization, other structural analyses, diffusion simulations, and a number of other calculations that require global electrostatic properties. The PB theory is approximate and, as a result, has several well-known limitations which can affect its accuracy, particularly for strongly charged systems or high salt concentrations [7,11]. Despite these limitations, PB methods are still very important and popular for biomolecular structural analysis, modeling, and simulation. Several software packages have been developed that solve the Poisson-Boltzmann equations to evaluate energies, potentials, and other solvation properties. The most significant (based on user base and citations) of these include CHARMM [12], AMBER [13], DelPhi [14], Jaguar [15], Zap [16], MIBPB [17], and APBS [18]. However the APBS and associated software package PDB2PQR has served a large community of ∼27,000 users by creating web servers linked from the APBS website [19] that support preparation of biomolecular structures (see Section 2) and a fast finite-difference solution of the Poisson-Boltzmann equation (see Section 3.1) that are further augmented with a set of analysis tools. Most APBS electrostatics calculations follow the general workflow outlined in Figure 1. An even broader  Figure 1: Workflow for biomolecular electrostatics calculations using the APBS-PDB2PQR software suite. This workflow is supported by the APBS tool suite with specific support by PDB2PQR for preparing biomolecular structures (see Section 2) and support by APBS for performing electrostatics calculations (see Section 3).
range of features and more flexible configuration is available when APBS and PDB2PQR are run from the command line on Linux, Mac, and Windows platforms, and which can be run locally or through web services provided by the NBCR-developed Opal toolkit [20]. This toolkit allows for the computing load for processor intensive scientific applications to be shifted to remote computing resources such as those provided by the National Biomedical Computation Resource (NCBR). Finally, APBS can run through other molecular simulation programs such as AMBER [13], CHARMM [12], NAMD [21], Rosetta [22] and TINKER [23]. General support for integration of APBS with 3 rd -party programs is provided by the iAPBS library [24,25].
This article provides an overview of the new capabilities in the APBS software and its associated tools since their release in 2001 [18,[26][27][28]. In particular new solutions to the PB equation have been developed and incorporated into APBS: a fully analytical model based on simple spherical geometries that treats full mutual polarization accurately, known as PBAM (Poisson-Boltzmann Analytical Model) [29,30] as well as its extension to a semi-analytical PBE solver PB-SAM (Poisson-Boltzmann Semi-Analytical Model) that treats arbitrary shaped dielectric boundaries [31,32]. APBS-PDB2PQR also now includes an optimized boundary element solver, a geometry-based geometric flow solvation model, a graph theory based algorithm for determining pK a values, and an improved web-based visualization tool for viewing electrostatics. We describe these new approaches in the remainder of the paper, and more detailed documentation for these tools is available on the APBS website [19]. As discussed in the following sections, APBS-PDB2PQR offers a wide range of features for users across all levels of computational biology expertise.

Preparing biomolecular structures
Electrostatics calculations begin with specification of the molecule structure and parameters for the charge and size of the constituent atoms. Constituent atoms are generally grouped in types with charge and size values specified by atom type in a variety of force field files developed for implicit solvent calculations [1]. APBS incorporates this information into calculations via the "PQR" format. PQR is a file format of unknown origins used by several software packages such as MEAD [33] and AutoDock [34]. The PQR file simply replaces the temperature and occupancy columns of a PDB flat file [35] with the per-atom charge (Q) and radius (R). There are much more elegant ways to implement the PQR functionality through more modern extensible file formats such as mmCIF [36] or PDBML [37]; however, the simple PDB format is still one of the most widely used formats and, therefore, continued use of the PQR format supports broad compatibility across biomolecular modeling tools and workflows.
The PDB2PQR software is part of the APBS suite that was developed to assist with the conversion of PDB files to PQR format [38,39]. In particular, PDB2PQR automatically sets up, executes, and optimizes the structure for Poisson-Boltzmann electrostatics calculations, outputting a PQR file that can be used with APBS or other modeling software. Some of the key steps in PDB2PQR are described below.
Repairing missing heavy atoms Within PDB2PQR, the PDB file is examined to see if there are missing heavy (non-hydrogen) atoms. Missing heavy atoms can be rebuilt using standard amino acid topologies in conjunction with existing atomic coordinates to determine new positions for the missing heavy atoms. A debump option performs very limited minimization of sidechain χ angles to reduce steric clashes between rebuilt and existing atoms [38].
Optimizing titration states Amino acid titration states are important determinants of biomolecular (particularly enzymatic) function and can be used to assess functional activity and identify active sites. The APBS-PDB2PQR system contains several methods for this analysis.
• Empirical methods. PDB2PQR provides an empirical model (PROPKA [40]) that uses a heuristic method to compute pK a perturbations due to desolvation, hydrogen bonding, and charge-charge interactions. PROPKA is included with PDB2PQR. The empirical PROPKA method has surprising accuracy for fast evaluation of protein pK a values [41].
• Implicit solvent methods. PDB2PQR also contains two methods for using implicit solvent (Poisson-Boltzmann) models for predicting residue titration states. The first method uses Metropolis Monte Carlo to calculate titration curves and pK a values (PDB2PKA); however, sampling issues can be a major problem with Monte Carlo methods when searching over the O 2 N titration states of N titratable residues. The second method is a new polynomial-time algorithm for the optimization of discrete states in macromolecular systems [42]. The paper by Purvine et al. [42] describes the performance of the new graph cut method compared to existing approaches. While the new deterministic method offers advantages for large systems, the traditional stochastic approach is often sufficient for small to moderate systems. This method transforms interaction energies between titratable groups into a graphical flow network. The polynomial-time O (N 4 ) behavior makes it possible to rigorously evaluate titration states for much larger proteins than Monte Carlo methods.
Adding missing hydrogens The majority of PDB entries do not include hydrogen positions. Given a titration state assignment, PDB2PQR uses Monte Carlo sampling to position hydrogen atoms and optimize the global hydrogen-bonding network in the structure [43]. Newly added hydrogen atoms are checked for steric conflicts and optimized via the debumping procedure discussed above.
Assigning charge and radius parameters Given the titration state, atomic charges (for ρ) and radii (for and V i ) are assigned based on the chosen force field. PDB2PQR currently supports charge/radii force fields from AMBER99 [44], CHARMM22 [45], PARSE [46], PEOE PB [47], Swanson et al. [48], and Tan et al. [49]. Many of these force fields only provide parameters for amino acid biomolecules. The PEOE approach [47] provides algorithms to parameterize ligands. However, we welcome contributions for other biomolecular force fields, particularly for lipids and amino acids.

Solving the Poisson-Boltzmann and related solvation equations
The APBS software was designed from the ground up using modern design principles to ensure its ability to interface with other computational packages and evolve as methods and applications change over time. APBS input files contain several keywords that are generic with respect to the type of calculation being performed; these are described in Appendix A.1. The remainder of this section describes the specific solvers available for electrostatic calculations, also described in more detail on the APBS website [19].

Finite difference and finite element solvers
The original version of APBS was based on two key libraries from the Holst research group. FEtk is a general-purpose multi-level adaptive finite element library [50,51]. Adaptive finite element methods can resolve extremely fine features of a complex system (like biomolecules) while solving the associated equations over large problem domain. For example, FEtk has been used to solve electrostatic and diffusion equations over six orders of magnitude in length scale [52]. The finite difference PMG solver [51,53] trades speed and efficiency for the highaccuracy and high-detail solutions of the finite element FEtk library. However, many APBS users need only a relatively coarse-grained solution of φ for their visualization or simulation applications. Therefore, most APBS users employ the Holst group's finite difference gridbased PMG solver for biomolecular electrostatics calculations. Appendix A.2 describes some of the common configuration options for finite difference and finite element calculations in APBS.

Geometric flow
Several recent papers have described our work on a geometric flow formulation of Poissonbased implicit solvent models [54][55][56][57][58]. The components of this geometric model are described in previous publications [56,57]. The geometric flow approach couples the polar and nonpolar components of the implicit solvent model with two primary advantages over existing methods. First, this coupling eliminates the need for an ad hoc geometric definition for the solute-solvent boundary. In particular, the solute-solvent interface is optimized as part of the geometric flow calculation. Second, the optimization of this boundary ensures selfconsistent calculation of polar and nonpolar energetic contributions (using the same surface definitions, etc.), thereby reducing confusion and the likelihood of user error. Additional information about the geometric flow implementation in APBS is provided in Appendix A.3. This equation is solved in APBS using a finite difference method.

Boundary element methods
Boundary element methods offer the ability to focus numerical effort on a much smaller region of the problem domain: the interface between the molecule and the solvent. APBS now includes a treecode accelerated boundary integral PB solver (TABI-PB) developed by Geng and Krasny to solve a linearized version of the PB equation (Eq. 1) [59]. A discussion of the relative merits between finite difference/element and boundary integral PB methods are provided in Geng et al. and Li et al. [59,60]. In this method, two coupled integral equations defined on the solute-solvent boundary define a mathematical relationship between the electrostatic surface potential and its normal derivative with a set of integral kernels consisting of Coulomb and screened Coulomb potentials with their normal derivatives [61]. The boundary element method requires a surface triangulation, generated by a program such as MSMS [62] or NanoShaper [63], on which to discretize the integral equations. A Cartesian particle-cluster treecode is used to compute matrix-vector products and reduce the computational cost of this dense system from O(N 2 ) to O(N log N ) for N points on the discretized molecular surface [60,61]. A comparison of electrostatic potential with PDB ID 3app from the APBS multi-grid solution method and the new TABI-PB solver are shown in Figure A1. Additional information about the boundary element method and its implementation in APBS is provided in Appendix A.4.

Analytical and semi-analytical methods
Numerical solution methods tend to be computationally intensive, which has led to the adoption of analytical approaches for solvation calculations such as generalized Born [65] and the approaches developed by Head-Gordon and implemented in APBS. In particular, the Poisson-Boltzmann Analytical Method (PB-AM), was developed by Lotan and Head-Gordon in 2006 [29]. PB-AM produces a fully analytical solution to the linearized PB equation for multiple macromolecules, represented as coarse-grained low-dielectric spheres. This spherical domain enables the use of a multipole expansion to represent charge-charge interactions and higher order cavity polarization effects. The interactions can then be used to compute physical properties such as interaction energies, forces, and torques. An example of this approximation, and the resulting electrostatic potentials from PB-AM, are shown in Figure 3. The Poisson Boltzmann Semi-Analytical method (PB-SAM) is a modification of PB-AM that incorporates the use of boundary integrals into its formalism to represent a complex molecular domain as a collection of overlapping low dielectric spherical cavities [31]. PB-SAM produces a semi-analytical solution to the linearized PB equation for multiple macromolecules in a screened environment. This semi-analytical method provides a better representation of the molecular boundary when compared to PB-AM, while maintaining computational efficiency. An example of this approximation, and the resulting electrostatic potentials from PB-SAM, are shown in Figure 3.
Because it is fully analytical, PB-AM can be used for model validation as well as for repre- senting systems that are relatively spherical in nature, such as globular proteins and colloids. PB-SAM, on the other hand has a much more detailed representation of the molecular surface and can therefore be used for many systems that other APBS (numerical) methods are currently used for. Through APBS, both programs can be used to compute the electrostatic potential at any point in space, report energies, forces, and torques of a system of macromolecules, and simulate a system using a BD scheme [66]. Additional details about these methods and their use in APBS are presented in Appendix A.5. In addition to these advantages, the strengths and weaknesses of PB-AM and PB-SAM over existing methods are discussion in Yap et al. [31] and Lotan et al. [29].

Visualization
One of the primary uses of the APBS tools is to generate electrostatic potentials for use in biomolecular visualization software. These packages offer both the ability to visualize APBS results as well as a graphical interface for setting up the calculation. Several of these software packages are thick clients that run from users' computers, including PyMOL [67,68], VMD [64,69], PMV [70,71], and Chimera [72,73]. We have also worked with the developers of Jmol [74,75] and 3Dmol.js [76,77] to provide web-based setup and visualization of APBS-PDB2PQR calculations and related workflows. APBS integration with Jmol has been described previously [78]. 3Dmol.js is a molecular viewer that offers the performance of a desktop application and convenience of a web-based viewer which broadens accessibility for all users. As part of the integration with 3Dmol.js, we implemented additional enhancements, including extending our output file formats and creating a customized user interface. Data from the APBS output file is used to generate surfaces, apply color schemes, and display different molecular styles such as cartoons and spheres. An example of the 3Dmol.js interface is shown in Figure 4. Examples of 3Dmol.js visualization options are shown in Figure 5.

Other applications
Besides visualization and the processes described in Section 2, there are a number of other applications where APBS can be used. For example, during the past four years, the APBS-PDB2PQR software has been used in the post-simulation energetic analyses of molecular dynamics trajectories [79], understanding protein-nanoparticle interactions [80,81], understanding nucleic acid-ion interactions [82], biomolecular docking [83] and ligand binding [84], developing new coarse-grained protein models [85], setting up membrane protein simulations [86], etc. APBS also plays a key role in PIPSA for protein surface electrostatics analysis [87] as well as BrownDye [88] and SDA [89] for simulation of protein-protein association kinetics through Brownian dynamics. As discussed above with PB-SAM, another application area for implicit solvent methods is in the evaluation of biomolecular kinetics where implicit solvent models are generally used to provide solvation forces (or energies) for performing (or analyzing) discrete molecular or continuum diffusion simulations with APBS in both of these areas [79,[89][90][91][92][93][94].

The future of APBS
To help with modularity and to facilitate extensibility, APBS was based on an object-oriented programming paradigm with strict ANSI-C compliance. This "Clean OO C" [95] has enabled the long-term sustainability of APBS by combining an object-oriented design framework with the portability and speed of C. This object-oriented design framework has made it relatively straightforward to extend APBS functionality and incorporate new features.
The Clean OO C design has served APBS well for the past 17 years. This design still forms the basis for many modules of APBS and PDB2PQR. Additionally, we have begun to explore a framework for integrating components that exploits the common workflows used by APBS-PDB2PQR users and maps naturally to cloud-based resources. Our vision for APBS is to build the infrastructure that can enable our users to implement their own models and methods so that they can run on a common system. Our goal is to have a well-designed, well-tested and well-documented industry-grade framework that will integrate the APBS-PDB2PQR capabilities and make straightforward the incorporation of new features and new workflows.

A Appendix
These appendices provide basic information about configuring APBS electrostatics calculations. Rather than duplicating the user manual on the APBS website [19], we have only focused on input file keywords that directly relate to the setup and execution of solvation calculations.

A.1 General keywords for implicit solvent calculations
This section contains information about the general keywords used to configure implicit solvent calculations in the ELEC section of APBS input files. These keywords are used in all of the solver types described in this paper: • ion charge charge conc conc radius radius : This line can appear multiple times to specify the ionic species in the calculation. This specifies the following terms in Eq. 1: charge gives q i in units of electrons, conc gives c i in units of M, and radius specifies the ionic radius (inÅ) used to calculate V i .
• mol id : Specify the ID of the molecule on which the calculations are to be performed. This ID is determined by READ statements which specify the molecules to import.
• pdie and sdie: These keywords specify the dielectric coefficient values for the biomolecular interior (pdie) and bulk solvent (sdie). A typical value for sdie is 78.54; values for pdie are much more variable and often range from 2 to 40.
• temp temp : The temperature (in K) for the calculation. A typical value is 298 K.
Additionally, READ statements in APBS input files are used to load molecule information, parameter sets, and finite element meshes. More detailed information about these and other commands can be found on the APBS website [19].

A.2 Finite element and finite difference calculations in APBS
The finite difference and finite element methods used by APBS have been described extensively in previous publications [18,[26][27][28]53]; this section focuses on the configuration and use of these methods in APBS.

A.2.1 Finite difference calculation configuration
APBS users have several ways to invoke the PMG finite difference solver [53] capabilities of APBS through keywords in the APBS input file ELEC block. The different types of finite difference calculations include: • mg-manual: This specifies a manually configured multigrid finite difference calculation in APBS.
• mg-auto: This specifies an automatically configured multigrid finite difference calculation in APBS, using focusing [96] to increase the resolution of the calculation in areas of interest.
• mg-para: This specifies a parallel version of the multigrid finite difference calculating, using parallel focusing to increase the resolution of the calculation in areas of interest [18].
These different types of calculations have several keywords described in detail on the APBS website [19]. Some of the most important settings are described below.
Boundary conditions Boundary conditions must be imposed on the exterior of the finite difference calculation domain. In general, biomolecular electrostatics calculations use Dirichlet boundary conditions, setting the value of the potential to an asymptotically correct approximation of the true solution. There are several forms of these boundary conditions available in APBS with approximately equal accuracy given a sufficiently large calculation domain (see below). The only boundary condition which is not recommended for typical calculations is the zero-potential Dirichlet condition.
Grid dimensions, center, and spacing Finite difference calculations in APBS are performed in rectangular domains. The key aspects of this domain include its length L i and number of grid points n i in each direction i. These parameters are related to the spacing h i of the finite difference grid by h i = L i /(n i − 1). Grid spacings below 0.5Å are recommended for most calculations. The number of grid points is specified by the dime keyword. For mg-manual calculations, the grid spacing can be specified by either the grid or the glen keywords, which specify the grid spacing or length, respectively. For mg-auto or mg-para calculations, the grid spacing is determined by the cglen and fglen keywords, which specify the lengths of the coarse and fine grids for the focusing calculations. Grid lengths should extend sufficient distance away from the biomolecule so that the chosen boundary condition is accurate. In general, setting the length of the coarsest grid to approximately 1.5 times the size of the biomolecule gives reasonable results. However, as a best practice, it is important to ensure that the calculated quantities of interest do not change significantly with grid spacing or grid length. The center of the finite difference grid can be specified by the gcent command for mg-manual calculations or by the cgcent and fgcent keywords for the coarse and fine grids (respectively) in mg-auto or mg-para focusing calculations. These keywords can be used to specify absolute grid centers (in Cartesian coordinates) or relative centers based on molecule location. Because of errors associated with charge discretization, it is generally a good idea to keep the positions of molecules on a grid fixed for all calculations. For example, a binding calculation for rigid molecules should keep all molecules in the same positions on the grid.
Charge discretization As mentioned above, atomic charge distributions are often modeled as a collection of Dirac delta functions or other basis functions with extremely small spatial support. However, finite difference calculations performed on grids with finite spacing, requiring charges to be mapped across several grid points. This mapping creates significant dependence of the electrostatic potential on the grid spacing which is why it is always important to use the same grid setup for all parts of an electrostatic calculation. The chgm keyword controls the interpolation scheme used for charge distributions and includes the following types of discretization schemes: • spl0: Traditional trilinear interpolation (linear splines). The charge is mapped onto the nearest-neighbor grid points. Resulting potentials are very sensitive to grid spacing, length, and position.
• spl2: Cubic B-spline discretization as described by Im et al. [97]. The charge is mapped onto the nearest-and next-nearest-neighbor grid points. Resulting potentials are somewhat less sensitive (than spl0) to grid spacing, length, and position; however, this discretization can often require smaller grid spacings for accurate representation of charge positions.
• spl4: Quintic B-spline discretization as described by Schnieders et al. [98]. Similar to spl2, except the charge/multipole is additionally mapped to include next-next-nearest neighbors (125 grid points receive charge density). This discretization results in less sensitivity to grid spacing and position but nearly always requires smaller grid spacings for accurate representation of charge positions.
Surface definition APBS provides several difference surface definitions through the srfm keyword: • mol: The dielectric coefficient is defined based on a molecular surface definition. The problem domain is divided into two spaces. The "free volume" space is defined by the union of solvent-sized spheres (size determined by the srad keyword) which do not overlap with biomolecular atoms. This free volume is assigned bulk solvent dielectric values. The complement of this space is assigned biomolecular dielectric values. With a non-zero solvent radius (srad), this choice of coefficient corresponds to the traditional definition used for PB calculations. When the solvent radius is set to zero, this corresponds to a van der Waals surface definition. The ion-accessibility coefficient is defined by an "inflated" van der Waals model. Specifically, the radius of each biomolecular atom is increased by the radius of the ion species (as specified with the ion keyword). The problem domain is then divided into two spaces. The space inside the union of these inflated atomic spheres is assigned an ion-accessibility value of 0; the complement space is assigned bulk ion accessibility values.
• smol: The dielectric and ion-accessibility coefficients are defined as for mol (see above). However, they are then "smoothed" by a 9-point harmonic averaging to somewhat reduce sensitivity to the grid setup [99].
• spl2: The dielectric and ion-accessibility coefficients are defined by a cubic-spline surface [97]. The width of the dielectric interface is controlled by the swin parameter. These spline-based surface definitions are very stable with respect to grid parameters and therefore ideal for calculating forces. However, they require substantial reparameterization of the force field [100].
• spl4: The dielectric and ion-accessibility coefficients are defined by a 7 th -order polynomial. This surface definition has characteristics similar to spl2, but provides higher order continuity necessary for stable force calculations with atomic multipole force fields (up to quadrupole) [98].

A.2.2 Finite elements calculation configuration
Users invoke the FEtk finite element solver [50] in APBS by including the fe-manual keyword in the ELEC section of the input file. Many aspects of the finite element configuration closely follow the finite difference options described above. This section only describes the configuration options which are unique to finite element calculations in APBS. Finite element calculations begin with an initial mesh. This mesh can be imported from an external file via the usemesh keyword or generated by APBS. APBS generates the initial mesh from a very coarse 8-tetrahedron mesh of size domainLength which is then refined uniformly or selectively at the molecular surface and charge locations, based on the value of the akeyPRE keyword. This initial refinement occurs until the mesh has targetNum vertices or has reached targetRes edge length (inÅ).
As described previously [26,27], APBS uses FEtk in a solve-estimate-refine iteration which involves the following steps: 1. Solve the problem with the current finite element mesh.
2. Estimate the error in the solution as a function of position on the mesh. The method of error estimation is determined by the ekey keyword which can have the values: • simp: Per-simplex error threshold; simplices with error above this limit are flagged for refinement.
• global: Global (whole domain) error limit; flag enough simplices for refinement to reduce the global error below this limit.
• frac: The specified fraction of the simplices with the highest amount of error are flagged for refinement.
3. Adaptively refine the mesh to reduce the error using the error metric described by ekey.
This iteration is repeated until a target error level etol is reached or a maximum number of solve-estimate refine iterations (maxsolve) or vertices (maxvert) is reached.

A.3 Geometric flow calculations in APBS
This section contains additional information about the geometric flow equation implementation in APBS introduced in Section 3.2. The geometric flow methods used by APBS have been described extensively in previous publications [54][55][56][57]101]; this section focuses on the configuration and use of these methods in APBS.

A.3.1 Geometric flow calculation configuration
Users invoke the geometric flow solver in APBS by including the geoflow-auto keyword in the ELEC section of the input file. Because the geometric flow solver is based on finite difference solvers, many of the keywords for this section are similar to those described in Section A.2.1. Three additional parameters are needed for geometric flow calculations to specify how nonpolar solvation is linked to the polar implicit solvent models: • gamma tension : Specify the surface tension of the solvent in units of kJ mol −1Å−2 . Based on Daily et al. [57], a recommended value for small molecules is 0.431 kJ mol −1 A −2 .
• press pressure : Specify the internal pressure of the solvent in units of kJ mol −1Å−3 . Based on Daily et al. [57], a recommended value for small molecules is 0.104 kJ mol −1 A −3 .
• vdwdisp bool : Indicate whether van der Waals interactions should be included in the geometric flow calculation through bool (1 = include, 0 = exclude). If these interactions are included, then a force field with van der Waals terms must be included through a READ statement in the APBS input file.

A.4 Boundary element method implementation
This appendix provides additional information about the boundary element method introduced in Section 3.3.

A.4.1 Boundary element method background
This section provides additional background on the TABI-PB boundary element solver [59], introduced in Section 3.3. As described earlier, this method involves solving two coupled integral equations defined on the solute-solvent boundary which define a mathematical relationship between the electrostatic surface potential and its normal derivative with a set of integral kernels consisting of Coulomb and screened Coulomb potentials with their normal derivatives. The boundary element method requires a surface triangulation, generated by a program such as MSMS [62] or NanoShaper [63], on which to discretize the integral equations. Figure A1 shows different types of surface discretizations as well as example electrostatic potential output. The coupled second kind integral equations employed by TABI-PB for calculating the surface potential φ and its normal derivative [61] are: where ε = ε m /ε s , the ratio of the dielectric constant in the solute region and the dielectric constant in the solvent region. The integral kernels K 1 , K 2 , K 3 , K 4 are defined in Eq. 3.
where G 0 and G κ are the Coulomb and screened Coulomb potentials defined as: The normal derivatives of the potential kernels G are: ν n (y) ∂ xn G (x, y) , for the three spatial components n of the normal direction. Additionally, the source terms S 1 and S 2 in Eq. 2 are: where N c is the number of atoms in the solute molecule, and q k is the charge of the kth atom. Note that S 1 is a linear superposition of the point charge electrostatic potentials, and S 2 is a linear superposition of the normal derivatives of the potentials. Given a surface triangularization with N elements -where x i and A i are the centroid and area, respectively, of the ith triangle -the integral equations are discretized as: The omission of the j = i term in the summation avoids the singularity of the kernels at that point. Note that the right hand sides of these equations consist of sums of products of kernels and the surface potential or its normal derivative. These are the analogues to the N -body potential in the treecode, with the surface potential or its normal derivatives playing the role of the charges. In the discretized form, the total electrostatic energy of solvation is given by Eq. 8: where q k is the charge on the kth atom of the solute molecule, and x k is its position. The matrix-vector products involve evaluation of the integral kernels over the surface elements. These evaluations effectively take the form of an N -body potential: a sum over a set of N positions of products between a kernel and a "charge" at each position. In this case, the locations of the N particles are the centroids of the surface triangularization elements. A Cartesian particle-cluster treecode is used to compute matrix-vector products and reduce the computational cost of this dense system from O(N 2 ) to O(N log N ) for N points on the discretized molecular surface [60]. In particular, to rapidly evaluate the N -body potential at the N particle locations, the treecode subdivides the particles into a tree-like hierarchical structure of clusters. At each location, the potential contribution from nearby particles is computed by direct sum, while for well-separated particle cluster interactions, a Taylor approximation about the center of the cluster is used to evaluate the contribution. The Taylor coefficients are calculated through recurrence relations. The resulting linear system is then solved with GMRES iteration [102].
Because the integral equations are defined on the molecular boundary, the singular charges are handled analytically and do not introduce the same issues present in grid-based schemes. The integral equations also rigorously enforce the interface conditions on the surface, and the boundary condition at infinity is exactly satisfied. Thus, the boundary integral formulation can potentially be superior to other methods for investigating electrostatic potential on the boundary.

A.4.2 Boundary element calculation configuration
APBS users can invoke TABI-PB with the bem-manual flag in the ELEC section of the input file. Major options include: • tree order order : An integer indicating the order of the Taylor expansion for determining treecode coefficients. Higher values of order will result in a more accurate -but more expensive -calculations. A typical choice for this parameter is 3.
• tree n0 number : The maximum number of particles allowable in a leaf of the treecode (clusters in the last level of the tree). A typical choice for this parameter is 500.
• mac criterion : Multipole acceptance criterion specifies the distance ratio at which the Taylor expansion is used. In general, a higher value of criterion will result in a more accurate but more expensive computation; while a lower value causes more direct summations and forces the particle-cluster interaction to descend to a finer cluster level. A typical choice for this parameter is 0.8. • outdata flag : Type of output data file generated; 0 = APBS OpenDX format [103] and 1 = ParaView format [104].
Additional information about parameter settings is provided via the APBS website [19]. TABI-PB produces output including the potential and normal derivative of potential for every element and vertex of the triangularization, as well as the electrostatic solvation energy. Examples of electrostatic surface potential on the protein 1a63 are shown in Figure 2 by using MSMS and NanoShaper.

A.5 Analytical and semi-analytical method implementations
This appendix provides additional information about the analytical and semi-analytic methods [29,30] introduced in Section 3.4.

A.5.1 Analytical method (PB-AM) background
The solution to the PB-AM model is represented as a system of linear equations: where A represents a vector of the effective multipole expansion of the charge distributions of each molecule, E is a vector of the fixed charge distribution of all molecules, Γ is a dielectric boundary-crossing operator, ∆ is a cavity polarization operator, and T is an operator that transforms the multipole expansion from the global (lab) coordinates to a local coordinate frame. The unknown A determined using the Gauss-Seidel iterative method and can then be used to compute physical properties such as interaction energies, forces, and torques. The interaction energy for molecule i, (Ω (i) ) is given in Eq. 10.
where s is the dielectric constant of the solvent and M, N denotes the inner product. When energy is computed, forces follow as: By definition, the torque on a charge in the molecule is the cross product of its position relative to the center of mass of the molecule with the force it experiences. The total torque on the molecule is a linear combination of the torque on all charges of the molecule, as illustrated in Eq. 12.
where α H , α = x, y, z, is a coefficient vector for each of the charges in the molecule, M i is the number of charges in molecule i, q (i) J is the magnitude of the j th charge, and p is its position in spherical coordinates. For more details on the PB-AM derivation, see Lotan and Head-Gordon [29].

A.5.2 Semi-analytical method (PM-SAM) background
The derivation details of PB-SAM have been reported previously [31,32], with the main points being summarized in this section. The electrostatic potential (φ r ) of the system at any point r is governed by the linearized form of the PB equation: where κ is the inverse Debye length. Eq. 13 is a linearization of Eq. 1 for βq i φ 1. For the case of spherical cavities, we can solve Eq. 13 by dividing the system into inner sphere and outer sphere regions, and enforcing a set of boundary conditions that stipulate the continuity of the electrostatic potential and the electrostatic field at the surface of each sphere. The electrostatic potential outside molecule (I) is described by: where h(r) is an effective surface charge that can be transformed into the unknown multipole expansion H (I,k) with inside molecule I and sphere k. In a similar manner, the interior potential is given as where N

(I)
C is the number of charges in molecule I, q α is the magnitude of the α-th charge, r is its position in spherical coordinates, and f (r) is a reactive surface charge that can be transformed into the unknown multipole expansion F (I,k) . The reactive multipole and the effective multipole, H (I,k) , are given as: n,m (θ , φ )dr (16) n,m (θ , φ )dr (17) where Y n,m is the spherical harmonics, Y n,m is the complete conjugate, and a (I,k) is the radius of sphere k of molecule I. These multipole expansions can be iteratively solved using: where W F (I,k) and W H (I,k) are scaled multipoles computed from fixed charges and polarization charges from other spheres. I (I,k) E,n,m is a matrix of the surface integrals over the exposed surface: I n,m (θ , φ ) sin θ dθ dφ (20) Using the above formalism, physical properties of the system, such as interaction energy, forces and torques can also be computed. The interaction energy of each molecule, (Ω (i) ), is the product of the molecule's total charge distribution (from fixed and polarization charges) with the potential due to external sources. This is computed as the inner product between the molecule's multipole expansion, (H (I,k) ), and the multipole expansions of the other molecules in the system, (LHN (I,k) ) as follows: As in the analytical PB-AM method, the torque on a charge in the molecule is the cross product of its position relative to the center of mass of the molecule with the force it experiences. For a charge at position P about the center of mass c (I) for molecule I, the torque is given by the cross product of its position r (I,k) P with respect to the center of mass and the force on that charge f P . We can re-express r (I,k) P as the sum of vectors from the center of molecule I to the center of sphere k (c (I,k) ) and from the center of sphere k to point P (r (I,k) P ). The total torque on molecule I is then given by Eq. 23.
where α = x, y, z. For the derivation of the PB-SAM solver please see the previous publications [31,32] .

A.5.3 PB-AM and PB-SAM configuration in APBS
PB-AM and PB-SAM have been fully integrated into APBS, and is invoked using the keyword pbam-auto or pbsam-auto in the ELEC section of an APBS input file. Major options include: • runname name : Desired name to be used for outputs of each run.
• pbc length : Size of the periodic simulation/calculation domain.
• ntraj number : Number of Brownian Dynamics simulations to run.
• term type value mol : Allows the user to indicate conditions for the termination of each BD trajectory. The following values of type are allowed: time time : A limit on the total simulation time.
x or y or y or z or r and >= or <= : Represents the approach of two molecules to a certain distance r or certain region of space given by x or y or y. The operators >= and <= represent the corresponding inequalities.
The parameter mol is the molecular index that this condition applies. mol should be 0 for time and for a termination condition of x, y or z, the molecule index that this termination condition applies to.
• xyz idx fpath : Molecule index idx and file path fpath for the molecule starting configurations. A starting configuration is needed for each molecule and each trajectory. Therefore, jf there are m molecules and ntraj n trajectories, then the input file must contain m × n xyz entries.
• tolsp val : Modify the coarseness of the molecular description. val is the distance (inÅ) beyond the solvent-excluded surface that the coarse-grained representation extends. Increasing values of val leads to fewer coarse-grained spheres, faster calculation times, but less accurate solutions. Typical values for val are between 1 and 5Å.
The commands (keywords) not included in this list are used to specify system conditions, such as temperature and salt concentration. These parameters are similar to those found in the ELEC section of a usual APBS run and are documented on the Contributions portion of the APBS website [19]. Additional information about parameter settings is provided via the APBS website [19]. Examples of the electrostatic potentials produced from PB-AM and PB-SAM are shown in Figure 3.