I want to be able to use ReaxFF but define some bonds that cannot be broken. Initially I thought about using a atom_style hybrid charge bonds but I am not sure, has anyone done anything similar in the past and could point me in the right direction?
Any help is appreciated.
Reaxff computes bonds impliticly, so there should not be any topological, explicit bonds.
Thank you for your reply, I understand that ReaxFF works via bond distance but is there no way of ensuring certain bonds don’t break?
Yes, there is, but you must recognize that you cannot replace ReaxFF interactions with explicit bonds, since that would cause all kinds of inconsistencies and thus lead to (partially?) bogus results.
Thus the simplest way would be to use the fix restrain command — LAMMPS documentation
If you need to restrain more bonds than just a few, it gets complicated, so I won’t even discuss it unless you have a good explanation why you need that, how you would determine the interaction strength and how you plan to justify that model.
Thank you for your comments, I have had a rethink. The reasoning behind my initial request was due to this part of the initial ReaxFF paper:
‘During the equilibration simulations, the C−O and H−O bond parameters were switched off to prevent reactions from occurring during the equilibration of the system’ - https://doi.org/10.1021/jp709896w
And I had assumed that the way they did this was by inducing strict bonds during equilibration. Alternatively, do you know how I can achieve this, I noticed a similar question posted earlier and the answer pointed toward there being new bonds being introduced.
In science it is always a bad idea to assume. Instead you should contact the authors of the paper to explain their procedure in more detail. Changing bond parameters is one thing (they are the 3rd block of the ReaxFF parameters), adding new explicit bonds is another.
Yes, you can introduce explicit bonds, but that is rather complex, as I have already explained and you should only introduce then as an additional interaction.
To explain: in traditional classical force fields, bonded interactions are parameterized without considering the non-bonded interactions for a given pair of atoms. That makes it rather straightforward to infer force constants for (harmonic) bond potentials (or angles, etc.). However, that requires that in the computation of non-bonded interactions, those pairs of atoms are excluded. In LAMMPS that is controlled with the special_bonds command. The default setting is to exclude all 1-2, 1-3, and 1-4 pairs. This is straightforward to apply to pairwise additive non-bonded interactions like Lennard-Jones. However, with manybody potentials (like ReaxFF) you are not only excluding a specific pair of non-bonded interactions, but since the neighborlist is also used for charge equilibration, collecting list of potential bond, hydrogen bonds, angle, and dihedral interactions, applying those exclusions will seriously taint the ReaxFF model and exclude contributions that must not be excluded.
The way around this limitation would be to use
special_bonds lj/coul 1.0 1.0 1.0. That will recover the original neighbor list without consideration for excluding bonded interactions, but that also means that now any explicit bonds will be applied on top of the corresponding ReaxFF interactions. For the purpose of prohibiting reactions, it may work, but - as already mentioned - it will be crucial to derive meaningful parameters and then you will need to strip off those additional interactions at some point and will have to have some additional (short) period of re-equilibration. But that latter step would also be required if you modify the force field parameters.
Thank you for your input, I appreciate it. Thank you for explaining why my question would not yield good results.