Computational strategies for modelling excited states in organometallic chemistry ^†

Recent decades have seen the emergence of a plethora of quantum mechanical methods to simulate and predict the properties and energies of excited states, with applications extending across a variety of chemical systems. To the untrained eye, such diversity might give the impression that any particular problem can be solved equally in many ways. However, navigating judiciously this complex array of methods requires an understanding of the theoretical basis and capabilities of each method in order to critically assess their suitability for the problem at hand. One must be able to tell alpacas from llamas, and this is not done by their hums.

In this chapter, we aim to cover the most significant aspects of electronic structure levels of theory (Section 1) and accompanying excited state methods (Section 2), offering a hierarchical classification whenever possible. Applications and current practical achievements in the context of excited state calculations of organometallic systems are discussed in Section 3. Careful attention is given to the approximations made and their impact on accuracy when treating organometallics. In doing so, we avoid diving into extensive mathematical formalism, but appropriate references are provided for interested readers. This chapter, therefore, serves as an introductory guide or roadmap to these methodologies. The reason we begin by discussing the various models of electronic structure theory is that all electronic and chemical properties—including excited states—are based on them. Therefore, the accuracy and quality of the computed properties are fundamentally dependent on the chosen models.

Quantum mechanics—or more generally, quantum field theory (QFT)—provides us a way of describing physical systems at the subatomic level, in principle, in a complete and exact way, provided we neglect gravitational interactions. This optimistic picture is, however, destroyed in many ways by mathematical or practical impossibilities. This forces us to approximate our model in different, and quite clever, ways. Our treatment of the description of chemical systems here will consider what is currently available as computational resources implemented in standard software packages, thus limiting the scope of this text. We will concentrate on a particle-based theory (i.e. ignoring the quantum field picture) which initially neglects any relativistic effects (we will come back to these in Section 1.8). In other words, we will use the Schrödinger equation, and assume the accompanying standard collection of postulates as the foundation of our model.¹  The electronic spin—essentially a relativistic property—will be introduced in an ad hoc manner, as it is indispensable for the correct description of atoms and molecules.

Since the HF method models the wavefunction as one Slater determinant, the variational principle optimises the energy subjected to this restriction. Such methods are called single-determinant methods, and the reader is referred to ref. 6 for a detailed description of the HF method.

The most straightforward way to introduce electronic correlation in the HF framework is by using configuration interaction (CI) methods. As discussed above, the ground state HF wavefunction is in the form of a single Slater determinant. Such a wavefunction consists of a set of lower energy occupied spin–orbitals and a set of higher energy virtual (i.e. unoccupied) spin–orbitals. By exchanging occupied and virtual orbitals in such a determinant, one can construct excited Slater determinants. Note that “excited” here refers to the form of the determinant, not to excited electronic states of the system. If only one such exchange is made, we call it a singles (S) determinant. Analogously, doubles (D), triples (T), quadruples (Q), and so on, determinants can be constructed depending on whether two, three, four, etc. exchanges are made, up to the number N of electrons in the system. The size of the basis set utilised to expand the orbitals determines the number of virtual orbitals present in the wavefunction, and thus, the number of possible excited determinants.

One can easily perceive that a full CI calculation is unfeasible but for the smallest systems. In fact, given 2K spin–orbitals for a system of N electrons, there are $( 2 K N )$ possible Slater determinants. This notation represents the number of ways to choose N spin–orbitals out of 2K without considering the order. Essentially, it calculates how many different ways we can arrange N electrons into 2K available spin–orbitals. As both N and 2K increase, the number of possible configurations increases rapidly, in a combinatorial explosion.

A common approximation for chemical purposes is to consider only the excited determinants arising from valence electrons, known as the frozen core approximation. The absolute error of introducing such a procedure is not small. However, since core electrons respond weakly to chemical changes, this leads to good relative energies due to error cancelation.

In addition to the frozen core approximation, CI calculations are typically performed using a truncated expansion, which limits the inclusion of excited determinants based on their order. Including only single excitations (CIS) does not improve the HF wavefunction because integrals involving singly excited determinants and the HF reference are zero because of orthogonality. Thus, the smallest improvement can be achieved by considering doubly excited determinants, either alone (CID) or with singles (CISD). The latter offers a slight improvement because the overlap between singles and doubles is non-zero.

Further improvements can be made by including triply (CISDT) and quadruply (CISDTQ) excited determinants. The latter often approaches sufficiently close the full CI energy for many systems, but it becomes extremely costly for large molecules and basis sets. While singly excited determinants contribute minimally to the electronic energy of a system, they are crucial for calculating molecular properties, such as the ones we will consider in Section 2.

In the CI method described above, the spin–orbitals composing each Φ _i Slater determinant are those obtained as solutions to the HF equations. Further variational freedom can be introduced to the various truncated CI methods by allowing the coefficients of each orbital expansion in the basis set to be optimised variationally alongside the CI coefficients. This approach, known as the multiconfiguration self-consistent-field (MCSCF) method, is clearly more expensive than standard CI methods, but the wavefunction converges better and faster to the variational limit.

An MCSCF wavefunction can more accurately describe systems that require different electronic configurations for a balanced representation, such as organometallic molecules with heavier metals, where spin–orbit coupling mixes different spin configurations into the electronic states. In essence, MCSCF methods are multideterminant, as they do not rely on a single HF determinant as the reference for excited determinants in the expansion.

One challenge with MCSCF methods is the selection of appropriate configurations to include in the calculations. A frequently used approach is the complete active space self-consistent-field (CASSCF) method, which partitions the molecular orbitals into active and inactive spaces. A full CI expansion is then performed in the active space, including all configurations in the MCSCF calculation. The active space typically consists of a few of the highest energy occupied orbitals and lowest energy virtual orbitals from a HF reference. Further partitioning of the active space into domains with different excitation schemes is possible, known as the restricted active space self-consistent-field (RASSCF) method.

These CI methods all consider excitations from a single HF reference determinant. However, one can also consider excitations from the determinants contained in an MCSCF wavefunction. For example, CISD can be extended to include all singles and doubles excitations from all determinants contained in an MCSCF wavefunction. Such methods are termed multireference configuration interaction (MRCI). Of all the approaches considered so far, excluding the exact full CI, MRCI offers the highest quality wavefunction, but it comes at an extremely high computational cost.

One significant problem with any truncated CI method, including the analogous MCSCF and MRCI approaches, is their lack of size consistency and extensivity. Size consistency refers to a method’s ability to correctly describe the energy of a collection of non-interacting particles. A size consistent method should predict that the energy of N non-interacting particles is exactly N times the energy of a single particle. Size extensivity, on the other hand, refers to systems of interacting particles and requires that the method correctly scales the energy in proportion to the number of particles.

Although the full CI is both size consistent and extensive, as an exact theory should be, this is not the case for any truncated CI or MCSCF method. Consequently, these methods are not suitable for describing macroscopic effects, both due to computational cost and size inconsistency, or for large molecular systems, as their accuracy deteriorates with an increasing number of particles (size extensivity). For a detailed explanation of this, see ref. 6.

Truncated CI methods offer a conceptually straightforward approach to correcting the HF wavefunction and can be extended to multireference descriptions. However, their lack of size consistency and extensivity restricts their utility. This is because it is often desired to compare energies and properties between chemical systems that differ in number of atoms, for example, when following a reaction mechanism or the dissociation of a molecule. Since the full CI equation provides an exact description, it is desirable to devise approximations that, like the full CI method, do not suffer from these limitations.

This goal is indeed achieved through coupled-pair theories, with the gold standard being coupled-cluster (CC) (also known as coupled-pair many-electron theory). All CC methods are both size consistent and extensive. However, they are not variational, meaning that CC energies can sometimes be lower than the exact energy, and there is no definitive upper or lower bound against which to measure their accuracy.

Additional pair theory methods, such as the independent electron pair approximation (IEPA) and the coupled electron pair approximation (CEPA)—a simplification of the CC approach—also exist but are not widely applicable. One reason is their relatively crude approximation and the fact that their energies are not invariant under unitary transformation of degenerate orbitals. More details on these methods can be found in ref. 6.

As it happens with the CI methods, one can extend the CC methods to multireference wavefunctions by using, for example, an MCSCF reference to construct an MRCC model. However, due to the complex nature of the CC expansion, this approach is less straightforward than MRCI. CC methods are considered the most accurate and advanced of the post-HF approaches, and many other levels of theory, including density functional theory (DFT) methods, are usually benchmarked against them.

CI and CC methods essentially rely on different ansätze to expand the wavefunction and then approximate the exact solution in a systematic way. The full expansions for both approaches are exact in the non-relativistic BO approximation framework. However, the approximations arising from the CC methods, unlike those from CI methods, are size consistent and extensive, though not variational. Another class of methods for introducing correlation corrections to IPM wavefunctions, which share these same properties with CC methods, is based on many-body perturbation theory (MBPT). We will briefly cover the basic concept of this general theory, provide an understanding of the most popular version currently in use—namely, Møller–Plesset perturbation theory (MPPT)—, and discuss some of the multi-reference extensions to this approach. In doing so, we will refrain from entering the rather intricate mathematical aspects, and again, the interested reader can find these details in the cited sources.

The first-order energy correction in MPPT only corrects the $H 0$ energy, yielding the HF energy. Therefore, the first correction for correlation energy occurs on the second-order energy correction of the MPPT theory. If we designate the MPPT equations restricted to corrections up to order n by MPn, we have MP2 as the lowest level correction. The second order correction to the energy can be written in terms of the reference wavefunction and doubly excited determinants thereof.^6,7  The first-order correction to the wave-function also contains only contributions of the doubly excited determinants. Since knowledge of this correction allows for the calculation of the third-order energy correction, MP3 is readily available and is not much more expensive than MP2.

Higher-order corrections to energy require higher than first-order corrections to the wavefunction. Knowledge of the second-order correction to the wavefunction allows for MP4 and MP5 energy corrections, which include contributions of up to quadruply excited determinants. A similar pattern to CC and CI calculations arises here: more and more configurations are introduced to increase the accuracy of the model. While MP4 has a computational cost comparable to CISD, MP5 is already too costly for large systems.

Although MPPT methods offer an interesting alternative, their convergence is highly dependent on the reference wavefunction. If HF description of the system is not adequate, the convergence of MPPT energies can be erratic, with MP3 sometimes yielding poorer results than MP2, for example. A more detailed discussion about these convergence issues can be found in ref. 7–9.

For organometallic systems, the presence of heavy atoms often necessitates a multireference wavefunction as a starting point, posing particular challenges for MPPT methods due to their convergence issues. As with CI and CC methods, MPPT can be extended to the multireference scenario by employing an MCSCF wavefunction as the reference. The implementation of such methods is complicated by the choice of partition for the Hamiltonian. When an MCSCF wavefunction is used as the reference, one obtains multireference perturbation MRPTn methods, where n denotes the order of perturbation correction applied. If the reference is a CASSCF wavefunction, the methods are termed CASPTn. Only CASPT2 and CASPT3 have been widely implemented so far.

Another approach based on a CASSCF reference is the N-electron valence state perturbation theory (NEVPTn) class of methods.^10–12  These differ from CASPTn methods in the choice of reference Hamiltonian and perturbation. A detailed treatment of these methods is beyond the scope of this text, but a comparison and perspective on these approaches can be found in ref. 13.

Before we conclude our journey through the methods for correlation correction to the HF model, it is worth calling attention to explicitly correlated methods. A wavefunction built from a Slater determinant is based on one-particle orbitals. However, the correlation effects arise from the two particle terms in the electronic Hamiltonian. Thus, the methods discussed above suffer from a fundamental incompatibility with the electronic correlation phenomena: Slater determinants fail to model the exact wavefunction at short interelectronic distances. A solution is to explicitly consider terms dependent on the distances of electron pairs. Such approaches are motivated by Hylleraas’ solution to the helium atom,¹⁴  and are called R12 or F12 methods. They depend on the use of Gaussian basis functions for two electrons, known as Gaussian geminals, and include extremely costly many-electron integrals. This poses a big challenge for the simulation of big systems, but efforts are being made to improve their scalability. Given the complexity of these methods and the current scarce use for bigger molecular systems, we do not enter into the details, but refer the reader to two excellent reviews on these methods.^15,16  R12/F12 approaches are currently implemented for MP, CI and CC wavefunction ansätze.

All approaches to the electronic problem discussed so far have involved a molecular orbital (MO) picture of the wavefunction. In fact, the HF method—the fundamental MO model—has only been improved by applying more complicated wavefunction expansions through CI and CC ansätze and/or by using perturbative methods to include electronic correlation. In the MO model, a molecule is viewed as composed by nuclei with a particular configuration (within the BO approximation) and molecular orbitals (and electronic density) delocalised throughout its entirety.

This picture is markedly different from the language chemists—particularly organic chemists—typically use to describe molecules and chemical reactions. Chemists often think about chemical bonds, their formation and cleavage, and the composition of molecules in terms of functional units or groups, relying fundamentally on a localised description of electrons. This perspective is preconised by the century-old concept of the chemical bond first introduced by Lewis, whereupon a pair of atoms is connected either covalently by an electron pair shared between only those atoms, or ionically by the electrostatic interaction arising from the localisation of the bonding electron pair in one of the two atomic centres.^17–19 

Valence bond (VB) methods, slightly older than MO methods, are much closer to this chemical description of molecular systems and are currently implemented in many quantum chemistry packages. The reason MO methods are more popular is that they are easier to implement and computationally cheaper to apply than VB methods. In the context of organometallic and transition metal systems, the application of VB methods is relatively rare. However, some recent studies have demonstrated their potential use in challenging systems such as metal–methyl (M–CH₃) complexes,²⁰  coinage metal halides,²¹  gold complexes featuring interanion coinage bonds,²²  and metal–metal bonded systems.^23–25  Therefore, it is worth dedicating a few paragraphs to VB methods here, as they will become more and more approachable with increasing computational capacity and new means to improve the scalability of computations. Our goal is to provide a broad overview of the theoretical aspects in modern VB theory and the practical differences compared to MO methods.

When describing the HF model above, we argued that the IPM wavefunction ansatz in the form of a Slater determinant correctly accounts for the antisymmetric property of fermionic systems, but this alone does not ensure its correct form. Since the wavefunction is a solution to the electronic Hamiltonian, it must be a representation for the symmetric groups to which the Hamiltonian belongs. This includes both point group symmetry for the nuclear configuration and permutation group symmetry for the electrons as independent particles. The theory that explains the necessity of the latter is well-established, but its complexity meant that it was not a primary concern for many scientists developing these early electronic structure methods.^26–29 

The key point here is that a properly constructed VB wavefunction is superior to a Slater determinant (the core of the HF method), mainly because it is a representation of the S _N permutation group. The second characteristic that makes VB wavefunctions a better representation is the departure from the double occupancy of orbitals. Let us examine these properties in more detail.

In current applications of CVBT, the HAOs are expanded into an atomic orbital basis set, and the expansion coefficients, as well as the different configurations (as in eqn (36)) are optimised using an SCF procedure. This method is termed VBSCF. If the HAO coefficients for each independent structure in an equation like eqn (36) are allowed to be independently optimised, a breathing orbital VB (BOVB) method is being applied. While VBSCF improves the description compared to HF by including “static correlation” only, BOVB acts as a multireference method, thus also including “dynamic correlation”, making it superior to VBSCF.

In the CVBT approach, the HAOs are essentially localised atomic orbitals. Another possibility is to use Coulson–Fischer HAOs. In this model, each HAO is allowed to possess a delocalisation “tail” by mixing a small amount of the neighbouring atoms’ HAOs. The wavefunction using these orbitals is then written as an HL covalent-like structure. This is the Coulson–Fischer (CF) ansatz. Notice that, by doing this, the orbitals are not orthogonal, which partially accounts for the high computational cost of such VB methods. For a molecule A–B, this would be eqn (33). The HAO coefficients are then optimised using an SCF procedure.

It is generally observed that one or a few of the spin-coupling terms dominate the SCVB description. Indeed, a few years earlier, Goddard et al. independently developed a similar approach, called generalized valence bond (GVB) theory.^34–38  Although completely equivalent to SCVB in its full form, which recently led VB practitioners to unify their names into a single framework known as spin-coupled generalized valence bond (SCGVB),³⁹  this wavefunction can be simplified in both its spatial and spin components.³⁹  An interesting application of SCGVB wavefunctions is in the development of the interference energy analysis (IEA) method,^40–43  which investigates the nature of the chemical bond as a manifestation of the quantum mechanical interference effect.

If the spin part is constructed so that covalent bonds are described by two electrons occupying two orbitals, one in each of the bonding atoms, coupled in a singlet pair, this spin function is called a perfect pairing spin function.⁴⁴  By enforcing strong orthogonality amongst orbitals from distinct electron pairs (but not those coupled in a bond or lone pair) and allowing core orbitals to be doubly occupied, we obtain what is called a perfect pairing GVB (GVB-PP) wavefunction. Although more computationally tractable than the full SCGVB wavefunction, such an approximation can lead to significant errors for some systems and must be used cautiously.

All VB methods can also be improved using CI, leading to methods such as VBCI, CASVB, and others available in various software packages. Although the implementation differs, the general idea is similar to that discussed for HF wavefunctions, allowing correlation corrections to be systematically added. A VB wavefunction as a starting point, however, is a much better reference than HF. In fact, the entire discussion about “static” and “dynamic” correlation can be viewed as an unfortunate artifact of considering the HF wavefunction as a reference. It can be shown that SCGVB wavefunctions account for what is defined as static correlation by respecting the correct permutational symmetry and spin eigenfunction conditions. Thus, HF lacks this type of correlation because it is simply an incorrect ansatz. What remains, usually termed dynamic correlation, is the true correlation effect, not included in the SCGVB approach because it is, like HF, an IPM for the many-electron system. A better definition than eqn (11) for the correlation energy is then obtained by replacing the HF energy E_HF/CBS by the corresponding full SCGVB energy E_SCGVB/CBS. More on this interesting discussion can be found in ref. 45–47 and more on VB methods can be found in ref. 25, 43, 48 and 49.

Up to this point, we have discussed wavefunction methods exclusively. There is, however, a class of methods that is conceptually different but has proven to be exceptionally effective for many computational tasks at a fraction of the cost. Instead of solving the Schrödinger equation for the full wavefunction of a system and its corresponding energies, these methods exploit the fact that there is a univocal correspondence between the energy of a quantum system and its one-electron density. This framework is known as density functional theory (DFT), and its foundation is based on the Hohenberg–Kohn theorems. An introduction to DFT can be found in ref. 50. A more detailed treatment is provided in ref. 51 whilst a more rigorous mathematical treatment can be found in ref. 52.

The Hohenberg–Kohn theorems form the foundation of DFT. The first of the two theorems states that the ground-state of a system (and therefore its energy) is a univocal functional of the density ρ( r ). The second theorem establishes a variational principle for the ground-state energy. According to this principle, given any approximate electron density for a system, the calculated energy through the density functional will be higher than the exact energy, unless the density is exactly the ground-state density ρ₀. These theorems can also be extended to excited states.

Since the electronic density depends only on the three spatial coordinates of a single electron, regardless of the system’s size, DFT offers a promising, exact, and remarkably simple and inexpensive alternative compared to wavefunction methods, which depend on 3N coordinates. The main issue is that the exact functional connecting the electronic density to the energy is still unknown. Therefore, only approximate methods are available within this framework.

The XC part of the KS functional holds all the approximations of the method, including electronic correlation. This is where most of the current efforts in the development of DFT methods are focused. The task of the user is to judiciously choose the appropriate DFT functional for the task at hand. The number of available functionals is bewildering, and discussing their intricacies in detail is beyond the scope of this chapter. Instead, we aim to provide a rough classification of these functionals and some conceptual foundations. We will also discuss some applications in the context of organometallic systems in Section 3.

Building on the Taylor expansion mentioned earlier, one can further improve the GGA approach by including the Laplacian of the density, ∇² ρ( r ). Functionals that incorporate this additional term, although not directly constructed by the simple addition of the Laplacian of the density, are known as meta-GGA functionals. These functionals typically perform better than LDA and GGA functionals for molecular systems, and more details on their implementation can be found in ref. 52.

Another class of functional that generally outperforms GGAs, and often surpasses meta-GGAs, are the hybrid or hyper-GGA functionals. These functionals leverage the fact that the exchange part of the XC energy can be computed exactly for a Slater determinant using the HF approach. One can thus use this parcel of the energy and only approximate the correlation part in the functional. Although this partition is a fictitious mathematical construct—the XC term is inherently a single entity—a careful partitioning scheme can yield effective XC functionals by mixing the HF exchange energy with LDA, GGA, or even meta-GGA XC energies. In fact, one of the most popular DFT functionals, B3LYP, is an example of this approach.

Two of the major inaccuracies of DFT calculations are the self-interaction error and the absence of long-range correlation effects, responsible for London’s dispersion forces.^53,54  Although the omission of long-range correlation is an inherent limitation of DFT, it can now be easily addressed by applying several established dispersion correction methods. Prominent among these is Grimme’s D3 method⁵⁵ —which can be applied using either the initial damping function, D3(0),⁵⁶  or the updated Becke–Johnson (BJ) damping function—or its more recent version D4.⁵⁷ 

An even higher level of improvement is achieved with double hybrid functionals. These functionals not only use information from occupied orbitals in the form of the HF exchange energy but also incorporate information from all the KS orbitals. Specifically, they include an MP2 component of the energy, obtained by applying the formalism described in Section 1.5 to the KS Slater determinant. These methods have a higher computational cost than the previous ones because of the extra computations required for the MP2 energy and the convergence difficulties associated with such a perturbative approach, as discussed in Section 1.5.

When treating excited-states via time-dependent DFT (TDDFT) methods (see Section 2.4), it has been observed that many functionals fail to accurately describe charge transfer and Rydberg states. This issue arises partly because DFT tends to artificially favour delocalised systems over localised ones. Additionally, the long-range behaviour of the electronic energy is generally misrepresented due to an imbalance between the density-based Coulomb (J[ρ( r )]) and exchange (HF part in hybrid functionals) components and the approximate XC potential. A solution to this problem is brought by the class of range-separated functionals. In this strategy, the electron–electron Coulomb operator for the exchange energy (only) is split into short-range and long-range components, using different approaches for each. Typically, a density exchange functional is used for the short-range part, while the exact HF part is used for the long-range part. The two domains are interpolated using an ω parameter, which can be tuned for specific systems by applying Koopman’s theorem.⁵⁸ 

Fig. 1 presents a schematic hierarchy of DFT functional classes, commonly referred to as Jacob’s ladder.⁵⁹ 

All of the theoretical treatment above assumed a non-relativistic scenario. This is, of course, an approximation. The fundamental problem with the non-relativistic treatment of a quantum system is that the Schrödinger equation is not invariant under a Lorentz transformation of space-time coordinates—a fundamental requirement of special relativity. The Dirac equation can be seen as a generalisation of the Schrödinger equation that does respect this property. It also accounts for the fact that the mass of a particle increases with increasing velocity. This implies that an electron orbiting a heavier nucleus will have a higher velocity, and therefore a higher mass, compared to one orbiting a lighter nucleus. The outcome is either a contraction or an expansion of the orbitals.

Another consequence of the Dirac equation is that the total angular momentum of an electronic system is no longer defined solely by the orbital angular momentum (L), but by a combination of this with the spin angular momentum (S), giving rise to a J term. In other words, the L and S operators do not commute with the Hamiltonian anymore, but J does. This is the spin–orbit coupling phenomenon. As a result, orbitals (and the entire system) do not have well-defined spin states, and the pure s, p, etc. non-relativistic orbitals now mix. In the presence of an external electromagnetic field, other relativistic effects also arise.

Additionally, the use of a Coulomb potential to describe the electric interactions is no longer valid in the relativistic realm since such interactions are instantaneous and thus violate the universal speed limit of light. The correct treatment of electric interactions would have to start from a quantum electrodynamics perspective, but modifications of the Coulomb potential can be made through approximation schemes, such as the inclusion of a Breit operator.

Here, we point out that relativistic corrections become increasingly important for systems containing heavier atoms. Since many organometallic systems contain metals beyond the fourth row of the periodic table, it is crucial to consider these effects. The simplest way to account for relativistic effects is by including relativistic effective core potentials (ECPs) to the basis set. Examples of such basis sets include the Los Alamos National Laboratory 2 double zeta (LANL2DZ)⁶⁰  and the Weigend and Ahlrichs def2 family of basis sets.⁶¹ 

Another relatively straightforward approach is to approximate a relativistic Hamiltonian using methods such as the zeroth-order regular approximation (ZORA) or the first-order regular approximation (FORA).^62–64  Both methods include spin–orbit coupling and scalar relativistic effects, which account for the mass dependency on velocity and the relativistic movement of electrons, known as the Darwin correction.

Perturbative operators have also been devised to include corrections to the non-relativistic solution.⁷  Since this strategy results in complicated equations, it is sometimes easier to directly solve the Dirac equation. Using spinors (relativistic spin–orbitals) in a Slater determinant, one can write a relativistic ansatz like the HF wavefunction for the Dirac equation, and this is known as the Dirac–Fock (DF) method.

The so-called two-component methods are also available. The reason for this name is that a relativistic wavefunction has four components, which can be separated into one pair called the large component, and another pair called the small component, which depend on one another. Two-component methods use unitary transformations to separate these components, making them (quasi-)independent and allowing them to be solved separately. Examples of this approach include the exact two-component (X2C)⁷  and Douglas–Kroll–Hess (DKH) methods.^65,66  A notable application of these methods for describing the spectroscopic properties of organometallic systems can be found in the ReSpect (Relativistic Spectroscopy) program, accessible at https://www.respectprogram.org/.

Having discussed the various levels of theory for solving the static electronic problem, we are now positioned to address the problem of excited states and their properties. These methods are sometimes referred to as dynamic methods because many are derived from the TDSE (eqn (1)). However, this nomenclature should not be confused with dynamics simulations such as molecular dynamics or quantum (adiabatic or non-adiabatic) dynamics.

Our focus here will be on methods that allow us to calculate energies and molecular properties resulting from the interaction of UV and visible light with matter, thus altering their vibronic (electronic coupled with vibrational) states. The class of time-dependent methods (and much of the theory discussed here) however, is much broader. With proper formulation, they can account for many types of time-dependent phenomena arising from magnetic and electric fields interacting with matter.

The case of electronic time-dependent methods is complicated by the fact that a general theory should consider the relativistic and quantum nature of both the photons and the molecular systems. Such a treatment falls within the realm of quantum electrodynamics or, more generally, QFT. Since routine methods at this level are not easily available, we will use an approximate, semi-classical theory. For exact electronic methods, we started from a non-relativistic perspective, although it is possible to correct for, or approximate many, relativistic effects, as discussed in Section 1.8. We will continue with this assumption for the electronic wavefunction.

Additionally, we will treat light classically, as oscillations of an electromagnetic field. Since the magnetic component is much smaller than the electric component, we will disregard its effect and use the dipole approximation (see below). As usual, we will be working within the BO approximation, although these methods can be expanded to account for non-Born–Oppenheimer effects.

Conceptually, the simplest way to consider excited states is by directly constructing their wavefunctions and proceeding analogously to the ground state calculations discussed above. If the excited state under consideration has a different (spatial and spin) symmetry than the ground state, this approach is straightforward. However, if it has the same symmetry of the ground state, the SCF or variational procedure will almost invariably return the ground state, as it is the lowest energy state of that symmetry. There are special techniques to obtain these excited states directly, but they are difficult and possess several technical problems.⁶⁷ 

It is often the case that excited states cannot be properly described by a single reference wavefunction, thus necessitating the use of MCSCF wavefunctions. Correlation can then be included by the various methods described above. Excitation energies can be calculated by taking the difference between the excited and ground state energies directly (ΔSCF approach). Molecular properties can also be calculated from various energy derivative techniques directly applied to the appropriate wavefunctions.

HF, MPn and CC type wavefunctions are not well-suited for calculating excited state wavefunctions in such a way because of their unsuitable references. We have already seen that CI methods straightforwardly give excited state wavefunctions as eigenstates of the CI matrix with higher energy eigenvalues. In fact, in most software packages, when HF is selected for excited states, what actually happens is that the CIS matrix is diagonalised to yield the excited states. However, care must be taken when choosing the CI level, as many excited states need two or more excited determinants for proper description. Even low-lying excited states of simple organic molecules may need at least CISDTQ wavefunctions for accurate description;^7,68  this is even more crucial for organometallic compounds.

Additionally, dynamic correlation can only be included in singly excited states by considering at least doubly excited determinants or by including perturbative corrections to the excited state wavefunctions. Thus, CIS excited states are described with a quality similar to HF ground states. This demanding level is often prohibitive except for small organic molecules. The methodology described in the next sections will offer more robust and cost-effective options for calculating excited state properties and wavefunctions. Nevertheless, CASPT2 or NEVPT2 models are deemed as gold standards for benchmarking excited state calculations whenever possible. It is worth noting that, in general, the higher the excited states considered, the higher the error.

A series of formalisms have been developed for obtaining dynamic properties of molecular systems from the TDSE. The following methods will be considered: (1) time-dependent perturbation theory (TDPT), which will then lead to the most popular excited state method, namely, linear response theory (LRT), sometimes also called random phase approximation (RPA); (2) variational methods, to obtain time-dependent wavefunctions and properties, which form the basis for time-dependent HF (TDHF) and time-dependent DFT (TDDFT) methods; (3) propagator methods, which circumvent the need to calculate time-dependent wavefunctions; (4) equation-of-motion (EOM) methods; and (5) the algebraic diagrammatic construction (ADC) family of methods.

For each formalism, it becomes evident that certain representations of the TDSE are more suitable than others, and therefore, we will define them here. As usual, we will not dive into the mathematical intricacies of the methods but focus on the essentials to give a superficial understanding of the concepts and critically assess each option.

Arguably, the most popular theory for treating time-dependent quantum phenomena is perturbation theory. We have encountered its principles when discussing MBPT above (Section 1.5). One can outline a time-dependent formalism in a remarkably similar way to the static case. We start by partitioning the total time-dependent Hamiltonian in two parts according to $H= H 0 +λ H 1$ ⁠. Here, $H 0$ is a time-independent Hamiltonian for which we know at least an approximate solution for its time evolution. $H 1$ is regarded as a perturbation, i.e. the effect of $H 1$ in the evolution of the system is much smaller than the effect of $H 0$ ⁠. The parameter λ is introduced as an ordering parameter, similar to its role in MBPT, and will ultimately be set to unity.

By adequately considering different forms of perturbation and using the basic perturbation formulas outlined above, various spectroscopic techniques can be simulated, and the wavefunction dynamics and other properties can be obtained. To illustrate, let us consider one-photon absorption and emission spectroscopies. We will provide a sketch of the theory, and more details can be found elsewhere.^5,68  We assume the BO approximation and describe our system as composed of two states: the initial state with Hamiltonian $H a$ ⁠, which is typically the ground state but not necessarily, and the final excited state after absorption of a photon with Hamiltonian $H b$ ⁠. In the dipolar approximation, the interaction of matter with light is considered to be independent of position, and the potential due to the light is $−μ⋅ϵ(t)$ ⁠, where μ is the molecular dipole moment and $ϵ(t)$ is the electric component of the radiation.

To apply the polarizability-based approach above, one only needs to know the set of ground and excited states wavefunctions. We recall once again that such an approach carries all the usual problems related to the choice of the reference wavefunction, and a complete set of Ψ _n is unattainable, demanding the use of an approximate basis set. LRT (or RPA) approaches are available in most software packages applied to a variety of electronic models such as HF, DFT, VB, CI, and CASSCF. LRT applied to CC methods are possible (e.g. LR-CC2), but the expressions derived in this section are not suitable for non-variational wavefunctions. The poles of the polarizability matrix are still a measure of excitation energies in such cases, but the expression must be formulated in terms of the derivative of the action Lagrangian (shown below). LR-CC formulations suffer from many problems, including lack of size extensivity and consistency. Notice that this procedure is primarily meant to obtain frequencies and intensities of transitions (properties), rather than attaining to calculate the wavefunctions and excited state dynamics, although some qualitative characterisation of the wavefunctions is possible through the consideration of the transition dipole moments.

The TDPT formalism above relies on an imperfect reference, in the sense that no exact solution to the unperturbed wavefunction exists. A variational theory of the TDSE is, in principle, exact and offers a more reliable way of approximating the time-dependent wavefunction and property calculations. Many possible formulations exist depending on the choice of parameters and approximation schemes. Many variational approaches are so-called local in time. This means that the wavefunction at a time t′ = t + Δt is variationally optimised taking the wavefunction at t as a reference. The reference can only be exact when t = 0, which is the ground state stationary wavefunction Ψ(x,0). Variational formulations that are global in time are possible through numerical methods.

We have come across the definition of an evolution operator or propagator above, as an operator which maps Ψ( x ,t = 0) → Ψ( x,t). After defining our perturbed Hamiltonian for light–matter interaction, we arrived at an expression which maps a wavefunction Ψ _i ( x ,−∞) into the first order correction to the wavefunction Ψ ⁽¹⁾( x,t) (eqn (63)). The quantity multiplying Ψ _i ( x ,−∞) is, therefore, a propagator. When discussing LRT, we have seen that this quantity leads to the correlation function, whose Fourier transform can be interpreted as the absorption spectrum. Furthermore, we have seen that this is equivalent to the response of the average polarizability in the frequency domain 〈α〉 _ω plus constants. From this, we can see that there is a close connection between the properties of a propagator and time-dependent molecular properties such as the polarizability, which leads to excitation energies. Writing the TDSE in the Heisenberg picture, one can in fact work directly with propagators, without any need to consider expansions of the time-dependent wavefunction. If we Taylor-expand the expectation value of the dipole moment operator 〈Ψ|μ|Ψ〉 in terms of the perturbation H₁ as defined in eqn (61), we arrive at a propagator-based expansion. It can be shown that the first order term is the frequency dependent polarizability which we are already familiar with. Propagators can be used in a broader way for a variety of time-dependent formulations, going beyond perturbation theory or variational approximations and even considering QFT formulations. For this, one considers the calculation of the general propagators as solution to their equation-of-motion (e.g. eqn (54)), thus avoiding the costly calculation of corrections to the wavefunction. In this context it is common to call the propagators Green’s functions, since they are essentially the solutions to differential equations involving a linear differential operator.⁸⁷  Such a treatment is beyond the scope of this text, but ref. 5, 7 and 69 discuss some propagator approaches in more detail.

Recall that, when discussing ΔSCF methods, we argued that CC wavefunctions are not suitable due to their HF reference. The EOM technique can be applied to CC wavefunctions with a suitable expansion of the excitation operators $O n †$ for a given cluster exponential operator. This is generally done by a similarity transform of the Hamiltonian such that it commutes with the approximate $O n †$ operator. To perform this, a CIS matrix is used as an initial guess for the singles and doubles space Hamiltonian. Currently, EOM-CCSD is the highest level of this approach since reliable implementations for EOM-CCSD(T) or EOM-CCSDT have not been developed.

EOM-CC methods are costly, but an alternative similarity transformed EOM (STEOM) approach exists.⁸⁸  In this version, a second similarity transform to the Hamiltonian further compresses the space to the CIS space only, thus decreasing the cost. Multireference formulations of this approach are also available. A word of caution regarding these methods is that they are not size consistent or extensive.

Except for the CC method discussed above, when applied to other reference wavefunctions at the usual level of approximation considered here, EOM, propagator, perturbative and variational methods all lead to the same working equations for discussing electronic excitations. For example, with a one-determinant closed-shell reference function and a single excitation manifold, all methods yield the standard LRT (or RPA) equations.⁶⁹ 

The last class of methods to be discussed is the so called algebraic diagrammatic construction (ADC) methods, amongst which the most popular approach is ADC(2).⁸⁹  These methods can effectively include correlation to the excited states through perturbative methods at a considerably lower cost than the CI or CC approaches discussed above. The origin of ADC methods lies in the propagators or Green’s function approach briefly discussed above, but current developments rely on the Liouvillian formalism similar to EOM methods.

Originally, for electronic transitions, the ADC method involved expanding the polarizability propagator in the perturbative regime in terms of diagrammatic perturbation theory. The MP partitioning of the Hamiltonian (Section 1.5) is used to construct the wavefunctions, and the perturbative order n to which this is carried out determines the ADC(n) level.

Because of this, such methods include correlation, but are not suitable for multireference systems, of which organometallic systems are an example. Only excitation properties and energies are available through this formulation, with no excited state wavefunction obtained. However, by using representations of the propagators (known as intermediate states representation), akin to what is obtained by the EOM approach, the ADC matrices can be built and used to calculate approximate wavefunctions. These representations, from ADC(2) onward, are also correlated references. ADC(2) offers good quality excitation spectra, comparable to the EOM-CCSD approach, but with a cheaper and easier implementation. For more details, see ref. 90.

One limitation of ADC lies in its single-reference framework, which restricts its applicability to systems exhibiting open-shell or multireference characters in their ground or excited electronic states. A comprehensive multireference ADC (MR-ADC) formulation has been proposed by Sokolov,⁹¹  achieved by integrating the Liouvillian formalism with MRPT. This formulation extends the conventional ADC theory to accommodate multiconfigurational reference wavefunctions. MR-ADC methods are comparable in computational expense to low-order multireference perturbation theories and are capable of calculating various electronic spectra, including UV-Vis, X-ray absorption, and UV/X-ray photoelectron spectra. MR-ADC calculations were, however, restricted to small molecules due to their inefficient spin–orbital implementation. Recent advancements by de Moura and Sokolov,⁹²  incorporating spin adaptation, density fitting, and automated code generation, promise to significantly enhance the efficiency of MR-ADC, potentially enabling calculations for systems with over 1500 molecular orbitals. A newly reported application of this improved MR-ADC methodology in organometallic chemistry involves the analysis of X-ray photoelectron spectra of Fe(CO)₅ and its photodissociation products (Fe(CO)₄, Fe(CO)₃) after excitation with 266 nm light (Fig. 2).⁹³  These calculations indicated that core-hole screening, spin–orbit coupling, and ligand-field splitting are all crucial factors in accurately reproducing the experimentally observed chemical shifts in the transient Fe 3p XPS spectra of iron carbonyl complexes.

The popularity of TDDFT methods for excited states cannot be overstated. Their ease of use even by non-experts, wide availability, and capacity of accounting for a major part of correlation effects at a fractional cost compared to post-HF methods make for a highly attractive method. This became clear in a 2010 survey, where, according to the statistics of the papers published from 2007 to 2010 in the field of excited states, the use of different methods is distributed as follows: TDDFT: 50%, single-configuration methods: 23%, and multiconfigurational methods: 26%.^94,95  Another, slightly more recent, survey shows a similar distribution with a more refined classification.⁹⁶  Other important reviews in the field address a range of methods and exemplary organometallic systems such as: (1) interpretation and accuracy of electronic absorption spectra, treatment of long-range charge transfer states and spin–orbit coupling effects, environmental effects, and dynamical and QM/MM methods;⁹⁷  (2) binary metal carbonyl photodissociation, non-adiabatic relaxation, Jahn–Teller and pseudo-Jahn–Teller effects, photoisomerization of transition metal complexes, and coupled cluster response theory for electronic spectroscopy;⁹⁸  (3) excited state dynamics of transition metal complexes, from the wavefunction to the dynamical level.⁹⁹ 

Despite the advantages and the remarkable success of DFT, we have repeatedly emphasized in this text the many flaws of this method. We have highlighted the necessity of benchmarking and carefully assessing the quality of the results obtained against higher quality computational methods and experiment. This is not a trivial task, and much is yet to be done, but the following examples demonstrate progress in this area. Latouche et al. ¹⁰⁰  benchmark DFT functionals and ECP basis sets against experimental data for absorption energies of Pt(ii) and Ir(iii) complexes (Fig. 3a). Paranthaman et al. ¹⁰¹  investigate the performance of DFT and relativistic ECP for Ru-based organometallic complexes (Fig. 3b). Garino & Salassa¹⁰²  review DFT and TDDFT methods for a variety of photochemically active organometallic compounds. Harvey¹⁰³  critically evaluates the use of DFT to treat transition metal chemistry in general. Niehaus et al. ¹⁰⁴  investigate the TDDFT description of charge transfer excited states, using a series of 17 pseudo-square planar platinum(ii) and pseudo-octahedral iridium(iii) complexes (Fig. 3c) that are known to feature quite different localization characteristics ranging from ligand-centred (LC) to metal-to-ligand charge transfer (MLCT) transitions. Maschietto et al. ¹⁰⁵  call attention to the existence of spurious low-lying excited states from TDDFT and offer a methodology to address the problem.

Two recent examples illustrate interesting contrasts between TDDFT and multiconfiguration methods. In the first, Costa et al. ¹⁰⁶  compare TDDFT and CASSCF/MS-CASPT2 applied to the description of excited states for the complex (CH₃)ReO₃ (Fig. 3d). Although both methods yield similar energies, a small difference in the characterisation of the second transition exists. This is a charge transfer from C and O to Re, but the participation of C, as described by CASSCF/MS-CASPT2, is much smaller than described by TDDFT. In their studies, the authors utilised the LB94 functional, which tries to correct for the wrong DFT asymptotic behaviour, but is found to be the worst amongst the tested functionals in describing the excitation energies.¹⁰⁶  In the second study, Escudero et al. ¹⁰⁷  compared RASPT2/RASSCF to range-separated/hybrid DFT methods in excited states of the Ru(ii) bipyridyl complex trans(Cl)-Ru(bpy)Cl₂(CO)₂ (bpy = bypyridyl, Fig. 3e). The authors discuss the difficulty in partitioning the RAS subspaces, and benchmark various types of TDDFT functionals including implicit solvent models. It is found that none of the functionals can optimally describe all the excited states simultaneously. However, the hybrid M06, B3LYP, and PBE0 functionals seem to be the best compromise to obtain a balanced description of the excited states of trans(Cl)-Ru(bpy)Cl₂(CO)₂, when comparing with the experimental spectrum.¹⁰⁷ 

Below, we will consider further applications and more critical assessment of DFT and multireference methods in different contexts. The intention is to highlight some of the complexities, successes, and challenges in different fields.

Spin-crossover (SCO) materials^108,109  are a fascinating area of research in organometallic chemistry, particularly in the context of first-row transition metals, where the electronic spin state can be reversibly switched between high-spin and low-spin configurations by external stimuli. Accurate prediction and understanding of SCO behaviour is crucial for the development of novel materials with potential applications in sensors and displays. Computational methods, in particular TDDFT and CASSCF/CASPT2, play a central role in these investigations. For example, Phung et al. ¹¹⁰  have shown that a combined CASPT2/CC approach, using high-quality CASPT2 with extensive correlation-consistent basis sets for valence correlation and low-cost CCSD(T) calculations with minimal basis sets, is efficient for spin-state energetics of a series of iron complexes (Fig. 4a), with errors estimated to be around 2 kcal mol⁻¹ in favour of the high spin states. Reimann & Kaupp¹¹¹  used another composite method, the CASPT2+δMRCI method, to calculate the SCO energy. For a series of [Fe(He)₆]ⁿ⁺ test complexes, the approach reproduces the full MRCISD+Q/CBS results with an accuracy of less than 0.04 eV. Cirera, Via-Nadal, & Ruiz¹¹²  presented a systematic study on the performance of different density functional methods to study SCO on first-row transition metal complexes (Fig. 4b). Amongst the tested functionals, the hybrid meta-GGA functional TPSSh with a triple-ζ basis set containing polarization functions for all atoms gives the best results for different metals and oxidation states, and its performance in predicting the correct ground state and energy window for the occurrence of SCO is quite satisfactory. Radoń¹¹³  also performed advanced benchmark studies showing, for example, that NEVPT2 performs worse than CASPT2, and that the double hybrid B2PLYP-D3 is able to provide a balanced description of spin state energies for the four studied iron complexes simultaneously (Fig. 4c). Such examples illustrate the continuous refinement of computational methods to better capture the complex phenomena underlying SCO processes.

Transitioning from SCO materials to photodevices, computational studies have provided important insights for the development of materials with optimal light emission properties. In the field of phosphorescent materials, which are central to the next generation of organic light emitting diodes (OLEDs), understanding the photophysics and relative position of the lowest energy triplet excitation, T₁, is crucial and has therefore been recently reviewed in various contexts. For example, Powell¹¹⁴  has reviewed theories of phosphorescence in cyclometallated complexes, highlighting the crucial role of both scalar relativistic effects and spin–orbit coupling in accurately modelling the phosphorescent properties. These effects, incorporated via TDDFT, provide quantitatively accurate predictions of radiative decay rates, which are crucial for the development of effective photodevices. Kumar & Escudero,¹¹⁵  in turn, used domain-based local pair natural orbital CC Theory (DLPNO-CCSD(T)) to calculate the phosphorescence energies of a variety of Pt(ii) complexes (Fig. 4d) with potential application as components in phosphorescent organic light-emitting diodes (PhOLEDs). This approach, which included a rigorous test of various relativistic effects and a comparison of several DFT functionals, revealed that M06HF performed the best. Finally, Ludowieg et al. ¹¹⁶  focused on the calculation of magnetic transition dipole moments and rotational strengths using the relativistic two-component TDDFT with ZORA. This approach was applied to evaluate the phosphorescence dissymmetry factors and lifetimes of iridium complexes based on N-heterocyclic carbenes and platinum helicenes (Fig. 4e), showing excellent agreement with experimental data.

Organometallic compounds have traditionally been celebrated for their catalytic ability, yet their potential in the realm of chemical biology, including for bioimaging applications, is increasingly being recognised.¹¹⁷  Despite initial reservations about their stability and cytotoxicity under physiological conditions, many compounds are proving to be invaluable due to their unique physicochemical properties such as robust chemical stability, structural diversity, and distinct photo- and electrochemical behaviours.

One of the most noteworthy developments in the bioimaging area is the use of BODIPY (4,4-difluoro-4-bora-3a,4a-diaza-s-indacene) fluorophores (Fig. 5a),¹¹⁸  renowned for their high brightness and exceptional chemical and photochemical stability. Recent advancements have seen BODIPY derivatives being engineered to serve as multimodal imaging probes, theranostic agents, and sensitisers for photodynamic therapy. The adaptability of BODIPY synthesis allows for the creation of metal-based BODIPY derivatives tailored for medical applications, merging the luminescent properties of BODIPY with the therapeutic or diagnostic functions of metal complexes. This synthesis strategy has led to the design of compounds that not only act as potent cytotoxic agents but also facilitate real-time tracking of their biodistribution and mechanism of action in vivo, enhancing the understanding of their biological interactions.

From a computational standpoint, accurately modelling the intricate structures of organometallic compounds and their interactions within biological systems is vital for predicting their performance and refining their design. Recent studies, such as those by Feldt & Brown,¹¹⁹  have utilised local CC methods to compute vertical excitation energies for a benchmark set of seventeen BODIPY/aza-BODIPY molecules, comparing these computational predictions with experimental results (Fig. 5b). These studies demonstrated that methods like DLPNO-STEOM-CCSD offer excellent correlation with experiments, making it one of the most accurate single-reference methods available. Similarly, Momeni & Brown¹²⁰  (Fig. 5b) conducted benchmark studies using TDDFT, revealing that the approach faces challenges related to handling differential electron correlation and the impacts of multireference character and double excitations, which contribute to discrepancies in computational accuracy.

Recently, in collaboration with the groups of Braunschweig & da Silva Jr., we have developed a novel method for the late-stage electrochemical diselenation of BODIPYs (Fig. 5c),¹²¹  revealing red-shifted absorption and tuneable colour emission with substantial Stokes shifts in selenium-containing derivatives. Photophysical analyses, supported by TDDFT and DLPNO-STEOM-CCSD computations, showed that DLPNO-STEOM-CCSD aligns closer with experimental results than TDDFT, although TDDFT better captures energy differences due to cancellation of errors and lack of size consistency of the CC approach. The selenium-containing BODIPYs notably demonstrate selective staining of lipid droplets with distinct fluorescence, highlighting their potential in bioimaging applications.

Postils, Ruipérez & Casanova,¹²²  in turn, rationalised the electronic structure of BODIPYs by combining a variety of quantum chemical methods and computational tools. They showed that BODIPYs have a mild open-shell character, which explains the usual failure of TDDFT methods. By comparing a large number of methods, the authors concluded that the S₀ → T₁ transition energies are significantly improved by using methods with double excitation effects, including post-HF, and single-reference correlation methods such as ADC(2). Finally, Berraud-Pache et al. ¹²³  used DLPNO-STEOM-CCSD to calculate the lowest vertical excitation energies of more than 50 BODIPY molecules (Fig. 5d). The method worked remarkably well, providing an accuracy of about 0.06 eV compared to the experimental data. These examples show that the DLPNO-STEOM-CCSD method is appropriate for studying the photophysical properties of medium to large organic compounds. The results also emphasise the complexity of computational approaches and highlight the need for continued refinement of these methods to improve prediction accuracy.

Nanofabrication and metallisation of organic thin films have diverse applications in fields such as energy harvesting, electronics, and sensing. Two central experimental techniques in this domain include photo-assisted chemical vapour deposition (PACVD)^124–126  and focused electron beam induced deposition (FEBID).¹²⁷  PACVD utilises photoinduced reactions to initiate the decomposition of organometallic precursors, offering a viable method for metallising organic films at or near room temperature. Conversely, FEBID enables the creation of 3D metallic structures with sub-10 nm dimensions on both planar and nonplanar surfaces. In this process, metal-containing precursor molecules are physisorbed on a substrate and degraded by a focused high-energy electron beam to yield a metallic deposit.

A critical aspect in the effectiveness of these nanofabrication techniques is the selection of suitable organometallic precursors. From a theoretical perspective, understanding the photoabsorption spectra and the dynamics of excited states following ligand dissociation is essential. Theoretical studies, particularly employing TDDFT, have been instrumental in elucidating these aspects. For instance, Walker, McElwee-White et al. ^124,125  have extensively studied the photoabsorption spectra of ruthenium carbonyl complexes featuring either η⁴ (Fig. 6a) or η³ (Fig. 6b) ligands. Their studies utilised hybrid functionals to characterise excited states, assigned, e.g. as metal-to-ligand charge transfer or ligand-field states. Similarly, Zlatar, Allan, & Fedor¹²⁷  have utilised TDDFT calculations integrated with ZORA to analyse tetrakis(trifluorophosphine)platinum(0), Pt(PF₃)₄, a model precursor for FEBID (Fig. 6c). Their studies not only validated the computed electronically excited states up to 13 eV through comparison with electron energy loss spectra, but also explored potential energy curves for these states. They demonstrated that the lowest excited states are either directly dissociative or proceed via conical intersections, facilitating dissociation pathways. Their findings underscored the predominance of the neutral dissociation channel over dissociative electron attachment for Pt(PF₃)₄, a conclusion that may also extend to other FEBID precursors.

Artificial photosynthesis is a technology that replicates natural photosynthesis by using sunlight to split water into oxygen, protons, and electrons; these electrons can then be used to reduce protons or carbon dioxide (CO₂) into energy-dense fuels like hydrogen or hydrocarbons. It aims to provide a sustainable energy solution by converting solar energy into chemical fuels, addressing the rising global energy demand and costs associated with finite conventional fuels.^128,129 

The integration of computer-aided design with synthesis and spectroscopic characterisation has significantly advanced the development of artificial photosynthetic systems. The employment of semiconductor electrodes that are modified through the covalent attachment of molecular dyes enhances visible light absorption. Computational modelling has proven indispensable in formulating sensitisers and anchoring groups robust enough to operate under aqueous and oxidative conditions, facilitating rapid electron transfer at the interfaces. Wavefunction-based methods, including EOM-CCSD, CASSCF, CASPT2, NEVPT2, and MRCI, are increasingly being utilized in this field. For instance, Head-Gordon et al. ¹³⁰  computationally investigated an iron(ii) polypyridine electrocatalyst (Fig. 7a) for the conversion of CO₂ to carbon monoxide, a complex synthesized by Chang, Long et al. ¹³¹  They specifically used CASSCF/NEVPT2 calculations to obtain static and dynamic correlations, determining the energy differences between isomers of a doubly reduced intermediate in the reaction.¹³⁰  Additionally, List, Kongsted et al. ¹³²  examined the efficacy of various single-reference methods for assessing pigment–protein complexes, specifically in predicting the relative site energies and transition moments of the Q bands in the bacteriochlorophyll a (BChl a) pigments of the Fenna–Matthew–Olson (FMO) complex, comparing these against a hybrid DFT/MRCI approach (Fig. 7b).

Naturally, the computational demands of these methods restrict their use primarily to benchmark smaller molecular models. As a more cost-effective alternative, TDDFT is extensively utilised in artificial photosynthesis and the calculation of absorption spectra of dyes and nanoclusters. To enhance the accuracy, especially for charge-transfer excitations, sophisticated functionals such as range-separated hybrid functionals have been developed, showing considerable promise in improving predictive performance in the field.

In this chapter, we have explored a comprehensive array of quantum mechanical methods for modelling excited states in organometallic chemistry. The landscape of computational strategies is vast and nuanced, requiring a thorough understanding of each method’s theoretical foundations and practical applications. From fundamental wavefunction methods like Hartree–Fock, configuration interaction, and coupled-cluster to the more computationally efficient density functional theory and its time-dependent variants, each approach has its strengths and limitations. We discussed the intricacies of perturbative methods, multireference approaches, and the importance of considering relativistic effects, especially for systems involving heavy atoms. Valence bond methods, while less common in organometallic contexts, provide an alternative perspective aligned with the chemist’s traditional view of chemical bonding. The chapter also highlighted the challenges and advances in excited state calculations, emphasising the need for appropriate method selection to balance computational cost and accuracy. Techniques such as the ΔSCF approach, time-dependent perturbation theory, and the equation-of-motion methods offer robust frameworks for studying dynamic properties and interactions of molecular systems with light. Furthermore, we discussed several applications of these theoretical approaches in the field of organometallic chemistry. Particular focus was given to the theoretical description of organometallics in various contexts, including spin-crossover materials and photodevices, bioimaging, nanofabrication, and artificial photosynthesis. Overall, we have provided a detailed roadmap for understanding and applying computational methods to model excited states in organometallic systems. As computational resources and algorithms continue to advance, the accuracy and applicability of these methods will only improve, opening new frontiers in chemical research and industrial applications.

This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/W52461X/1]. We also acknowledge the University of Kent for their financial support. The authors are grateful for partial financial support by the European Commission through the RADON project (GA 872494) within the H2020-MSCA-RISE-2019 call. This chapter is based on work from the COST Action CA20129 – Multiscale Irradiation and Chemistry Driven Processes and Related Technologies (MultIChem), supported by COST (European Cooperation in Science and Technology).