What are the types of bond orders?

Question

I am studying the interaction of a metal with a nanostructure. As we can not trust graphical software about the bond formation, I decided to analyze them using the software Multiwfn.

It happens that there is a Zoo of several theories about how to calculate the bond orders:

• Mayer bond order analysis (Chem. Phys. Lett, 97, 270 (1983))
• Standard multi-center bond order (Struct. Chem., 1, 423 (1990))
• Multi-center bond order in natural atomic orbital (NAO) basis
• Wiberg bond order analysis in Löwdin orthogonalized basis (Tetrahedron, 24, 1083 (1968), J. Mol. Struct. (THEOCHEM), 870, 1 (2008))
• Mulliken bond order analysis
• Orbital occupancy-perturbed Mayer bond order (J. Chem. Theory Comput., 8, 908 (2012))
• Fuzzy bond order (Chem. Phys. Lett., 383, 368 (2004))
• Laplacian bond order (J. Phys. Chem. A, 117, 3100 (2013))
• Intrinsic bond strength index (J. Phys. Chem. A, 124, 1850 (2020))
• AV1245 index (approximate multi-center bond order for large rings) (Phys. Chem. Chem. Phys., 18, 11839 (2016))
• DDEC6 atomic population analysis (RSC Adv. 7, 45552 (2017))

What method/theory do you suggest to be used?

Answer 1

Wiberg (1968)

Let's start with the "classic" bond order paper by Ken Wiberg (born in 1927 and still alive!).

The Wiberg Bond Index (WBI) between fragments A and B of AB is calculated as follows:

\begin{equation} \tag{1} W_{AB} \equiv \sum_{\mu \in A}\sum_{\nu \in B}D_{\mu \nu}^2~ , \end{equation}

where $D$ is the following density matrix:

\begin{equation} \tag{2} D_{\mu \nu} \equiv \sum_{i\in \alpha }C_{\mu i}C_{\nu i} + \sum_{i\in \beta }C_{\mu i}C_{\nu i}~, \end{equation}

where $\mu$ and $\nu$ are atomic orbitals for fragments A and B respectively, and $\alpha$ and $\beta$ denote spins, and $C$ is the matrix of atomic orbital coefficients in the LCAO formalism.

Cons:

• It is old and primitive. Easily you could take into account the overlap matrix $S$ and the spin density matrix $Q$ to get the Mayer Bond Index (MBI) or an analogous WBI.
• It relies on wavefunctions, which Walter Kohn has said is not legitimate for large systems.
• While this bond order can easily be decomposed into $\sigma$, $\pi$ and $\delta$ bonding components for linear molecules, some more work may be necessary for more general molecules.

Pros:

• Despite the MBI being objectively more sophisticated, a detailed 2008 study of both MBI and WBI concluded with a recommendation to use WBI and cited two reasons:
  • It was quite stable with respect to basis set size (see Fig. 1) and this concluding quote "For the description of bond orders the preference is given to WBI calculated with NAO, because these indices are stable enough with respect to basis set variation."
  • Their WBI results were close to what we might expect: "The valences calculated using the WBI values taking into account both the covalent and ionic components of the chemical bonds are close to the formal atomic valences of main-group elements."
• It's very simple, easy to calculate, and therefore likely available in most software.

Answer 2

Laplacian Bond Order

This method is an extension of the QTAIM (Quantum Theory of Atoms In Molecules) concept of using the Laplacian of the electron density $\nabla^2\rho$ to characterize bonding. Standard QTAIM just uses the value of $\nabla^2\rho$ at bond critical points (minimum electron density along the bond, max along orthogonal directions), but these values are unintuitive and often don't accurately describe the strength or character (covalent/noncovalent) of the bond.

To make this idea more accurate and the values more in line with chemical intuition, the Laplacian Bond Order uses a weighted integral of the Laplacian:

$$\text{LBO}_{AB}=-10\int_{\nabla^2\rho<0}w_A(r)w_B(r)\nabla^2\rho(r)dr$$

The weight for any given atom is 1 close to it and decays smoothly to 0 as it approaches another atom. At any point in space, the weight of all atoms has to add up to 1. This weight function could be arbitrary (provided it meets these criteria), but LBO uses the so called Becke weighting function.

Integrating over the negative region of the Laplacian is done to specifically capture covalent bonding. The $-10$ was chosen ad hoc to make the magnitude fit chemical intuition about covalent bonds.

Pros:

• Matches chemical intuition for a wide range of systems.
• Well correlated with bond dissociation energy (at least for organic molecules and diatomics) and bond frequency (at least for CO in metal complexes).
• Tracks well with structural changes (e.g. weaker bonds in strained rings, weakening during double bond rotation or diatomics dissociation).

Cons:

• Basis set dependent, should probably use at least aug-cc-pvdz.
• Can only capture covalent interactions (though this can be an advantage in some cases).
• Currently only implemented for two body interactions (though in principle this shouldn't be difficult to extend)

Answer 3

How to test bond order methods for chemical consistency

Some of the bond order methods that claim to work do not give consistent results across different SZ values of a spin multiplet or different levels of theory (i.e., basis sets and exchange-correlation functionals). Consequently, there have been several published methods claimed to compute bond orders that do not work. Tests of this kind were performed in the article referenced below for the oxygen molecule in different spin states (singlet, triplet, and quintet) at different bond lengths and for different exchange-correlation approaches (e.g., DFT, CCSD, SAC-CI) and basis sets. For a fixed bond length and spin multiplet (e.g., triplet spin state at 200 pm bond length), the SZ=S and SZ=0 molecules are almost energy degenerate (except for a tiny spin-orbit coupling energy), have nearly identical electron density distributions, and therefore should have similar computed bond orders. However, many of the existing bond order methods fail this simple test for chemical consistency, because they give huge differences (1.5-2 bonds difference in some cases) between nearly chemical equivalent states.

see this article (open access): T. A. Manz, “Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders,” RSC Advances, 7 (2017) 45552-45581 (open access) DOI: 10.1039/c7ra07400j

Additionally, one should require the bond order to work well for both molecular and periodic materials, including materials with localized bonding electrons (e.g., insulators or small molecules) as well as ones with highly delocalized bonding electrons (e.g., metallic conductors). A new bond order definition satisfying those properties was introduced in the article cited above.

As another test of bond order definitions, one could start with the smallest molecules (i.e., diatomic molecules) and ask for which bond order definitions have worked well enough to compute bond orders for a large number of diatomic bond orders. Only one quantum-mechanically computed bond order definition has ever been applied in a published systematic study of bond orders for a large number of diatomic molecules. See the article below:

T. Chen and T. A. Manz, “Bond orders of the diatomic molecules,” RSC Advances, 9 (2019) 17072-17092 (open access) DOI:10.1039/c9ra00974d

Answer 4

It depends on what you want from it. If you're an organic chemist, you'd be mostly interested in the type of bond (single, double, triple, aromatic), and for that the bond order is most relevant (but tricky to compute). I sent a paper away earlier this year which deals with the choice of bond order analysis, and will post a link to it when it will be online. But if done properly, you can see e.g. the change from aromatic bond in benzene to triple bond in acetylene if you dissociate benzene into its three components $\ce{H-C≡C-H}$.

If you're however interested in the interaction of a metal with a nanostructure, then probably you're more interested in the energy. In 2007, I published this paper on the metal-ligand bonding in metallocenes (Inorg. Chim. Acta 2007, 360, 179),using the energy decomposition analysis (EDA) as available in the ADF program, using fragments. I.e. you prepare two fragments: (i) metal, (ii) ligands, and compute the change in energy if you bring them together. This will work as well for your system of metal/nanostructure, by preparing a fragment of the metal, and one for the nanostructure. In a followup paper on the metallocenes we thus showed that V(0) binds very strongly to a super[6]cyclophane cage (belonging to "Predicting new molecules by quantum chemical methods" special issue).

Of course you could also compute the intrinsic force constant for that particular bond, with all kind of if's and but's. But that would be opening up a whole new discussion.

The paper I mentioned above can be read here.