I am simulating how atomic oxygen reacts with kapton, for that I am using the reaxff potential. In my simulation box I have Kapton (my polymer) at the bottom of the cell (between 0-50 A in z direction) the supercell is empty between 50 and 400A in z direction. During the simulation I bombard Kapton by introducing atoms of oxygen at 5eV and direct them towards the kapton surface.
I wanna study the chemical modification of Kapton layer by layer. For that I am defining differnt regions with steps of 9 A between 0 and 45 A in z direction. Then I determine the moleuclar composition of each layer using fix reaxff/species. The problem is that when defining my regions I cut some of the molecules, so in the final file I donāt see the entire molecule but just the pieces that are inside each region. What I want to do is to include in each region the entire molecule without cutting it.
I know it is possible to do it in a classic system with chunk/atom and group dynamic all include molecule. But as I am using ReaxFF I donāt have that option (no molecude IDs). Someone knows how to solve that?
Here is how I define the different groups:
# 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_$(v_z_min) block INF INF INF INF ${z_min} $(v_z_min+ v_step) units box
group kapton_layer_$(v_z_min) dynamic all region region_layer_$(v_z_min) every 500
fix layer_$(v_z_min) kapton_layer_$(v_z_min) reaxff/species 1000 1 1000 kapton_molecules_$(v_z_min)_$(v_z_min+step).out
variable z_min equal $(v_z_min + v_step)
if "${z_min} < ${z_max}" then "jump SELF loop_layers"
I was able to save the IDs of the molecules using the next commands :
region region_layer_$(v_z_min) block INF INF INF INF ${z_min} $(v_z_min+ v_step) units box
group kapton_layer_$(v_z_min) dynamic all region region_layer_$(v_z_min) every 1
fix layer_$(v_z_min) kapton_layer_$(v_z_min) reaxff/species 10 1 10 kapton_$(v_z_min)_$(v_z_min+v_step).out cutoff * * 0.4
dump dumpMolID_$(v_z_min) all custom 10 molecule_ids_$(v_z_min).dat id type f_layer_$(v_z_min)
The problem is that in the output I have more than 73 IDs even at the beginning of the simulation. But I only have 3 types of molecules in my kapton (before it reacts with the oxygen). So i donāt think this is a good method to do what I want.
When describing the different regions, some molecules are divided in different regions. What I want is that each group of atoms include all the molecules not only the part that is present in each region. To be clear, what I want is lammps to add to every atom group, the whole molecule even if a part of the molecule is in another region.
I donāt know if that is possible with ReaxFF? Someone has an idea?
molecule IDs != molecule types. That is irrelevant from whether you using ReaxFF or not.
It may be possible, but youāll have to write your own C++ code based on the existing code.
Or you do the whole thing in post-processing. Then you may be better off using fix reaxff/bonds and also do the whole topology/molecule classification on your own.
This is becoming to narrow a use case and to specific to be something that would make sense to implement in LAMMPS as a generic feature.
I knew about this conceptual issue but I was thinking about correcitng it in post processing.
I am going to calculate the molecular composition of all my kapton (working only with one region), compute the center of mass of each molecule. Then Iāll write a python code that reads all this information and will count the molecules in each region. This way I will not cut the molecules.
Given the information being outputted by fix reaxff/species it SHOULD be possible to use compute chunk/atom ... f_FIX_NAME ... to āchunkā atoms into molecules.
Once chunking is successful, you can calculate each chunkās CoM (and thus each āmoleculeā's CoM) using compute com/chunk. You can also write a custom āelement variableā and sum across a chunk:
variable elem atom 100^(type-1) # types start from 1
compute molformula all reduce/chunk COMPUTE_NAME sum v_elem
Then a chunk with (for example) a molformula of 50401 is a molecule with 1 atom of type 1, 4 atoms of type 2, and 5 atoms of type 3.
You can then use compute chunk/spread/atom so each atom āknowsā what its molecular formula is, and so on.
