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
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
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.