Counting an atom

LAMMPS LAMMPS Beginners
1

Hi everybody,

I need to count oxygen atoms passing through a virtual sheet in a water channel. As I understood there may be possible to use compute reduce/region command. The reason that I told a virtual sheet is to recognize passing forward or backward atoms. Anyway, could you please write any suggestion on using compute reduce/region command for computing this atom type?
Thanks

2

There are different ways in LAMMPS to count atoms. If you want to use compute reduce/region you would need to use the “sum” keyword. There is no “count” option, but you can easily just sum a property that is the same for all atoms (e.g. define a group with oxygen atoms, sum the mass and then divide by the oxygen mass to get the count).
A different approach would be to use atom style variables and code the conditions to count inside the variable since those will turn logical expressions into 1/0 numbers.

If you want to differentiate between forward and backward moving atoms, then the situation is a bit more complex. For the human eye it is rather straightforward to decide when a particle is moving through a barrier, but that is because the brain will not only look at the current state but rather interpolates and extrapolates the motion of particles. When you want to do the same analysis from a MD simulation software, it gets a bit trickier since you usually only have the information about the current state.

Thus it is possible to tell if a particle is inside a given region but that may not be exactly how you would count when viewing:

  • due to discretization you can only “see” the system at specific discrete points in time, so you cannot exactly monitor the point where a particle is passing through a specific plane
  • due to (fast) local vibrations/librations you cannot use the current velocity as a reliable indicator of direction in condensed fluids
  • since particles are moving at different velocities you may count slow moving particles more often than faster moving ones
  • if the region is too thin, you may be missing fast moving particles; if it is too thick particles get counted multiple times.

So the best way to count atoms passing through is actually in post processing where you can write a custom script that stores some extra properties with atoms that can help. E.g.

  • a flag which keeps of whether an atom has been counted as passing through the virtual barrier already in one direction and which will be reset if it is moving the other direction
  • a counter where you keep track when the “passing through” state of an atom has last been changed

Then you can rather reliably count how many atoms are moving forward or backward by counting only atoms that have changed their state and kept it for a minimum number of steps.

For analysis during the simulation, you are limited in what you can do as explained above. One possible approach would be to not use the current velocity as a direction indicator but the average velocity over the last N steps. Example (to be added to the “melt” example):

variable outfreq index 50
fix vel all ave/atom 1 10 ${outfreq} vx
variable fwd atom (f_vel>0.0)*(x>8.0)*(x<8.5)
variable rev atom (f_vel<0.0)*(x>8.0)*(x<8.5)
compute count all reduce sum v_fwd v_rev
fix count all print ${outfreq} "Forward: $(c_count[1])  Backwards: $(c_count[2])"

The atom style variables fwd and rev will be 1 if all 3 conditions are true (average velocity over last 10 step positive/negative and within 8.0 and 8.5 sigma in x coordinate) and compute reduce will thus count atoms conforming to those conditions.

1 Like
3

I think Axel’s solution above is a good one. It is effectively looking at the time integral of velocity of particles that are now in the control region, labelling them as left-moving and right moving. A similar result could be obtained by setting up two regions “red” and “green” separated by the virtual sheet, possibly with a gap in between. Atoms are dynamically assigned color or type according to which region they are in. Green atoms in the red region are switched to red and vice versa. Every time a green atom in the red region is recolored red, that would count as one atom moving from red to green. There are multiple ways to program this using LAMMPS commands e.g. using the group command with dynamic keyword.

1 Like
4

Great thanks everybody for your fruitful replies.

Regarding @athomps comment, I read some papers and authors exactly followed your procedure to calculate spatial dependent physical quantities. Hereby, it is really interesting for me how do you establish your explained procedure. Would you mind letting me some example or guidelines to follow this?

I’m eagerly looking forward to hearing from you