"bond atom missing" when using the bond/react command for a reversible reaction

LAMMPS LAMMPS General Discussion
1

Hi there,

I am encountering a “bond atom missing” error when using the bond/react command for a reversible reaction. I will describe my problem in as much detail as possible. Any help on the issues would be greatly appreciated!

As shown in the diagram below, there are two structures, A and B. I want to achieve a probability-based polymerization reaction between A and B within a specific distance range, leading to the removal of a water molecule. Please find my input file attached below.

image
image1116×556 143 KB

My pre and post template rules are as follows:

The T3(C)I11-1 atom and the T9(N)I56-7 atom initially form a bond. At the same time, the T5(O)I15-6 atom breaks its bond with the T3(C)I11-1 atom, and the T10(H)I59-9 atom and the T10(H)I60-8 atom each break their bonds with the T9(N)I56-7 atom. The T5(O)I15-6 atom then forms a bond with both the T10(H)I59-9 atom and the T10(H)I60-8 atom, resulting in the formation of a water molecule.

I have completed the writing of the pre and post templates and have successfully implemented the desired functionality. Following the instructions in the manual, I swapped the positions of the pre and post templates to enable reversible reactions. While the reversible reactions now occur as expected, there is a persistent instability issue with frequent “bond atom missing” errors. As shown below:

image
image1134×1125 302 KB

I have observed the trajectory files and found that most of the errors occur at the junction between the A monomer and the B monomer, specifically involving the T3(C)I11-1 atom and the T9(N)I56-7 atom. Additionally, there are also instances where bond atom missing errors occur with the generated water molecules.

I have tried various methods to resolve this issue. For example, I increased the Nevery value for the reverse reaction in an attempt to stagger the evaluation times of the forward and reverse reactions. I also experimented with adjusting the Rmin and Rmax values for the reverse reaction from 1 and 3 to 2 and 3 respectively. However, these attempts did not yield satisfactory results, and I am uncertain about how to proceed in resolving this problem.

Anyway, thank you all so much for your time and help.

Zilin

in.lmp (1.6 KB)
output (1.0 MB)
rxn1_stp1.map (223 Bytes)
rxn1_stp1_post.data_template (1.5 KB)
rxn1_stp1_pre.data_template (1.5 KB)
rxn1_stp1_re.map (223 Bytes)

2

The “Bond atom missing” error is a consequence of the domain decomposition in LAMMPS. In order to handle interactions that cross the subdomains, each subdomain has a “buffer” of atoms copied from the surrounding subdomains. The size of this buffer is determined by the largest cutoff radius plus the neighbor list skin. This distance is usually long enough for common force fields, if a reasonable cutoff is chosen. Your cutoff of 10 angstrom is rather short. More common is a choice of cutoff between 12 and 14 angstrom these days. That will increase the computational cost but lower the error from cutting of the Lennard-Jones interactions.

When using fix bond/react you are changing interactions drastically and that may lead to (temporarily) large forces and thus bonds stretched very wide and then it can happen that one of the bond atoms is too far away to be included in the “buffer” when the neighbor lists are recomputed.

The straightforward way to address this would be to either increase the non-bonded cutoff, of if that is not an option, change the communication cutoff with comm_modify. This will only increase the “buffer” size without adding more atoms to the neighbor list and non-bonded force computation.

3

A pairwise LJ cutoff of 10 angstroms is not a problem for ‘fix bond/react’ when the reaction cutoff distances are 3 angstroms, as in this example. However, ‘fix bond/react’ is known to not play well with constrained dynamics, such as ‘fix shake.’ Please try disabling ‘fix shake’ to eliminate the stability issue.

4

Thank you. I tried a cutoff of 14, but it still didn’t resolve my error.

5

Thank you, sir. I tried commenting out the “fix shake” command, but I still encountered a “bond atom missing” error. I also tried disabling the reverse reaction, which resolved the error. So, it seems to be an issue with the reversibility of the reaction.

Is it not sufficient to simply switch the order of pre and post steps for reversible reactions? Are there any other considerations I should be aware of?

6

sorry for the delay. I am unable to reproduce your issue. What version of LAMMPS are you using?

7

My LAMMPS version is 23Jun2022_update3.

Is there anything specific to consider when dealing with reversible reactions? I’m thinking that at certain times, simultaneous occurrence of forward and reverse reactions might lead to the bond atom missing.

8

Due the active development of fix bond/react, there have been several fixes since June 2022 (more than a year ago). I recommend updating to see if the issue persists.

There should not be simulation stability issues associated with reversing the reaction, in general. Issues could arise if you have reactions that create or delete atoms, or if molecules are fragmented into many small molecules. However, in these cases, the reverse reaction should simply not be able to occur, rather than causing simulation stability issues.

9

no, simultaneous reactions involving the same atoms are prevented

10

Thank you, sir. After updating to the latest version, the error still persists, which could indicate an issue with my molecular template. However, I am surprised to observe that with the same input file, the newer version of LAMMPS has reduced the number of reactions from 29 to 17. Additionally, overall computational efficiency has greatly improved. Could this also be a result of code updates?

11

can you add a link to your output file (and any input files, if you changed anything)? I will take another look. Yes, there were multiple updates for efficiency, and reactions may not occur on the same timesteps as they did for previous versions.

12

Thank you, sir. Here is my output file.
in.lmp (1.7 KB)
log.lammps (1.8 MB)

13

I was able to reproduce the unstable simulation. With your previous settings, the instability did not occur until more than one million steps, so I adjusted the settings so that reactions happen much more frequently (such that the error occurred within two to three thousand steps). The issue was with the water molecule coming off, because it was not well stabilized during the reaction. Adjusting the ‘xmax’ argument of the ‘stabilization’ keyword, or perhaps the stabilization time (‘stabilize_steps’ keyword), can generally fix this sort of issue. In this case, increasing xmax to 0.3 seems to have solved the issue. Please see attached for input file and output file.
in.lmp (1.9 KB)
out.lammps (8.9 MB)

1 Like
14

Thank you very much, sir. I tried your code, and everything is working fine now. However, I found that it is because of the line “fix 4 bond_react_MASTER_group temp/rescale 1 310 310 10 1” that the previous error was resolved. When I removed this line, the previous error occurred again. So, does it mean that for a less stable system, we need to use temp/rescale to prevent temperature fluctuations?

15

Great to hear everything is working for you. It is not always necessary to use temp/rescale, but yes it can help with stability and make your dynamics more physically realistic, such that small molecules do not end up with too much kinetic energy.