# Response to EiC Nature Comp. Sci. and Reviewer #3

## Reviewer's comments on revision

> In Table 6, why are the computational domains chosen as cubics for all proteins for APBS? It is unfair for APBS in terms of efficiency.

We added results with APBS to satisfy the reviewers' request to compare with other trusted software in the community, as more evidence of the _correctness_ of the computational results. The results provided make no claims or comparisons with regards to _efficiency_. If another choice of computational domain for APBS can make these runs more efficient it will make no difference to the results we present and the claims we make.

> In Table 7, solvation free energy comparison shows very large differences (i.e., up to 1.9%) between MIBPB and the proposed method for some molecules. This is a series problem. The convergence of MIBPB is known in the literature (see for example: Nguyen et al. "Accurate, robust, and reliable calculations of Poisson–Boltzmann binding energies." Journal of computational chemistry 38.13 (2017): 941-948.). As shown in Figure 3 of this reference, the averaged relative absolute error of the electrostatic solvation free energies for all the 153 molecules with mesh size refinements from 1.1 to 0.2 ̊A is under 0.3%. I suggest the authors to carry out the same calculations at Nguyen et al. to find out the averaged relative absolute errors of the present method and plot them against those of MIBPB. Please report the numbers of elements at all resolutions.
>
> As noticed by Fenley and coworkers (Influence of Grid Spacing in Poisson-Boltzmann Equation Binding Energy Estimation, J Chem Theory Comput. 2013 August 13; 9(8): 3677–3685), the calculations of Poisson–Boltzmann binding energies is a challenging task. The authors need to produce reliable Poisson–Boltzmann binding energies as shown by Nguyen et al. (see Figure 6 of the above-mentioned JCC reference) for the challenging task proposed Fenley and coworkers. This test will reveal the level of performance of the present method.

The difference in the solvation energy computed with MIBPB versus Bempp is explained similarly to the differnce with APBS, which we detail in the results section. Quoting from the paper:
"A discrepancy is expected because Bempp-Exafmm and APBS use different geometrical representations for the molecular surface, while the atomic charges in the molecule are mapped to grid points through interpolation in APBS."
Just as APBS, MIBPB is a volumetric grid-based solver, while Bempp is a boundary integral solver.

To clarify, these third-party solvers discretize the Poisson-Boltzmann partial differential equation using finite difference or finite element methods in a discretized volumetric domain. The molecular boundary is approximately represented in the volumetric mesh, and so are atomic positions. This is especially cumbersome in APBS, where the rectangular mesh generates a staircase representation of the molecular surface, and the charges need to be smeared out to the mesh. Bempp, in contrast, solves the equivalent _boundary integral equation_ on a surface triangulation at the molecular interface: the molecular structure is represented with greater fidelity than in volumetric meshes and the atomic positions are exact. While MIBPB is more careful than APBS in this sense, by using a matched-interface technique on the surface, and a regularized form of the equation to avoid smearing out the charges, the molecular geometry will still be slightly different with respect to Bempp. Boundary conditions are also different in both approaches. While MIBPB and APBS impose Dirichlet boundary conditions at the surface of the box enclosing the molecule (with MIBPB not including any dielectric jump in the calculation of the boundary conditions), boundary integral methods (as implemented in Bempp) satisfy boundary conditions at infinity exactly. In sum: even if representing the same physical model, MIBPB and our software solve _different (discretized) equations_ with different boundary conditions on a different geometry, so they cannot be expected to give exactly the same results.

Due to the same reasons, [Geng, Krasny 2013](https://doi.org/10.1016/j.jcp.2013.03.056) observed a 1% discrepancy in the solvation energy when comparing between TABI and APBS using 1A63 (Table 4). Even among grid-based codes, discrepancies of such magnitude have been reported. For example, Table IX in [Geng et al. 2007](https://doi.org/10.1063/1.2768064) reported a 2% difference between MIBPB and APBS in results for protein 1SH1, using a grid spacing of 0.25A.

In Appendix B, we compare our results with Bempp using the finest mesh to the results with MIBPB using a grid spacing of 0.25A (corresponding to a fine grid). The comparison is meant to show agreement, rather than compare the convergence rate of each software, so it is not clear how Nguyen et al.’s study would help in this regard, as well as why we should carry out a similar experiment, given that we have reported Bempp's convergence result in the main body of the paper.

In fact, the observed order of convergence of MIBPB is reported in [Chen et al. 2011](https://doi.org/10.1002/jcc.21646). The error was calculated based on a "designed" analytical solution in [Geng et al. 2007](https://doi.org/10.1063/1.2768064). The result from [Nguyen et al.](https://doi.org/10.1002/jcc.24757) only looks at the "error" with respect to the solution from the finest grid—this is not in fact a measure of the _error_, but just the difference between the solutions obtained with two meshes at different resolutions. The problem with defining the error as the solution difference with the finest mesh is that the latter still contains an error w.r.t. the exact solution. The correct approach is to use Richardson extrapolation to estimate the solution with an infinitely refined mesh. Fig. 3 of Nguyen et al. shows a 0.3% difference of the solution with MIBPB _compared with itself_ using a finer mesh. This is merely a grid-independence result (grids are refined enough that they give almost the same solution). The observation cannot be used to make any claim about how the results with this software should compare to the results with another software, especially if they use different formulations (differential equation versus integral equation) of the physical model, use different discretization schemes, and impose different boundary conditions.

The reviewer is correct in suggesting that binding energies are challenging, as they are small differences between larger solvation energy calculations. In fact, [Harris et al. 2013](https://doi.org/10.1021/ct300765w) (article co-authored by M. Fenley mentioned by the reviewer) report problems in some complexes with a grid spacing as low as 0.3A. However, we don't see how binding energy calculations add to the point we make in this paper: both finite difference and boundary element methods are known to work well in this application, and we are showing evidence of the code correctness. If one were to see a high error in a binding energy calculation (with either finite differences of boundary elements), the solution would be to use a higher mesh density.