Matter Modeling Stack Exchange is a question and answer site for materials modelers and data scientists. It only takes a minute to sign up.
Sign up to join this community
How does one introduce the possibility of breaking a bond in a molecule during MD simulation?
I only found the cases when we just introduce the harmonic potential $U_{ij} = \frac{1}{2}k(r_i - r_j)^2$ for atoms in a molecule or introduce an additional equation restricting the length of a bond during simulaton.
Nevertheless, it is still not clear for me how we can take into account the breaking bond result?
The quadratic potential is the simplest possible model for a bond. You can derive it by considering the Taylor expansion of the potential around the natural bond length
$V(r - r_0) = V(r_0) + \frac{d V(r_0)}{d r} (r - r_0) + \frac{1}{2} \frac{d^2 V}{dr^2} (r - r_0)^2$
The constant term can be set to 0 since it does not contribute to the force and just sets the 0 point of your energy scale. The linear term is 0, because you are at a stationary point [1].
That leaves you with $V(r - r_0) = \frac{1}{2} \frac{d^2 V}{dr^2} (r - r_0)^2$
If we require the second derivative to be constant, we have recovered the harmonic potential. As @Charlie Crown illustrated, this force resulting from this potential does not go to 0 at infinity, while the Morse potential does. You can of course take a polynomial of higher order than two, but not every order is suitable. A third order polynomial results in a potential that (typically) goes to negative infinity at large $r$, so instead a quartic potential is some times used. It has the advantage of being slightly "wider" than the quadratic one. That said, a completely unrelated reason why none of these can simulate bond breaking/formation is that the implementation requires explicit declaration of which atoms should interact via the stretching potential.
Still, at large $r$ both differ significantly from the Morse potential. Why then is the Morse potential not used? The restoring force for large $r$ is very low in case of the Morse potential, hence it takes longer for the bond length to return to the equilibrium position. The quadratic potential describes the potential well for displacements close to equilibrium and for moderate temperatures, this is the part of the potential you care about.
Obviously that still leaves the question of how to simulate bond breaking in a force field. ReaxFF assumes that the bond order of a pair of atoms can be determined from the interatomic distance alone.
(qualitative recreation from [2])
The sigma, pi and double pi bonds contribute increasingly to the overall bond order (max individual bond order is 1) as the atoms get closer together. For simplicity I am leaving out the corrections made to the overall bond order necessitated by overcoordination. The bond stretching potential takes the form of a modified Morse potential
$E_{Bond} = -D_e \cdot BO_{ij} \cdot \exp(p \cdot (1 - BO_{ij}^p))$
where $p$ is a bond specific parameter[2].
References:
[1]: Frank Jensen, Introduction to Computational Chemistry Chap. 2
[2]: J. Phys. Chem. A 2001, 105, 9396-9409
The potential you showed is the most common form of bond, the harmonic potential a.k.a. Hookes law.
If you have ever broken a bond in QM (calculated a dissociation curve), you know it is a bit tricky, you need to use "unrestricted" settings, meaning, that a given pair of electrons does not need share the same orbital. As an added point for accuracy at the risk of distraction... saying you need "QM" is a bit hand wavey, although, not wrong. For real accuracy in QM, you need more than one Slater determinant to approximate the wave function. Now this is somewhat besides the point, but it lets me show a picture hastily drawn in powerpoint, so don't judge...
As you can see, the potential has a minimum, and at that minimum, and for a small ways out from it, as you stretch or compress the bond from equilibrium (the minimum is the equilibrium energy/bond length) the form of the potential is pretty close to quadratic. SPOILER ALERT: this is why the potential you showed in your question is used. It is quadratic, and describes the bond energy fairly well, BUT, only close to the minimum. As you get further away the bond energy clearly is not well modeled by a harmonic potential (quadratic)
Thus, if you try to break a bond with in classical MD, you will do a terrible job of accurately modelling it. QM can't even do it well without using un-Restricted or Restricted open shell approaches!
As I alluded to in my comment "non-classical" methods can be used, for instance ab-initio MD which incorporated QM, thus allowing the bond to break. Also, ReaxFF which is closer to a classical FF uses bond-orders rather than actual bonds to describe molecules. I am not experienced with either of these, so I simply mention them and supply these two links for your further reading if you are interested:
ab-initio ab-initio Wiki
reaxFF ReaxFF Wiki
You can do something with Lennard-Jones potentials instead of harmonic ones. I know this has been done for Coarse-Grained MD simulations using Go-Martini models as in https://pubs.acs.org/doi/abs/10.1021/acs.jctc.6b00986
Now this needs more testing, I have done some force-probe simulations with it but nothing too serious.
