Bin size with respect to DEM simulation

voro · February 7, 2023, 12:57am

Hi all,

I notice some non-physical behavior in my DEM jamming simulation and I wonder if it had anything to do with the bin size. For reference, my granular system (surrounded by periodic walls) is shown below and I compress it beyond the RCP volume fraction by using fix deform to shrink in the direction along the red and green plane (x and y). The top and bottom (z) along the direction perpendicular to compression is just periodic.

When compressed beyond the RCP threshold, I get a very high pressure, but the kinetic energy does not dissipate even with a tangential friction contact when let to equilibrate. I instead get oscillatory behavior that I don’t have an explanation for… I expect friction to dissipate this energy and the granular velocities to eventually die down to very low values (Kinetic energy of individual granules ~ < 1e-10J). When I use viscous damping, the kinetic energy momentarily goes down but shoots back up when viscosity is removed so it doesn’t help.

I wonder if it has anything to do with the number of bins… My granules have radii between 0.0045m and 0.005m and the length of the channel is 0.05m x 0.05m x 1m, therefore there are slightly more than 5 granules per length of the box in x and y direction, and 100 granules in the z-direction.

I get the following bins: binsize = 0.0149999 → bins = 4 4 67 when I use neighbor 0.02 bin. I noticed that I get the following info:

The ave neighs/atom is really large and it should ideally be close to 6 and not 106… But when I don’t use neighbor 0.02 bin, I get a smaller binsize and more bins, and a more realistic neighs/atoms of 5.3:

If I build the neighbor list regularly, say like 10 timesteps, will the inclusion of the neighbor command just mean more computation time but would give the same result? Also, how does binsize affect in a serial run?

At the end of a sufficiently large timestep for the case with neighbor command, I get total kinetic energy as 7.0970525J and when its not being used I get 7.095222J so it looks like I get nominally the same result. I am worried if I’m not there are some force pairs that are going undetected resulting in the lack of decrease of kinetic energy?

I’ve attached the kinetic energy and pressure as a function of timesteps to show what I get but I don’t reach near 0 kinetic energy even after playing around with simulation parameters

My input file is as follows, and I’ve attached my data file:

units si
atom_style  	sphere
boundary 	pp pp pp
newton  	off
comm_modify 	vel yes
read_data	n_0p63_periodic_mid.sim
[n_0p63_periodic_mid.sim|attachment](upload://uiyoCBKz1QPQVOquHk0itnLhBWJ.sim) (934.1 KB)

[n_0p63_periodic_mid.sim|attachment](upload://uiyoCBKz1QPQVOquHk0itnLhBWJ.sim) (934.1 KB)


log 			n_0p63_periodic.lammps


#neighbor 	0.02 bin
neigh_modify 	delay 0

pair_style      gran/hertz/history 8.424908425E+10 NULL 0 0 0.05 0

pair_coeff	* *

timestep 	1e-7
fix  		1 all nve/sphere
fix 11 all viscous 5


region 		atom_region block 0 0.05 0 0.05 0 1
region 		atom_region2 block 0.005 0.045 0.005 0.045 0.05 0.8
group 		atoms region atom_region
group		atoms2 region atom_region2


compute        peratom atoms stress/atom NULL pair
compute        p atoms reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
compute        p1 atoms reduce sum c_peratom[4] c_peratom[5] c_peratom[6]
variable       press equal -(c_p[1]+c_p[2]+c_p[3])/(3*0.05*0.05*1)#((y[896]-1e+5)-(y[895]+1e+5))*((z[893]-1e+5)-(z[894]+1e+5)))
variable       stress_ratio equal c_p[1]/c_p[2]


compute        peratom2 atoms2 stress/atom NULL pair
compute        p2 atoms2 reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]
compute        p12 atoms2 reduce sum c_peratom2[4] c_peratom2[5] c_peratom2[6]
variable       press2 equal -(c_p2[1]+c_p2[2]+c_p2[3])/(3*0.04*0.04*0.75)#((y[896]-1e+5)-(y[895]+1e+5))*((z[893]-1e+5)-(z[894]+1e+5)))
variable       stress_ratio2 equal c_p2[1]/c_p2[2]


compute  	rot_kin all erotate/sphere # rotational kinetic energy of the spherical particles

variable myvfrac equal 0.001575126176137/vol#0.005040275060683/vol
thermo_style 	custom step atoms ke c_rot_kin v_press v_press2 c_p[1] c_p[2] c_p[3] c_p2[1] c_p2[2] c_p2[3] vol v_myvfrac
thermo  	1000
thermo_modify 	lost ignore norm no # atom number non-constant
dump myDump all atom 10000 n_0p63_periodic.atom






run  		200000 # final longer run

unfix 11

fix 111 all viscous 0

run 	        0

fix 1111 all deform 1 x vel -0.01 y vel -0.01

run 		400000

unfix 1111

run 2000000









write_data      n_0p63_periodic.sim

srtee · February 7, 2023, 1:30am

Your simulation sounds very complex and combines many non-standard features:

Box aspect ratio very far from 1
Unusual pair style (self-coded?)
Non-constant box size
Viscous dynamics
On particles with internal DoFs

Any one of the above factors (or more than one!) could be causing you problems. In order:

your box aspect ratio may be exciting unnatural phonons because of the high periodicity in x and y
your pair style may be buggy or not behave properly in some regimes (after all even the LJ potential will explode if atoms come too close)
your changing box size might be setting off shock waves (or might be interacting with subtle bugs in your pair style’s handling of ghost atoms)
viscous damping might not be correctly set up (and if it is, it might be correctly putting energy into your system, instead of removing energy, as it sometimes would if connected to a heat bath of positive temperature)
you may not be measuring temperature correctly from your particles’ internal DoFs.

And this is not an exhaustive list of possible problems. I am not giving you this list as a checklist. You will probably get stuck if you focus on investigating or eliminating the specific issues I’ve mentioned, in the same way that you will probably get stuck if you focus on average neighbors per atom (which, for all we know, might not be the issue, or might just be the symptom rather than the source of problems – here you are falling prey to the “streetlight effect”).

Rather, I am trying to emphasize that there is a lot going on in your simulation, and the surest way to proceed is to simplify. In particular, you need to eliminate as many unknowns as possible. For example, I don’t know what viscous damping is supposed to do in your system, but I do know that many MD simulations use energy-conserving dynamics, so if you can pare back your simulation to a point where energy is supposed to be conserved and energy is conserved then you have something you can be confident does work.

You can then add complexities back in one by one – and when you see that something does break you will be much more confident that you should focus on that particular aspect.

Also, because of everything piled on in your simulation, you are at the point where frankly you know your system better than I or anyone else on this forum does. You are far out of the realm of bulk SPC/E water or OPLS hydrocarbons or CHARMM proteins or anything else where a non-expert would be able to make useful specific comments. That is another motivation to pare things back.

Finally, when coding up a custom MD method, it is pivotal for you to have a base case using standard methods, so that you can point to the differences between the base case and your custom method and say something like “See, this custom method correctly predicts the radius of gyration of the average broppling widget, which is pivotal to the correct refobulation of broppled widgets in the bropple industry”. Have you got your base case set up? That can be an invaluable source of “ground truth” against which you can debug your new method.

voro · February 7, 2023, 2:18am

Thanks for your response @srtee. I am using the default gran/hertz/history in LAMMPS and not my custom potential (yet!). I see what you are saying and have tried without viscous damping where only friction is the dissipative source but I have gotten the same result.

I was worried about shockwaves too but my velocities are relatively small (but not as small for it to be considered stationary) with respect to the internal pressure… Actually, my goal is to induce a shock after confinement hence the long shock tube like geometry. I will try your suggestions for an aspect ratio that is more cubic than the current geometry and see.

Is it possible for you to comment on the bin size question I had in my original discussion as I am unsure about the neighbor command? Will the same pair wise neighbor contributions to the force calculation show up with the neighbor bin command as would if I didn’t use it? I see very small differences in results (for example using kinetic energy as a metric) between the two cases but I would like an expert opinion.

For granular systems, the cutoff distance is the sum of the radii so even if the neighbor bin command generates more neighbors, I expect for both to give me the same results if I generate a neighbor list frequently.

Many thanks!

srtee · February 7, 2023, 3:04am

The documentation for neighbor is here: neighbor command — LAMMPS documentation with details of the algorithm in the LAMMPS paper (https://www.sciencedirect.com/science/article/pii/S0010465521002836?via%3Dihub), section 3.3.

I do not know what is the correct number of average neighbors each particle should have. You have the data you need to calculate it yourself, either by postprocessing a snapshot and counting how many particles are within the (cutoff + skin distance) bins of each particle, or on pen and paper using the preset density of the sample.

Note that pair cutoffs are applied after neighbor binning. Hence, tests of the sort if r_ij < r_cutoff still appear in the source code for pair force calculations, because not all neighbor pairs are within the pair cutoff (which is how MD codes get away with less frequent neighbor list updates).

The most straightforward way for neighboring to fail is if a particle crosses from outside the skin to inside the pair cutoff before the neighbor list is updated, leading to some force interactions being improperly excluded.

Since you are able to run scripts with various settings for neigh_modify, you are in a position to observe and verify for yourself if that happens or doesn’t happen on various settings. (You may find the Anna Karenina principle helpful – all correct settings of neigh_modify will give correct results alike, while most wrong settings of neigh_modify will give results that are wrong in their own individual way.)

I am not trying to be difficult, but I do not know what you will observe when you run LAMMPS and I do not want to give you false information that would lead you astray.

Any fully periodic system has zero confinement.

voro · February 8, 2023, 12:28am

Thanks for the info @srtee. Yes, I will do some postprocessing to verify neighbors.

I had a question about this statement and its validity to a granular system. By zero confinement do you mean 0 pressure? I would expect there to be a non-zero pressure if compressed beyond the critical jamming volume fraction because granules are not lost and would be contained in a smaller volume, and are thus bound to interact.

srtee · February 8, 2023, 1:02am

This is far out of topic for the LAMMPS forum, so you will have to do your own reading and thinking to consider whether I am right about this (or wrong! I am just an internet stranger, not your supervisor or colleague or mentor):

Consider one cubic centimeter of water in a 1cm by 1cm by 1cm box. Ordinary hydrostatics tells us that its pressure equalizes and it exerts as much pressure on the sides of the box as it does on the bottom, no more, no less (the amount of water is small enough that we can neglect barometric differences between the top and bottom of the water).

Now let’s say I divide this water up into some incredibly large number of 1nm by 1nm by 1nm imaginary cubes. Again, ordinary hydrostatics says the pressure exerted by each imaginary block of water on each of its imaginary boundaries is exactly the same (this time, up to thermal fluctuations, since each block of water is now small enough that we will start noticing these things).

Now let’s say I divide this water up, again, into some even larger number of 0.1nm by 0.1nm by 1nm imaginary long prisms. Once again, ordinary hydrostatics says the pressure exerted by each imaginary block of water on each of its imaginary boundaries is exactly the same. (Otherwise, one can imagine constructing a perpetual motion machine using flexible nanofluidic cells of changeable aspect ratio …)

So how can the aspect ratio of a simulation box, in a fully periodic simulation, ever change the isotropic pressure one would expect from hydrostatic equilibrium – unless the periodic repetition is now so severe that the simulation is no longer physically reliable for representing real bulk fluids?

I think the answer is it can’t. But again, you have all the tools you need to work this out for yourself. You can actually measure the pressure in a molecular simulation in different directions. You can work out for yourself whether you are seeing any anisotropy, and whether it is real or an artifact (for example, whether it persists if you double the box cross-section or halve it). As a reviewer I would absolutely demand these kinds of tests if someone wrote a paper depending on artificially short periodicity to simulate nanoconfinement, especially since it is not difficult to make an actually nano-confined simulation system with walls (which is my simulation bread-and-butter).