Density functional theory.

Density functional theory (DFT) finds increasing use in applications related to biological systems. Advancements in methodology and implementations have reached a point where predicted properties of reasonable to high quality can be obtained. Thus, DFT studies can complement experimental investigations, or even venture with some confidence into experimentally unexplored territory. In the present contribution, we provide an overview of the properties that can be calculated with DFT, such as geometries, energies, reaction mechanisms, and spectroscopic properties. A wide range of spectroscopic parameters is nowadays accessible with DFT, including quantities related to infrared and optical spectra, X-ray absorption and Mössbauer, as well as all of the magnetic properties connected with electron paramagnetic resonance spectroscopy except relaxation times. We highlight each of these fields of application with selected examples from the recent literature and comment on the capabilities and limitations of current methods.

Introduction

The present contribution is devoted to the use of density functional theory (DFT) in bioinorganic chemistry and more specifically in the modeling of structures, properties, and processes related to photosynthesis. DFT has been established as a valuable research tool because it can serve either to validate the conclusions that have been reached from the analysis of the experiments or to distinguish between those possibilities that were left open. The calculation of a wide range of molecular properties with DFT allows a close connection between theory and experiment and often leads to important clues about the geometric, electronic, and spectroscopic properties of the systems being studied. Here, we will first introduce briefly the general theoretical principles that constitute the basis of the DFT approach. Our priority in this paper is to describe what can be achieved through the practical use of DFT; therefore, we will then focus on the properties that can be computed and mention selected applications of DFT for molecular property calculations, drawn principally from examples relevant to photosynthetic research. We will comment not only on the strengths but also on the technical pitfalls and the current limitations of the technique, discussing the performance of DFT and the foreseeable achievements in the near future.

Theoretical background

To appreciate the special place of DFT in the modern arsenal of quantum chemical methods, it is useful first to have a look into the more traditional wavefunction-based approaches. These attempt to provide approximate solutions to the Schrödinger equation, the fundamental equation of quantum mechanics that describes any given chemical system. The most fundamental of these approaches originates from the pioneering work of Hartree and Fock in the 1920s (Szabo and Ostlund 1989). The HF method assumes that the exact N-body wavefunction of the system can be approximated by a single Slater determinant of N spin-orbitals. By invoking the variational principle, one can derive a set of N-coupled equations for the N spin orbitals. Solution of these equations yields the Hartree–Fock wavefunction and energy of the system, which are upper-bound approximations of the exact ones. The main shortcoming of the HF method is that it treats electrons as if they were moving independently of each other; in other words, it neglects electron correlation. For this reason, the efficiency and simplicity of the HF method are offset by poor performance for systems of relevance to bioinorganic chemistry. Thus, HF is now principally used merely as a starting point for more elaborate “post-HF” ab initio quantum chemical approaches, such as coupled cluster or configuration interaction methods, which provide different ways of recovering the correlation missing from HF and approximating the exact wavefunction. Unfortunately, post-HF methods usually present difficulties in their application to bioinorganic and biological systems, and their cost is currently still prohibitive for molecules containing more than about 20 atoms.

Density functional theory attempts to address both the inaccuracy of HF and the high computational demands of post-HF methods by replacing the many-body electronic wavefunction with the electronic density as the basic quantity (Koch and Holthausen 2000; Parr and Yang 1989). Whereas the wavefunction of an N electron system is dependent on 3N variables (three spatial variables for each of the N electrons), the density is a function of only three variables and is a simpler quantity to deal with both conceptually and practically, while electron correlation is included in an indirect way from the outset. Modern DFT rests on two theorems by Hohenberg and Kohn (1964). The first theorem states that the ground-state electron density uniquely determines the electronic wavefunction and hence all ground-state properties of an electronic system. The second theorem establishes that the energy of an electron distribution can be described as a functional of the electron density, and this functional is a minimum for the ground-state density. Thus, the problem of solving the many-body Schrödinger equation is bypassed, and now the objective becomes to minimize a density functional. Note, however, that although the Hohenberg–Kohn theorems assure us that the density functional is a universal quantity; they do not specify its form.

In practice, the common current realization of DFT is through the Kohn–Sham (KS) approach (Kohn and Sham 1965a). The KS method is operationally a variant of the HF approach, on the basis of the construction of a noninteracting system yielding the same density as the original problem. Noninteracting systems are relatively easy to solve because the wavefunction can be exactly represented as a Slater determinant of orbitals, in this setting often referred to as a Kohn–Sham determinant. The form of the kinetic energy functional of such a system is known exactly and the only unknown term is the exchange–correlation functional. Here lies the major problem of DFT: the exact functionals for exchange and correlation are not known except for the free electron gas. However, many approximations exist which permit the calculation of molecular properties at various levels of accuracy.

The most fundamental and simplest approximation is the local-density approximation (LDA), in which the energy depends only on the density at the point where the functional is evaluated (Kohn and Sham 1965b). LDA, which in essence assumes that the density corresponds to that of an homogeneous electron gas, proved to be an improvement over HF. While LDA remains a major workhorse in solid state physics, its success in chemistry is at best moderate due to its strong tendency for overbinding. The first real breakthrough came with the creation of functionals belonging to the so-called generalized gradient approximation (GGA) that incorporates a dependence not only on the electron density but also on its gradient, thus being able to better describe the inhomogeneous nature of molecular densities. GGA functionals such as BP86 (Becke 1988) or PBE (Perdew et al. 1996) can be implemented efficiently and yield good results, particularly for structural parameters, but are often less accurate for other properties. The next major step in the development of DFT was the introduction of hybrid functionals, which mix GGA with exact Hartree–Fock exchange (Becke 1993). Nowadays, hybrid DFT with the use of the B3LYP functional (Becke 1988; Lee et al. 1988) is the dominant choice for the treatment of transition metal containing molecules (Siegbahn 2003). This method has shown good performance for a truly wide variety of chemical systems and properties, although specific limitations and failures have also been identified.

More recent theoretical and methodological developments include the “meta-GGA” functionals, which extend the GGA corrections to higher derivatives, and the “double hybrid” functionals (Grimme 2006a, b; Neese et al. 2007a), which contain not only a fraction of exact exchange but also a fraction of orbital-dependent nonlocal correlation energy estimated at the level of second-order many-body perturbation theory. These new functionals, such as TPSSh (Staroverov et al. 2003) and B2PLYP (Grimme 2006a, b), respectively, yield improved energetics and spectroscopic properties, and will likely see more use in the future as their performance and range of applicability is established.

Properties and applications

Geometries

Optimizing the geometry of the species under investigation is the first step in most theoretical studies. Geometries predicted by DFT tend to be quite reliable and the optimized structures usually agree closely with X-ray diffraction (XRD) or extended X-ray absorption fine structure (EXAFS) data. From our experience, the achievable accuracy for short and strong metal-ligand bonds is excellent, whereas intra-ligand bonds are predicted typically within 2 pm of experiment. Weaker metal-ligand bonds are usually overestimated by up to 5 pm (Neese 2006a, b). A reasonable choice of basis set has to be made, although this condition does not pose particularly stringent requirements since the structures predicted by all DFT methods generally converge quickly with basis set size, thus making geometry optimization rather economical. Basis sets of valence triple-ζ quality plus polarization are usually enough to get almost converged results for geometries; however, results obtained with smaller basis sets should be viewed with caution. An extended study of the performance of various modern functionals and basis sets for the geometries of all first-, second-, and third-row transition metals has recently appeared (Bühl et al. 2008). Weak interactions are not satisfactorily treated with current density functionals owing to the wrong asymptotic behavior of the exchange-correlation potential, but this deficiency can be overcome to some extent by inclusion of functional-specific empirical dispersion corrections (Grimme 2006a, b).

Concerning the choice of method, the differences between density functionals are usually small for structural parameters making the choice of functional not critical for the success of a geometry optimization. GGA functionals provide good geometries and are sometimes even better than hybrid functionals, which also tend to be more expensive (Neese 2006a, 2008a). The computational efficiency of GGA in practical applications stems from the density fitting approximation (Baerends et al. 1973; Vahtras et al. 1993; Eichkorn et al. 1997) that is implemented in many quantum chemistry programs and significantly speeds up GGA calculations. This allows for fast optimizations, an important advantage especially when many different probable structures have to be considered. This has been the case, for example, in a recent computational study focusing on possible models of the oxygen evolving complex (OEC) in photosystem II (PSII), which used Mn4O topologies derived from accurate polarized EXAFS measurements (Yano et al. 2005). The great number of possible protein ligation patterns and the additional potential for a multitude of protonation and hydration states (Fig. 1) creates the need for efficient geometry optimizations which can be performed with GGA functionals such as BP86. Once optimized structures have been obtained, other molecular properties can be evaluated using a potentially more accurate hybrid functional (Zein et al. 2008a). Exploring many structural alternatives and their corresponding spectroscopic properties in this way is an important step in cross-validating theory and experiment, forming the basis for further elaboration toward more realistic models.

Fig. 1

Optimized geometry of an OEC model constructed on top of a polarized EXAFS topology for the Mn4O5Ca cluster; side-chain and water ligation shown are one out of many possibilities (Zein et al. 2008a)

Despite the overall good performance of GGA functionals, it is still likely that for certain systems high accuracy can be achieved only with hybrid functionals. In this case, the obvious choice has traditionally been the B3LYP functional. More recent studies, however, have accumulated evidence that the hybrid PBE0 and TPSSh functionals are superior performers for systems within the field of inorganic and bioinorganic chemistry (Bühl et al. 2008; Jensen 2008), the latter yielding improved energies as well. The particularly promising performance of TPSSh has been attributed in part to the use of 10% exact exchange, a value half-way between GGA and B3LYP (20%). It should be noted at this point that the computational disadvantage of hybrid functionals mentioned earlier will likely be diminished with the arrival of new state-of-the-art and potentially linear-scaling procedures such as the ‘chain of spheres’ (COSX) approximation to HF exchange (Neese et al. 2008).

Energetics and reaction mechanisms

Locating transition state structures is a more complicated task for the researcher, but in many ways it is computationally the same as optimizing a geometry; the difference is simply that the target now is not a minimum on the potential energy surface but rather a saddle point. Once this stationary point is found and its energy is computed, one gains immediate access to energy barriers and is therefore able to study reaction mechanisms. However, if this effort is to have any real value, the calculated relative energies must be reasonably accurate. A great number of studies over the years have converged to the conclusion that energetic predictions with the B3LYP functional tend to be systematically more accurate and reliable than GGA functionals. Hence, this hybrid functional is widely used for predicting and/or elucidating the major features of various mechanisms in bioinorganic chemistry (Siegbahn 2006b). The estimated error of B3LYP is approximately 2–3 kcal mol−1 for the G2 reference set of molecules, but this figure is probably too optimistic when one is faced with electronically complicated open-shell transition metal systems like those encountered in almost any bioinorganic setting. It is safe to say that there is no consensus regarding the optimal choice of method when one considers additionally the prediction of energies for electronically distinct states of the same species, such as those arising from different electronic configurations of a metal center, from a different distribution of oxidation states within a metal cluster, or even from the interplay between metal-centered and ligand-centered redox processes. When these factors come into play, the error margin can easily exceed by far the optimistic range mentioned earlier.

Nevertheless, even if the estimated errors may be already too large for quantitative predictions in cases of small activation energies such as those observed during the catalytic cycle of the OEC (Sproviero et al. 2007), the simulation of reaction pathways is a fundamentally important application of DFT. A representative example that stands out in the field of photosynthesis research is the systematic work that has been focused on elucidating mechanistic aspects in the catalytic cycle of OEC (Lundberg and Siegbahn 2004; Siegbahn 2006a, 2008a, b; Sproviero et al. 2008a, b). This line of work demonstrates that DFT calculations can offer significant input to mechanistic investigations, sometimes revealing possibilities that were not previously considered. It should be kept in mind, however, that a reaction mechanism predicted by DFT cannot be validated on the basis of computed energies alone, especially when the structure of the principal component is itself debatable. All such efforts should attempt to combine and incorporate many lines of evidence, taking into account additional criteria such as the spectroscopic properties of the putative intermediates.

Vibrational frequencies

Closely connected in research practice to the procedure of structural optimization is the calculation of vibrational frequencies. They are used not only for simulating infrared (IR) or Raman spectra but also for characterizing the nature of stationary points as minima or transition states. Moreover, the information obtained from such a calculation is used to compute statistical thermodynamic corrections to the electronic energy and thus to make direct comparisons with experimentally determined free energies. It is well established that the predicted harmonic frequencies with GGA functionals such as BP86 and PBE typically agree well with measured vibrational fundamentals if basis sets of polarized triple-ζ quality are used (Murray et al. 1992; Sosa et al. 1992; Stratmann et al. 1997). It has been shown, however, that this excellent accuracy is a result of error cancelation between the neglect of anharmonicity and systematic errors in the prediction of the correct harmonic frequencies (Neugebauer and Hess 2003). We note that the effects of anharmonicity are practically impossible to be computed with DFT for large systems of interest to biology. Intensities of IR as well as Raman modes can, however, be obtained straightforwardly.

Theoretical studies on a model of the oxygen evolving complex of PS II have demonstrated how computed vibrational frequencies can provide valuable feedback for the interpretation of experimental data. Specifically, calculations by Gascon et al. (2007) suggested that the vibrational modes of carboxylate groups ligated to manganese ions of the OEC might be insensitive to changes in the formal oxidation states of the ions because of electron delocalization within the cluster. At the same time, it was shown that the charge rearrangement associated with the S-state transitions in the OEC might induce shifts in the vibrational frequency of carboxylate groups that do not function as direct ligands to the manganese ions. These theoretical results imply that the vibrational frequency shifts observed in experimental FTIR measurements do not necessarily have to be interpreted as reflecting changes in the first coordination sphere of the Mn cluster, thus providing ways to reconcile the perceived discrepancies between FTIR and XRD data (Sproviero et al. 2008b).

Optical spectra

Density functional theory is restricted from its foundations to ground states only; therefore, the calculation of excited states and their properties has to be approached indirectly. This is achieved using time-dependent linear response theory, in which one studies the frequency dependence of a time-dependent electric field perturbation, the poles of which provide excitation energies. Thus, time-dependent DFT (TD-DFT) calculations yield the transition energy rather than the total energy of the excited state, which therefore is never explicitly calculated (Bauernschmitt and Ahlrichs 1996; Casida et al. 1998; Stratmann et al. 1998). It should be noted that the TD-DFT approach allows also for a full determination of the central quantities involved in the calculation of both absorption and circular dichroism (CD) spectra. It is also possible to predict magnetic circular dichroism (MCD) spectra through TD-DFT calculations (Seth et al. 2004, 2005; Seth and Ziegler 2006), although ab initio multireference approaches are preferred in this respect since they explicitly cover the correct physics involved (Ganyushin and Neese 2008).

Optical spectra predicted by TD-DFT with the use of either the BP86 or B3LYP functionals may occasionally be of acceptable quality (Fiedler et al. 2005; Jackson et al. 2005; Schenker et al. 2005; Stich et al. 2005) even though many problematic cases and a multitude of artifacts plague this methodology (Grapperhaus et al. 2001; Neese 2008a). TD-DFT problems arise principally from the shortcomings of current functionals and include the erroneous treatment of states with ionic or charge transfer character, poor prediction of highly excited Rydberg states, entirely missing states when double excitations are involved, and the inability to obtain the correct multiplet structures for open-shell systems. In general, one has to apply TD-DFT calculations with utmost caution and it is imperative to seek critical feedback from experimental data. With this provision, TD-DFT can be a useful interpretative tool, as was recently demonstrated by Sun et al. (2007) in their study of the P700 system found in the reaction center (Fig. 2) of photosystem I (PSI). The authors used TD-DFT in conjunction with the statistical average of different orbital potentials (SAOP) model (Gritsenko et al. 1999) to examine the excitation processes in the pair of chlorophylls that comprise P700. The detailed analysis of the individual excitations in terms of molecular orbital contributions and transition dipole moments revealed that, despite the apparent symmetric disposition of its two branches of cofactors, the P700 pair is intrinsically excited in an asymmetric fashion. On the basis of the TD-DFT results the authors were further able to establish connections with the experimentally observed asymmetric electron transfer process in PSI and propose a charge separation mechanism for P700 (Sun et al. 2007).

Fig. 2

A view of the electron-transfer chain in the reaction center of photosystem I. Chlorophyll pairs are arranged in two symmetric branches that diverge at P700 and reconverge at the iron–sulfur cluster. TD-DFT calculations have probed the nature of the excitation at the P700 pair

X-ray absorption spectroscopy

X-ray absorption spectroscopy (XAS) is a powerful probe of the electronic and geometric structure of metal sites in inorganic and biological systems since it provides valuable information on the oxidation state, geometry, and, in some cases, spin state of the metal centre (Roe et al. 1984; Westre et al. 1997). The shape, position, and intensity of absorption peaks in the X-ray absorption near-edge structure (XANES) of the metal result from core electron excitations to valence orbitals below the ionization threshold and carry information on the oxidation state, coordination, and character of the bonding with the ligands. As with optical spectra, TD-DFT can be used for the computation of metal or ligand pre-edge features, by allowing excitations into the virtual orbital space only out of localized core-holes (Ray et al. 2007; DeBeer George et al. 2008a). Although absolute transition energies are not predicted accurately, this simple and effective protocol yields relative transition energies for a series of related complexes or for a sequence of transitions to within a few tenths of an electron volt (DeBeer George et al. 2008a; Neese 2008a).

Applications to inorganic molecules have shown that good agreement between theory and experiment can be achieved with TD-DFT calculations using the BP86 functional in conjunction with a large decontracted doubly polarized triple-ζ basis set for the metal (Neese 2002) and polarized triple-ζ basis set for the remaining atoms (Ray et al. 2007; DeBeer George et al. 2008a, b). With relevance to potential catalytic intermediates involved in the water oxidation chemistry of PSII, Yano et al. (2007) have successfully correlated TD-DFT and experimental pre-edge spectra (1s to 3d excitations) of mononuclear Mn(V) nitrido and oxo compounds. More recently, Jaszewski et al. (2008) performed TD-DFT calculations of Mn core excitations in a series of Mn complexes with nitrogen and oxygen donor ligands. Excitations were allowed not only from 1s but also from 2p orbitals, yielding results that could be compared with 1s2p resonant inelastic X-ray scattering (RIXS) studies. The computed values at the BP86/TZP level were found to agree well with the experimental correlation between Mn oxidation state and the Mn K-edge and L-edge energies, confirming that TD-DFT is a robust method for analysis of XAS features. It remains to be seen how this approach extends to larger clusters such as the OEC.

Mössbauer spectroscopy

Mössbauer spectroscopy is an invaluable spectroscopic technique in bioinorganic chemistry, since it is able to probe selectively the charge and spin distribution around iron centers (Gütlich et al. 1978; see also, the contribution by Krebs and Bollinger in the present issue). The combination of DFT calculations with 57Fe-Mössbauer spectroscopy has emerged as a particularly fruitful strategy for the study of the ground-state properties of iron-containing enzymes (Schünemann and Winkler 2000; Gütlich and Ensling 1999). In the zero-applied magnetic field, the two main quantities that are extracted for a given iron site are the quadrupole splitting (ΔE Q) and the isomer shift (δ). Both quantities are related to the total electron density and are sensitive reporters of the spin state, valence state, and covalency of iron sites. The estimation of ΔE Q requires the calculation of the electric gradient field at the iron nucleus, which can be done with basis sets of sufficient flexibility in the core region (Neese 2002). Many studies at the B3LYP level have demonstrated that the sign and the magnitude of ΔE Q is predicted accurately, although absolute errors ranging from 0.3 to 1.00 mm s−1 are not uncommon (Berry et al. 2008; Godbout et al. 1999; Han et al. 2006; Salzmann et al. 1999; Sinnecker et al. 2005). Moreover, it has been shown that the computed ΔE Q values react fairly sensitively to details of the surrounding, such as counter ions.

The isomer shift is known from basic principles to be directly proportional to the electron density at the iron nucleus. Thus, it can be determined to good accuracy (often better than 0.1 mm s−1) from ground-state DFT calculations using a suitable method-specific calibration procedure on the basis of a linear correlation between the calculated electron density at the nucleus versus the measured δ (Han et al. 2006; Liu et al. 2003; Neese 2002; Sinnecker et al. 2005; Zhang et al. 2002). The predicted δ is in good agreement with experimental data when using B3LYP-calibrated curves in combination with basis sets that are flexible in the core region and extensively polarized (Berry et al. 2008; Schoneboom et al. 2005; Sinnecker et al. 2005).

Exchange couplings

In the case of bioinorganic systems which contain two or more interacting open-shell magnetic ions, the interaction is typically described in terms of the phenomenological Heisenberg–Dirac–van Vleck Hamiltonian. Thus, the main problem from the theoretical point of view becomes the evaluation of the exchange coupling constants (J) that measure the “strength” of the supposed interactions between local spins. Such systems are presently handled in the DFT framework by the broken symmetry (BS) approach, which gives access to exchange coupling constants, geometries, and total energies (Noodleman 1981). Experience indicates that hybrid functionals such as B3LYP may be slightly more accurate than GGAs for the prediction of exchange coupling constants. The finer details on the procedure are a subject of ongoing controversy, but among the different formalisms to extract the J values from separate high-spin and BS calculations, Yamaguchi’s method appears to be most suitable since it correctly reproduces the limit of both weak and strong interaction (Yamaguchi et al. 1986). It is worth emphasizing that the BS method provides excellent electron densities owing to the variational adjustment of the ionic and neutral components of the wavefunction (Neese 2004). Therefore, this approach should be able to predict geometries that faithfully reflect those of the true low-spin states. On the other hand, the spin density remains unphysical and thus for the prediction of magnetic properties based on the BS-DFT approach, it is mandatory to use spin-projection techniques (Mouesca et al. 1995; Sinnecker et al. 2004). Several computational studies of biomimetic oxomanganese complexes have been dedicated to the prediction of J values and valuable correlations between theory and experiment were found on the basis of BS-DFT calculations (Sinnecker et al. 2004, 2006).

On extension to oligonuclear systems, complications in the application of BS-DFT might arise due to the inherent indeterminacy in the values of the exchange coupling parameters. In a recent contribution (Pantazis et al. 2009), we investigate the magnetic properties of a tetramanganese complex bearing resemblance to the OEC of PSII (Fig. 3). Our results reveal that the absolute values of the exchange coupling constants J are not a safe criterion for comparing theory and experiment owing to their indeterminacy when more than a few interactions among the metals exist. Instead, one should use the J values computed with BS-DFT to extract the actual energies of the magnetic levels by diagonalizing the Hamiltonian. These energy levels can subsequently be used for constructing magnetic susceptibility curves and form a more physically meaningful basis for comparison with experimental data. Following this approach, the energetic levels computed with the TPSSh hybrid meta-GGA functional are found to agree well with experiment despite discrepancies in the fitted exchange coupling constants. Similar observations were made by Cauchy et al. (2008) in their study of a pentanuclear iron complex. The authors point out that many different sets of J values can reproduce the experimental data and proceed to exact diagonalization of the Hamiltonian and construction of a theoretical magnetic susceptibility curve to make comparisons to experiment. This approach clearly emerges as the only credible way of studying magnetic interactions with BS-DFT in oligonuclear clusters similar to the oxygen evolving complex in PSII (Pantazis et al. 2009).

Fig. 3

The tetranuclear manganese complex [Mn4O6(bipyridine)6]4+ and magnetic susceptibility curves constructed from BS-DFT results with various functionals. A direct comparison of computed and experimentally fitted exchange coupling constants is not meaningful for such systems owing to the indeterminacy of the exchange parameters

EPR spectroscopy

Electron paramagnetic resonance (EPR) spectra are parameterized in terms of an effective spin Hamiltonian (SH) which contains adjustable numerical parameters that are fitted to experiments. These SH parameters are the g-tensor, the zero-field splitting (ZFS), and the hyperfine coupling (HFC). The accuracy of EPR parameter calculations with DFT is somewhat variable. For organic radicals and biradicals (including amino acid radicals) usually good results are obtained for the g-tensor, the hyperfine and quadrupole coupling and also for the ZFS (Neese 2008b). In all DFT investigations of EPR parameters specifically developed basis sets with extra flexibility in the core region such as Barone’s EPR-II and EPR-III (Barone 1997) or the CP(PPP) basis sets (Neese 2002) should be employed. As regards the choice of functional, it is by now established that hybrid functionals are more accurate than GGA functionals (Neese 2008a). For transition metal complexes, the situation turns out to be more complicated. The g-values are usually underestimated by standard functionals, and errors of a factor of two are not uncommon. The performance of different density functionals is similar although hybrid functionals like B3LYP tend again to be slightly more accurate than GGAs like BP86 (Neese 2001a). The modeling of ZFS parameters with DFT is particularly difficult owing to the complicated spin dependence of this property (Neese 2006b). For transition metal complexes, it was shown that DFT predicts the ZFS parameter with the correct sign but tends to underestimate its magnitude, often by a factor of 2 (Neese 2003). Meanwhile, a certain number of applications have demonstrated the usefulness of ab initio treatments for the calculation of the ZFS (Ganyushin and Neese 2006; Neese et al. 2007b).

The intricacies in the application of DFT in this area are highlighted in a detailed evaluation of DFT performance for the prediction of ZFS in Mn(II) coordination complexes (Zein et al. 2008b). The study revealed that regardless of whether the spin–orbit coupling (SOC) part of the ZFS was estimated with the Pederson–Khanna or the quasi-restricted orbitals approach, accounting for the spin–spin (SS) interaction always improves the results. The physical necessity of accounting for the SS interaction is shown from its 30% contribution to the axial D parameters. In general, the calculations were found to overestimate systematically the experimental D values by 60%. The authors call attention to the fact that the signs of calculated axial ZFS parameters are unreliable once E/D > 0.2. Calculated D and E/D values were found to be highly sensitive to small structural changes; disconcertingly, the use of optimized geometries was found to lead to a significant deterioration of theoretical predictions relative to experimental XRD geometries. A subsequent study (Zein and Neese 2008) showed that using the coupled-perturbed spin–orbit coupling (CP-SOC) approach (Neese 2007) together with hybrid DFT functionals leads to a slope of the correlation line between experimental and calculated D values that is essentially unity, provided that the direct SS interaction is properly included in the treatment.

For the case of the hyperfine coupling to the metal, DFT performance is not entirely satisfactory (Munzarova and Kaupp 1999; Munzarova et al. 2000). Since this property involves three contributions (Fermi contact, spin–dipolar, and spin–orbit coupling) which feature different physical mechanisms, it is difficult to calculate all of them simultaneously with quantitative accuracy. Ligand HFCs are easier to compute but, again, results are less accurate than for organic radicals, and errors of 30% must be tolerated (Neese 2001b). Kossmann et al. (2007) investigated the performance of modern DFT functionals for the prediction of molecular hyperfine couplings in extended test calculations for a series of small radicals and transition metal complexes. It was shown that for the prediction of metal and ligand HFCs, TPSS is better than BP86, but more importantly, that the hybrid variant TPSSh is significantly superior to TPSS and probably even better than the “de facto standard” B3LYP functional. The double-hybrid B2PLYP functional also affords accurate predictions, particularly for HFCs of metal nuclei, but the existence of outliers suggests that this method may lack stability. The reliable performance of the TPSSh functional has since received additional confirmation in our recent study (Pantazis et al. 2009) aimed at the analysis of hyperfine coupling parameters in tetramanganese models of the OEC.

A topical application of DFT for the determination of EPR parameters, highlighting the capabilities of the broken-symmetry approach, was the study of the trapped-valence, exchange-coupled [MnIIIMnIV(μ-O)2(μ-OAc)DTNE]2+ complex (Sinnecker et al. 2006). Interestingly, BP86-optimized geometries were better than those obtained from B3LYP; however, B3LYP yielded exchange coupling constants in excellent agreement with experiment. The coupled perturbed Kohn–Sham equations were employed for the g-tensor calculations, and a strategy for the computation of g-tensor site values was presented that provided single-site g-tensors in good agreement with the expectations for the respective Mn formal oxidation states. Spin projection gave the g-tensor of the coupled manganese complex in good agreement with the experimental results. Small values were found for the nuclear quadrupole splitting of 55Mn. Hyperfine tensors were furthermore calculated and spin-projected. 14N and 1H ligand hyperfine data were found to compare well with experiment. 55Mn HFCs were qualitatively in line with experimental results, tracing the source of anisotropy to the MnIII center. However, isotropic 55Mn HFCs were distinctly underestimated. The authors indicated that this deficiency is systematic in character and does not originate from the broken symmetry approach. Similar deviations were found between theory and experiment for DFT calculations on mononuclear Mn complexes, suggesting that the use of a universal scaling factor of approximately 1.5 might be appropriate.

Summary and perspectives

Density functional theory methods have already been established as a valuable research tool both in independent applications and as a complement of experimental investigations. In favorable cases, the calculated properties are sufficiently accurate to discriminate between structural alternatives for reaction intermediates or other species that are not amenable to experimental structure elucidation. DFT appears generally reliable for geometries, vibrational frequencies, and total energies, having over wavefunction-based methods the advantage of quick convergence to the basis set limit. DFT appears to be quite successful for the prediction of molecular properties as well, since a number of spectroscopic properties of interest to the bioinorganic community can be predicted with good accuracy. Hybrid functionals are in most cases better performers, with the TPSSh functional emerging as a potential new standard. There are still cases, however, where quantitative accuracy may be difficult to achieve, especially for the prediction of EPR parameters or optical spectra, necessitating a cautious and critical approach from the part of the researcher.

It is important for both practitioners of DFT and the nontechnical audience of DFT studies to keep in mind that errors do arise and they can be significant. Despite the enormous advances in density functional implementations and the sufficiently documented accuracy of results for many applications, there is no systematic way of improving DFT or converging its results to the “correct” answer, in contrast to some of the traditional wavefunction-based methods. Moreover, the success of a particular functional in one setting does not guarantee its performance in a different one. Therefore, to enhance their credibility, DFT applications must include some form of validation or estimation of the error range on the basis of careful comparison between calculated and measured observables.

A final point of interest is that DFT studies of bioinorganic systems have usually employed simplified models in vacuo. Therefore, the issue of modeling the interaction of the active site with the protein environment and the solvent comes into play (Noodleman and Han 2006; Noodleman et al. 2004, Schoneboom et al. 2005). A realistic and computationally feasible modeling of these effects can be achieved at present by combining the DFT treatment of the active site with a classical force-field description of the surrounding protein. This is the concept behind quantum mechanics/molecular mechanics (QM/MM) approaches (Senn and Thiel 2007), which are discussed by Batista and coworkers in the present issue. In a broader theoretical context, many issues can be identified that warrant further developments. We anticipate that in the future we will witness developments regarding functionals that provide a consistent treatment of exact exchange, improvements in the treatment of electronic relaxation and excited states, and a more proper treatment of magnetic and relativistic effects. A longer term target is certainly the reliable, consistent and efficient treatment of system dynamics or of very large systems.