The strength of a chemical bond can be defined as the change in enthalpy when a bond is homolytically broken into two radicals. Bond dissociation energy is thus the measure of the strength of a bond and, while it is temperature dependent, it can be calculated by DFT or ab initio methods.
In this tutorial we’re going to calculate the C-H bond dissociation energy (BDE) in the methane (CH4) molecule, one of the strongest aliphatic bonds in organic chemistry. To that end, we start from the geometry optimization and frequencies calculation of CH4 and the radicals product of dissociation in equation 1.
CH4 —> CH3. + H. (equation 1)
For the sake of reproducibility, please know that I’m using the wB97XD/cc-pVDZ level of theory, but a suitable reference or benchmark should be consulted to choose the ideal level of theory for your molecules under study. Bear also in mind that the multiplicity changes from singlet to doublet in the radical products. Since one of our fragments is just a Hydrogen atom optimization and frequency analysis becomes superfluous, but if we’re analyzing a molecular fragment then those steps cannot be omitted.
Now, since the definition of BDE is of enthalpic nature, the not-corrected electronic energies could be a good approximation, but conceptually wrong. We should look into the thermochemistry section of the output, and thus this is why the freq keyword is paramount to these calculations, more specifically to the line labeled as “Sum of electronic and thermal Enthalpies=”. So, for the calculations made at 0 K we have the following values:
| Fragment | H (0K) [Hartrees] |
| CH4 | -40.453173 |
| CH3. | -39.78826 |
| H. | -0.501881036 |
| ΣH(prods) – ΣH(reactants) | 0.163031964 |
The resulting value of 0.163 a.u. translates to 102.304 kcal/mol (or 428.039 kJ/mol) while the reported value at 298 K is 103.011 kcal/mol (431 kJ/mol), this is a very good agreement for such a strong bond, but the agreement will vary with the sizes of the fragments and the level of theory.
Now, suppose you want to calculate BDE at a given temperature, say 300 K, you can just type the keyword temperature=300 in the route section to set the value of T at which the thermochemistry analysis is performed after the frequencies calculations. With this value the obtained values shift in the following way:
| Fragment | H (300K) [Hartrees] |
| CH4 | -40.453148 |
| CH3. | -39.788232 |
| H. | -0.501881036 |
| ΣH300(prods) – ΣH300(reactants) | 0.163034964 |
This seemingly similar value corresponds to 102.306 kcal/mol (428.048 kJ/mol), a slight change but nonetheless the T value can be incorporated into these calculations.
Now, these enthalpy values, usually reported in kcal/mol (or alternatively in kJ/mol) refer to the breakage of a given bond per mole of the reagent under study, but the energy of a single bond should be reported in eV by using the conversion factor of 23.060 kcal/mol (96.485 kJ/mol).
Dr. Joaquin Barroso's Blog 
Dear Dr. Joaquin,
In my opinion, the energy, the enthapy and the Gibbs free enthapy of hydrogen atom at 0K is -0.500000 a. u without any any calculation.
It is the definition of the Hartree unit of energy.
At 0K, E(H.) = H(H.) = G(H.) = -1/2=-0.5 Hartree
Sincerely,
S. M. Mekelleche
You’re right, however that is true at the exact calculation, these were made under the DFT scheme with a given basis set that are not the exact solutions to Schrödinger’s equation for the H atom. Your appreciation is fundamentally correct. What I’m showing here is an approximation that can be extrapolated to any bond breaking.
Thanks for comment and congrats on having such a great insight!
Would love to see the G16 input code for this.
Well the input is quite straightforward, just use the opt and freq (optionally freq=noraman to speed up the calculation) keywords and search for the corresponding lines in the output as pointed out in the post.
Thank you for reading
Dear Dr. Joaquim,
The chose of the level of theory is certainly necessary, but how to do it “correctly”?
Using several combinations of basis-set/pseudopotentials/functionals and comparing with experimental value, I think it is not the best approach. For learning, it is ok but for production and bigger systems, kind of we lose the best characteristic of simulations that it is to be able to predict such interesting properties.
Sincerely,
Camps
Dear Camps,
Your question is actually THE question. Knowing the way each functional was created, calibrated, and tested, helps a lot to discern what are the reaches of each. A thorough literature search is always a good idea to find out which functionals have been successful in describing what kind of properties. There are many benchmarks out there. Also, if you’re trying to describe something not fully covered in the literature, one may try to find a benchmark that covers a related property that you could unambiguously relate to the one of interest to your own research.
Jaguar by Schrödinger INC. had a module to suggest the best level of theory for your needs, based on your system and computational resources; I think we should have an app for that, and I’m sure there has to be an AI tool out there for making such suggestions.
Thank you for reading and commenting.
Thank you Dr Barroso for this info. There is a small typo error in converting kcal/mol to kJ/mol in the first calculation.
Apologies it’s my mistake, there is no typo. I am reading using my phone and the numbers were blurred for sometime before clearing.
Thank you for noticing it! There was indeed a typo, more like a conversion error on my part, but it was already taken care of.
Thanks for reading!
Did you mean E of reactants? I see reagents.
Thank you for noticing it, I’ll fix it right away.