Local Molecular Dynamics with Coulombic Interactions

We propose a local, $\mathcal{O}(N)$ molecular dynamics algorithm for the simulation of charged systems. The long ranged Coulomb potential is generated by a propagating electric field that obeys modified Maxwell equations. On coupling the electrodynamic equations to an external thermostat we show that the algorithm produces an effective Coulomb potential between particles. On annealing the electrodynamic degrees of freedom the field configuration converges to a solution of the Poisson equation much like the electronic degrees of freedom approach the ground state in ab initio molecular dynamics.

Figure 1

Time for the field integration (□), field-particle couplings (△), and total time (◇) in a system of size $L=30a$ as a function of $N$ on an single AMD Athlon CPU. Also shown are results for the same system treated with a Fourier-based method using a charge interpolation of the same order (filled symbols, ■ corresponds to the total reciprocal work).

Figure 2

Static structure factor $S(q)=⟨s(\mathbf{q})s(-\mathbf{q})⟩,s(\mathbf{q})=(1/\sqrt{N})\sum _i{e}_i\exp (i{\mathbf{r}}_i\mathbf{q})$ of a charged symmetric electrolyte $(L=15a)$; three densities $n=0.05a^{-3}$, $n=0.1a^{-3}$, $n=0.2a^{-3}$, and $T=ϵ/4\pi k_B$ using the present algorithm in mode I (□), $c=1a/\tau$; transverse electric field and particles thermostated to $T=ϵ/4\pi k_B$, $\overline{v}/c\approx 0.3$ using Eqs. (4, 6). For ■, the electric field is damped, mode II, $c=32a/\tau$, $\overline{v}/c\approx 0.01$. Also shown is a corresponding Monte Carlo simulation (△). Solid lines: Debye theory. Inset: average instantaneous force for a pair of charges, $L=8a$ for modes I and II. Solid curve: $f=-dV/dr$, where $V=-1/4\pi r-r^2/6{ϵ}_0{L}^3$ from Ewald summation [20]. Lattice artifacts removed using dynamic subtraction, $\mu =0.9/a$.

Figure 3

Velocity autocorrelation functions for $L=10a$, 50 particles, mode I. Intrinsic damping ${\gamma }_1=0.1{\tau }^{-1}$, dotted line: $\exp (-0.1t/\tau )$. Solid lines from left to right: correlations for $c=0.32,0.35,0.41,0.5,0.7,1,3.2,10a/\tau$. Inset: initial decay rate ${\omega }_0$ from exponential fits.