Hello,
There are a couple old emails in the archive asking about how to specify a random molecule from a group. However, all of the answers only work if the IDs are sequentially numbered. I have a system where I am assigning a group in a defined geometric region, meaning that the IDs are not sequential.
I am wondering if there is a work-around that will allow me to select a random molecule within my pre-defined group.
Thanks,
Hello,
There are a couple old emails in the archive asking about how to specify a
random molecule from a group. However, all of the answers only work if the
IDs are sequentially numbered. I have a system where I am assigning a group
in a defined geometric region, meaning that the IDs are not sequential.
you should be able to convert the non-sequential numbering into a
sequential numbering via compute chunk/atom using the "compress yes"
option, and then take it from there.
axel.
Hi Axel,
Thanks for your reply. I followed your advice but it seems that when I try to use the random command to find a number within the range of the IDs in my group, the original IDs are used instead. The commands I am using are:
variable xpos equal lx
variable xpos2 equal lx/2
variable ypos equal ly
variable ypos2 equal ly/2
variable zpos equal lz
variable zpos2 equal lz/2
region rRad cylinder z {xpos2} {ypos2} 10.0 {zpos2} {zpos}
group gRad region rRad
group gRad include molecule
variable naRad equal count(gRad)
variable nmRad equal v_naRad/3
#----- Select a random atom ----
compute chRad all chunk/atom molecule region rRad nchunk once ids once compress yes
variable randomN equal round(random(1,{nmRad},2987659))
variable nPKA equal {randomN}
variable nPKA equal {randomN}
group gPKA molecule {nPKA}
variable xpos equal xcm(gPKA,x)
variable ypos equal xcm(gPKA,y)
variable zpos equal xcm(gPKA,z)
print ‘The PKA position is (x,y,z) {xpos} {ypos} ${zpos}’
The corresponding output for the region command is:
region rRad cylinder z 22.4951839757269 23.3776637294881 10.0 22.0406766029625 44.0813532059251
The corresponding output for the last line is:
The PKA position is (x,y,z) 26.7948178203765 2.46710537966203 26.1828706773078
As you can see, the y position is outside of the predefined ‘rRad’ region above. I guess I don’t understand how ‘compress yes’ is reseting the IDs of my group without conflicting with the IDs of the remaining particles. Any additional insight you could provide would be appreciated.
Best,
Micah
Hi all,
Sorry to pester about this. Maybe there is a better way to accomplish my goals than to convert a chunk ID to a group ID so that I can set the velocity of the group. Basically, I am trying to determine the effect of exciting an atom in a water ice solid, and I wanted to see how the result changed as a function of distance from a free surface. I could go find the specific molecule ID from the dump files, but I thought that selecting a random molecule would be more efficient. If this would require changes to the lammps code then it probably isn’t worth it.
However, if there is a relatively simple way to define a region in a solid, randomly select a molecule in that region, and then assign that molecule a given velocity, I would really appreciate any suggestions on how to accomplish it.
Thanks
I am still trying to find a solution that will allow me to select a random molecule within a specified group/region and assign this molecule a velocity. Working from Axel’s previous solution, I used
compute chRad all chunk/atom molecule region rRad nchunk once ids once compress yes
which gives me a set of sequentially numbered chunk IDs for the molecules contained in the region ‘rRad’. I can also find the number of molecules in this region and, using the random number generator, specify a random chunk ID and set this equal to a per-atom variable.
# – Define a region in which to excite a molecule –
region rRad cylinder z lx/2 ly/2 10.0 lz/2 lz
group gRad region rRad
group gRad include molecule
variable naRad equal count(gRad)
variable nmRad equal v_naRad/3
variable randomN equal round(random(1,${nmRad},2987659))
_variable nPKA equal {randomN}**_ <i>**variable chunkID atom c_chRad[{nPKA}]
However, when I try to assign this atom-style variable to a group
group gPKA variable chunkID
I receive the error “ERROR: Compute used in variable between runs is not current”. When I searched for an explanation for this error I learned that “There is no mechanism in LAMMPS to use a per-atom compute variable to create a group.” See: http://lammps.sandia.gov/threads/msg24099.html
So, I am still left with the problem of assigning a random molecule to a velocity. If anyone has suggestions I would appreciate them.
Thanks,
Micah
I am still trying to find a solution that will allow me to select a random
molecule within a specified group/region and assign this molecule a
velocity. Working from Axel's previous solution, I usedcompute chRad all chunk/atom molecule region rRad nchunk once ids once
compress yeswhich gives me a set of sequentially numbered chunk IDs for the molecules
contained in the region 'rRad'. I can also find the number of molecules in
this region and, using the random number generator, specify a random chunk
ID and set this equal to a per-atom variable.# -- Define a region in which to excite a molecule --
region rRad cylinder z lx/2 ly/2 10.0 lz/2 lz
group gRad region rRad
group gRad include molecule
variable naRad equal count(gRad)
variable nmRad equal v_naRad/3
variable randomN equal round(random(1,\{nmRad\},2987659\)\) variable nPKA equal {randomN}
variable chunkID atom c_chRad[${nPKA}]
this all looks needlessly convoluted and is in part outright wrong
(definition of region rRad).
However, when I try to assign this atom-style variable to a group
group gPKA variable chunkID
I receive the error "ERROR: Compute used in variable between runs is not
current". When I searched for an explanation for this error I learned that
"There is no mechanism in LAMMPS to use a per-atom compute variable to
create a group." See: http://lammps.sandia.gov/threads/msg24099.htmlSo, I am still left with the problem of assigning a random molecule to a
velocity. If anyone has suggestions I would appreciate them.
this message is a several years old and the comment about dynamic
groups is no longer true. in fact, that is how you can make it work.
region rRad cylinder z \(lx/2\) (ly/2) 10.0 \(lz/2\) (lz)
group gRad region rRad
group gRad include molecule
compute cRad gRad chunk/atom molecule compress yes
variable nmRad equal count(gRad)/3
variable rMol equal round(random(1,\{nmRad\},2987659\)\)
variable pMol atom "c\_cRad==v\_rMol"
group gPKA dynamic all var pMol
print "group gPKA has (count(gPKA)) atoms now"
run 0 post no
print "group gPKA has $(count(gPKA)) atoms now"
velocity gPKA set NULL NULL 0.15
some notes:
- since you already narrowed down the group to the region, no need to
use the region as additional qualifier for chunk/atom. that will
exclude atoms from molecules that are partially inside the region and
thus will lead to incorrect velocity assignments
- the pMol variable is how you best select atoms with atom style variables
- the run 0 statement triggers the update of the dynamic group using
the atom-style variable based on comparing the chunk-id to the random
number
from here on you should be all set. if you want to select another
molecule, you better delete all used groups, computes and variables
and start over.
axel.