Better Than Counting: Density Profiles from Force Sampling

Calculating one-body density profiles in equilibrium via particle-based simulation methods involves counting of events of particle occurrences at (histogram-resolved) space points. Here, we investigate an alternative method based on a histogram of the local force density. Via an exact sum rule, the density profile is obtained with a simple spatial integration. The method circumvents the inherent ideal gas fluctuations. We have tested the method in Monte Carlo, Brownian dynamics, and molecular dynamics simulations. The results carry a statistical uncertainty smaller than that of the standard counting method, reducing therefore the computation time.

Figure 1

(a) Density profiles obtained with MC simulations for different numbers of MCS, as indicated. The bin size is $\mathrm{\Delta }x/\sigma =0.01$, $N=25$, and $L/\sigma =10$. The top panels show $\rho (x)$ obtained via the traditional counting method (black solid lines). The density profiles obtained via force sampling are represented in the bottom panels (blue-dashed lines). The inset in the top panel with $10^{11}$ MCS is a close view of both methods in the vicinity of the density peak. Only one-fifth of the simulation box $x/\sigma \in [0,2]$ corresponding to one density peak is represented. (b) Logarithmic plots of the sampling error $\mathrm{\Delta }$ as a function of (i) the number of Monte Carlo steps in MC simulations (top), (ii) the simulation time $t/\tau$ in BD simulations (middle), and MD simulations (bottom). In BD, the time is measured in units of $\tau ={\sigma }^2\gamma /\epsilon$, with $\gamma =1$ the friction coefficient. In MD, $\tau =\sigma \sqrt{m/\epsilon }$, with $m=1$ the mass of the particles. Data obtained via counting (black circles) and via force sampling (blue squares).

Figure 2

(a) Density profiles obtained with MC simulations using counting (left) and force sampling (right). The number of MCS is $10^5$ (top panels) and $10^7$ (bottom panels), as indicated. The bin size is $\mathrm{\Delta }x/\sigma =0.01$. The particles are in equilibrium in an external parabolic potential. The insets are close views of the region where the density vanishes. Using force sampling, the density in this region reaches an artificial negative value of $\sim -2\times {10}^{-2}$ for $10^5$ MCS and $\sim -4\times {10}^{-3}$ for $10^7$ MCS. (b) Second ${\rho }^{(2)}$ (black lines) and third ${\rho }^{(3)}$ (red lines) numerical derivatives (centered difference) of $\rho (x)$ with respect to $x$, ${\rho }^{(n)}(x)={\partial }^n\rho (x)/\partial x^n$. Data obtained via counting (left) and force sampling (right) in a MC simulation with $4\times {10}^{11}$ MCS.

Figure 3

Sampling error $\mathrm{\Delta }$ as a function of the bin size. Data obtained via counting (black circles) and force sampling (blue squares) using MC simulations with $10^9$ MCS. The “true” equilibrium profile used to compute $\mathrm{\Delta }$ is approximated by the average profile of both methods after $10^{12}$ MCS.

Figure 4

(a) Contour plots of the equilibrium 2D density profile in the external potential $V_{\mathrm{ext}}(\mathbf{r})=V_0\sin (2\pi n_{w}x/L)\sin (2\pi n_{w}y/L)$, with $V_0/\epsilon =1$ and $n_w=5$. Data obtained with $10^6$ MCS using counting (left) and force sampling (right). The bins are squares of side length $0.025\sigma$. The simulation box is a square of side length $10\sigma$. Only the central region of the box is shown. (b) Density profile vs $x$ at constant $y/\sigma =5$ obtained via counting (left) and force sampling (right).