I’m having some problems with stabilization keyword in fix bond/react command. Any help on the issues would be greatly appreciated!
The following is described in the official documentation，
Fix bond/react creates and maintains two important dynamic groups of atoms when using the stabilization keyword. The first group contains all atoms currently involved in a reaction; this group is automatically time-integrated by an internally-created nve/limit integrator. The second group contains all atoms currently not involved in a reaction. This group should be controlled by a thermostat in order to time integrate the system. The name of this group of non-reacting atoms is created by appending ‘_REACT’ to the group-ID argument of the stabilization keyword, as shown in the second example above.
In my computational system, I deposited a lot of small molecules down a designated area of the box (using fix deposit command). Then, I used fix bond/react to make these small molecules bond to each other when deposited on the lowermost plane. I defined the molecular template as shown below，
But I found that no matter whether the thermostat is performed on the second group (‘_REACT’), the whole system is running normally without error.
My input script is as follows，
group addatoms type 2 3 4 #Atom type in deposited molecule
region slab block 50 300 50 250 40 45 # Molecules are deposited down from this region
region mobile block INF INF INF INF 2 INF
group mobile region mobile
molecule dimer CHP_charge_resp # molecule template for fix deposit command
molecule mol1 rxn1_stp1_pre.data_template # molecule template for bond/react command
molecule mol2 rxn1_stp1_post.data_templatem # molecule template for bond/react command
fix 1 addatoms nve
fix rxns all bond/react stabilization yes statted_grp .03 react rxn1_stp1 all 1 0.0 3 mol1 mol2 rxn1_stp1.map
fix 2 mobile nve
fix 3 mobile langevin 475 475 100 587283
fix 4 addatoms deposit 187 1 5000 12345 region slab near 2.0 mol dimer vz -0.002 -0.002
fix 5 addatoms wall/reflect zhi EDGE
In this input script, I did not use a specific thermostat to control the second group (‘_REACT’), but there was no abnormality in the running results of the system.
I don’t know if this understanding is correct: when the second group is already controlled by another thermostat, there is no need to perform a thermostat on the second group (‘_REACT’) alone.
Anyway, thank you all so much for your time and help.
Yes, this is correct, and the input you provided looks fine. Essentially, you are performing your own reaction stabilization. For the default stabilization, the very effective (but less physical) ‘nve/limit’ was chosen to locally relax high-energy reactions, but if using Langevin and a short reaction distance works for your system, that is a great alternative.
Sir, thank you for your reply, it solves a problem that has been bothering me for a few days.
I still have a problem with the stabilization keyword. According to the lammps manual, there is a description of the stabilization keywords as follows,
1: “When using reaction stabilization, you should generally not have a separate thermostat that acts on the “all” group.”
I don’t understand what the above sentence means. My understanding is that when the stabilization keyword is used, the system will create two dynamic groups. We must apply the same thermostat to both groups? If so, what does “all” mean here?
2：" The group-ID set using the stabilization keyword can be an existing static group or a previously-unused group-ID. It cannot be specified as “all”."
fix rxns all bond/react stabilization yes statted_grp .03 react rxn1_stp1 all 1 0.0 3 mol1 mol2
In the above input file, my group id is “all”. So I don’t quite understand the expression of the above sentence.
Thanks for bringing this up, and for carefully reading the docs! Actually, in your script, reacting atoms are being time integrated by multiple fixes: both by the ‘nve/limit’ fix created internally by ‘fix bond/react’ for stabilization, and by the your global nve+Langevin thermostat. So, you should be using bond/react’s internally-created group for nve+Langevin, e.g., ‘statted_grp_REACT’. However, in your case, you probably want to use the existing ‘mobile’ group instead of creating a new group, which will default to including all atoms. In other words:
fix rxns mobile bond/react stabilization yes mobile .03 react rxn1_stp1 all 1 0.0 3 mol1 mol2 rxn1_stp1.map
fix 2 mobile_REACT nve
fix 3 mobile_REACT langevin 475 475 100 58728
I think this question is mostly answered above. To clarify, the “group-ID set using the stabilization keyword” refers to ‘statted_grp’ in your original script.
Sir, please forgive me for one more thing I don’t quite understand.
The group-ID you mentioned above is ‘static_ grp’. But there are two other all group-ID in the script above. I don’t quite understand the difference between the two all here.
My understanding is that it is possible to control the group ID here to limit the location of reactions in the box?
Referring the notation in the documentation, only atoms in the group react-group-ID (the second all) are considered for that particular reaction. The first all is ignored, as is the case with some other LAMMPS fixes.
I did not use a specific thermostat to control the second group (‘_REACT ’), but there was no abnormality in the running results of the system.
You said that I’m performing my own reaction stabilization. Indeed, everything is running fine, but when I need to read the restart file I get an error：
ERROR on proc 3: Bond atoms 57558 57588 missing
Similarly, the error disappeared when I comment out fix 1 addatoms nve in the restart input file.
Therefore, I guess the problem comes from When using reaction stabilization, you should generally not have a separate thermostat that acts on the “all” group.
What I don’t understand is why the error occurs only during the restart process. (my LAMMPS version is 23Jun2022_update3)
I also tried what you suggested to me，
But the same problem exists.
In my system, the target group using bond/react has to be controlled by another thermostat (fix nve). So I have to give up reaction stabilization?
The way I suggested above is correct. The system is time integrated using ‘fix nve’ and the Langevin thermostat, and reacting atoms are stabilized (using nve/limit). If your system is unstable using these commands, please include a simple example that reproduces the problem with all files needed to run it.
