"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?