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