Questions about REAXFF package commands

Hello,

I have a few questions regarding the REAXFF user package commands, which I address in the following bullets:

For atom_style charge the molecule id of the reax/c/bonds file is set to zero everywhere. Using atom_style full does indeed assign a molecule id (including the possibility to output it via dump mol), however this value refers to initial molecules and stays static throughout the simulation, even though reactions occur. I.e.:

atom_style charge:

id type nb id_1…id_nb mol bo_1…bo_nb abo nlp q

183 3 1 184 0 1.584 1.584 2.000 -0.057

184 3 1 183 0 1.584 1.584 2.000 -0.055

215 3 1 216 0 1.584 1.584 2.000 -0.092

atom_style full:

id type nb id_1…id_nb mol bo_1…bo_nb abo nlp q

21 1 4 22 24 25 23 5 0.967 0.967 0.967 0.967 3.951 0.000 -0.794

28 3 2 29 58 7 1.584 1.003 2.589 1.121 -0.009

35 3 1 34 10 1.584 1.584 2.000 -0.043

On the contrary, using the reax/c/species command with the attribute position outputs a different molecule id which changes dynamically for both cases of atom_style. I.e.:

ID Atom_Count Type Ave_q CoM_x CoM_y CoM_z

1 5 CH4 0.08069482 0.68288650 0.19484946 0.39072814

2 5 CH4 0.07244117 0.88381249 0.28605048 0.65662160

3 9 CH4O4 0.01159357 0.65641404 0.02605985 0.81482721

4 5 CH4 0.12252271 0.62514317 0.93977286 0.56528468

5 7 CH4O2 0.03767128 0.45542362 0.28864767 0.45142155

6 2 O2 -0.01736948 0.94474667 0.66583586 0.68254349

Thus, my question is why doesn’t the molecule id of the dump and reax/c/bonds commands take the values of the reax/c/species position ones, which can change with time? To my understanding the current implementation pretty much limits the usage of all compute chunk commands with the compute chunk/atom molecule, since information do not correspond to the realistic values.

  • Another question I have is about the reax/c/species command resetting the neighboring criteria based on the values of Nevery Nrepeat Nfreq. Inside the fix_reaxc_species.cpp there is a comment that says “Neighbor lists must stay unchanged during averaging of bonds, but may be updated when no averaging is performed”. For reliable results I believe one should choose a maximum of Nevery or Nrepeat=10 or less, since this will result in update every 10 steps, delay 0 steps, check no, which -I think- may still be dangerous for certain cases.
    However, then a Nfreq of 10 should be chosen or else reactions may be lost in between the averaging, which eventually will result in huge files for big simulations. Please correct me if my thinking is wrong here.
    Additionally, I would like to ask: if I put the command neigh_modify every 10 delay 0 check yes after a reax/c/species 10 100 1000 command for example, is it guaranteed that I will get unrealistic results? (For a simple and small case I tested I got correct results)

  • And finally, I just wanted to note that when using the reax/c/species with the position attribute, the resulting position file is empty if there is a minimize command before.

Thank you very much for your help in advance,

Efstratios

Hello,

I have a few questions regarding the REAXFF user package commands, which I address in the following bullets:

For atom_style charge the molecule id of the reax/c/bonds file is set to zero everywhere. Using atom_style full does indeed assign a molecule id (including the possibility to output it via dump mol), however this value refers to initial molecules and stays static throughout the simulation, even though reactions occur. I.e.:

atom_style charge:

id type nb id_1…id_nb mol bo_1…bo_nb abo nlp q

183 3 1 184 0 1.584 1.584 2.000 -0.057

184 3 1 183 0 1.584 1.584 2.000 -0.055

215 3 1 216 0 1.584 1.584 2.000 -0.092

atom_style full:

id type nb id_1…id_nb mol bo_1…bo_nb abo nlp q

21 1 4 22 24 25 23 5 0.967 0.967 0.967 0.967 3.951 0.000 -0.794

28 3 2 29 58 7 1.584 1.003 2.589 1.121 -0.009

35 3 1 34 10 1.584 1.584 2.000 -0.043

On the contrary, using the reax/c/species command with the attribute position outputs a different molecule id which changes dynamically for both cases of atom_style. I.e.:

ID Atom_Count Type Ave_q CoM_x CoM_y CoM_z

1 5 CH4 0.08069482 0.68288650 0.19484946 0.39072814

2 5 CH4 0.07244117 0.88381249 0.28605048 0.65662160

3 9 CH4O4 0.01159357 0.65641404 0.02605985 0.81482721

4 5 CH4 0.12252271 0.62514317 0.93977286 0.56528468

5 7 CH4O2 0.03767128 0.45542362 0.28864767 0.45142155

6 2 O2 -0.01736948 0.94474667 0.66583586 0.68254349

Thus, my question is why doesn’t the molecule id of the dump and reax/c/bonds commands take the values of the reax/c/species position ones, which can change with time? To my understanding the current implementation pretty much limits the usage of all compute chunk commands with the compute chunk/atom molecule, since information do not correspond to the realistic values.

This can only be explained by looking at the history of the USER-REAXC package. This was initially a standalone code written in C (not C++), that has then be adapted and integrated into LAMMPS. So it didn’t know anything about molecule IDs and molecule IDs in LAMMPS predate the use of reactive force fields in LAMMPS. Molecule IDs are just a way to partition the system, so I am not certain that making those dynamical would be a good idea. It may break features that expect molecule ID to remain unchanged (like atom IDs) during a run.
Rather than using the explicit molecule ID in compute chunk/atom, you can use the per-atom vector that is produced by fix reax/c/species and that would allow you to do the partitioning that you desire. this is the way of performing this kind of operation consistent with how LAMMPS is designed.

  • Another question I have is about the reax/c/species command resetting the neighboring criteria based on the values of Nevery Nrepeat Nfreq. Inside the fix_reaxc_species.cpp there is a comment that says “Neighbor lists must stay unchanged during averaging of bonds, but may be updated when no averaging is performed”. For reliable results I believe one should choose a maximum of Nevery or Nrepeat=10 or less, since this will result in update every 10 steps, delay 0 steps, check no, which -I think- may still be dangerous for certain cases.
    However, then a Nfreq of 10 should be chosen or else reactions may be lost in between the averaging, which eventually will result in huge files for big simulations. Please correct me if my thinking is wrong here.
    Additionally, I would like to ask: if I put the command neigh_modify every 10 delay 0 check yes after a reax/c/species 10 100 1000 command for example, is it guaranteed that I will get unrealistic results? (For a simple and small case I tested I got correct results)

sorry, cannot help you here. this needs some explanation from somebody familiar with using USER-REAXC regularly for research.

  • And finally, I just wanted to note that when using the reax/c/species with the position attribute, the resulting position file is empty if there is a minimize command before.

this i have to answer with a (bad) joke:

Doctor, doctor! It hurts when I do this.
Well, don’t do it, then!

Axel.

Dear Axel,

Thank you very much for your quick reply! I had not thought of the fix reax/c/species atom vector, which seems to be the solution to my issue. I definitely agree with you on my third question, however I only point it out in case anyone has this issue in the future. Using a restart file after the minimization withdraws this problem.

Thank you again,
Efstratios

PS. I would really appreciate if someone could ask my question regarding the reax/c/species neighbour list issue (2nd question).