# Simulating breaking bonds in molecular dynamics

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?

• You cannot break bonds in classical MD. See ReaxFF. Commented May 17, 2020 at 10:22
• Ab-initio molecular dynamics (AIMD), ie using a quantum chemical/quantum mechanical description of the ensemble does it
– Greg
Commented Nov 23, 2020 at 23:06

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

• +1. Great effort went into this. Now if only there's a way to "tag" Stack Exchange user Frank Jensen to show him that you cited his book! Commented May 17, 2020 at 21:32

## Very short answer: No, classical molecular dynamics cannot break bonds.

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!

## There is also the subtle point that the quadratic we use for bonds will never actually break. The energy will increase to infinity without "breaking".

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

• That's a REALLY good picture, especially considering it was done in PowerPoint! I didn't even know PowerPoint could make such decent figures! Commented May 17, 2020 at 19:15
• Being a poor grad student for several years looks to have paid off :) Commented May 17, 2020 at 19:18
• One quibble - it's not simply a matter with QM of open shell calculations. To really do bond breaking, you need more than one Slater determinant, e.g. multi-reference methods, etc. Commented May 18, 2020 at 0:50
• @JamesFlash One thing neither of us discussed, is that bond breaking does not effect only the bond energy/forces. As you break a bond, the angles and torsions those two atoms are invovled in need to be accounted for (eventually turned off). Also, you have to start switching to applying non-bonded electrostatic and vdW forces between them. This also means, how do you figure out the partial charges on each atoms? reaxFF addresses all of these issues, often with empirical parameters. It is good to keep all these other issues in mind though. Commented May 18, 2020 at 22:44
• Thanks @CharlieCrown, your comments advanced my understaning of the process Commented May 19, 2020 at 9:36

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.

• Lennard-Jones potentials describe only van der Waals effects, i.e. not really bond breaking. Commented May 23, 2020 at 9:47