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.

  • 2
    $\begingroup$ To calculate a dissociation energy, don't you just need the energy at the equilibrium, and then the energies of the constituents into which the system dissociates? The problem of near-degenercies, for example near the dissociation of the N$_2$ molecule, means one needs a multi-reference method to calculate the potential energy surface, but a single-reference method should be fine for the N$_2$ dissociation energy right? A single-reference method is fine for the equilibrium of N$_2$ and and fine for the N atoms into which N$_2$ dissociates. What molecule do they want to simulate? $\endgroup$ Commented May 23, 2020 at 1:29
  • 2
    $\begingroup$ The simplest case of an aryl R$_2$P-H would be something with two aromatic rings and P-H? Minimum of 12 carbon atoms + 24H + P + H, meaning 38 atoms? I think it would not hurt to ask the experimental colleague how many atoms, but I have a rough idea of how to answer. $\endgroup$ Commented May 23, 2020 at 2:03
  • 1
    $\begingroup$ How exact of an answer do they want? $\endgroup$ Commented May 24, 2020 at 20:15
  • 3
    $\begingroup$ If this is the case, more simple DFT methods can be used. Unless I am interpreting the question incorrectly. However, a DLPNO CCSD(T) cc-pVTZ combination is probably best for that system size. $\endgroup$ Commented May 25, 2020 at 6:12
  • 1
    $\begingroup$ @Greg It is definitely overkill. Personally I like using B3LYP 6-311++G**. My comment was also made before OP added the DFT part of the question. $\endgroup$ Commented May 29, 2020 at 2:14

3 Answers 3


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.
  • 1
    $\begingroup$ Great advice, especially on CCSD(T) flavors $\endgroup$
    – Greg
    Commented Jun 21, 2020 at 14:06
  • $\begingroup$ Are there any benchmark papers that support the claim "(the DLPNO-CCSD(T) method of) ORCA is extremely slow (compared to PNO-LCCSD(T) and LNO-CCSD(T))"? If yes, what's the applicability range of this statement (for example, is it only true for QZ basis sets and above)? $\endgroup$
    – wzkchem5
    Commented Jan 27, 2023 at 12:47
  • $\begingroup$ @wzkchem5 I didn't say that "the DLPNO-CCSD(T) method of ORCA is extremely slow", I simply said that "ORCA is extremely slow" compared to Molpro and MRCC. I think that you and I both know that papers comparing timings in different programs are extremely rare, but MOLPRO is known for speed (especially for single-CPU calculations). I haven't used ORCA since early 2019 though, and a lot could have changed in the last 4 years. $\endgroup$ Commented Jan 27, 2023 at 17:56
  • $\begingroup$ I agree that Molpro is very fast for multireference calculations, especially MRCISD. And MRCC is of course very fast for high order coupled cluster. I guess you probably use high order coupled cluster and multireference methods much more than the average user does, which may give you the impression that ORCA is extremely slow. To how much extent this carries over to local CCSD(T) is not so clear, since ORCA put much more effort in optimizing DLPNO methods compared to multireference correlation and high order coupled cluster methods. $\endgroup$
    – wzkchem5
    Commented Jan 27, 2023 at 19:48
  • $\begingroup$ @wzkchem5 In February 2019 my ORCA calculation on (BN)2-pyrene in aug-cc-pVDZ said "TOTAL RUN TIME: 90 days 9 hours 16 minutes 24 seconds 410 msec" to do EOM-CCSD whereas in MOLPRO it took 7.5 hours to do CCSD(T) in cc-pVDZ, and 204 hours to do CCSD(T) in cc-pVTZ. The SCF in ORCA took 10613.407 sec in aug-cc-pVDZ and MOLPRO took 19.82 sec in cc-pVDZ and 472.66 sec in cc-pVTZ (all are wall times, on the same machine). $\endgroup$ Commented Jan 27, 2023 at 20:09

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.


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

enter image description here

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


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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .