Issues with fix bond/react not creating the specified bond

Hello,
I have run into some odd behaviour with the bond/react fix that I am struggling to understand. For reference, I am using the latest stable version of LAMMPS, lammps-2Aug2023. My input files, data files, and log file can be found here: Sign in to your account .

I am trying to create some crosslinks in this polymer system. I specify that the reaction should be initiated when two carbons get within 3.2 angstroms. The prereact.mol and postreact.mol templates specify that two hydrogens should break their bonds with the carbons, and that the carbons should form a bond between them.

The reactions get triggered, and the hydrogens get deleted as expected. However, no bond (bond type 4 in the files I posted) is formed between the carbons. I have tried a number of things and I have found that:

  • I can make a bond between the hydrogens.
  • Changing the initiator atoms to be the hydrogens makes no difference
  • Creating an angle between the carbons and a hydrogen works fine as well
  • If I convert the hydrogens to another atom type, these are properly converted (so the reaction is definitely taking place)
  • If I set the reset_mol_ids keyword, two different molecules will become one molecule despite there being no explicit bond between them.

I’ve compared my inputs with the examples provided with LAMMPS, and I am finding no significant difference in how I’ve set up the inputs. There are no errors, and only one warning about the overall charge not being neutral. The only think I can think of is that there may be some sort of cutoff that may be somehow excluding the bonds (though I could not immediately find this in the documentation), but even running an NPT after the reaction run does not lead to the eventual appearance of the desired bonds.

If anyone (though this may be most pertinent for @jrgissing ) could give me any insight into what I have potentially done wrong, it would be greatly appreciated!

Many thanks,

Christophe

looks like your edge atoms to too close to your bond-forming atoms. adding one more atom attached to your current edge atoms should fix the issue. you are getting away with this without bond/react complaining because none of your reacting atoms change atom type (should they?)

Thank you for that bit of advice. That did indeed work. Just for my own understanding then:

  • There is some kind of distance criterion regarding the proximity of edge atoms to the initiator atoms? I hadn’t spotted this in the documentation.
  • As far as I know the atom type shouldn’t change, though I have not yet double-checked this. This was mostly just a prototype to see what the procedure needs to be before building the final model. Is this change of atom type somehow also a mandatory criterion?

Yes, this part of the bond/react docs describes the issue:

“The pre-reacted molecule template should contain as few atoms as possible while still completely describing the topology of all atoms affected by the reaction (which includes all atoms that change atom type or connectivity, and all bonds that change bond type). For example, if the force field contains dihedrals, the pre-reacted template should contain any atom within three bonds of reacting atoms.”

Usually, either the atom type or the total number of bonds changes for reacting atoms. These are not mandatory criterion, they are just what are currently used for error checking. In your case, I may need to add another more stringent type of error check.

Ah of course, I had misinterpreted what you had said. It’s not the distance per se, it’s the number of atoms away. I think I had misunderstood/misread this part of the documentation then; in hindsight this should really have been obvious. Thanks for the clarification.

One last question: I was able to make an angle even when the bond was not made (I believe in the case above it was a C-C-H angle that I had added just as a test between the initiator atoms). Should this have been possible? It seemed strange to me at the time that a bond would not form, but an angle could.

No, an angle should not form when the corresponding bonds do not form, under normal operation. In the current case, an error should have been thrown in the first place, so there could be some strange behavior.

Hi,
I also encountered this problem as well when I run a condensation reaction. Once the reaction occurred, a water molecule will be deleted and a C-O bond should form. In my simulation, a water molecule was removed. I also found the charge of reaction atoms change and 2 molecules become 1 molecule after using “reset_atoms mol all”, but no new C-O bond was observed. There’s no warning or error message about “edge atoms to too close to your bond-forming atoms”.

If anyone could give me any advise or suggestion, it would be greatly appreciated.
Thanks in advance,
Zongyun

This problem was resolved. Can you resolve your problem the same way using the advice above? If not, please post a new issue that includes details of your reaction, preferably including files for a small system that reproduces your problem.

Sorry for late reply, I solved this problem by adding enough atoms. Thanks for your kind help!

Sincerely,
Zongyun