Mismatch in reaction number with fix bond/react

LAMMPS LAMMPS General Discussion
1

I am exploring reactions with colloidal particles. The idea is to have rigid-body particles with 5 non-bonded sites, and a total 3 different types of beads. The system is designed to have 2 types of non bonded beads, where interactions between one bead type cause a reaction that activates the other bead type. Reactions are made through a fix bond/react command, which changes the bead type. Specifically, the bead types can be thought of as 1) type 1: the 2 sites that have a constant interaction energy that initiate reactions (blue), 2) type 2/4: 3 sites that react from type 2 to type 4 based on interactions by type 1 (purple before reacting, red after reacting), and 3) type 3: one bead at the center of mass (gray). For now, I am trying to have reactions occur whenever beads of type 1 are within a cutoff distance (rcut), and be undone by the reverse reaction when they separate by further than rcut. This idea is expressed visually in Figure 1.

Simple_Example.pdf (90.8 KB)
Figure 1: Example of a reaction. Colloidal particles consist of 4 atom types, type 1 (blue), type 2 (purple), type 3 (gray), and type 4 (red). Reactions occur when the distance between type 1 atoms (r) is less than a cutoff (rcut), switching type 2 to type 4.

To test out this system, I created a simple system with 5 colloidal particles, and ran a short simulation. The idea is that 2 pairs of particles are close enough to start a rection, causing 1 favorable type4-type4 interaction, while the type2-type2 interactions and type2-type4 interactions are negligible. The initial system (before any reactions occur) is illustrated in Figure 2.

5body_example.tga (1.3 MB)
Figure 2: Simple system before reaction. Type 1: blue, type 2: purple, type 3: gray, and type 4 red. Note no type 4 particles are present because the simulation has not started.

I performed this test with 1 timestep on a variety of type1-type1 and type2-type2 distances, got the results I expected, and started to look at larger systems. However, when I ran the larger systems, I noticed that I systematically got less reactions than expected. To understand these results, I turned back to the simple, 5-particle system and ran longer simulations, consisting of 1000 steps.

When I run these simulations, I can track how many reactions I expect to occur, based on how many particles are within rcut (should_react.txt), and how many reactions actually occurred, based on how many type 4 particles are in the system (reacted.txt). For a significant fraction of the simulation, less reactions occur than expected (Figure 3A). To understand this better, I examined the distance of the 2 pairs of type 1 beads that start off within the cutoff distance, and compared this to when they actually separate by rcut and when the back reaction actually occurs in the simulation (Figure 3B). Interestingly, the distances do not seem to be realated to when the back reaction occurs.

Reaction.pdf (435.8 KB)
Figure 3: Comparison between expected and actual interactions. A) Number of reacted particles (measured as the number of type 4 particles in the system) compared to the number that should have reacted (based on the number of type 1 beads within rcut). The number of type 4 particles is less than expected for over half the simulation. B) Distance as a function of time for the 2 pairs of type 1 beads that start off within the cutoff (blue/red). This is compared to rcut (black) and the time at which the beads switch from type 2 to type 4 (green). The reaction occurs before expected, even though the distance is much smaller than the cutoff.

I was wondering if anyone had any thoughts or guidance on what might be going on and why the fix bond/react does not seem to be working as I would expect in this case? (I have been using LAMMPS: 23 Jun 2022 - Update 4) My code for the simple, 5 particle simulation can be found here: bond_react_test.tar - Google Drive

2

Is it the bond-forming or bond-breaking reaction that is not behaving as expected? If bond-forming, try using the ‘molecule inter’ keyword.

3

Follow-up: I see in your attached files that you are already making good use of the ‘molecule’ keyword. One issue might be that, even for non-stabilized reactions, more than one reaction involving the same atoms is prevented from occurring at the same time. This restriction is not necessary in your case. Please find a (not thoroughly tested) version of bond/react that removes this restriction.

As a side note, I am unsure what you are going for with this set of reactions, but it seems like the connectivity of type-2 vs type-4 patchy particles might be ambiguous, at any given time. For example, if three particles ‘bond’ together to form three type-4s, one could debond, leaving a type-2 bonded to a type-4.
modified_fix_bond_react.zip (39.9 KB)

4

Thank you so much for your help. While expected “molecule inter” to be the most appropriate in this case, I have also removed that keyword to allow all reactions, and the results remain identical.

When you say “more than one reaction involving the same atoms is prevented from occurring at the same time”, do you mean that multiple reactions with the same atoms are not allowed at the same timestep? Thus, even in the default version of the code, you can react in timestep 1, and then react again in timestep 2 to further extend the chain, allowing for a model of “polymerization.”

If so, that does not seem to be the primary issue here. I did test your new code on my model, but the problem is not resolved (though it may help in the future one the immediate problem is resolved). The current problem seems to be that the back reaction seems to occur in instances where it should not occur (there may also be issues with the forward reaction not occurring when it should occur because it would be expected at the following timestep, but let’s start with the back reaction). For example, here, I have plotted the pairwise type 1 distances, which should be the criteria for interacting (Figure 4). Upon close analysis of these distances, it is clear that A and C and B and D remain interacting for the majority of the simulation, yet the backwards reaction happens quickly in the simulation. This back reaction occurs despite no pairwise interactions being formed / broken between that point in the simulation and the simulation start. (Note the probability of both reactions is set to 1.0 in this case, so they should occur whenever particles are within the appropriate cutoff).

Distances.pdf (1.2 MB)

Figure 4: Pairwise distances for intermolecular type 1 beads for each of the 5 molecules in our simulation. This is compared to rcut1 and rcut2, the distance cutoffs (black) for the forward reaction (~0.4) and back reaction (4) respectively, and the time at which the back reaction occurs, switching the beads from type 2 to type 4 (green). The back reaction occurs before expected, despite no beads crossing either distance cutoff near this timepoint.

You also made a very going point about reactions occurring type-2 with type-2 but not type-2 with type-4, as both should be possible. I added an extra reaction call to allow for such interactions to occur. However, this somehow further increases the frequency at which back-reactions occur, which I did not change in my code (Figure 5). Further, the back reaction occurs with an increased frequency, allowing for multiple forward and backwards interactions that is even more perplexing.

Reaction_double.pdf (484.6 KB)

Figure 5: Comparison between expected and actual interactions when adjusting for type2-type4 interactions. A) Number of reacted particles (measured as the number of type 4 particles in the system) compared to the number that should have reacted (based on the number of type 1 beads within rcut). The number of type 4 particles is less than expected for over half the simulation. B) Distance as a function of time for the 2 pairs of type 1 beads that start off within the cutoff (blue/red). This is compared to rcut (black) and the time at which the beads switch from type 2 to type 4 (green). The reaction occurs before expected, even though the distance is much smaller than the cutoff.

This leads to a few questions I had: How does the logic of fix bond/react work if multiple reactions are possible? What is the logic behind only allowing for 1 fix bond/react call? What if multiple fix calls could be allowed, allowing for ordering of reactions within a single timestep? That may help in this case.

Thank you so much for your time.

5

Okay I think I understand. I simplified your system to include just two patchy particles, and observed the back reaction even without dynamics. The issue is that the far blue atom (which are actually type 3 in your files) is within the distance cutoff of 4. Switching the cutoff to 3 results in no backreaction for this simplified system. This type of thing is more prone to happen when using bond/react without actually creating any bonds, as you do here (in some similar cases, it helps to create dummy bonds that do not actually affect the dynamics).

6

When you say “more than one reaction involving the same atoms is prevented from occurring at the same time”, do you mean that multiple reactions with the same atoms are not allowed at the same timestep?

Yes

Thus, even in the default version of the code, you can react in timestep 1, and then react again in timestep 2 to further extend the chain, allowing for a model of “polymerization.”

Correct

How does the logic of fix bond/react work if multiple reactions are possible?

It chooses randomly between eligible reactions.

What is the logic behind only allowing for 1 fix bond/react call? What if multiple fix calls could be allowed, allowing for ordering of reactions within a single timestep?

Consolidating all reactions into one bond/react command in itself did not remove any functionality, and was required for things such as preventing multiple atoms from participating in more than reaction on the same timestep and not having reactions necessarily happen in a certain order. You are using bond/react is a creative way, but these are usually desired features and were much more difficult to implement than not enforcing these requirements. Good news is ‘undoing’ them can be relatively straightforward; for example, the modified code I sent above consisted of commenting out two lines. A modification could also be made to enforce an ordering or ‘priority’ of one reaction over another. This is a relatively niche use case: although this sort of feature has been brought up to me once before, the user was able to achieve what he wanted without it. Please let me know if the same is the case for you.

7

Thank you so much for your help and detailed replies. I incorporated dummy bonds into the system, and that now seems to give the desired results on the small system. Thanks again!