Stabilization keyword in fix bond/react command

[wzlclhy]

Hi there,

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.

Zilin

[jrgissing]

@wzlclhy

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.

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.

[wzlclhy]

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.

The above is my confusion.

Finally, thank you again for your generous help,

Zilin

[jrgissing]

@wzlclhy

1: “When using reaction stabilization, you should generally not have a separate thermostat that acts on the “all” group.”

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

@wzlclhy

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”."

In the above input file, my group id is “all ”. So I don’t quite understand the expression of the above sentence.

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.