Density profile by stress/cartesian

Hello all,

I am performing simulations to measure the osmotic pressure of a system of hard spheres. The setup looks like:

where the red planes represent semipermeable walls that act only on the red particles. Blue particles are free to move in all of the simulation cell. Both particle types have the same size (range of hard-core repulsion) but red ones are pictured larger for clarity.

I use compute stress/cartesian to get the pressure profile in the x-direction. However, I noticed that the density profile calculated by it shows an anomalous behaviour near the periodic boundaries. As an example, the density profile of the above pictured system, with lx=120 and a bin size of 5 for the compute stress/cartesian gives:

The “bump” in the middle of the plot is the expected consequence of osmotic pressure. The problem is encountered near the system boundary, where you can see that on the left side there is a sudden increase in density, balanced out by a sudden decrease on the right. It looks like some particles belonging to the right-most bin are considered a part of the left-most one, leading to this imbalance.

This behavior is consistently present in all the examples I run. I made sure that the bin center positions are placed correctly, in the sense that the left-most and right-most bins start and finish at the edge of the simulation box respectively.

Any ideas are welcome,
Christos

I guess this has to do with the sampling along x. Can you plot the number of particles per bin?
You can also try to apply an offset along x, which will create 25 bins:

variable bindelta equal lx/24
compute  layers  all chunk/atom bin/1d x $(xlo-v_bindelta/2.) ${bindelta} units box

I tried your suggestion. By using compute chunk/atom bin/1d I was able to recover the expected density profile even when using 24 bins (without applying an offset). Bin positions are exactly the same as those reported by compute stress/cartesian, so it looks like a problem limited to this specific compute. It also affects the calculation of the kinetic part of the pressure tensor, so it cannot be ignored by using a different compute just for the density profile.

I recently noticed this as well. I think it is a bug that needs to be fixed.

@chrisp This issue should be fixed by Collected small changes and fixes by stanmoore1 · Pull Request #4445 · lammps/lammps · GitHub.

2 Likes

I think there may be a similar problem with the stress calculation and periodic boundary conditions, will try to debug when I get a minute.

Another bug should be fixed by Collected small changes and fixes by stanmoore1 · Pull Request #4445 · lammps/lammps · GitHub. This was in the stress contour integration when using periodic boundary conditions.