I have a question that I know that is not directly related to Lammps, but it would be great if I could get some answers about it here.
I would like to compute the local volume fraction of beads in an biaxial confinement (a cylinder with the PBC along its longitudinal (z) direction). For the histogram, local number density, or pair distribution of beads/particle, each particle only contributes to the frequency in one bin. However, the volume of each bead, depending on it size, can contribute to the local volume fraction in more than one bins; for instance, a particle of diameter 10 nm contributes to the local volume fraction in 20 bins of width 0.5nm along the z axis.
The PBC is not important for the histogram or number density of particles along the z direction since the center of each particle can only be in one of bins. Is the same true for the local volume fraction along z direction? My guess is that it is not true since the volume of the particles in bins near one end of the cylinder is not preserve unless we add the rest of their volumes to the bins near the other end of cylinder.
First off. In matters of science, you should never “guess”.
Advancement in science happens by making a hypothesis and then trying to prove or disprove it (by experiment or compute simulation or mathematical proof or similar).
While analysis using “chunks” in LAMMPS will only consider point particles and thus the assignment with bins is unique, there are other methods conceivable to do the analysis you want to do.
Those have not been programmed into LAMMPS, mind you (but would not be too difficult now that we have a generic grid communication available). To quantify how much of a sub-volume is “used”, there are essentially two approaches: 1) you compute the “density” usually used for “soft” particles 2) you compute the “occupancy” usually used for “hard” particles. In both cases you overlay your simulation domain with a grid, then loop over the particles and determine which grid points are within the radius from the center of a particle and then compute a contribution to that grid point. in case of density it would be a distance dependent contribution (e.g. assuming a gaussian distribution), in case of occupancy, it would be assigning a value of 1 (after all points were initialized to 0).
Then you can compute the contribution to a given volume by combining the data from the grid points describing that volume.
Periodic boundaries need to be applied in just the same way as for computing interactions.
To come back to my original statement. With my suggestion it should be straightforward for you to write some analysis script/program that can do this procedure. If you are familiar with Tcl scripting, you may also consider using VMD, since it has a “volmap” tool that can do both of these methods (and then some) but you will likely need additional Tcl scripting to get all the properties set for your system in VMD so that the available tool can work. VMD is assuming a CHARMM force field biomolecular simulation by default and expects data to conform to that unless explicitly set differently.
P.S. If you go with VMD’s volmap tool (which is arguably the quickest approach) please use one of the 1.9.4 alpha builds.
Versions of VMD up to 1.9.3 contain a bug in volmap when processing atom selections that vary between frames (e.g. because atoms move in or out of a specified region or other dynamic criterion is used to select). This bug is fixed in all the currently available 1.9.4 alpha builds. All VMD/volmap versions work fine if the selection is constant (i.e. a fixed list of atoms is processed over the whole trajectory).