Comparison of three-dimensional Poisson solution methods for particle-based simulation and inhomogeneous dielectrics

Particle-based simulation represents a powerful approach to modeling physical systems in electronics, molecular biology, and chemical physics. Accounting for the interactions occurring among charged particles requires an accurate and efficient solution of Poisson's equation. For a system of discrete charges with inhomogeneous dielectrics, i.e., a system with discontinuities in the permittivity, the boundary element method (BEM) is frequently adopted. It provides the solution of Poisson's equation, accounting for polarization effects due to the discontinuity in the permittivity by computing the induced charges at the dielectric boundaries. In this framework, the total electrostatic potential is then found by superimposing the elemental contributions from both source and induced charges. In this paper, we present a comparison between two BEMs to solve a boundary-integral formulation of Poisson's equation, with emphasis on the BEMs' suitability for particle-based simulations in terms of solution accuracy and computation speed. The two approaches are the collocation and qualocation methods. Collocation is implemented following the induced-charge computation method of D. Boda et al. [J. Chem. Phys. 125, 034901 (2006)]. The qualocation method is described by J. Tausch et al. [IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 20, 1398 (2001)]. These approaches are studied using both flat and curved surface elements to discretize the dielectric boundary, using two challenging test cases: a dielectric sphere embedded in a different dielectric medium and a toy model of an ion channel. Earlier comparisons of the two BEM approaches did not address curved surface elements or semiatomistic models of ion channels. Our results support the earlier findings that for flat-element calculations, qualocation is always significantly more accurate than collocation. On the other hand, when the dielectric boundary is discretized with curved surface elements, the two methods are essentially equivalent; i.e., they have comparable accuracies for the same number of elements. We find that ions in water—charges embedded in a high-dielectric medium—are harder to compute accurately than charges in a low-dielectric medium.

Figure 1

(a) Sphere test case (a). A dielectric sphere of radius 5 Å is embedded in a different dielectric medium. An elementary point charge is set 4 Å off sphere center. (b, c) Ion-channel test case. The model is obtained by rotating (b) by 180${}^{\circ }$ around the axis. The water and the membrane feature different dielectric constants (${\epsilon }_W=80$ and ${\epsilon }_M=2$, respectively). The trajectories used for comparison tests are noted in (b) with the labels T1 (along the channel axis), T2 (3 Å off the channel axis), and T3 (radial direction from the channel axis at the center of the channel). (c) 3D model of the channel. The filled (blue) circle at the right represents the ion moving on the channel axis (red line).

Figure 2

Sphere test cases with flat tiles. Reaction potential and its relative error along the sphere diameter that passes through the source charge (a, b) for the high-dielectric sphere and (c, d) for the low-dielectric sphere. The sphere goes from $z=-5$ Å to $z=5$ Å. The source charge is located at $z=4$ Å. In this the subsequent figures, for each data curve, legends give the kind of tiles (flat or curved), the BEM implementation (ICC or QUAL), and the number of tiles and subtiles per tile used to discretize the dielectric boundary (e.g., $364\times 100$ stands for 364 tiles, 100 subtiles for each tile).

Figure 3

Sphere test cases, flat tiles. Gauss's Law consistency check (a) for the high-dielectric sphere and (b) for the low-dielectric sphere. Graphs show the total induced charge on the boundary, expressed in elementary charge units ($h_{\text{tot}}/\left|e\right|$), as a function of the number of tiles used to discretize the dielectric boundary. The number of subtiles per tile is varied.

Figure 4

Sphere test cases with curved tiles. Reaction potential and its relative error along the sphere diameter that contains the source charge (a, b) for the high-dielectric sphere and (c, d) for the low-dielectric sphere. Results for both ICC and QUAL were obtained using different numbers of tiles. The sphere goes from $z=-5$ Å to $z=5$ Å. The source charge is located at $z=4$ Å.

Figure 5

Sphere test cases with curved tiles. Gauss's Law consistency check (a) for the high-dielectric sphere and (b) for the low-dielectric sphere. Graphs show $h_{\text{tot}}$ as a function of the number of tiles used to discretize the dielectric boundary. The number of subtiles per tile is indicated for each curve.

Figure 6

Ion-channel test case for trajectory T1. Reaction potential and $h_{\text{tot}}$ as a function of the position of a cation that moves along the channel axis (T1). (a, b) Flat and (c, d) curved tiles are used. Results for both ICC and QUAL were obtained using different numbers of tiles and subtiles.

Figure 7

Ion-channel test case for trajectory T2. Reaction potential and $h_{\text{tot}}$ as a function of the position of a cation that moves 3 Å off the channel axis (T2). (a, b) Flat and (c, d) curved tiles. Results for both ICC and QUAL were obtained using different numbers of tiles and subtiles.

Figure 8

Ion-channel test case for trajectory T3. Reaction potential and $h_{\text{tot}}$ as a function of the position of a cation that moves radially from the channel axis at the center of the channel (T3). (a, b) Flat and (c, d) curved tiles. Results for both ICC and QUAL were obtained using different numbers of tiles and subtiles.