ERROR : compute chunk/atom molecule for non-molecular system

Hello everyone,

I am new to using the compute chunk/atom command in LAMMPS and would appreciate some guidance.

I am simulating the bombardment of Kapton with atomic oxygen (AO) using the ReaxFF potential to study how AO interacts with Kapton. In my setup, Kapton is positioned at the bottom of the simulation cell, and I have extended the z-dimension to allow the addition of atomic oxygen for bombardment.

During the simulation, I observe the formation of various molecules such as OH, H₂O, O₂, carbon chains, CO, CO₂, and others. I now want to identify and count all the molecules generated during the reaction using LAMMPS, rather than relying solely on visualization.

Based on my understanding of the compute chunk/atom documentation, this command can assign an ID to each molecule. To achieve my goal, I included the following lines in my script:

Temperature Profile Across z-Direction

compute Temp_cc1 all chunk/atom bin/1d z lower 2
compute tempChunk all temp/chunk Temp_cc1
fix avgTemp all ave/time 1000 1 1000 c_tempChunk file temp_profile.dat

Tracking Desorbed Gases

region region_desorbtion block INF INF INF INF 55 198 units box
group At_desorbed dynamic all region region_desorbtion every 500

compute mol_cc1 At_desorbed chunk/atom molecule
compute countMol At_desorbed count/chunk mol_cc1
compute idMol At_desorbed property/atom mol

fix saveMolCount all ave/time 1000 1 1000 c_countMol file desorption_gases.dat
fix saveMolID all ave/time 1000 1 1000 c_idMol file desorption_mol_IDs.dat

However, I encounter the following error:

ERROR: Compute chunk/atom molecule for non-molecular system (src/compute_chunk_atom.cpp:312)
Last command: compute mol_cc1 At_desorbed chunk/atom molecule

I am unsure if I am using compute chunk/atom molecule correctly. Could someone help me understand what I might be doing wrong? Any advice would be greatly appreciated!

Thank you!

Nope. It works the other way around, it uses the molecule ID information to group atoms into chunks.

You are not understanding what compute chunk/atom does and how this relates to your simulation.

If you ask compute chunk/atom to create chunks for “molecules” it does so using the molecule-IDs of individual atoms. That is an (arbitrary) parameter is can be applied to atoms for atom styles that have bonds, angles, dihedrals etc. Those static assignments. However, that is not what you are doing. In a ReaxFF model, there are no static molecules which is exactly why you are using it. Unless for molecular atom styles, the bond topology is recomputed dynamically in every step based on the computed bond order. This does not involve molecule IDs and if you are using atom style charge (which is typical for ReaxFF simulations), then there is no room to store that information in the atom style.

If you want to record which kinds of molecules are forming during your simulation, you need to do the same as ReaxFF does internally using the same kind of information. This is done, for example, by fix reaxff/species. It can also assign molecule ID values to those detected species and stores this in a per-atom vector and thus you can output them in a dump file.

1 Like

Thank you very much for your answer and for your time. It is running now.

I am still having an issue. I need the atom groups to be dynamic as the atoms contained in each region will change as the simulation evolves. But reaxff/species only works with static groups.

To overcome this issue I am changing the atom groups after each 10000 steps. Something tells me that this might not be a good idea; Could you please have the time to tell if it is a good idea?

(I erased somer $ characters so the text is readble in this forum)

------------------------ Desorbed Gases -------------------------

region region_desorbtion block INF INF INF INF 55 198 units box
group At_desorbed region region_desorbtion
fix 1 At_desorbed reaxff/species 500 2 1000 desorbed_gases.out

----------------- Kapton Composition Layer by Layer -----------------

variable z_min equal 0
variable z_max equal 45
variable step equal 9

label loop_layers
region region_layer_{z_min} block INF INF INF INF {z_min} {z_min + step} units box
group kapton_layer_{z_min} region region_layer_{z_min}
fix layer_{z_min} kapton_layer_{z_min} reaxff/species 500 2 1000 kapton_molecules_{z_min}_{z_min + step}.out
variable z_min equal {z_min} + {step} if "{z_min} < ${z_max}" then “jump SELF loop_layers”

------------------------ Bombardment + Chunk -------------------------

Adding AO

variable xlo equal 5.0
variable xhi equal 190
variable ylo equal 4
variable yhi equal 22
variable zhi equal 133
variable zlo equal 123
variable vxlo equal 0.0
variable vxhi equal 0.0
variable vylo equal 0.0
variable vyhi equal 0.0
variable vzhi equal -0.077
variable vzlo equal -0.07765721
region region-AO block {xlo} {xhi} {ylo} {yhi} {zlo} {zhi}
fix AO_added all deposit 800 3 50 2 region region-AO id max vx {vxlo} {vxhi} vy {vylo} {vyhi} vz {vzlo} {vzhi}

variable tot_cycles equal 200
variable cycle equal 0

label bombardment_loop
group At_desorbed region region_desorbtion

variable z_min equal 0
variable z_max equal 45
variable step equal 9

label loop_group
group kapton_layer_{z_min} region region_layer_{z_min} variable z_min equal {z_min} + {step} if "{z_min} < ${z_max}" then “jump SELF loop_group”

variable cycle equal {cycle} + 1 run 10000 write_data data.after_erosion if "{cycle} < ${tot_cycles}" then “jump SELF bombardment_loop”

It is still not very readable. Please see this post for explanations on how to correctly quote text in this forum. It has “please read this first” in the subject for a reason.

That is not correct. It does work with dynamic groups, but then you must not use averaging (i.e. nrepeat must be 1, not 2) as that would result in bogus averages if the group members change.

That is for you to figure out. But it doesn’t look like you are doing it correctly. I don’t see you deleting a group. If you don’t delete it first, the group command will just add members to the existing group, but not delete. This is documented behavior, of course.

1 Like

Again, thank you very much, it works if I remove the averaging.

But was that information written in the fix reaxff/species documentation or anywhere else (I am not talking about the group command, only the dynamic option)? I read the following documentation and I didn’t find where it was written that without the averaging the fix reaxff/species command works with a dynamic group. Is it somewhere else? I am asking so I can read that documentation too.

fix reaxff/species command — LAMMPS documentation

For these kind of things, I usually look it up straight in the source code since I know what to look for. I’ll add a note to the documentation page under “Restrictions”.

1 Like

Okey, thank you for your help