A Synchronized Stereo and Plenoptic Visual Odometry Dataset

We present a new dataset to evaluate monocular, stereo, and plenoptic camera based visual odometry algorithms. The dataset comprises a set of synchronized image sequences recorded by a micro lens array (MLA) based plenoptic camera and a stereo camera system. For this, the stereo cameras and the plenoptic camera were assembled on a common hand-held platform. All sequences are recorded in a very large loop, where beginning and end show the same scene. Therefore, the tracking accuracy of a visual odometry algorithm can be measured from the drift between beginning and end of the sequence. For both, the plenoptic camera and the stereo system, we supply full intrinsic camera models, as well as vignetting data. The dataset consists of 11 sequences which were recorded in challenging indoor and outdoor scenarios. We present, by way of example, the results achieved by state-of-the-art algorithms.

While there are various datasets available for traditional types of cameras -monocular [18,19,20], stereo [21,22,23], and RGB-D [24,25] -there are no public datasets available for plenoptic camera based VO. The few existing algorithms were only tested on very simple and short sequences.
In this paper, we present a versatile and challenging dataset which supplies synchronizes light field and stereo images. These sequences were recorded based on a hand-held platform on which a micro lens array (MLA) based plenoptic camera and a stereo camera system is mounted.
Of course, there exist larger dataset for monocular and stereo cameras. However, the goal of this dataset is to evaluate the versatility of plenoptic VO algorithms and rank these algorithms with respect to methods based on traditional cameras (i.e. monocular and stereo).
The entire synchronized stereo and plenoptic VO dataset is available at: https://www.hs-karlsruhe.de/odometry-data/

Outline
In Section 3 we present the platform which we used to record the image sequences. Afterwards, in Section 4, we describe the entire structure of the dataset. This includes the calibration of all cameras as well as the ground truth data and the suggested evaluation metrics. Section 5 shows the results of existing algorithms, which were obtained for the proposed dataset. Furthermore, some 3D reconstructions calculated by [17] are shown to get an impression of the recorded sequences. Section 6 mentions some limitations, which should be taken into account, when rating the results of different algorithms against each other.

Data Acquisition Setup
To be able to compare plenoptic with stereo or monocular algorithms, a hand-held platform was developed. On this platform a MLA based plenoptic camera and a stereo camera system are assembled. This platform is shown in Figure 1. The stereo camera system is based on two monochrome, global shutter, industrial grade cameras by IDS Imaging Development Systems GmbH (model: UI-3241LE-M-GL). On both cameras, a lens from Lensation GmbH (model: BM4018S118) with 4 mm focal length is mounted. Furthermore, the stereo system has a baseline distance of 100 mm.
The utilized plenoptic camera is a R5 by Raytrix GmbH. This camera is based on a xiQ sensor from Ximea GmbH (model: MG042CG-CM-TG). The xiQ is also a global shutter, industrial grade sensor. To achieve a suitable trade-off between a wide field of view (FOV) and a high angular resolution of the captured light field, a main lens with focal length f L = 16 mm from Kowa (model: LM16HC) was mounted on the camera. To receive synchronized image sequences from all three cameras, one camera of the stereo systems runs in master mode and generates a signal which, in turn, triggers the other two cameras (second camera of the stereo system and plenoptic camera) running in slave mode. To record the data, all three cameras are connected to a single laptop. Using this platform, we are able to record synchronized image sequences for all three cameras running at the maximum image resolution and 8 bit quantization with frame rates of more than 30 fps.
Because of the small FOV, the images of the plenoptic camera in particular are affected by motion blur. To keep the motion blur in an acceptable range, for all cameras the exposure time is upper bounded at 8 ms. Below this boundary, the automatic exposure adjustment of the respective camera controls the exposure time 3 .

The Dataset
Based on the platform presented in Section 3, we recorded a synchronized dataset for the quantitative comparison of plenoptic and stereo VO systems. The inspiration for this dataset was taken from the Benchmark for monocular VO presented by Engel et al. [19].
Especially for long, large scale trajectories, it is impossible to obtain reference measurements which are accurate enough to serve as ground truth. Hence, guided by the idea of [19], we perform all sequences in the dataset as a single very large loop, for which beginning and end of the sequence capture the same scene.
For all recorded sequences we supply geometric camera parameters and vignetting data for both the stereo camera system and the plenoptic camera. Furthermore, for each sequence, ground truth data is obtained from the loop closure between beginning and end.

Geometric Camera Calibration
For all cameras, we perform geometric calibration based on a 3D target, which is presented below in Section 4.1.3. Aside from the camera models for the plenoptic camera and the stereo system, we supply the entire sets of images as well as the marker positions detected in the images, which were used for calibration. This way, one can later test new camera models and calibration approaches on the basis of this data.

Plenoptic Camera Model
The projection model which is used for the plenoptic camera in this dataset, is the one proposed in [16] and visualized in Figure 2. In this model, the main lens of the plenoptic camera is a thin lens, while the micro lenses in the MLA are pinholes. In [16] it was shown, that this model forms, in fact, the equivalent to a virtual camera array, as it is shown in Figure 2(b), where each micro image represents the image of a small pinhole camera with a very narrow field of view.
Using this model, one obtains the coordinates of a point x p in a virtual camera from the camera coordinates x C of a 3D object point as follows: (1) Here, p M L defines the center of the virtual camera, or projected micro lens. The center p M L is calculated from the center of the real micro lens c M L as follows: The parameter z C0 defines the distance from the real main lens of the plenoptic camera to the virtual camera array (see Fig. 2(b)).
Furthermore, a point x M L in a real micro image can be calculated from the corresponding point x p in the respective virtual camera as follows: The micro image point x M L , given in eq. (4), is an image point relative to its micro lens center c M L . Hence, corresponding raw image coordinates x R , which are unique for each single point in the entire raw image recorded by the plenoptic camera, are defined: The common way to obtain the centers of the micro lenses in the MLA is to estimate them based on a recorded white image [26]. However, these centers, in fact, do not represent the micro lens centers c M L , but instead the corresponding micro image centers c I . As it is shown in [17], the micro lens center c M L and the micro image centers c I have the following relationship: To correct for lens distortions, a distortion model is applied to the raw image coordinates x R . Hence, the following connection between the distorted coordinates x Rd , which, in fact, are the coordinates of the image recorded by the camera, and the undistorted coordinates x R is defined: Here, the distortion terms, ∆x dist and ∆y dist , consist of a radial symmetric as well as a tangential distortion component and are defined as follows: For the plenoptic camera, we found two radial symmetric parameters (A 0 , and A 1 ) to be sufficient to model the distortion. The micro image centers are detected also on distorted raw image coordinates c Id and therefore have to be corrected by the same distortion model. Hence, the corrected micro image centers c I will not be arranged on a regular hexagonal grid anymore, but will slightly deviate from this grid.
So far, the coordinates x Rd were defined in metric dimension and relative to the optical axis. Hence, they still have to be transformed into pixel coordinates x Rd as follows: Here, s defines the size of a pixel and c = [c x , c y ] T is the so-called principal point. As already mentioned, the micro image centers c Id can be estimated from a recorded white image. Furthermore, the pixel size s can be taken directly from the sensor specifications. All other parameter, have to be estimated in a geometric calibration. These parameters are: • main lens focal length: f L • distance between main lens and MLA: b L0 • distance between MLA and sensor:

Stereo Camera Model
For the two monocular cameras in the stereo setup we define the pinhole camera model as given in eq. (12). In contrast to the plenoptic camera model, we define different focal lengths (f x , f y ) in xand y-direction. Thereby, we are able to consider rectangular, instead of squared, sensor pixels. Lens distortion is applied on normalized image coordinates as given in eq. (13).
Due to the much larger FOV, rather than two, we have to consider three parameters for radial symmetric distortion. Furthermore, the effect of tangential distortion is negligible.
The variable r = x 2 + y 2 defines the distance to the principal point on the sensor in normalized image coordinates. In addition to the intrinsic parameters, the orientation of the slave (right) camera with respect to the master (left) camera is defined by a rigid body transformation G(ξ M S ) ∈ SE(3), which is represented by the respective tangent space element ξ M S ∈ se(3).

Calibration Approach
For both systems, the plenoptic camera and the stereo camera, the model parameters are estimated from a set of images in a full bundle adjustment. For this purpose we use the 3D calibration target shown in Figure 3(a). In the bundle adjustment, all parameters of the camera models, the extrinsic orientations of the single images as well as the 3D coordinates of the calibration markers are estimated. The micro images of the plenoptic camera generally do not cover a complete marker point. Therefore, the marker points cannot be detected reliably in the micro images. Hence, we calculate from each raw image, recorded by the plenoptic camera, the corresponding totally focused image. Afterwards, the marker points are detected in the totally focused image and are projected back to the micro images in the respective raw image. This procedure was described already in [16].

Vignetting Correction
Especially for direct VO approaches, vignetting has a negative effect on the performance. The model parameters are estimated based on photometric measurements and therefore the measured intensity of a point must be independent of its location on the sensor. Indirect methods are more robust to vignetting, as extracted feature points generally rely on corners in the images, which are invariant to absolute intensity changes.
While in a monocular camera vignetting generally results in a continuously increasing attenuation of pixel intensities from the image center towards the boundaries, in a plenoptic camera further vignetting effects are present for each single micro lens.
Mathematically we can describe the vignetting as follows: The function I(x) is the observed intensity value measured by the sensor, while B(x) is the irradiance image which represents the scene in a photometrically correct way. The vignetting V (x) defines a pixel-wise attenuation, with V (x) ∈ [0, 1]. We use the notations τ for the exposure time and I for the sensor noise. In this simplified model, we considered the image sensor to have a linear transfer characteristic which is, in fact, not the case for a real sensor. While the nonparametric vignetting compensation based on white images, recorded with a white balance filter, is commonly applied to plenoptic cameras, we correct the vignetting of the two cameras in the stereo system in the same way. For each camera we recorded a set of 10 white images and calculated an average attenuation image from this set. Figure 3(b) shows the used white balance filter. For the two cameras of the stereo system, we additionally filtered the attenuation images using a Gaussian kernel. This cannot be done for the attenuation image of the plenoptic camera, as in this case, the MLA produces quite high frequent components in the attenuation image which must be preserved (see Figure 4).
In [19] a different method is described, where the attenuation map is calculated based on a sequence of images capturing a white wall. However, the method [19] is quite time  consuming and, in our experience, error prone 4 . For the method [19], one has to capture a sequence of hundreds of images and then run the estimation for up to one hour. The attenuation map based on the white balance filter, by contrast, is obtained in just a few seconds. Furthermore, we want to apply comparable calibrations to both systems; the plenoptic camera as well as the stereo cameras. Figure 4 shows the vignetting for the plenoptic camera. As one can see, vignetting is visible in each individual micro image. Furthermore, outer image regions are attenuated stronger than the image center. This is due to the influence of the main lens. In addition, there are some small irregularities visible in the map resulting from defect micro lenses and dirt on the MLA. Figure 5 shows the attenuation maps estimated for the left camera of the stereo system. Figure 5(a) shows the result using the white balance filter, while Figure 5(b) shows the one obtained from the method of Engel et al. [19]. From Figure 5(c) one can clearly see that for the white balance filter the attenuation is sightly stronger at the sensor boundaries than for the method in [19]. However, the deviation is quite small and furthermore, it is difficult to evaluate which map describes the vignetting of the camera in a more accurate way.

Ground Truth and Evaluation Metric
It is almost impossible to obtain ground truth trajectories for long and large-scale sequences recorded by hand-held cameras. We decided to obtain ground truth data for our dataset in a similar way as suggested by Engel et al. [19], where the accuracy of a VO algorithm is evaluated based on a single, larger loop closure. Each trajectory in the dataset starts with a winding sequence while capturing a nearby object. This starting sequence is followed by the actual trajectory which finally leads back to the starting point in a large loop, followed by a short, winding, finishing sequence.
Using the winding sequence at the beginning and at the end, we are able to register both segments to each other using a standard structure from motion (SfM) approach. These registered segments then can be used as ground truth information. In contrast to monocular datasets, we also want ground truth data for the absolute scale of the trajectory. The absolute scale of the trajectory is obtained from stereo images. Since the stereo cameras have a much larger stereo baseline than the micro images in the plenoptic camera, the observed scale is accurate enough to serve as reference for the plenoptic sequences.
While one can basically use any SfM algorithm to register the start and end segment to each other, we use a modified version of ORB-SLAM2 (stereo) [9]. Instead of selecting keyframes, we build up a frame-wise pose graph which is optimized in a global bundle adjustment. Figure 6 shows, by way of example, the registered segments for the beginning and the end of the sequence in blue and red respectively, overlaid with the point cloud calculated by [17].
Using the registered ground truth data, based on each recorded trajectory two similarity transformations T gt s and T gt e (∈ Sim (3)) with respect to the start and the end segment of the sequence can be calculated, as given in eqs. (17) and (18).
T gt e := arg min The vectors p i ∈ R 3 are the estimated points of the trajectory while p gt i ∈ R 3 are the respective points of the ground truth. S and E define the sets of indices of the start end and segment respectively. One may notice that in eqs. (17) and (18) actually the homogeneous representations of 3D points p i and p gt i have to be used. Using the similarity transformations T gt s and T gt e we define evaluation metrics similar to [19]. From the two transformations the accumulated drift T drift ∈ Sim(3) from the start to the end of the trajectory can be calculated as follows: From T drift we can directly extract the scale drift e s , the rotational drift e r , and the translation drift e t := t . The rotational drift e r is defined by the rotation angle around the Euler axis w → w corresponding to the rotation matrix R ∈ SO (3): The mapping w → w defines the mapping of the vector w ∈ R 3 to the skew-symmetric matrix w ∈ so (3): For an easier interpretation of the scale drift e s := max{e s , e −1 s } is defined. The drift metrics e s , e r , and e t define more or less independent quality measures. Looking at just one of these values offers only a very limited insight into the overall quality of the estimated trajectory. Furthermore, e t , as it is defined here, is proportional to the absolute scale of the estimated trajectory and therefore is not meaningful at all without considering the absolute scale. In [19] the alignment error e align , as a more meaningful combined metric, is defined: The parameter N is the number of points (frames) in the complete trajectory. In comparison to e t , e align is always scaled with respect to the ground truth and implicitly incorporates all drifts e s , e r , e t in a single number. All previously defined metrics consider only the relative drift from the beginning to the end of the trajectory, but not the error of the absolute scale. We define the absolute scale difference d s of the front and end segment as another metric: (a) sample images from the sequence of the left camera in the stereo camera system (b) sample images from the sequence of the plenoptic camera Similar to the scale drift e s , we define d s := max{d s , d −1 s }. Of course, d s must be considered only for plenoptic and stereo algorithms and not for monocular approaches. Furthermore, d s has significance only in combination with the scale drift e s : To obtain reliable ground truth data, all sequences start and end in a scene showing objects in a distance of several meters, which are easy to track. However, for these nearby objects it is easier to estimate the correct scale. Hence, to consider only the scale drift e s or the absolute scale d s might be misleading. In combination with the alignment error e align , these values become more meaningful since the alignment error would reflect large scale drifts along the trajectory.
Following the scheme described above, we recorded a set of 11 sequences in versatile environments. The recorded scenes range from large scales to small scales, from man-made environments to environments with abundant vegetation. The sequences capture moving objects like pedestrians, bikes or cars. The sequences also cover difficult and changing lighting conditions due to shadows, moving clouds, and automatic exposure adjustment. The path lengths of the performed trajectories range from 25 m to 274 m. Figure 7 shows a set of sample images extracted from a single sequence. The corresponding trajectory is shown in Figure 8. Figure 6 visualizes the registered start and end segments which are used to calculate the metrics. As one can see, the monocular images of the stereo system ( Fig. 7(a)) have a much wider FOV than the images of the plenoptic camera ( Fig. 7(b)). In Figure 7(b), the images of the plenoptic camera seem to be a bit blurred. This is not, in fact, the case and is only due to the multiple projections of a point in neighboring micro images. Due to the narrower FOV and the higher number of pixels on the sensor, images of the plenoptic camera actually have a much higher spatial resolution than those from the monocular cameras.  and 3D point cloud estimated by [17]. Figure 9 shows a magnified subsection of the trajectory and the point cloud of Figure 8(b). This subsection shows the beginning and end of the sequence. One can clearly see the drift accumulated over the complete sequence, resulting in the same scene being reconstructed twice in slightly different locations.

Exemplary Results
By way of example, this section shows results which were obtained for different algorithms based on the presented dataset. The tested algorithms are: monocular: • DSO [7] • ORB-SLAM2 [6,9] stereo: • ORB-SLAM2 [9] plenoptic: • SPO [17] We also ran LSD-SLAM [5] on the dataset. Though, the algorithms failed on most of the sequences or resulted in extremely high drift metrics. For none of the algorithms did we enforce real time processing. For the algorithms which include a full SLAM framework (ORB-SLAM2 and LSD-SLAM), large scale loop closure detection and relocalization was disabled. The implementations of DSO and LSD-SLAM are not able to handle the high image resolution of 1.3 megapixel of the monocular images. Thus, for these algorithms, the image resolution is reduced to 960 pixel × 720 pixel. Both versions of ORB-SLAM2 run at the full image resolution of 1280 pixel × 1024 pixel. Figure 10 shows the results for all algorithms stated above, which were obtained based on the presented dataset. Obviously, no absolute scale error d s can be measured for the monocular algorithms. Depending on the implementation a VO algorithm either signals Fig. 9. Example of the accumulated drift of [17] from the beginning to the end of a sequence. Due to the drift in the trajectory, the same scene is reconstructed twice at different locations. The green line represents the camera trajectory estimated by the algorithm.
a tracking failure or results in an abnormally high tracking error. In Figure 10 we chose appropriate graph limits. All values at the upper graph border signify that the respective algorithm either failed, and therefore no metric could be measured, or that the measured metric lies above the upper graph limit.
Furthermore, Figure 11 shows, by way of example, The 3D point clouds, for some of the sequences, calculated by [17]. These point clouds are supposed to give an impression about the recorded sequences. More results can be found in [17].
From Figure 10 one can see that there are particular sequences for which the plenoptic camera based approach (SPO [17]) perform worse than the other algorithms. These sequences are on one side indoor sequences (e.g. #7 and #8) which show corridors and staircases with lots of white walls. For these scenes monocular and stereo approaches benefit from the wider field of view. On the other side SPO fails in an outdoor sequence (#5), where a van is driving trough the scene. Here, the monocular and stereo approaches again benefit from the wider field of view, while large areas of the light field image are covered by the driving car.

Known Limitations
While the two monocular cameras of the stereo system both have a monochromatic sensor, the plenoptic camera has a RGB sensor. Hence, even though all image sensors have a similar pixel size, a pixel of the plenoptic camera captures only approximately a third of the light energy compared to a pixel of the monocular cameras 5 , when we, in fact, assume that all other parameters are similar. For the plenoptic camera the F-number is predefined by construction, due to the aperture of the micro lenses. This F-number is higher than the one of the lenses used for the two monocular cameras. For this reason, the monocular cameras gather even more light energy on the same sensor area compared to the plenoptic  [17] DSO [7] ORB2 (mono) [9] ORB2 (stereo) [9] (b) scale drift SPO [17] DSO [7] ORB2 (mono) [9] ORB2 (stereo) [9] (c) alignment error (e) (f) Fig. 11. Examples of point clouds reconstructed by [17]. The lengths of all 11 recorded trajectories range from 25 m to 274 m.
camera. To compensate for these two issues, a hardware-sided amplification of 6 dB was set for the plenoptic camera. Thus, for all cameras the exposure times are within the same order of magnitude, although they do not match exactly. Due to the amplification, the images of plenoptic cameras will contain more noise when compared with the monocular cameras of the stereo system. For the stereo camera system, it is important that both cameras run synchronized and with the same exposure time. For this reason, the automatically calculated exposure time of the master camera must be used to set the exposure time of the slave camera. It can happen that if the exposure time changes, an image pair is captured for which the two cameras had slightly different exposure times.
Currently, there exists no plenoptic camera based VO algorithm which performs loop closures. Therefore, the loop closure ground truth, calculated on the basis of the stereo images, is also used as ground truth for the plenoptic camera. With respect to the plenoptic camera, the ground truth might be slightly inaccurate due to the slightly different positions of the master camera of the stereo system and the plenoptic camera. The superior way would be to calculate a second ground truth on the basis of the plenoptic images.
Due to the reason that we use different sensors and lenses which have different properties, and the fact that the cameras see the scene from sightly different perspectives, one has to keep in mind that even though we are able to perform quantitative evaluations based on the presented dataset, these quantities are only valid up to a certain degree. However, the dataset helps to emphasize the strength of VO based on a certain sensor with respect to the other sensors. The results presented in this paper as well as [17] especially show that plenoptic camera based VO offer a promising alternative to approaches based on traditional sensors.


extfieldsofstudy:
Computer Science
Mathematics

text:
Promises and Caveats of Uplink IoT Ultra-Dense Networks

In this paper, by means of simulations, we evaluate the uplink (UL) performance of an Internet of Things (IoT) capable ultra-dense network (UDN) in terms of the coverage probability and the density of reliably working user equipments (UEs). From our study, we show the benefits and challenges that UL IoT UDNs will bring about in the future. In more detail, for a low-reliability criterion, such as achieving a UL signal-to-interference-plus-noise ratio (SINR) above 0 dB, the density of reliably working UEs grows quickly with the network densification, showing the potential of UL IoT UDNs. In contrast, for a high-reliability criterion, such as achieving a UL SINR above 10 dB, the density of reliably working UEs remains to be low in UDNs due to excessive inter-cell interference, which should be considered when operating UL IoT UDNs. Moreover, considering the existence of a non-zero antenna height difference between base stations (BSs) and UEs, the density of reliably working UEs could even decrease as we deploy more BSs. This calls for the usage of sophisticated interference management schemes and/or beam steering/shaping technologies in UL IoT UDNs.

I. INTRODUCTION
Recent years have seen rapid advancement in the development and deployment of Internet of Things (IoT) networks, which can be attributed to the increasing communication and sensing capabilities combined with the falling prices of IoT devices [1]. For example, base stations (BSs) can be equipped with the latest IoT technologies, such as the new generation of machine type communications [2], to collect data from gas, water, and power meters via uplink (UL) transmissions. In practice, such BSs can take the form of both terrestrial and aerial ones [3].
This poses, however, a challenge to the wireless industry, which must offer an increasing volume of reliable traffic in a profitable and energy efficient manner, especially for the UL communications. In this context, the orthogonal deployment of ultra-dense (UD) small cell networks (SCNs), or simply ultra-dense networks (UDNs), have been selected as one of the workhorse for network enhancement in the fourth-generation (4G) and fifth-generation (5G) networks developed by the 3rd Generation Partnership Project (3GPP) [2]. Here, orthogonal deployment means that UDNs and macrocell networks operate on different frequency spectrum, which simplifies network management due to no inter-tier interference.
The performance analysis of IoT UDNs is, however, particularly challenging for the UL because i) UDNs are fundamentally different from the current sparse/dense networks [4] and ii) the UL power control mechanism operates according to the random user equipment (UE) positions in the network, which is quite different from the constant power setting in the downlink (DL) [5].
In this paper, by means of simulations, we evaluate the network performance of UL IoT UDNs in terms of the coverage probability and the density of reliably working UEs. The main findings of this paper are as follows: • We find that for a low-reliability criterion, such as achieving a UL signal-to-interference-plus-noise ratio (SINR) above γ = 0 dB, the density of reliably working UEs quickly grows with the network densification, showing the benefits of UL IoT UDNs. In contrast, for a highreliability criterion, such as achieving a UL SINR above γ = 10 dB, the density of reliably working UEs remains low in UDNs due to excessive inter-cell interference, which should be considered when operating UL IoT UDNs. • We find that due to the existence of a non-zero antenna height difference between BSs and UEs, the density of reliably working UEs could even decrease as we deploy more BSs in a UL IoT UDN. This calls for the usage of sophisticated interference management schemes and/or beam steering/shaping technologies in UL IoT UDNs. • We find that the correlated shadow fading allows a BS with a lower environmental fading factor to provide connection to a larger number of UEs. Thus, its theoretical analysis is an open problem for further study. • We find that the optimized hexagonal-like BS deployment can improve network performance for relatively sparse networks, but not for UDNs. Thus, its theoretical study is not urgent.

II. DISCUSSION ON THE ASSUMPTIONS OF UL IOT UDNS
In this section, we discuss several important assumptions in UL IoT UDNs.

A. BS Deployments
In general, to study the performance of a UL IoT network, two types of BS deployments can be found in the literature, i.e., the hexagonal and the random deployments as shown in Fig. 1.
In Fig. 1, BSs are represented by markers "x" and cell coverage areas are outlined by solid lines. Note that the hexagonal BS deployment leads to an upper-bound performance because BSs are evenly distributed in the network scenario,  and thus very strong interference due to close BS proximity is precluded [6]. In contrast, the random BS deployment reflects a more realistic network deployment with more challenging interference conditions [2,7]. For completeness, in our following performance evaluation, we will consider both BS deployments.

B. Antenna Heights
In the performance analysis of the conventional sparse or dense cellular networks, the antenna height difference between BSs and UEs is usually ignored due to the dominance of the horizontal distance. However, with a much shorter distance between a BS and its served UEs in an UDN, such antenna height difference becomes non-negligible [7]. The performance impact of such antenna height difference between BSs and UEs on the DL UDNs has been investigated in [8]. More specifically, the existence of a non-zero antenna height difference between BSs and UEs gives rise to a non-zero cap on the minimum distance between a BS and its served UEs, and thus a cap on the received signal power strength. Thus, and although each inter-cell interference power strength is subject to the same cap, the aggregate inter-cell interference power will overwhelm the signal power in an UDN due to the sheer number of strong interferers. Consequently, the antenna height difference between BSs and UEs should be considered in the performance evaluation of UL IoT UDNs.

C. Line-of-Sight Transmissions
A much shorter distance between a BS and its served UEs in an UDN implies a higher probability of line-of-sight (LoS) transmissions. The performance impact of such LoS transmissions on the DL of UDNs has been shown to be significant in [9,10]. Generally speaking, LoS transmissions are more helpful to enhance the received signal strength than non-line-of-sight (NLoS) transmissions. However, after a certain level of BS densification, not only the signal power, but also the inter-cell interference power will significantly grow due to the emergence of LoS interfering paths. Thus, the transition of a large number of interfering paths from NLoS to LoS will overwhelm the signal power in an UDN. Consequently, the probabilistic LoS transmissions should also be considered in the performance evaluation of UL IoT UDNs.

D. BS Idle Mode
Considering the surplus of BSs in UDNs, a BS can be put to sleep when there is no active UE connected to it, which is referred as the BS idle mode capability (IMC) [11,12]. As a result, the active UEs' SINR can benefit from i) a BS selection diversity gain, i.e., each UE has a plurality of BSs to select its server from, and ii) a tight control of inter-cell interference, as effective UL inter-cell interference only comes from active UEs served by active neighboring BSs. The surplus of BSs together with the IMC can be seen as a powerful tool that can mitigate the interference problems presented in the previous subsections. However, it should be noted that switching off BSs has a negative impact on the number of IoT UEs that can concurrently transmit.

III. SYSTEM MODEL
For a certain time-frequency resource block (RB), we consider a UL IoT network with BSs deployed on a plane according to the hexagonal deployment or the random deployment, as shown in Fig. 1. For both BS deployments, the density of BSs is denoted by λ BSs/km 2 . Furthermore, we consider a homogeneous Poisson point process (HPPP) Φ to characterize the random deployment.
Active UEs are assumed to be distributed following an HPPP with a density of ρ UEs/km 2 . Here, we only consider active UEs in the network because non-active UEs do not trigger data transmission. Note that the total number of UEs, e.g., phones, gateways, sensors, tags, etc., in a UL IoT network should be much higher than the number of the active UEs. However, we believe that in a certain time-frequency RB, the active UEs with non-zero data traffic demands should still be not many. For example, a typical density of active UEs is around 300 UEs/km 2 in 5G [2].

A. Channel Model
The two-dimensional (2D) distance between a BS and a UE is denoted by r. Moreover, the absolute antenna height difference between a BS and a UE is denoted by L. Thus, the 3D distance between a BS and a UE can be expressed as (1) Note that the value of L is in the order of several meters [13]. Following [10], we adopt a general path loss model, where the path loss ζ (w) is a multi-piece function of w written as where each piece ζ n (w) , n ∈ {1, 2, . . . , N } is modeled as where • ζ L n (w) and ζ NL n (w) , n ∈ {1, 2, . . . , N } are the n-th piece path loss functions for the LoS and the NLoS cases, respectively, • A L n and A NL n are the path losses at a reference 3D distance w = 1 for the LoS and the NLoS cases, respectively, • α L n and α NL n are the path loss exponents for the LoS and the NLoS cases, respectively, and • Pr L n (w) is the n-th piece LoS probability function that a transmitter and a receiver separated by a 3D distance w has an LoS path, which is assumed to be a monotonically decreasing function with respect to w. Existing measurement studies have confirmed this assumption [13]. Moreover, we assume that each BS/UE is equipped with an isotropic antenna, and that the multi-path fading between a BS and a UE is independently identical distributed (i.i.d.) Rayleigh distributed [9,10,14].

B. User Association Strategy
We assume a practical user association strategy (UAS), in which each UE is connected to the BS giving the maximum average received signal strength [10,14]. Such UAS can be formulated by whereR b (w) denotes the average received signal strength from BS b and the UE of interest, separated by a distance of w. Assuming a constant BS transmission power,R b (w) can be equivalently evaluated by ζ (w) defined in (2). As a special case to show our numerical results in the simulation section, we consider a practical two-piece path loss function and a two-piece exponential LoS probability function, defined by the 3GPP [13]. More specifically, in (2) we use N = 2, ζ L 1 (w) = ζ L 2 (w) = A L w −α L , ζ NL 1 (w) = ζ NL 2 (w) = A NL w −α NL , Pr L 1 (w) = 1 − 5 exp (−R 1 /w), and Pr L 2 (w) = 5 exp (−w/R 2 ), where R 1 = 156 m, R 2 = 30 m, and d 1 = R1 ln 10 = 67.75 m [13]. For clarity, this path loss case is referred to as the 3GPP Case hereafter.

C. BS Activation Model
As discussed in Subsection II-D, a BS will enter an idle mode if there is no UE connected to it. Thus, the set of active BSs is determined by the UAS. Since UEs are randomly and uniformly distributed in the network and given the adopted UAS strategy, we can assume that the active BSs also follow an HPPP distributionΦ [11], the density of which is denoted byλ BSs/km 2 , whereλ ≤ λ andλ ≤ ρ. Note thatλ also characterizes the density of active UEs because no collision exists in the centralized cellular IoT UDNs.
For illustration purposes, considering a single-slope path loss model and a nearest-BS UAS,λ can be calculated as [11] where an empirical value of 3.5 was suggested for q in [11] 1 .

D. UL Power Control Model
The UE power, denoted by P UE , is subject to semi-static power control (PC) in practice. In this paper, we adopt the fractional path loss compensation (FPC) scheme standardized in 4G [13], which can be modeled as where P 0 is the target received power in dBm on each RB at the BS, η ∈ (0, 1] is the FPC compensation factor, and N RB is the number of RBs in the frequency domain.

E. The Coverage Probability
Based on this system model, we can define the coverage probability that the typical UE's UL SINR is above a designated threshold γ as where the UL SINR is calculated by where b o denotes the serving BS of the typical UE, P UE bo is the UE transmission power given by (6), w bo is the distance from the typical UE to its serving BS b o , h is the Rayleigh channel gain modeled as an exponentially distributed random variable (RV) with a mean of one as mentioned above, P U N is the additive white Gaussian noise (AWGN) power at the serving BS b o , and I U agg is the UL aggregate interference.

F. The Density of Reliably Working UEs
Based on the definitions of the active BS density in Subsection III-C and the coverage probability in Subsection III-E, we can further define a density of reliably working UEs that can operate above a target UL SINR threshold γ as where the active BS densityλ measures the maximum density of UEs that can simultaneously transmit, and the coverage probability p cov (λ, γ) scales downλ, giving the density of reliably working UEs. The larger the UL SINR threshold γ, the higher the reliability of the IoT communications, and thus the less the UEs that can simultaneously achieve such reliability.

G. More Refined Assumptions
In performance analysis, the multi-path fading is usually modeled as Rayleigh fading for simplicity. However, in the 3GPP, a more practical model based on generalized Rician fading is widely adopted for LoS transmissions [15]. In such model, the K factor in dB scale (the ratio between the power in the direct path and the power in the other scattered paths) is modeled as K[dB] = 13 − 0.03w, where w is defined in (1). More specifically, let h L denote the multi-path fading power for LoS transmissions. Then, for the 3GPP model of Rician fading, h L follows a non-central chi-squared distribution with its PDF given by [16] where K is the distance-dependent value discussed above and I 0 (·) is the 0-th order modified Bessel function of the first kind [16]. Moreover, the shadow fading is also usually not considered or simply modeled as i.i.d. RVs in performance analysis. However, in the 3GPP, a more practical correlated shadow fading is often used [13,15,17], where the shadow fading in dB unit is modeled as zero-mean Gaussian RV [13]. More specifically, the shadow fading coefficient in dB unit between BS b and UE u is formulated as [13] where τ is the correlation coefficient of shadow fading, S UE Shad . In [13], it is suggested that τ = 0.5 and σ Shad = 10 dB.
Considering the distance-dependent Rician fading for LoS transmissions and the correlated shadow fading, we can upgrade the 3GPP Case to an advanced 3GPP Case. In the next section, we will present simulation results of UL IoT UDNs for both the 3GPP Case and the Advanced 3GPP Case. It should be noted that for the Advanced 3GPP Case, shadow fading should be considered in the computation ofR b (w), i.e.,R b (w) should be evaluated by ζ (w) × 10 S bu 10 in (4).

A. Performance Results of the 3GPP Case
In Figs. 2 and 3, we plot the performance results of the 3GPP Case for γ = 0 dB and γ = 10 dB, respectively. Here, we only consider realistic networks with the random deployment of BSs. From these two figures, we can draw the following observations: • Figs. 2a and 3a show that the active BS densityλ monotonically increases with the network densification, and it is bounded by ρ = 300 UEs/km 2 . Such results are in line with the analytical results in (5) [11,12]. However, the density of reliably working UEsρ, i.e.,λp cov (λ, γ) defined in (9), does not necessarily grow as the BS density λ increases. This is because p cov (λ, γ) is a non-monotone function with respect to λ, which will be explained in the following. • When the BS density λ is around λ ∈ 10 −1 , 70 BSs/km 2 , the network is noise-limited, and thus p cov (λ, γ) increases with λ as the network is lightened up with coverage and the signal power strength benefits form LoS transmissions. • When the BS density λ is around λ ∈ [70, 400] BSs/km 2 , p cov (λ, γ) decreases with λ. This is due to the transition of a large number of interfering paths from NLoS to LoS, which accelerates the growth of the aggregate intercell interference. Such performance behavior has been reported in [10] for the DL and [18] for the UL. • When λ ∈ 400, 10 4 BSs/km 2 , p cov (λ, γ) continuously increases thanks to the BS IMC [12], i.e., the signal power continues increasing with the network densification due to the BS diversity gain, while the aggregate interference power becomes constant, asλ is bounded by ρ. However, as shown in Figs. 2b and 3b, the antenna height difference L between BSs and UEs has a large impact on p cov (λ, γ) because a non-zero L places a bound on the signal power strength, which degrades the coverage performance. In more detail, when λ = 10 4 BSs/km 2 and γ = 0 dB, the coverage probability with L = 8.5 m loses 13 % compared to that with L = 0 m. Such performance degradation further enlarges to 32 % when λ = 10 4 BSs/km 2 and γ = 10 dB, showing a much less chance of UE working reliably above 10 dB.    densification, showing the benefits of UL IoT UDNs. In contrast, for a high-reliability criterion, such as surpassing a UL SINR threshold of γ = 10 dB, the density of reliably working UEsρ does not exhibit a satisfactory performance even in UDNs, e.g., we merely getρ < 50 UEs/km 2 when λ = 10 3 BSs/km 2 . The situation only improves when the BS IMC fully kicks in, e.g., λ > 10 3 BSs/km 2 . -Considering the existence of a non-zero L, the density of reliably working UEsρ could even decrease as we deploy more BSs (see Fig. 3c, λ ∈ [200, 600] BSs/km 2 ). This calls for the usage of sophisticated interference management schemes [19] in UL IoT UDNs. Another solution to mitigate such strong inter-cell interference is beam steering/shaping using multi-antenna technologies [2].

B. Performance Results of the Advanced 3GPP Case
In Figs. 4 and 5, we plot the performance results of the Advanced 3GPP Case for γ = 0 dB and γ = 10 dB, respectively. To make our study more complete, here we consider both the random and the hexagonal deployments of BSs. From these two figures, we can see that the previous conclusions are qualitatively correct, which indicates that it is not urgent to investigate Rician fading and/or correlated shadow fading in the context of UDNs. However, there are two new observations that are worth mentioning: • From Figs. 2a/3a and Figs 4a/5a, we can see that the active BS densityλ of the Advanced 3GPP Case is smaller than that of the 3GPP Case. This means that (5) is no longer accurate to characterizeλ for the Advanced 3GPP Case. This is because the correlated shadow fading allows a BS with a lower environmental fading factor, i.e., S BS b , to attract more UEs than other BSs with higher values of S BS b . Its theoretical analysis is an open problem for further study.
• The hexagonal deployment of BSs can improve network performance for relatively sparse networks (e.g., λ < 10 2 BSs/km 2 ), but not for UDNs (e.g., λ > 10 3 BSs/km 2 ). Such conclusion indicates that it is not urgent to investigate the performance of UDNs with the hexagonal deployment of BSs, which has been a longstanding open problem for decades [5].

V. CONCLUSION
We presented simulation results to evaluate the network performance of UL IoT UDNs. From our study, we can see that for a low-reliability criterion, the density of reliably working UEs grows quickly with the network densification. However, in our journey to realize a more reliable UL IoT UDNs, we should be aware of several caveats: • First, for a high-reliability criterion, the density of reliably working UEs remains low in UDNs due to excessive inter-cell interference, which should be considered when operating UL IoT UDNs.
• Second, due to the existence of a non-zero antenna height difference between BSs and UEs, the density of reliably working UEs could even decrease as we deploy more BSs. This calls for further study of UL IoT UDNs. • Third, well-planned hexagonal-like BS deployments can improve network performance for relatively sparse networks, but not for UDNs, showing that alternative solutions other than BS position optimization should be considered in the future UL IoT UDNs.


extfieldsofstudy:
Medicine
Biology
Computer Science
Chemistry

text:
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.