Help with fix bond/react for Model Construction

Dear lammps users,

I would like to use the fix bond/react command to construct a molecular model. However, I am currently confused about how to implement it correctly for my specific system. I would greatly appreciate any guidance or suggestions you could provide.

Below is a description of my model:
I want to construct a linear polymer (LJ unit, FENE model, no angle) using fix bond/react. To simplify the system, I initially pack many short molecules of the form type 2–1–1 into the simulation box. My goal is to use fix bond/react during the simulation to form a linear polymer, where type 2 beads form the main chain and type 1 beads act as side chains.

The intended structure looks like this:
2–2–2–2–2–2–2–2–2
| | | | | | | | |
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
I’m not sure how to properly define the templates and reactions for this process in fix bond/react, and would greatly appreciate any advice or example input files that could help.


blue bead is type 2 bead, orange bead is type 1 bead

@jrgissing Thank you.

I would start simple. Have you tried building a template that simply forms a bond between two monomers?

Thank you for your suggestions. I think I’d like to approach this by organizing it into five steps. From monomer to dimer, trimer, tetramer, and eventually extended chains. Thank you once again for your helpful suggestion.

Sounds like a good plan. Here are some resources related to bond/react:
Examples that come with LAMMPS: examples/PACKAGES/reaction
REACTER website: https://www.reacter.org/tutorial
Detailed tutorial for atomistic simulation: Reactive Molecular Dynamics -
docs: fix bond/react command — LAMMPS documentation

1 Like

I’m not familiar with fene potential, but unless you do have dihedral terms it might be easier to use fix bond/create instead of bond/react.

Thanks for your answer. I use a simplified coarse-grained FENE model that excludes dihedral interactions.

As said, then you might want to have a look at fix bond/create, which does not need templates.

I have constructed the template file and successfully achieved step-by-step polymerization — from monomer to dimer, and from dimer to trimer. However, I am now puzzled about how to simulate radical polymerization more realistically. In radical polymerization, chain growth can occur in a stochastic manner, where, for example, a dimer can react with a trimer, or two trimers can couple. But in fix bond/react, the polymerization proceeds strictly step by step according to predefined templates. How can I simulate such random reactions occurring simultaneously?

I think each of those would be considered a different reaction and need it’s own template (until you’ve reached the necessary size of the templates / “edge”).
E.g.
Reaction 1: 2 Monomers → Dimer
Reaction 2: Monomer + Dimer → Trimer
Reaction 3: Monomer + Trimer
Reaction 4: 2 Dimers
Reaction 5: Dimer + Trimer
Reaction 6: 2 Trimers

I think after that growth should continue with these reactions as you only specify bonds and angles. And the rest can be behind the “edge” particles.
I am not sure if the sequence of the reaction definitions play a role role here. You could play around with the probability parameter, so that not all reactions happen instantly, but the system has time to evolve.

as @wolfgang mentioned, ‘edge atoms’ are atoms that connect to connect to the rest of a polymer chain (or other larger molecule). All of the examples shipped with LAMMPS involve polymerization and include edge atoms, so those could be a good place to start.

I really appreciate your help. I will try it again.

Thank you very much for your valuable suggestion. I will try it again.

One second thought, there are also building/preprocessing tools available. E.g. EMC (Welcome) or Moltemplate (lammps/tools/moltemplate at develop · lammps/lammps · GitHub) among others might be worth looking at.