Randomly select a certain number of atoms (more than 1)

LAMMPS LAMMPS General Discussion
1

Dear Lammps users,

I am trying to select a certain number of atoms as PKAs using the following commands:
group select empty
group PKAstore empty
variable numPKA equal count(PKAstore)
variable n loop 500
variable check atom grmask(Si,FREE)
label jump
variable select equal (random(1,atoms,4224):%.0f) variable a equal (v_check[v_select])

if “v_check[v_select] == 1” then &
“group PKAstore id {select}" if "{numPKA}<10” then &
“next n” &
“jump SELF jump”

variable n delete

However, I ran into the following error:
ERROR: Indexed per-atom vector in variable formula without atom map. (…/variable.cpp:4604)
Last command: variable a equal $(v_check[v_select])

I have looked up LAMMPS manual and added “atom_modify map yes id yes”in my input, yet the error still exists. I don’t know how to fix this atom mapping problem, and if there is a better way to select a certain number of atoms randomly within a specified group and region?

Many thanks for any suggestion!

Best wishes,
Dingyu

The whole in script is as follows:

variable T1 equal 300.0 # K
variable RDN equal 4928459
variable PROD equal 10000 #100ps test
variable dt1 equal 0.1
variable P1 equal 1.0
variable Tpka equal 55830.0

units real
boundary p p p
atom_style charge
atom_modify map yes id yes
read_restart PHPS50_slab_Eq300K.restart
pair_style reaxff lmp_control safezone 32 mincap 2000 minhbonds 1000
pair_coeff * * ffield N Si H

fix 1 all qeq/reaxff 1 0.00 10.00 1e-6 reaxff
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
balance 1.1 shift z 10 1.05
timestep ${dt1}
reset_timestep 0
region BOTTOM block INF INF INF INF INF -16.0 units box
region FREE block INF INF INF INF -16.0 23.0 units box
region SURF block INF INF INF INF 10.0 23.0 units box
group BOTTOM region BOTTOM
group FREE dynamic all region FREE
group SURF region SURF
group FREEstatic region FREE

group N type 1
group Si type 2
group H type 3
group FREESi intersect Si FREEstatic
group SURFSi intersect Si SURF

#randomly select a Si PKA
group select empty
group PKAstore empty
variable numPKA equal count(PKAstore)
variable n loop 500
variable check atom grmask(Si,FREE)
label jump
variable select equal (random(1,atoms,4224):%.0f) variable a equal (v_check[v_select])

if “v_check[v_select] == 1” then &
“group PKAstore id {select}" if "{numPKA}<10” then &
“next n” &
“jump SELF jump”

variable n delete

velocity BOTTOM set 0.0 0.0 0.0
fix FREEZE BOTTOM setforce 0.0 0.0 0.0

fix spec all reaxff/species 1 1 100 species_0_1000_ps.out element N Si H &
delete species.del masslimit 0 47
fix bond all reaxff/bonds 1000 bonds_0_1000_ps.reax

compute mytemp FREE temp
compute_modify mytemp dynamic/dof yes

thermo_style custom step time atoms temp c_mytemp pe ke evdwl ecoul etotal press xlo xhi ylo yhi zlo zhi cpuremain
thermo 1000
dump 1 all custom 10000 PHPS50_slab_rad.lammpstrj id type q x y z vx vy vz
dump_modify 1 sort id

#setup PKA
velocity select set 0.3388362 0.3388362 0.3388362 units box
fix 2 all nve
fix 3 all dt/reset 1 1.0E-3 1.0E-1 0.2 units box

restart 1000000 PHPS50_slab_rad.restart
run ${PROD} # 100ps

2

Hi @hdy2022,
Some comments here:

  1. Please help people read your post by enclosing your script extract between three backticks characters ``` (do not use apostrophes '). See more about formatting your post here. You can read the other useful general advices there.
  2. Your script is long and complicated. Try to restrict it to a simple case if you want meaningful help from other users.
  3. The error gives you the line at which the code fails: variable a equal $(v_check[v_select]) But this line does not appear elsewhere in your script extract and full script. Is this normal?
  4. I am not sure the command random(1, atoms, 4224) is legal. I think you should give it the variable natoms instead of atoms. Not sure this is related to your problem though.
3

Hi Germain,

Thank you very much for your time and reply. Sorry I was not being clear about my problem.
I tried to restrict my input script to a simple case. I used the following commands try to randomly select 10 Si atoms (Si atoms are in group Si) in the FREE region.

group           select empty
group           PKAstore empty
variable        numPKA equal count(PKAstore)
variable        n loop 500
variable        check atom grmask(Si,FREE)
label           jump
variable        select equal $(random(1,atoms,4224):%.0f)
variable        a equal $(v_check[v_select])

if "$a == 1" then &
    "group PKAstore id ${select}"
	if "${numPKA}<10" then &
	  "next n" &
	  "jump SELF jump"
	  
variable n delete

An atom style variable “check” is defined. if an atom is Si and it is in the FREE region, variable “check” should be 1. Then this atom is added to the group “PKAstore”. Keep selecting Si atoms until there is 10 atoms in the group “PKAstore”.

Running this script lead to the following error:

ERROR: Indexed per-atom vector in variable formula without atom map (../variable.cpp:4604)
Last command: variable        a equal $(v_check[v_select])

Thanks for pointing out. The command random(1, atoms, 4224) is legal in my case.

The entire input script is as follows:

variable        Tpka    equal 55830.0
units           real
boundary        p p p
atom_style      charge
atom_modify     map yes id yes
read_restart    PHPS50_slab_Eq300K.restart
pair_style      reaxff lmp_control safezone 32 mincap 2000 minhbonds 1000 
pair_coeff      * * ffield N Si H

fix             1 all qeq/reaxff 1 0.00 10.00 1e-6 reaxff
neighbor        2.0  bin
neigh_modify    every 10 delay 0 check yes
balance         1.1 shift z 10 1.05

timestep        0.1
reset_timestep  0

region          FREE   block INF INF INF INF -16.0 23.0 units box
group           Si type 2

#randomly select Si PKA
group           select empty
group           PKAstore empty
variable        numPKA equal count(PKAstore)
variable        n loop 500
variable        check atom grmask(Si,FREE)
label           jump
variable        select equal $(random(1,atoms,4224):%.0f)
variable        a equal $(v_check[v_select])

if "$a == 1" then &
    "group PKAstore id ${select}"
	if "${numPKA}<10" then &
	  "next n" &
	  "jump SELF jump"
	  
variable n delete

thermo_style    custom step time atoms temp pe ke evdwl ecoul etotal press
thermo          1000
dump            1 all custom 10000 PHPS50_slab_rad.lammpstrj id type q x y z vx vy vz
dump_modify     1 sort id

velocity        PKAstore scale ${Tpka}
fix             2 all nve
run             10000

Many thanks again!

Best wishes,
Dingyu

4

What is your LAMMPS version? Can you try with the most recent version?
Perhaps you are running into a bug that has already been fixed, since I cannot reproduce the error message with your atom_modify settings on the current development version.

5

Thank you very much! My version is 15 Jun 2023.

6 
atom_modify map yes
region box block 0 10 0 10 0 10
create_box 2 box
lattice fcc 0.8
create_atoms 1 box
set type 1 type/fraction 2 0.5 23435234
region          FREE   block INF INF INF INF -16.0 23.0 units box
group Si type 2
group select empty
variable check atom grmask(Si,FREE)
variable select equal $(random(1,atoms,4224):%.0f)
variable a equal $(v_check[v_select])

print "select = ${select}  a = ${a}"

Here is the most minimal input I could distill from your post and that only produces the error you report for me, if I comment out the atom_modify command. I tried with the 15 Jun 2023 version and the behavior is the same. Perhaps if you try atom_modify map array?

If this input also works for you, try to experiment with it and move over pieces from your input and see if you can find the minimum set of changes that takes this from no error to creating the error.

These kinds of issues can only be addressed if they can be reproduced…

7

Thank you very much! I tried atom_modify map array but I still got the error. It seems like atom_modify map array does not work before read_restart PHPS50_slab_Eq300K.restart command. Here is the output from log.lammps:

atom_style      charge
atom_modify     map array
read_restart    PHPS50_slab_Eq300K.restart
Reading restart file ...
  restart file = 15 Jun 2023, LAMMPS = 15 Jun 2023
WARNING: Restart file used different # of processors: 112 vs. 128 (../read_restart.cpp:661)
  restoring atom style charge from restart
  orthogonal box = (-20.950562 -20.950562 -25.950562) to (20.950562 20.950562 105)
  4 by 4 by 8 MPI processor grid

I cannot upload files here. Is there any way that I can send all the files to you?
Many thanks!

8

There is no need (and I am not interested in debugging files that need to run on over a hundred processors. that is just a massive waste of time).

For the time being, you should be able to work around this by converting your restart file to a data file and then use read_data.

9

You can try making the following modification to your source and recompile.
After this, you should be able to use atom_modify map before a read_restart and its setting will persist.

  Author: Axel Kohlmeyer <[email protected]>
  Date:   Wed Aug 30 23:51:25 2023 -0400
  
      make atom_modify map settings in restart file overridable
  
  diff --git a/src/read_restart.cpp b/src/read_restart.cpp
  index 0de36ebe1e..cb5754909b 100644
  --- a/src/read_restart.cpp
  +++ b/src/read_restart.cpp
  @@ -831,9 +831,12 @@ void ReadRestart::header()
       } else if (flag == ATOM_ID) {
         atom->tag_enable = read_int();
       } else if (flag == ATOM_MAP_STYLE) {
  -      atom->map_style = read_int();
  +      // we should be able to enable an atom map, even
  +      // if the simulation in the restart didn't use one
  +      int itmp = read_int();
  +      if (atom->map_user == Atom::MAP_NONE) atom->map_style = itmp;
       } else if (flag == ATOM_MAP_USER) {
  -      atom->map_user  = read_int();
  +      int itmp = read_int();  // ignored
       } else if (flag == ATOM_SORTFREQ) {
         atom->sortfreq = read_int();
       } else if (flag == ATOM_SORTBIN) {
10

Thank you very much! I will first try read_data see if it works. If needed, I will modify source and recompile Lammps.

11

Thank you very much! I convert restart file to a data file and it works fine.