How to calculate homolytic bond dissociation energies?

Question

An experimental colleague asked me how hard it would be to calculate homolytic bond-dissociation energies for different phosphonates which are involved in a Hydrophosphination. The compounds include aryl and alkyl phosphonates (R2P=O), and the oxygen can be replaced with BH3 (R2P-BH3).

Maybe I'm overthinking this but isn't it difficult to calculate BDE because at the dissociation limit the $\sigma$ and $\sigma^*$ become degenerate, and more than one electronic configuration can describe the system (multi-configurational), e.g. $(\sigma)^2$, $(\sigma^*)^2$? Therefore, the correct wavefunction should be sufficiently flexible to treat both configurations on the same footing, which single-references like Hartree-Fock, and DFT cannot. However, multi-configuration techniques like CASPT2/CASSCF are complicated and I'd probably not be able to easily instruct a colleague how to use them easily.

Besides CASPT2/CASSCF, what other methods for calculating homolytic BDEs are good and are fairly "black-box"?

They just want a trend, so the absolute value doesn't matter as long as the model is useful for getting relative trends correct. Blackbox methods like DFT are preferred.

Answer 1

Multi-reference or single-reference?

While it is appreciated that near dissociation there will be a near-degeneracy of orbitals, requiring a multi-reference treatment, fortunately we don't have to worry about this when calculating dissociation energies. For example, the N$_2%$ molecule has profound multi-reference character as you approach the dissociation limit, which is why single-reference methods like CCSD(T) would fail miserably in calculating the potential energy curve. However you do not need the potential energy curve here, you just need a good energy of the overall molecule at its single-reference equilibrium, and an "equally good" energy (hopefully with good error cancellation) for the single-reference fragments after dissociation. Single-reference methods tend to describe such small aryl and alkyl phosphonates quite well at equilibrium.

Should you use DFT?

DFT (every single flavor) is notorious for giving bad results for computations outside of what the functional was optimized for. B3LYP can be used as a black-box, and is cheaper than most wavefunction based methods, but B3LYP is not likely going to calculate a dissociation energy correctly to within about 3 kcal/mol (see the figure here). If your dissociation energies are just a few kcal/mol, then your noise will be roughly the same size as your signal, and there's no easy way to systematically improve a B3LYP energy. For these reasons, if the system has only dozens of atoms, I would recommend wavefunction based methods, and I would resort to DFT for something like this only when reaching 100s of atoms.

I asked a few times in the comments how many atoms we're dealing with, and after that all I have available to work with is that we have at most some (aryl)$_2$P-BH$_3$ which has about 30 atoms if aryl=phenyl and about 40 atoms if aryl=naphthyl, and about 100 atoms if the aryl group is derived from heptacene. On the largest end of this spectrum, local-correlation wavefunction methods are still within reach, and on the lower end of this spectrum it is not too challenging to do without local-correlation techniques (meaning calculations would be even more black-box and more accurate).

What wavefunction based method should you use?

CCSD(T) is the ultimate "gold standard", black-box, method for molecules where it would not be too costly (such as your aryl phosphonates with only a few dozen atoms). You can even use the cc-pVDZ basis set without too big of a challenge, and your CCSD(T) calculations will greatly benefit from error cancellation (the coupled-cluster and cc-pVDZ errors for the total bound molecule and for the fragments, will both have similar sources of error, which will be eliminated when calculating the energy difference, resulting in dissociation energies that are much more accurate than you would otherwise imagine possible).

Local-correlation coupled-cluster type techniques can be used if you need to deal with 100s of atoms:

• DLPNO-CCSD(T) in ORCA might be the most famous, but compared to compared to the options I will give next, ORCA is extremely slow and less black-box in my opinion.
• PNO-LCCSD(T) in MOLPRO is essentially the same type of thing but more accurate, and much faster (and you might agree that MOLPRO is easier to use). This is only in MOLPRO 2019 and later though (and MOLPRO is not free).
• LNO-CCSD(T) in MRCC is essentially the same type of thing again, but might be the best option because:
  • MRCC is free (unlike MOLPRO).
  • MRCC is very well-maintained (with frequent new releases, the latest being in 2020).
  • MRCC is much easier to use (in my experience) than ORCA, and often even MOLPRO.
  • The LNO-CCSD(T) method in MRCC is extremely well implemented: It's been used with a QZ basis set for a system with 1023 atoms and 44,712 AOs. I think it's unlikely for other free programs to compete in both accuracy, speed, and ease-of-use.

Answer 2

BDEs are calculated by the below procedure:

1. Calculate initial energy

2. Perform homolytic bond cleavage and separate fragments

3. Calculate energies of the fragments, add the energies together

4. Calculate BDE by comparing the fragment energies to the initial energy

The level of theory and basis set is dependent on how accurate you want the results to be. B3LYP 6-311++G** is typically a good compromise between accuracy and speed for this system size.

Answer 3

One thing the other answers haven't mentioned is the zero point and thermal corrections to the BDE.

As mentioned in the Wikipedia Morse potential article, a geometry optimization will take you, by definition to the equilibrium bond length ($r_e$). If you use this as the low energy state, you're calculating $D_e$.

The problem is that a quantum harmonic oscillator always has zero point energy even at 0K. So a better answer requires the zero-point vibrational energy. Usually a vibrational / frequency calculation will also produce the thermal corrections, e.g. to 300K since there will be some thermal occupation of higher vibrational levels.

I've seen a few papers from groups that specialize in accurate thermochemistry, and their procedures are roughly:

• Optimize the geometry, e.g. B3LYP-D3BJ / def2-TZVP
• Calculate vibrational levels (for ZVPE and thermal correction) at the same level
• Calculate the energies (bottom of the well) with DLPNO-CCSD(T) / def2-QZVP

Examples:

• Efficient DLPNO–CCSD(T)-Based Estimation of Formation Enthalpies …
• Generalized Prediction of Enthalpies of Formation Using DLPNO‐CCSD(T) Ab Initio Calculations…

You ned the optimized geometry and DLPNO-CCSD(T) energies for the molecule and fragments. In your case, where the fragments are also molecules (unlike the diatomic in the figures) you'll probably also want the ZPE and thermal correction of the fragments.

I won't debate the merits of the various CCSD(T) schemes - I think Nike's answer already did that.