ONETEP + TOSCAM: uniting dynamical mean ﬁeld theory and linear-scaling density functional theory

We introduce the uniﬁcation of dynamical mean ﬁeld theory (DMFT) and linear-scaling density functional theory (DFT), as recently implemented in ONETEP, a linear-scaling DFT package, and TOSCAM, a DMFT toolbox. This code can account for strongly correlated electronic behavior while simultaneously including the eﬀects of the environment, making it ideally suited for studying complex and heterogeneous systems that contain transition metals and lanthanides, such as metalloproteins. We systematically introduce the necessary formalism, which must account for the nonorthogonal basis set used by ONETEP. In order to demonstrate the capabilities of this code, we apply it to carbon monoxide-ligated iron porphyrin and explore the distinctly quantum-mechanical character of the iron 3 d electrons during the process of photodissociation.


Introduction
In the last few decades, density functional theory (DFT) has established itself as a key method in computational materials science. [1][2][3][4] Facilitated by exponentially increasing computing power, modern DFT codes are capable of routinely calculating the electronic structure of hundreds of atoms, opening the door to quantum-mechanical modeling of a vast landscape of systems of considerable scientific interest.
The range of computationally accessible systems has broadened even further with the advent of linear-scaling DFT codes (that is, codes whose computational cost scales linearly with the number of atoms in the system, rather than the cubic scaling of traditional methods).
ONETEP 5,6 is one such code, notable for its equivalence to plane-wave approaches due to the in situ optimization of its basis (a set of local Wannier-like orbitals). Its ability to routinely perform DFT calculations on systems containing thousands of atoms allows more detailed study of nanostructures, 7,8 defects, 9,10 and biological systems. [11][12][13][14] That said, DFT is not without its shortcomings. Many of these stem from its approximate treatment of exchange and correlation via an exchange-correlation (XC) functional. These shortcomings become especially evident in "strongly-correlated" systems, which typically contain transition element or rare-earth atoms whose 3d-or 4f -electron shells are partially filled. Electrons in these shells are in especially close proximity with one another, and their interaction is too pronounced to be adequately described by DFT, which can provide even qualitatively incorrect descriptions of the electronic structure. For example, DFT often yields magnetic moments inconsistent with experiment, 15 predicting some insulators to be metallic, 16,17 and yielding equilibrium volumes dramatically different to experiment. 18 DFT also fails to capture important dynamic properties that are enhanced by strong correlation, such as satellite peaks in photoemission spectra. 19,20 These cases motivate the need for more accurate theories. One such approach is dynamical mean field theory (DMFT), 21-25 a Green's function method that maps the lattice electron problem onto a single-impurity Anderson model with a self-consistency condition.
Local quantum fluctuations are fully taken into account, allowing DMFT to capture complex electronic behavior such as the intermediate three-peak states of the Mott transition, the transfer of spectral weight, and the finite lifetime of excitations. 22 Furthermore, it is possible to embed DMFT within a DFT framework, whereby only atoms with strongly-correlated electrons are treated at the DMFT level, while the rest of the system can be treated at the DFT level. 17 This is critical, as DMFT alone is prohibitively expensive for studying most realistic computational models.
In the past decade, numerous codes have been written to add DMFT functionality to existing DFT packages. These include EDMFTF 26,27 and DFTTools 28 on top of Wien2K, 29 EDMFTF 27 on top of VASP, [30][31][32] DCore 33 on top of Quantum Espresso 34 and OpenMX, 35,36 TOSCAM 37 on top of CASTEP, 38,39 Amulet 40 on top of Quantum Espresso 34 and Elk, 41 and ComDMFT 42 on top of FlapwMBPT. 43,44 Many of these make use of stand-alone libraries such as TRIQS, 45 ALPS, 46 iQIST, 47 or W2dynamics. 48 This paper introduces the implementation of TOSCAM (A TOolbox for Strongly Correlated Approaches to Molecules) on top of ONETEP, a linear-scaling DFT code. In contrast to the packages mentioned above, this approach uniquely enables us to perform DMFT calculations on large and aperiodic systems such as nanoparticles and metalloproteins. This code has already seen success: it has been used to explain the insulating M 1 phase of vanadium dioxide, 49 to demonstrate the importance of Hund's coupling in the binding energetics of myoglobin, 50,51 and to reveal the super-exchange mechanism in the di-Cu oxobridge of hemocyanin and tyrosinase. 52 But until now it has not been available to the scientific community at large. The DMFT module in ONETEP has been included in version 5.0, and TOSCAM is being released at <github link to accompany publication>. This paper presents an overview of this methodology, its implementation, and an example of its application.

The ONETEP framework
In the ONETEP implementation of linear-scaling DFT, we work with the single-particle density-matrix: where {φ α } are a set of localized non-orthogonal generalized Wannier functions (NGWFs) and K αβ is the density kernel. These orbitals are variationally optimized in situ during the energy-minimization carried out as part of the DFT calculation. 53 For the purposes of this optimization, the NGWFs are in turn expanded in terms of a systematic basis of psinc functions 54 -systematic, in the sense that the size of this basis is determined solely by a scalar parameter (a plane-wave kinetic energy cutoff determining the grid spacing) that can be increased until convergence is reached. In this scheme, a DFT calculation does not involve cyclically calculating the Kohn-Sham density and potential, but instead involves the direct minimization of the DFT energy with respect to both the density kernel and the NGWF expansion coefficients ( Figure 1). Due to the fact that the NGWFs are localized, the associated matrix algebra is sparse. Meanwhile, because the NGWFs are optimized in situ, the calculations are not prone to basis set incompleteness or superposition error, 55   G αβ (ω) (ω here may be ω + iη or iω n if operating in the finite-temperature Matsubara representation) and the self-energy Σ αβ , which are related via where µ is the chemical potential (fixed at the mid-point of the energies of the highest occupied and lowest unoccupied KS orbitals), and S αβ is the NGWF overlap matrix (that Treating most physical systems at the DMFT level would usually be prohibitively expensive. The DFT + DMFT scheme takes advantage of the fact that strong electronic correlation is often confined to identifiable localized subspaces (for instance, the 3d orbitals of a transition metal atom), with the remainder of the system having a delocalized, free-electron character. In such systems, the correlated subspaces can be treated at the DMFT level, while DFT alone should be sufficient everywhere else.

The Anderson impurity model
In order to efficiently find a self-consistent solution to equation 2, DMFT relies on mapping correlated subspaces to auxiliary Anderson impurity models (AIMs). The AIM is a simplified Hamiltonian that describes the interaction of a number of sites (known as impurity sites) with a bath of additional electronic levels: whereĤ bath describes the non-correlated behavior of the bath (parameterized by the hopping matrix ε ij ),Ĥ loc the impurity (parameterized by the impurity hopping t mm and the interaction HamiltonianĤ U ), andĤ mix the coupling between the two (parameterized by V mi ). This Hamiltonian is depicted pictorially in Fig. 2. The bath and impurity sites have a shared chemical potential µ, andĉ/f are the annihilation operators for the bath/impurity.
The convention throughout will be that Greek indices correspond to NGWFs, m and m to Hubbard subspaces and their corresponding impurity sites, and Latin indices to bath sites.
σ is the spin index.
The non-interacting Anderson model (i.e. H U = 0) has the Green's function where the full hopping matrix is of the block matrix form It follows that the (non-interacting) impurity Green's function -that is, the top-left-hand block of G 0 tot (ω) -simplifies to where is the so-called impurity hybridization function. This quantity is of particular importance because it encapsulates all of the contributions of the bath sites to the physics of the impurity sites; the AIM impurity Green's function is given by

A DMFT calculation
This subsection will walk through the steps in a standard DMFT calculation as performed in TOSCAM + ONETEP. It is important to note that DMFT typically invokes a mean field approximation across multiple correlated sites (hence dynamical "mean field" theory), an approach that only becomes exact in the limit of infinite coordination (or equivalently, dimensions). This is not the case in our following real-space approach, where correlated sites are treated as a (possibly multi-site) AIM. Hamiltonian, an estimate of the system self-energy (zero is a reasonable starting point), and a total Green's function (obtained via equation 2) are each projected onto the correlated subspaces. For instance, the local Green's function is given bỹ where W I mα = ϕ I m |φ α is the overlap of the NGWFs and the Hubbard projectors. In a similar manner one can obtain the projected self energyΣ I (ω) and the projected Kohn- We can determine the appropriate impurity hopping parameters t mm for the AIM by comparing the AIM impurity Green's function G imp (ω) and the local Green's function G I mm (ω): if these are to match in the high=frequency limit to order O(1/ω 2 ) then it follows that Meanwhile, in order to define {V mi } and {ε ij }, we define the local hybridization function for our physical system which is analogous to the definition of the impurity hybridization function (equation 7).
We choose the impurity model bath parameters such that the AIM hybridization function matches this local hybridization function as closely as possible. This is done by minimizing the distance function using a conjugate gradient (CG), Broyden-Fletcher-Goldfarb-Shanno (BFGS), or similar minimization algorithm. Here, ω c is a cutoff frequency and γ is a user-specified parameter that can allow for preferential weighting of agreement at low frequencies.
In order to complete the construction of the auxiliary AIM Hamiltonian we choose H U to be of the Slater-Kanamori form 59,60 This Hamiltonian is well-suited to capturing multiplet properties of low energy states. Now that we have defined ε, V , t, and H U , the mapping of a real system to an auxiliary AIM is complete. In theory, this mapping can be exact: as long as ∆ imp (ω) and∆ I (ω) match exactly, G imp (ω) andG I (ω) will also. Getting this mapping right is therefore of the utmost importance.

Solving the AIM
Having constructed the AIM Hamiltonian H AIM , the next step is to calculate the Green's function of the AIM (known as the impurity Green's function): where • is the thermodynamic average, which at zero temperature becomes ψ 0 | • |ψ 0 .
Resolving equation 12 is highly expensive, and becomes one of the most substantial computational barriers in a DMFT calculation. If there are m bath sites and n impurity orbitals, the Hilbert space of this problem scales as 4 m+n . (For a system containing a single transition metal there will be five impurity orbitals -one for each 3d orbital -and then typically six to eight bath sites.) This is far larger than any of the other matrix inversions that we need to calculate during the DMFT loop (for instance, G αβ is only as large as the number of Kohn-Sham orbitals, which in turn will be of the order of the number of electrons in the physical system -typically several thousand at most). There are a number of approaches for obtaining G imp , such as exact diagonalization (ED) and continuous time Monte Carlo algorithms. The calculations in this work employ ED via the Lanczos algorithm to evaluate equation 7, a process which is explained in detail in the Supplementary Information. Given a solution G imp (obtained via ED or otherwise), the impurity self-energy can then be obtained via where the non-interacting impurity Green's function is given by equation 6. Note that this operation is far less expensive than equation 7 because these matrices are only m × m in size.

Upfolding and double-counting
Having obtained the impurity Green's function Σ I imp for each AIM, the final step is to upfold this result to the complete physical system. Since the original DFT solution already contains the influence of the Coulomb interaction to some degree, double-counting becomes an issue.
A popular form of the correction is where n is the total occupancy of the subspace, n σ is the occupancy of the subspace for the spin channel σ, and with N being the number of orbitals spanning the correlated subspace (and recall that U = U − 2J). 61 This double-counting is derived by attempting to subtract the DFT contributions in an average way; U av is the average of the intra-and inter-orbital Coulomb parameters.
The self-energy is upfolded to the NGWF basis via -and with that, we are back where we started, having generated a new estimate of the self-energy Σ αβ for the full system.
To summarize, the scheme is as follows: 1. perform a DFT calculation to construct the system Hamiltonian 2. initialize the self-energy as Σ αβ (ω) = 0 3. obtain the Green's function for Note that if we only have one correlated site in our system, this mapping is exact, and the local lattice Green's function at step 9 will already match the impurity Green's function.
This is not the case for bulk systems. There, the mean field approximation that we adopt means that the self-energy of a correlated site is also inherited by the "bath" i.e. one would typically solve a single Anderson impurity problem but then in equation 16, the index I would run over all correlated sites. This means that after step 9 we must return to step 3, and repeat this loop until the local lattice and impurity Green's functions match.
Once the calculation is converged, we can extract system properties from the Green's function. For example, DFT+DMFT total energies can be calculated via where the first term is the DFT energy functional, ε iσ are the Kohn-Sham eigenenergies, Other quantities that can be extracted from the Green's function include the density of states and the optical absorption. One can also apply standard ONETEP analysis techniques to the electron density (such as natural bonding orbital analysis). These techniques will be demonstrated in Section 3.

Extensions
There are several possible extensions to the theory described thus far. These are not essential but often useful.

Enlarged AIM via cluster perturbation theory
If an AIM has too few bath sites at its disposal, it will be insufficiently flexible to fit a given

Self-consistency
For a system with a single correlated site, there is no feedback from the self energy to the hybridization function, and -provided the AIM is sufficiently representative -the DMFT algorithm will converge in a single step. In this case the algorithm is not a mean-field approximation, but exact. This scheme is shown in Fig. 3a.
However, there are a number of reasons why we may not be content with the resulting solution. Firstly, the total number of electrons in the system is related to the total retarded Green's function via where ρ αβ (ω) is the basis-resolved DMFT density matrix. There is no reason a priori why the Green's function, updated via the DMFT loop, should yield the same number of electrons as we started with -in fact, this is almost never the case. For this reason, charge conservation can optionally be enforced by adjusting µ so that µ −∞ ρ(ω) = N . This update is performed during each DMFT cycle, which means that our total Green's function (now adjusted by our altered µ) will not necessarily be consistent with the self energy -and consequently more than one DMFT loop will likely be required to iterate to self-consistency (Fig. 3b). We will refer to this as "charge-conserving" DMFT.
Finally, in the DFT formalism, the Hamiltonian is a functional of the density. It could be argued that if we are to be fully self-consistent, if ever the density changes the Hamiltonian should be updated accordingly. In this scheme, one iterates until Σ, H, and µ all converge ( Fig. 3c). This we will refer to as "fully self-consistent" DMFT. We use Pulay mixing 67,68 to update the Hamiltonian (via the density kernel) and the self-energy. Performing this doubleloop naturally makes the calculations much more expensive, but they remain feasible. This approach was taken in Refs. 65-70, for example.

Practical implementation
In our implementation, ONETEP and TOSCAM are responsible for separate sections of the DMFT loop, as shown in Figure 4. As the calculation proceeds, these two programs are called in alternation, with the entire procedure being driven by an overarching script.

Scaling
One of our primary considerations is how ONETEP+TOSCAM calculations scale. As discussed already, obtaining the Green's function of the AIM scales very poorly with the number of AIM sites. This is shown in Fig. 5a. We are not entirely in a position to dictate the number of AIM sites: a 3d correlated site is represented as a five-site impurity, and typically we need to include at least six bath sites to give the AIM sufficient flexibility to fit the hybridisation function.
To some extent, poor scaling can be overcome by efficient parallelization. Both ONETEP and TOSCAM employ hybrid MPI and OpenMP parallelization schemes. ONETEP's parallelization is highly optimized. Individual atoms are distributed across MPI threads, with lower-level computationally-intensive operations (including 3D FFT box operations, sparse

Iron porphyrin
To demonstrate the use of the ONETEP+TOSCAM interface, the second half of this paper presents some calculations on an archetypal strongly-correlated system: an iron porphyrin ring with imidazole and carbon monoxide as the axial ligands (FePImCO) shown in Fig. 6a, a toy model for the full carboxymyoglobin complex (Fig. 6b). By translating the carbon and there are unresolved questions surrounding the process of carbon monoxide photodissociation, as we shall discuss.

Computational details
All DFT calculations were performed using a modified copy of ONETEP. 5,58,73-76 Those modifications were subsequently included in ONETEP 5.0. All calculations used the PBE XC functional, 77 were spin-unpolarised, and had an energy cut-off of 908 eV. There were thirteen NGWFs on the iron atom, four on each carbon, nitrogen, and oxygen, and one on each hydrogen. All NGWFs had 6.6Å cut-off radii. Open boundary conditions were achieved using a padded cell and a spherical Coulomb cut-off. 78 Scalar relativistic pseudopotentials were used, generated in-house using OPIUM, 79-86 and the Hubbard projectors were constructed from the Kohn-Sham solutions for a lone iron pseudopotential. 58 The bound structure was taken from Ref. 87 Example input and output files can be found on Materials Cloud. 88

The quantum-mechanical state of the 3d iron subspace
A large effort (largely in the quantum chemistry community) has been made to correctly predict the spin state of Fe(II)P with (and without) a variety of axial ligands. These range from decades-old Hartree-Fock calculations to recent FCIQMC studies. [89][90][91][92][93][94] FePImCO is one of the simpler cases, with a singlet state universally predicted. Meanwhile, FePIm has proven to be more of a challenge. Experiment characterises FePIm as a quintet. Semi-local DFT wrongly predicts it to be a triplet (as shown in Fig. 7). DFT + U remedies this, 95  Fe-C distance (Å) This is plotted in Fig. 8a as a function of the Fe-C distance. The unbinding is plainly visible as a sudden step in the total occupancy, at the same distance that DFT predicted the lowto-high-spin crossover (refer back to Fig. 7). The effect of DMFT is especially pronounced at large Fe-C distances, where it drives the subspace occupancy towards the expected formal As a means of analysing the spin state of the iron atom during the dissociation process with DMFT, we construct the reduced density matrix where we take the partial trace of the low-lying eigenstates of the AIM over the bath degrees of freedom, leaving a mixed density operator for the impurity alone. It is then straightforward Unfortunately self-consistent DMFT calculations proved very difficult to converge beyond the low-to-high spin transition, so these results have been excluded throughout. Below this transition, the two methods qualitatively agree. (b) The effective spin S eff of the reduced density matrix, defined via Tr[Ŝ 2ρ ] =h 2 S eff (S eff + 1) for the DMFT calculations and as S eff = 1 2 (n ↑ 3d − n ↓ 3d ) for the broken-symmetry DFT+U +J calculations. (c) The decomposition of the reduced density matrix by spin state. The colours correspond to the respective weights of the different contributions; if a colour occupied all the vertical axis, it would mean that all eigenvectors of the density matrix are in that particular quantum sector.
to calculate the expectation value ofŜ 2 = i,jŜ i ·Ŝ j and extract the effective spin S eff (Fig. 8b). Here we can see that at large distances we approach the quintet S eff = 2. At small distances we are closer to the triplet value S eff = 1. Note that this does not mean that DMFT has failed to predict that FePImCO is a singlet. Rather, this result is compatible with (but does not confirm the existence of) a singlet forming across the Fe-CO bond. By limiting ourselves to the Fe subspace we cannot directly measure such a singlet.
For comparison, we also present the occupancy and spin of the 3d subspace as calculated by DFT+U +J, using the same U and J parameters as for the DFT+DMFT calculations. DFT+U (+J) is a widespread and computationally cheap correction to DFT for accounting for electronic correlation, and is equivalent to solving the DMFT impurity problem at the level of Hartree-Fock. 75,[96][97][98] We note that DFT+U +J recovers the correct quintet spin state in the dissociated limit, in line with previous DFT+U studies, 95 with a window where the triplet state is favored. Furthermore, the singlet state becomes unstable at a shorter Fe-C distance of around 2Å (see Fig. 8b). This is a common feature of DFT+U , 99,100 where the corrective +U potential reduces the hybridization between the correlated subspace and the ligand orbitals, thereby weakening the bond between them. Addressing this within a Hubbard model framework requires more sophisticated approaches such as inter-site terms 99 or applying corrective Hubbard potentials to ligand subspaces. 100 This is not a problem that DMFT suffers from.
To inspect the DMFT reduced density matrix in more detail, one can construct the spin-projectorP as the sum of the eigenstates |s of the operatorŜ 2 with eigenvalue S(S + 1). This allows us to evaluate the fraction of the reduced density matrix in singlet, doublet, triplet, and higher states via Tr[P SρPS ] for S = 0, 1 2 , 1 etc. Note, however, that this approach is incompatible with the CPT extension. The CPT extension involves solving an auxiliary AIM Hamiltonian that shares the same impurity Green's function as a larger AIM Hamiltonian, and consequently any quantities derived directly from the Green's function will be unaffected.
However, there is no such guarantee for the reduced density matrix, because the hybridization function of this auxiliary system does not necessarily match that of the physical system. To overcome this, the CPT extension was at first applied in order to obtain an approximate solution, but then removed for the final DMFT step. Typically this final step required the addition of an extra bath site so that the AIM acquired sufficient flexibility to fit the impurity hybridization function to the local hybridization function without the assistance of the CPT extension.
The decomposition of the reduced density matrix into spin sectors is displayed in Fig. 8c.
It reveals a large quintet state contribution in the limit of dissociation, but also that, regardless of Fe-C distance, many different spin sectors are important. This would be missed if we only examined S eff or only performed DFT.
Evidently, a multitude of states play an important role throughout CO-unbinding, and therefore the success of DFT + U and HF in predicting the quintet ground state must be for the wrong reasons, as neither go beyond the single-determinantal picture. (Note that HF is known to overly favour high-spin states. 101 ) It should be noted that the precise details of Fig. 8 are somewhat sensitive to various simulation parameters -most notably the definition of the Hubbard projectors -but qualitatively the results are expected to hold generally.

Photodissociation
The photodissociation mechanism of carboxymyoglobin is already relatively well understood.
Irradiation at 570 nm (2.18 eV) causes the excitation of electrons in the porphyrin ring into low lying singlet states with π/π * character (the so-called Q band). 102 The carbon monoxide ligand then dissociates within 50 fs, as the system adiabatically crosses to a repulsive antiback-bonding orbital. 103,104 There is a small (but not insignificant) predicted energy barrier of 0.08 eV between these two states, as calculated by B3LYP and TDDFT. 87 The porphyrin then undergoes the "intersystem crossing", a complicated, multi-step process which ultimately takes the dissociated system to its high-spin ground state.
Semi-local DFT captures this process qualitatively. The energies of the lowest unoccupied KS molecular orbitals as predicted via DFT are shown in Fig. 9. The Q band is present, and the pathway from the Q band to the anti-back-bonding orbital is clearly visible via their To compare the results of DMFT to these KS eigenenergies, the analogous quantity we must extract is the DOS. The DOS is given by the trace of the many-body density matrix The DMFT DOS is compared to the KS eigenenergies in Fig. 10. Qualitatively, they yield very similar results, but with DMFT (unlike DFT) we obtain the lifetime of the excitations.
To reveal the contribution of individual atoms (or groups of atoms) towards the DMFT where I denotes a subset of NGWFs typically belonging to atoms that are a particular element or part of a spatially distinct subsystem (e.g. all the NGWFs belonging to atoms in the porphyrin ring). One such LDOS is plotted in Fig. 11, along with isosurfaces of the spectral density at energies corresponding to the various peaks in the DOS. The Q-band π/π * orbitals and the Fe-CO back-and anti-back-bonding orbitals are all clearly identifiable.
Another important quantity that can be extracted from DMFT calculations is the optical spectrum. The theoretical optical absorption spectrum can be obtained within the linearresponse regime (that is, Kubo formalism) as π/π * Fe CO imidazole porphyrin σ d xy π/π * σ * Figure 11: Self-consistent DMFT density of states for carboxy-heme with a Fe-C distance of 2.06Å. The DOS is further decomposed into contributions from the iron atom, CO molecule, imidazole ligand and porphyrin ligand. Above, isosurfaces of ρ(r, ω peak ) have been plotted for each peak.
where Ω the simulation cell volume, f (ω) is the Fermi-Dirac distribution, ρ is the basisresolved spectral density, the i and j indices correspond to Cartesian directions, the velocity which includes the effect of non-local pseudopotentials V nl on the velocity operator matrix elements, and adopts the no-vertex-corrections approximation. 105 Optical spectra for heme are typically carried out in liquid or gas phases, and so are described by the isotropic part of the optical conductivity tensor The optical absorption spectra for carboxy-heme complexes as given by self-consistent DMFT are plotted in Fig. 12. These spectra are dominated by a feature at around 2 eV associated with π-π * transitions on the porphyrin ring -that is, the Q band.  Figure 12: Optical spectra of FePImCO calculated using self-consistent DMFT, going from ligated (dark) up to the point of dissociation (light). Also pictured are the Q-band peaks from experimental spectra of carboxymyoglobin. 106 photoexcitation of the anti-back-bonding orbital.

Conclusions
This paper has described how DMFT has been interfaced with linear-scaling DFT in the ONETEP+TOSCAM implementation. Crucially, for the purposes of simulating metalloproteins, this DFT + DMFT implementation does not compromise our ability to model thousands of atoms at the DFT level, opening up a new frontier for accurate simulation of complex and heterogeneous systems containing transition metals and lanthanides.
The ONETEP+TOSCAM interface will continue to be developed. In particular, work is underway to compute forces at the DFT+DMFT level (as has been achieved elsewhere 37,107 ), a GPU implementation of the ED solver will be incorporated, as well as a CTQMC solver (which will allow us to solve substantially larger AIMs.) Note that it is straightforward to add additional solvers due to the modularity of the code.
Calculations on the photodissociation of carboxymyoglobin showcased the kinds of results one can extract from such a DFT + DMFT calculation on a metalloprotein, including somesuch as the mixed quantum state of the iron 3d subspace and the finite lifetime of excitations -that are inaccessible via DFT. And while the calculations do not reveal any previously unknown physics, there is scope here to resolve some unanswered questions surrounding the photodissociation process. In particular, the remarkably fast rate of photodissociation (∼ 50 fs) is at odds with the gentle slope of the potential energy surface (discussed above) and the predicted barrier on the order of 0.1 eV (compared to the 0.028 eV zero-point energy of the Fe-C stretching mode). 87 Further study could investigate this apparent contradiction.

Supporting Information Available
The Supporting Information contains a brief description of the TOSCAM/ONETEP interface, as well as relevant links, and an introduction to the Lanczos algorithm and how it may be used to solve the Anderson impurity model. This information is available free of charge via the Internet at http://pubs.acs.org.