Problems about combining fix bond/react and fix gcmc


I want to realize a crosslinking process of a polymer chain based on two metathesis reactions with fixed the chemical potential of small molecules. In the data file, the atomtype of polymer is 1 and there are some of reaction site of the polymer chains with atomtype 2. The template of small molecules react with reaction sites is 3-4-3. In the fist step reaction, the reaction sites on polymer react with the small molecules:

1-2 + 3-4-3 = 1-4-3 + 2-3

where 2-3 is the small molecules produced after reaction. Here I want to fix the chemical potential of the two kinds of small molecules (3-4-3 and 2-3) so the fix gcmc would be performed. For dynamically group atoms in gcmc, I change the atomtype after reaction:

1-2 + 3-4-3 = 1-6-5 + 7-8

In the second step, the reaction sites reacted in the fist step will react with other reaction sites:

1-4-3 + 1-2= 1-4-1 + 2-3

For dynamically group atoms in gcmc, the atomtype would be changed:

1-6-5 + 1-2= 1-9-1 + 7-8

Therefore achieving cross-linking of polymer chains. The two steps are reversible.

Here is my lammps in. file:

#box information

units lj

boundary p p p

atom_style bond

#read data

read_data vitrimers.dat

molecule primer data_templates/primer.data_template

molecule byproduct data_templates/byproduct.data_template

molecule step1_unreacted data_templates/step1.unreacted.data_template

molecule step1_reacted data_templates/step1.reacted.data_template

molecule step2_unreacted data_templates/step2.unreacted.data_template

molecule step2_reacted data_templates/step2.reacted.data_template

#neighbor and velocity

neighbor 0.4 bin

neigh_modify every 1 one 100

velocity all create 1.0 123

#set LJ pair_style

pair_style lj/cut 2.5

pair_modify shift yes

pair_coeff * * 1.0 1.0 2.5

#set bond

bond_style harmonic

bond_coeff * 500.0 1.0

#reaction and gcmc

fix metathesis all bond/react stabilization yes nvt_grp .03 &

react step1_forward all 100 1.0 1.5 step1_unreacted step1_reacted map_files/step1_forward_map &

react step1_backward all 100 1.0 1.5 step1_reacted step1_unreacted map_files/step1_backward_map &

react step2_forward all 100 1.0 1.5 step2_unreacted step2_reacted map_files/step2_forward_map &

react step2_backward all 100 1.0 1.5 step2_reacted step2_unreacted map_files/step2_backward_map

run 50

variable type_primer atom “type==3 || type==4”

variable type_byproduct atom “type==7 || type==8”

group dynamic_primer dynamic all var type_primer every 1

group dynamic_byproduct dynamic all var type_byproduct every 1

variable Nprimer equal count(dynamic_primer)/3.

variable Nbyproduct equal count(dynamic_byproduct)/2.

fix gcmc_primer dynamic_primer gcmc 100 1 0 0 123 1.0 -7.0 0.5 mol primer group nvt_grp_REACT

fix gcmc_byproduct dynamic_byproduct gcmc 100 1 0 0 123 1.0 -7.0 0.5 mol byproduct group nvt_grp_REACT

fix 1 nvt_grp_REACT nvt temp 1.0 1.0 0.05

compute_modify thermo_temp dynamic/dof yes

compute_modify 1_temp dynamic/dof yes

thermo_style custom step atoms temp press vol epair & v_Nprimer v_Nbyproduct

run 1000000

where primer and byproduct is the molecule 3-4-3 and 2-3 before mentioned.

while after some steps, it will output an error message:

ERROR on proc 0: Bond atoms 12 73 missing on proc 0 at step 18151 (…/ntopo_bond_partial.cpp:64)
Last command: run 1000000

The message usually occured in a step that reaction happened. The lammps version is 7 Aug 2019.

I doubt the identity of some atoms may changed incorrectly in the bond/react or gcmc process.

Any ideas what’s going on?


Combining bond/react and gcmc is new territory, as far as i know, and would come along with a number of challenges. Both fixes execute a lot of internal logic/operations that would need to be communicated to each other. For example, bond/react generally assumes that other fixes are not creating or deleting atoms, or modifying its internally-created groups.

Frankly, i doubt compatibility between these two fixes will be acheived any time soon. However, if you send a simple example, we could start thinking about ways to make bond/react more friendly to other complex fixes

Hi Jacob,

Here is an example. where
./data_templates contains some template files for gcmc, primer.data_template, byproduct.data_template, and template files before and after bond/react (two steps).
./map_files contains mapping files for bond/react.
vitrimers.dat is initial configuration file, which contains a polymer chain with 10 atoms where 2 site atoms bonding on it, and 2 primer (3-4-3) molecules.
vitrimers_methathesis.test.lmp is lammps in file.

Actually I thought some questions concerning combining fix gcmc and bond/react. Here are something I did in the lammps in file:
Atoms increasing/deleting by fix gcmc would be added/deleted in group all. fix bond/react recognise possible reactions in the group all, maybe it can avoid some problems due to atom counting.
For fix gcmc, I used dynamic group by counting the atom type, but I not sure if I can use it like this. In addition, the increasing molecules are added in group nvt_grp_REACT for stabilization in fix bond/react.


vitrimers_test.tar.gz (2.65 KB)

In case this is still relevant to you, an upcoming improvement to bond/react that correctly updates molecule IDs, seems to improve the compatibility between bond/react and gcmc.

Your system now runs without issue (I tested up to 10000 time steps).

Currently available as a pull request on Github: