Monte Carlo simulation of electrostatic interactions in inhomogeneous dielectric media: Correct sampling for the local lattice simulation algorithm

TL;DR

This paper identifies and corrects a bias in a lattice Monte Carlo algorithm for simulating electrostatic interactions in inhomogeneous dielectrics, improving the physical accuracy of the simulations.

Contribution

The authors introduce a parallel tempering approach to fix unphysical biases in the existing lattice Monte Carlo algorithm for dielectric media.

Findings

Corrected the spurious attractive interactions in the original algorithm.

Validated the improved method on binary mixtures and polymer solutions.

Enhanced the physical realism of electrostatic simulations in inhomogeneous dielectrics.

Abstract

We present a lattice Monte Carlo algorithm based on the one originally proposed by Maggs and Rossetto for simulating electrostatic interactions in inhomogeneous dielectric media. The original algorithm is known to produce attractive interactions between particles of the same dielectric constant in the medium of different dielectric constant. We demonstrate that such interactions are spurious, caused by incorrectly biased statistical weight arising from particle motion during the Monte Carlo moves. We propose a simple parallel tempering algorithm that corrects this unphysical bias. The efficacy of our algorithm is tested on a simple binary mixture and on an uncharged polymer in a solvent, and applied to salt-doped polymer solutions.

Equations

$$\exp (-U_{\mathrm{el}}[\vec{D}])$$

$$⟨S⟩$$

$$\vec{D}(\vec{r})=⟨\vec{A}(\vec{r})⟩.$$

$$⟨\vec{A}(\vec{r})⟩$$

$$p_{\mathrm{swap}}$$

Full text

There is a persistent need for efficient computational methods for modeling dielectric responses of materials and electrostatic interactions in dielectric media. In spite of many advances in methodology development over the past few decades, the computational effort for directly simulating charge systems remains considerably demanding Cisneros et al. (2014). Specifically, accurate treatment of systems with spatially varying dielectric permittivity still poses a substantial challenge in both implicit Maggs and Everaers (2006); Shen and Freed (2002); Barros and Luijten (2014) and explicit Papazyan and Warshel (1997); Levy and Gallicchio (1998); Sagui and Darden (1999) solvent models. Continued progress in the study of electrostatic interactions in diverse systems and phenomena requires a fundamental breakthrough in the simulation techniques capable of efficiently treating dielectric inhomogeneities Messina (2009); Jadhao et al. (2012); Barros and Luijten (2014).

Maggs and Rossetto proposed a local lattice Monte Carlo (MC) algorithm based on solving Gauss’s law Maggs and Rossetto (2002), which scales with the particle number $N$ as $O(N)$. While the algorithm has proven to be strikingly efficient in computational performance Maggs and Rossetto (2002); Rottler and Maggs (2004); Duncan et al. (2005); Rottler and Maggs (2011), the original method is known to produce attractive interactions between uncharged bodies of one dielectric permittivity in the medium of different dielectric permittivity. Duncan et al. Duncan et al. (2006) argued that such attractive force is spurious, likely caused by an inaccurate statistical sampling; yet the nature and origin of the attractive force remain controversial Rottler and Maggs (2011).

The local simulation algorithm to solve Gauss’s law can be articulated by introducing the identity Nakamura and Wang (2013)

[TABLE]

Here, $U_{\mathrm{el}}[\vec{D}]=\int d\vec{r}\,\frac{[\vec{D}(\vec{r})]^{2}}{2\varepsilon_{0}\varepsilon(\vec{r})}$ is the electrostatic energy, with $\vec{D}(\vec{r})$ being the electric displacement and $\rho(\vec{r})$ the local charge density. $\varepsilon_{0}$ and $\varepsilon(\vec{r})$ are the vacuum permittivity and local dielectric function, respectively. The particle position is denoted collectively by $\vec{r}_{i}$ and includes all charged and uncharged species. $Z_{\mathrm{fluc}}=\int{\mathcal{D}}\vec{A}_{t}\{\prod_{r}\delta[\mathrm{div}\vec{A}_{t}]\}\exp(-U_{\mathrm{el}}[\vec{A_{t}}])$ where $\vec{A}_{t}(\vec{r})$= $\vec{A}(\vec{r})-\vec{D}(\vec{r})$ is the transverse vector field and div$\vec{D}(\vec{r})=\rho(\vec{r})$. Using Eq. (1), the statistical average of any property $S$ can be written as Nakamura and Wang (2013):

[TABLE]

Taking $S$ to be $\vec{A}(\vec{r})$, the electric displacement is calculated to be

[TABLE]

The local lattice MC algorithm involves statistical sampling of both the position of the species $\vec{r}_{i}$ and the auxiliary field variable $\vec{A}(\vec{r})$ in accordance with Eq. (2).

The auxiliary field $\vec{A}(\vec{r})$ is not to be confused with the physical electric displacement field $\vec{D}(\vec{r})$. Indeed, the electric displacement field satisfies both rot$\vec{D}(\vec{r})=0$ and div$\vec{D}(\vec{r})=\rho(\vec{r})$ in the absence of a time-dependent magnetic field, whereas div$\vec{A}(\vec{r})=\rho(\vec{r})$ is the only constraint for the auxiliary field. Since $\vec{A}(\vec{r})$ can be written as $\vec{A}(\vec{r})=\vec{D}(\vec{r})+\text{rot}\,\vec{Q}(\vec{r})$, where $\vec{Q}(\vec{r})$ is an arbitrary vector, the functional integration over the transverse vector can be cast into $Z_{\mathrm{fluc}}=\int{\mathcal{D}}\vec{Q}\exp\{-[\text{rot}\,\vec{Q}(\vec{r})]^{2}/[2\varepsilon_{0}\varepsilon(\vec{r})]\}$.

For homogeneous dielectric media, $Z_{\mathrm{fluc}}$ is constant and therefore thermodynamically inconsequential Maggs and Rossetto (2002). However, for systems where the dielectric function depends on the particle positions of the species, no simple treatment of $Z_{\mathrm{fluc}}$ exists. Therefore, $Z_{\mathrm{fluc}}$ is commonly neglected. Such a treatment is known to result in the attractive interactions between uncharged species, which was rationalized as arising from the dispersion interaction Maggs (2004). However, dispersion interaction has its origin in quantum mechanics London (1937) whereas the local simulation algorithm is just a convenient mathematical scheme to solve Gauss’s law (or equivalently Poisson’s equation). Accordingly, the scheme should only account for purely classical effects when the Hamiltonian includes only classical degrees of freedom and interactions. Therefore, when there is no charge in the system, we expect $D(\vec{r})=0$ from div$\vec{D}(\vec{r})=0$ and rot$\vec{D}(\vec{r})=0$, regardless of any dielectric inhomogeneity. Thus, the energy landscape of the system given by $U_{\mathrm{el}}[\vec{D}]$ should be uniform and independent of particle configurations, indicating no net electrostatic interaction between the uncharged species. In this Letter, we first demonstrate that, indeed, the local simulation algorithm does not produce any force in the systems containing no charge. We then present a simple, general method to avoid the undesirable spurious interactions in the local lattice MC simulations via parallel tempering.

We first consider the statistical average of the auxiliary field $\langle\vec{A}(\vec{r})\rangle$ for a system of uncharged species with their positions fixed. Eq. (2) then becomes

[TABLE]

Note that the factor $Z_{\mathrm{fluc}}$ does not appear in this equation because the positions of the species are fixed. When the particles are allowed to move, the value of $Z_{\mathrm{fluc}}$ changes with the particle configuration, and does not cancel out in performing the statistical average Eq. (2). Neglecting $Z_{\mathrm{fluc}}$ therefore generates unphysical bias during the simulation that involves particle moves.

To model the dielectric inhomogeneity, we adopt the same scheme for liquid mixtures as proposed by Maggs Maggs (2004): we assign the dielectric constants $\varepsilon_{A}$ and $\varepsilon_{B}$ to the lattice grids Maggs and Rossetto (2002); Nakamura and Wang (2013) occupied by lower-dielectric and higher-dielectric components, respectively. We perform Metropolis sampling using the Boltzmann weight with energy $A_{i}^{2}/(2\varepsilon_{0}\varepsilon_{i})$ on the i-th lattice site. The current of $A_{i}$ through a surface of the cube that bisects the lattice bond connecting $i$ and $j$ is expressed as $b_{ij}$. A group of the lattice variables $b_{ij}$ are sequentially updated on a lattice plaquette Nakamura and Wang (2013); Duan and Nakamura (2015) using the Boltzmann weight subject to the lattice version of Gauss’s law $\sum_{j}b_{ij}=q_{i}$. Here, $q_{i}$ is the charge placed on the lattice grid. The sizes of the simulation box are set to $L^{3}=(30u)^{3}$ for Fig. 1 and Fig. 3, $L^{3}=(10u)^{3}$ for Fig. 2, and $L^{3}=(20u)^{3}$ for Fig. 4, where the lattice unit $u=2.8$ Å . Periodic boundary conditions are used in all directions. Interested readers are referred to a more complete explanation of the update scheme for the auxiliary field $\vec{A}(\vec{r})$ in Refs. Nakamura and Wang (2013) and Duan and Nakamura (2015).

In Fig. 1, we consider a fixed stripe pattern of the two dielectric species with no charge. The uncharged components with $\varepsilon_{A}=6.3$ and $\varepsilon_{B}=80$ are distributed periodically. We performed statistical analysis by averaging on the order of 3000 samples; good statistical convergence was observed on the order of 100 samples. Note that the simulation for Eq. (4) results in $\vec{D}(\vec{r})=0$ throughout the simulation box, indicating no electrostatic interaction between the species. We also considered a non-periodic pattern but again obtained $\vec{D}(\vec{r})=0$. This implies that any spurious force between uncharged species when particles are allowed to move must be due to incorrect sampling as a result of neglecting $Z_{\mathrm{fluc}}$ in the MC moves. Furthermore, the inset shows the statistical average of the auxiliary field $\vec{A}(\vec{r})$ as a function of the MC steps between consecutive measurements in equilibrium. We find that a minimum of $10^{4}$ MC steps are required to obtain statistically independent samples to reach the correct average of $\vec{A}(\vec{r})$. This result also suggests that at least $10^{4}$ MC steps are required to equilibrate the MC simulation. Thus, combining a straightforward implementation of the local simulation algorithm at fixed particle positions that correctly samples the electrostatic interactions with other simulation methods for moving the particles such as molecular dynamics, poses formidable computational challenges. Existing hybrid methods essentially take $\vec{A}(\vec{r})$ to be $\vec{D}(\vec{r})$Rottler and Maggs (2004); Pasichnyk and Dunweg (2004); Rottler and Maggs (2011); Fahrenberger and Holm (2014); Fahrenberger et al. (2015), which is tantamount to neglecting the factor $Z_{\mathrm{fluc}}$.

We propose a scheme that uses parallel tempering to avoid the biasing errors introduced through the neglect of $Z_{\mathrm{fluc}}$ in the conventional Metropolis implementation of the lattice MC moves. We start by recalling the key ideas of the parallel tempering method as relevant to our current problem. We consider states $a$, $b$, $c$ and $d$ with the energies $E_{a}$, $E_{b}$, $E_{c}$, and $E_{d}$ that correspond to different particle configurations and temperatures $T_{0}$, $T_{1}$, $T_{2}$, and $T_{3}$. The corresponding canonical probability distribution for the state $\alpha$ at temperature $T_{i}$ is p_{i}(E_{\alpha})=Z_{i}^{-1}\exp\bigl{(}\frac{-E_{\alpha}}{T_{i}}\bigr{)}, where $Z_{i}$ is the normalization factor for temperature $T_{i}$. For simplicity, we set the Boltzmann constant $k_{b}$ to 1. We consider parallel tempering with the swaps of the states $a$ and $d$ between $T_{1}$ and $T_{2}$, and $b$ and $c$ between $T_{0}$ and $T_{3}$. The Boltzmann weight for these two simultaneous swaps is p_{\mathrm{swap}}=\exp{\bigl{[}-\frac{(E_{a}-E_{d})}{T_{2}}\bigr{]}}\exp{\bigl{[}-\frac{(E_{d}-E_{a})}{T_{1}}\bigr{]}}\exp{\bigl{[}-\frac{(E_{b}-E_{c})}{T_{0}}\bigr{]}} \exp{\bigl{[}-\frac{(E_{c}-E_{b})}{T_{3}}\bigr{]}}.

For subsequent discussions, we consider the special limiting case with $T_{2}\rightarrow\infty$ and $T_{3}\rightarrow\infty$. The Boltzmann weight for the swaps then becomes

[TABLE]

We will use Eq. (5) to deal with the contributions from the functional integration over the transverse vector in the local simulation method.

We write the dielectric constants for species A and B as $\varepsilon_{B}=\varepsilon_{A}+\eta\Delta\varepsilon$. The real system has $\eta=1$ but for parallel tempering, we allow the dielectric contrast parameter $\eta$ to take any value between [math] and $1$. Let $E(i,\eta,\vec{A})$ denote the energy of a microstate with particle configuration $i$, dielectric contrast $\eta$, and auxiliary field $\vec{A}$. In the local simulation algorithm, Eq. (2) defines $E(i,\eta,\vec{A})$ as $E(i,\eta,\vec{A})=\ln Z_{\mathrm{fluc}}(i,\eta)+\int d\vec{r}\vec{A}(\vec{r})^{2}/[2\varepsilon(i,\eta)]$. For simplicity, we scale the field $\vec{A}$ so that $\varepsilon_{0}=1$.

To make use of Eq. (5), we let the state be represented by the particle configuration $i$ and auxiliary field $\vec{A}$, and let $\eta$ play the role of the temperature $T$. We consider the Boltzmann weight for the following swaps:

[TABLE]

Since Eq.  (6) holds for arbitrary values of the vector field, we choose $\vec{B}(\vec{r})$ and $\vec{C}(\vec{r})$ such that $E(i,\eta,\vec{B})=\ln Z_{\mathrm{fluc}}(i,1)+\int d\vec{r}\,\vec{A}(\vec{r})^{2}/[2\varepsilon(i,\eta)]$ and $E(j,\eta,\vec{C})=\ln Z_{\mathrm{fluc}}(j,1)+\int d\vec{r}\,\vec{A}^{\prime}(\vec{r})^{2}/[2\varepsilon(i,\eta)]$. We thus arrive at

[TABLE]

Note that the unknown factor $Z_{\mathrm{fluc}}$ does not appear in Eq. (7). Equation (7) thus allows us to effect the correct MC transitions from the particle configuration $i$ to $j$ and field configuration $\vec{A}$ to $\vec{A}^{\prime}$ in the physical system with $\eta=1$ without knowledge of $Z_{\mathrm{fluc}}$!

To ensure good overlap of distribution between the different replicas in parallel tempering, differences in the energies of replicas should be small. To this end, we generate $\vec{A}_{n}^{\prime}$ from $\vec{A}_{n}$ by $k$ updates, using the statistical weight $\exp{\{-\sum_{n\in\mathrm{plaquette}}[(A_{n}+a_{n})^{2}-A_{n}^{2}]/[2\varepsilon(i,\eta)]\}}$ with a random number $a_{n}$. In the context of parallel tempering, the number of replicas increases with increasing $k$.

We first test our parallel tempering method on a binary mixture of particles A and B with $\varepsilon_{A}=10.5$ and $\varepsilon_{B}=80$. $\Delta\varepsilon$ is therefore $69.5$. Given the incompressible nature of the lattice model, the volume fractions of particles A and B are denoted by $\phi$ and $1-\phi$, respectively. We set $\phi=0.3$. We consider A-A pairs defined as when two neighboring particles A are separated by a distance of less than or equal to $\sqrt{3}u$ [Fig. 2].

In Fig. 2 we show the average number of neighboring A particles around a tagged A particle as a function of the number of updates. In the absence of any interactions, random mixing gives 7.8, as is confirmed for $\varepsilon_{A}=\varepsilon_{B}$. This is shown as the blue dashed line. The red dashed line for $\varepsilon_{A}\neq\varepsilon_{B}$ without parallel tempering is higher than the blue dashed line, indicating effective attraction between the like particles due to incorrect sampling. We note that this spurious attractive force is indeed obtained in Ref. Maggs (2004). With parallel tempering, for any value of $\eta$ (that is appreciably less than 1), the results converge to the blue dashed line with increasing number of updates of the auxiliary field, $k$.

As the second example, we consider an uncharged polymer of $N=200$ units with lower dielectric constant $\varepsilon_{A}=10.5$ immersed in a higher dielectric medium of $\varepsilon_{B}=80$ (mimicking water). To model the polymer conformation, we employ the bond fluctuation model of Shaffer Shaffer (1994) on a simple cubic lattice with nearest-neighbor distance $u$. The bond vectors of the consecutive monomers on the polymer are allowed to take the bond lengths, 1, $\sqrt{2}$, and $\sqrt{3}$ in the lattice unit $u$. An empty site corresponds to the coarse-grained density of water molecules. The details of the model were described in an earlier publication Nakamura and Wang (2013).

Figure 3 shows the radius of gyration $R_{g}$ for a single uncharged polymer with different values of $\varepsilon_{A}$. In the absence of dielectric contrast and other energetic interactions, the polymer is in an athermal solvent, and so the chain is expanded with $R_{g}$ given by the blue dashed line. If sampling is done correctly, the polymer conformation should be independent of the dielectric contrast. When we ignored $Z_{\mathrm{fluc}}$ using conventional MC moves in the case of $\varepsilon_{A}\neq\varepsilon_{B}$ (red dashed line), the chain adopts a globular structure due to the spurious attractive interaction. Our parallel tempering algorithm eliminates this artificial interaction; the colored symbols from our parallel tempering sampling converge to the same $R_{g}$ as in the case of $\varepsilon_{A}=\varepsilon_{B}$ for $k>8$ .

As a final example, we use our parallel tempering algorithm to study the effect of salts on the miscibility in an uncharged-polymer solution. We consider a solution of polymers with $N=50$ at volume fraction $\phi=0.125$ and assign the nearest neighbor interaction energy $\sigma_{\mathrm{pp}}=-0.5\,k_{b}T$ for the monomer-monomer interaction. For simplicity, the interaction energies for the monomer-solvent and solvent-solvent interactions are set to zero. The dielectric constant of the polymer and solvent are respectively $\varepsilon_{A}=10.5$ and $\varepsilon_{B}=80$. To avoid the volume effects of the salt ions, we allow an anion or cation to occupy the same lattice site as a solvent, but forbid more than one ions from occupying the same site. To calculate the change in the miscibility, we compute the structure $S(q)=\langle|\sum_{n}\exp(i\vec{q}\cdot\vec{r}_{n})[{\delta_{A}}^{n}-{\delta_{B}}^{n}-\langle{\delta_{A}}^{n}-{\delta_{B}}^{n}\rangle]|^{2}\rangle/L^{3}$ Sariban and Binder (1991), where the local concentration variable $\delta_{A}^{n}$ (or $\delta_{B}^{n}$) is 1 or 0 when the lattice site designated by $\vec{r}_{n}$ is occupied or unoccupied by a monomer (or solvent). Due to the periodic boundary condition, the wavevector is given by $\vec{q}=\frac{2\pi}{L}(n_{x},n_{y},n_{z})$, where $(n_{x},n_{y},n_{z})$ is a group of integers. We perform spherical averaging to obtain the structure factor as a function of the wavenumber $q$. As shown in Fig. 4, $S(q)$ at small $q$ increases with increasing salt concentration $c$. The results indicate that the polymer (green sphere) and solvent become less miscible with the addition of the salt ions (red and blue spheres), with the salt ions concentrated in the higher-dielectric-constant solvent regions (see the insets). This feature is consistent with predictions from previous coarse-grained mean-field theories Wang (2008); Nakamura et al. (2012); Nakamura and Wang (2012).

In summary, we have shown that the original algorithm by Maggs and Rossetto Maggs and Rossetto (2002) for treating configuration-dependent dielectric inhomogeneities introduces an unphysical bias in the sampling by neglecting the factor $Z_{\mathrm{fluc}}$ from the integration over the transverse auxiliary field. We have proposed a parallel tempering scheme that effectively corrects for this bias without requiring knowledge of $Z_{\mathrm{fluc}}$. The extra computational cost associated with generating the auxiliary field in the replica exchange is quite modest and the overall $O(N)$ scaling of the original algorithm is preserved. Our new algorithm therefore represents an efficient method for simulating a wide range of phenomena involving charge interactions in inhomogeneous, configuration-dependent dielectric media.

This work was supported by the National Natural Science Foundation of China (21474112 and 21404103) and the start-up funds of Michigan Technological University. We are grateful to Computing Center of Jilin Province for their essential support.