How to count particles with zero velocity


I model the flow of small granular particles through the packing of large granular particles that are freezed, and some small particles become retained between large particles, i.e. their velocity becomes equal to zero. So I want to get the number of the particles with zero velocity every few timesteps. It would be perfect to get this number in thermo output but I didn’t find the way to do it, please advice is ther some easy way to count them?

Currently I count them in a loop through all particles with +1 counter if velocity is less than 1.0^-10, like that (this is for vz velocity):

variable retainedAtoms equal 0
variable Natoms equal count(particlesGroup)
compute cvAtom particlesGroup property/atom vz
compute cvAtom_red particlesGroup reduce ave c_cvAtom
thermo_style custom step … c_cvAtom_red <- I had to add this to thermo output to get the “current” compute state
thermo 50000

run 0
label runLoopLabel
variable retainedAtomsOnStep equal 0
label startLoopLabel
variable n loop ${Natoms}
variable atomID equal n variable vAtom equal c_cvAtom[{atomID}]
print ‘Vz speed of atom n: {vAtom}’
if ‘{vAtom}<=1.0e-10' then 'variable retainedAtomsOnStep equal {retainedAtomsOnStep}+1’
next n
jump SELF startLoopLabel
print ‘Retained atoms on step: {retainedAtomsOnStep}' if '{retainedAtomsOnStep} > {retainedAtoms}' then 'variable retainedAtoms equal {retainedAtomsOnStep}’
print ‘Retained atoms: ${retainedAtoms}’
run 500000
jump SELF runLoopLabel

Another way that I can imagine is to get this information from dump file, but is there a way to do it in thermo output?

Nikita Tropin | Scientist

Take a look at the "math operators" section of the variable doc page.
If you write a atom-style variable formula with something like
(vx > -0.01) && (vx < 0.01) it will return a 1 if the Vx is near 0.0,
else a 0. If you then sum that variable over all atoms via compute
reduce, you would get an effective count of atoms with (near) 0.0
velocity, which you could print out with thermo stats. You could
replace vx with v_myV, where myV computes a scalar velocity from
all 3 components.