In condensed matter simulations, how is particle number density computed in practice?

Question

I have been reading a recent paper. In it, the authors performed molecular dynamics (MD) simulations of parallel-plate supercapacitors, in which liquid resides between the parallel-plate electrodes. To simplify the situation, let us suppose that the liquid between the electrodes is argon liquid.

The system has a "slab" geometry, so the authors are only interested in variations of the liquid structure along the $z$ direction. Thus, the authors compute the particle number densities averaged over $x$ and $y$: $\bar{n}_\alpha(z)$, where $\alpha$ is a solvent species. (That is, in my simplified example, $\alpha$ is argon -- an argon atom.) $\bar{n} _\alpha(z)$ has dimensions of $\frac{\text{number}}{\text{length}^3}$ or simply $\text{length}^{-3}$, I think.

The $xy$-plane is given by the inequalities $-x_0 < x < x_0$ and $-y_0 < y < y_0$. The area $A_0$ of the $xy$-plane is thus given by $A_0 = 4x_0y_0$.

So, the authors define the particle number density averaged over $x$ and $y$ as follows: $$\bar{n}_\alpha(z) = A_0^{-1} \int_{-x_0}^{x_0} \int_{-y_0}^{y_0} dx^\prime dy^\prime n_\alpha(x^\prime, y^\prime, z)$$ where $A_0 = 4x_0y_0$ and $n_\alpha(x, y, z)$ is the local number density of $\alpha$ at $(x, y, z)$.

Thus, $\bar{n}_\alpha(z)$ is simply proportional to $n_\alpha$ integrated over $x$ and $y$. But, my question is, what is $n_\alpha(x, y, z)$? How is $n_\alpha(x, y, z)$ determined in practice?

As far as the computer is concerned, the argon atoms are point particles; they are modeled as having zero volume (although they interact by Lennard-Jones interactions). So how is it possible to define a number density?

Do we simply "cut" the "slab" in "slices" along $z$ and then assign the particles to these slices? There might be 5 particles in the first $z$ slice, 10 in the second, 7 in the third, and so on. If I then divide 5, 10, and 7 by the volume of the respective slice, then I have a sort of number density, with units of $\frac{\text{number}}{\text{length}^3}$ or simply $\text{length}^{-3}$. But how do I now integrate this $n_\alpha(x^\prime, y^\prime, z)$ over $x$ and $y$? Do I have to additionally perform binning in the $x$ and $y$ directions?

Answer 1

Without seeing the paper, it's hard to know for sure, but the actual particle density probably takes the form

$$n_\alpha(x,y,z) = \sum_{i\in\text{ particles}} \delta^{(3)}(x_i, y_i, z_i)$$

When you integrate this over $x$ and $y$ and some small range $\Delta z$, you get the number of particles in the region you integrated over. So a computer wouldn't actually have to do an integral, it would just count the number of particles in the region. In other words, the simulation probably works with $\bar{n}_\alpha$ directly, not $n_\alpha$.

Answer 2

To calculate $\bar{n}_\alpha$ it's pretty much just what you said. You take the slice between, say $z=2.3$nm and $z=2.301$nm, and count the average number of atoms in it. Divide that number by the volume of the slice (cross-sectional area of the simulation box, multiplied by the slice thickness, i.e. 0.001nm). The answer you get is the number density at $z=2.3$nm

In practice: Each simulation snapshot, you write down the z-coordinate of each atom. As the simulation goes on, you get a larger and larger list of real numbers---all those z-coordinates. Now, plot those numbers in the form of a histogram. If you have a long enough simulation, you can make the bin size of your histogram very very small, so the histogram will look like a smooth curve. (Make sure you scale the histogram so that the integral under the curve is the total number of particles in the simulation divided by the cross-sectional area.)

You never have to explicitly bin or integrate over x and y, if all you need is $\bar{n}_\alpha$.

An alternate approach to calculating $\bar{n}_\alpha$---although it makes no sense to do it this way---is to calculate $n_\alpha$ first, then $\bar{n}_{\alpha}$ second. For the first step, you need to bin in the x,y,z directions---draw little cubes, count the average number of atoms in them, divide by volume. For the second step, you use the formula you cited to integrate $n_\alpha$ over x and y, then divide by cross-sectional area (or in simpler terms, take the mean value of $n_\alpha(x,y,z)$ as $x$ and $y$ vary but $z$ is fixed).

I think you may have gotten confused because the authors discuss the concept of averaging over $x$ and $y$, but you can and should calculate $\bar{n}_\alpha$ without actually explicitly doing that as a separate step.