Stress/atom computatuion to pressure in granular molecules

Dear all,

I am simulating granular molecules being first poured into a periodic box, and then compressed. I am trying to extract the pressure on each particle from the stress/atom compute in the following way:

compute 1 all stress/atom NULL

variable p atom “c_1[1]+c_1[2]+c_1[3]”

I am using the gran/hooke/history pairwise potential with fix rigid.

When I try this compute on granular spheres that are not molecules, I get what I would expect: the pressure on each sphere is 0 as they are free falling, and then is greater than or equal to zero when settled at the bottom of the box.

However, when I run this compute on granular spheres that are molecules (so that I only have 1 sphere in the molecule file), or when I try it on granular molecules that have more constituent spheres, I get nonzero pressures as the molecules are falling, and then when they are settled the per-atom pressures are both large positive and large negative values, and frequently nearly symmetric about zero.

I have tried changing the “keywords” in the compute, including fix, pair, etc. but nothing seems to change this behavior. I have also tried manually calculating the inertia tensor in the molecule file, which does not change anything either.

I would greatly appreciate any help!

Many thanks,

Leah Roth

I don’t see why a granular sphere that is a single particle

would give different answers depending on whether you also

flag it as a molecule or not.

If you have molecules that are composed of highly overlapped

spherical granular particles then the forces between them

are huge. If you try to sum those for pairwise force or stress

you can get round-off errors. So I would try turning

off pairwise forces between particles in the same molecule.

The neigh_modify exclude command has an option for

doing that.


Dear Steve,

Thank you for your reply!

None of my spheres are overlapping in these simulations, but when I use

neigh_modify exclude molecule mol1

where mol1 is my group of spherical molecules with the following molecule file:

1 atoms
4.18888 mass
1 0.0 0.0 0.0
1 1
1 2.0

I am still getting pressures that range from about -1e5 to 1e6 or more, in cgs units. Are the units important, or am I missing a normalization factor? I say that because when I look at the pressure distribution in the packing, the top layers have the least pressure (which ends up being negative), and the pressure increases with depth, as expected. I can provide images of the pressure heat map in the packing if that would be helpful. It seems almost as if there is some sort of offset - are the periodic boundary conditions important?

Again, thank you for your help!



A quick followup, I have just run the simulations with both side walls and with periodic boundaries, and while the molecules are pouring and settling, there are negative pressures. When all of the molecules have settled, the negative pressures disappear. This appears to be a dynamics / timestep problem, perhaps solved by averaging pressures over many timesteps?

Many thanks,


In general pressures can be either negative or positive, and they

will certainly fluctuate as particles settle. Granular force

fields are also dissipative (not conservative of energy),

so pressures will dampen over time.