Pick using id a random atom within a group

Dear All,

I need to pick an atom (atom id) randomly from the surface of a sphere which I have created as a group. Using the math_function random available in LAMMPS, I can only pick up a random id within the group all. I tried various approaches, but the only one that I could come up with was a “rejection method” using the if command and looping over all atoms till I chance upon an atom within the region of interest. The problem with this is, if the region of interest holds a very few atoms compared to the simulation box, a lot of time will be spent to satisfy the criteria.

If I need to write code to do this, what would be the best approach? My inclination would be to extend the variable command to handle something like:
variable myRandomAtm equal rand_atm(group)

I look forward to your suggestions.
Thanks in advance.

Manoj

Manoj Warrier
Bhabha Atomic Research Center,
Vizag, India

Hi Manoj,
May be my suggestions will be off because I misunderstood your system geometry, but wouldn’t it be possible to :
- option 1 : attribute from the start a different Id to the atoms of the surface (probably work if your system is a solid)
- option 2 : detect continuously which atoms are at the surface using regions (if you system is a fluid)
- option 3 : detect and remove atoms using an external home-made script (e.g. Python + MDAnalysis)
Best,
Simon

It think adding a “random_select_group()” and “random_select_region()” function to the variable command might be a useful addition.

There could be a workaround for that, assuming that you have no other use for the molecule ids.
You initially set all molecule ids to 0. Then you apply the reset_mol_ids single yes command on the group with the atoms of interest. Assuming there are no bonds, this will give those atoms a molecule ID in the range of 1 to N. Then you get a random number from 1-N and assign it to a variable, say selectid, and define an atom style variable based on mol==v_selectid. This will contain all zeros but for the atom with the selected molecule ID. This can be used to define a group. If
you need to access the value of that atoms atomID, you can use compute reduce max replace on the atom ids and this atom style variable.

Hi Simon.

Thank you for your time. A brief description of the system - It is a cubical solid crystal. I define two spheres using the region command and subtract (using the group command) the inner sphere from the outer thereby creating a group of atoms on the surface of a hollow sphere. Therefore I have a group of atoms, identified their IDs. Now I would like to carry out a series of simulations where I choose atoms randomly from this group. I do not know how to randomly pick any of the atoms belonging to the specific group within the LAMMPS script and need help/suggestions to do that. Regarding your options:
option 1: Even if I attribute a different ID and type to the atoms in the region of interest, how do I randomly choose from them within the LAMMPS script?
Option 2: Detected, how do I randomly choose from them within the LAMMPS script?
Option 3: Prefer adding functionality to LAMMPS … As a quick fix, yes. We are working on pre-processing this.

The simplest way for your process (in general, not just because there is no built in method). Would be to write a dump of just the atoms ID of the group, load that into a python script as a list, use the python shuffle function and write out those numbers to a file, one per line.
Then you can use a file style variable to read those randomized atom ids one at a time (just use the “next” command to get the next one) in your input.

Thank you. Proceeding with this approach seems the quickest way to go ahead with the simulations. I was trying to get my head around your earlier reply.

Regards
Manoj

#start with approximate center - cx,cy,cz, then find region with sphere at this center. this method reduces time needed to find a random atom

label loop
variable a loop 10
variable radius equal 0.1*$a
region mysph sphere (…syntax using cx,cy,cz, radius) units box
group sphere region mysph

variable nsph equal count(sphere)

if “${nsph} > 0” then “jump SELF break”

variable radius delete
region mysph delete
group sphere delete

next a
jump SELF loop
label break