I am simulating the dynamic evolution of a slab under radiation using ReaxFF. I want to delete the gas-phase species such as NH3, H2 produced during the radiation. By using the fix reaxff/species command in LAMMPS, I can delete all of the NH3 and H2 molecules. Is there any way that I only delete the gas-phase species when they have moved into the upper region of the simulation box (above the slab)? Currently with the fix reaxff/species delete species.del specieslist 2 H2 NH3 command, I noticed many NH3 molecules are deleted when they are still inside the slab, which is not reasonable. Do I have to modify the source code to realize this function?
Probably not. All fix commands require to specify a group ID. Often people use “all”, but you can use a different group and that usually means, that the fix command in question is only applied to that specific group of atoms. If the fix command supports it, that group may also be a dynamic group based on a region. A quick look at the source code of fix reaxff/species indicates that both prerequisites should be fulfilled and thus it should be possible to use it as suggested. For details, you need to have a closer look at the documentation for the group and the fix command.
Thank you very much! I defined a dynamic group UPPER for atoms in the upper region of the simulation box, and applied fix UPPER reaxff/species delete specieslist 2 H2 NH3 command only on the UPPER group instead of all. Then only H2 and NH3 in the upper region are deleted. However, for species that across the boundary of the UPPER region, for example, 2 H atoms of one NH3 molecule are in the UPPER region while the remaining N and H are not in the UPPER region. We want to delete the entire NH3 molecule. However, I suppose the command will delete the 2 H atoms in the UPPER region because they may be recognized as a H2 molecule. To avoid this problem, LAMMPS needs to know the molecule id for each species. Do you have any suggestion to deal with the deletion of species across the region boundary?
Please have a closer look at the group command. It can augment the group by including all atoms with the same molecule ID. For that you for that you could use fix reaxff/species as well, but you need to apply it to group all. So you would need to use two copies, one to get (and then assign) molecule IDs and then update the dynamic group for all included molecules and then use a second fix reaxff/species to delete the molecules that are at least partially inside the “to be deleted” region.
Thank you very much for your suggestion! However, I am not sure how to assign molecule IDs by using fix reaxff/species all command. I ran into the following error:
fix spec all reaxff/species 1 1 100 species.out element N Si H
group UPPER dynamic all region UPPER
dynamic group UPPER defined
group UPPER include molecule
ERROR: Group include molecule command requires atom attribute molecule (../group.cpp:373)
Last command: group UPPER include molecule
That is easy to solve. Please think about it for a bit, do some searching. This forum is no replacement for you doing the thinking and researching, so please do not post immediately if you run into an error and wait for somebody to do your job for you.
I used the atom_style full to run the simulation so it now has the mol_id property. Then I used the following command to access and assign molecule id determined by fix reaxff/species to each atom:
fix spec all reaxff/species 1 1 100 species.out element N Si H
variable molid atom f_spec[3]
set group all mol v_molid
However, I get the error:
ERROR on proc 8: Variable molid: Variable formula fix vector is accessed out-of-range (../variable.cpp:1769)
Last command: set group all mol v_molid
As I understand, fix reaxff/species output one global vector of length 2 and a per-atom vector. I have tested that f_spec[1] is the total number of molecules and f_spec[2] is the total number of distinct species. So I suppose f_spec[3] is the per-atom vector molecule ID. Please correct me if i am wrong.
Thank you very much for your reply. In fact I have tried variable molid atom f_spec, however, the following error message appears:
fix spec all reaxff/species 1 1 100 species.out element N Si H
variable molid atom f_spec
set group all mol v_molid
Setting atom values ...
ERROR: Variable molid: Fix global vector in atom-style variable formula (../variable.cpp:1810)
Last command: set group all mol v_molid
With this error message, I thought maybe I should use f_spec[3]. It seems LAMMPS takes f_spec as a global vector by default. I am very confused about how to use the per-atom vector calculated by fix reaxff/species?
Dear akohlmey,
Thank you very much! That does solve my problem. Another question I would like to ask is that besides run every, are there any other approach can set the molecule ID based on fix reaxff/species at each time step? Currently I am using the run every command as following. I delete species every 100 time steps. However, if I want to delete species at a lower frequency, i.e., every 1000 time steps, then fix specDel UPPER reaxff/species 1 1 1000 won’t work due to the run every 100 command.
fix spec all reaxff/species 1 1 100 species_all.out element N Si H &
position 1000 NSiH.pos
fix cache all ave/atom 1 1 100 f_spec
variable molid atom f_cache
group UPPER include molecule
fix specDel UPPER reaxff/species 1 1 100 species_upper.out element N Si H &
delete species_upper.del masslimit 0 48
run 1000000 every 100 "set group all mol v_molid"
Then do a run xxx every 1000. There is no point in updating the molecule IDs more often than you need them to be updated. In fact, all the other “every” commands can be changed to be called only every 1000 steps. That is just plain common sense.
Many thanks for your reply! I was thinking that initially all mol ID = 0. Then I first set mol ID at time step 1000. At time step 2000, some species will be deleted. I suppose when determining which species are in the group to be deleted, the mol ID at time step 1000 is taken. That’s the reason why I would like to set mol id more frequently but perform the deletion at a lower frequency (because mol IDs are very likely to be the same in 10 time steps). Please correct me if I am wrong.
Whether the values of the molecule IDs change is irrelevant in this case, what matters is that the same group of atoms will have the same molecule ID, and that is highly likely for the situation you describe.
Other than that, I can think of some rather convoluted way to reset the molecule IDs during a run without even using run every by using fix python/invoke and have it set up so it updates the molecule IDs right before they are used through the LAMMPS python interface.
