Computing the Coulomb interaction in inhomogeneous dielectric media via a local electrostatics lattice algorithm

The local approach to computing electrostatic interactions proposed by Maggs and adapted by Rottler and Pasichnyk for molecular-dynamics simulations is extended to situations where the dielectric background medium is inhomogeneous. We furthermore correct a problem of the original algorithm related to the correct treatment of the global dipole moment, provide an error estimate for the accuracy of the algorithm, and suggest a different form of the treatment of the self-energy problem. Our implementation is highly scalable on many cores, and we have validated and compared its performance against theoretical predictions and simulation data obtained by other algorithmic approaches.

Figure 1

Recursive scheme for the initial solution of the $E$ field. The average charge in $z$ plane is scaled and added to each node, following $E_z^{(n+1)}=E_z^{\left(n\right)}+q_{\mathrm{plane}}/\left(\varepsilon a^2\right)$. Then the charge $q_{\mathrm{plane}}$ is subtracted from each charge in the $z$ plane. Analog with $y$ lines and the single nodes in $x$ direction.

Figure 2

(a) Discretization of the currents, fields, and permittivities onto a lattice cell. (b) Interpolation of dielectric permittivity values on the lattice. $\varepsilon \left(\mathit{r}\right)$ has a position and a direction (blue arrow). The values for ${\varepsilon }_1$ and ${\varepsilon }_2$ are determined and the value on the connecting link is set to the average value. If the gradient is too large, the value is determined by forming the harmonic average.

Figure 3

Error estimate for the MEMD algorithm: the interpolation and the finite speed of light create numerical errors. Both depend on the lattice spacing. For affirmation, the errors from three simulations of different systems (see Fig. 4) are included, compared to high precision P3M force calculations.

Figure 4

Three example systems to determine the numerical error of the algorithm. Two polyelectrolytes in aqueous solution (a), a melting silica crystal (b), and an artificial system (c) with two oppositely charged walls and a surrounding cloud of charges.

Figure 5

(a) Dipole moments of a simple electrolyte system are compared for the MEMD and a P3M algorithm. The original MEMD keeps the dipole moment at zero, which corresponds to P3M boundary conditions of $\varepsilon =0$ at infinity. The corrected version shows the same behavior as P3M for metallic boundary conditions. (b) Two particles are dragged apart over several lengths of the simulation box via applying an external field. As the dipole moment of the system increases, these two boundary conditions also diverge quantitatively in the force calculation.

Figure 6

(a) Charge is placed in the center, with the dielectric constant linearly rising with distance from the center. This graph shows the absolute value $\sqrt{{\mathit{E}}^2}$ of the resulting field. As theoretically predicted, the system follows a $1/r^3$ behavior. (b) The relative rms error of the electric field. The graph shows that the error stays well below the ${10}^{-3}$ limit, and this relative error is smallest in the cell center.

Figure 7

(a) Two particles are simultaneously pulled through the simulation box (${\varepsilon }_W=80$) with two infinite dielectric walls (${\varepsilon }_C=2$). (b) The results are compared to the analytical prediction. The box center is located at $z=10$, and one wall at $z=0$ on the left of the graph.

Figure 8

(a) Dielectric sphere of ${\varepsilon }_C=2$ is placed in the center of the box and salt ions assemble in bulk water ${\varepsilon }_W=80$ around it. (b) The ion density as a function of the radial distance is analyzed and compared to simulations using the $\mathrm{ICC}☆$ algorithm and a Monte Carlo simulation.

Figure 9

For constant density, the algorithm scales linearly with the particle number from ${10}^3$ to $2\times {10}^6$ particles. For this system, the implementation is slower than the optimized P3M solver, but the inclusion of dielectric changes comes at no extra cost.