I would like to use bond/react. I have looked through the examples, which worked fine, but I’m running into an issue in my case.
My system consists of linear polymer chains. I have bonds, angles, and dihedrals in the system. I want to change and ultimately break a bond in the linear chain.
For this case, I thought that I need to specify 8 atoms in the molecule files and map_file. This way the edge atoms are 3 bonds away from the bond ID atoms.
However, when I specify this the simulation crashes with an error of “bond atoms missing”. The reason for this is that the potential energy becomes really high. However, there is no reason the potential energy should be so high. In fact, even when the template-ID(pre-reacted) and template-ID(post-reacted) are the same, there is a huge increase in the pair-wise potential energy and the simulation crashes.
I tried only specifying 6 atoms in the map_file and the template-ID molecule files. If there is no change in atom types in the template-ID(post-reacted), this simulation runs. The fix bond/react shows the output of cumulative reactions occurring. However, even though the template-ID(post-reacted) had specified that a bond had changed type, no bonds changed type as seen by a data file written after the fix bond/react run a while. Thus, even though the thermo output of fix bond/react showed reactions happening, the bond change specified in the templates did not occur.
I tested that there was no issue in the example case when the template-ID(pre-reacted) and template-ID(post-reacted) were the same. I think there perhaps is something I’m missing with my linear polymer case.
I would like some guidance on this issue. I’m happy to provide any additional info/files that would be helpful in addressing this issue.