TURBOMOLE: Today and Tomorrow

TURBOMOLE is a highly optimized software suite for large-scale quantum-chemical and materials science simulations of molecules, clusters, extended systems, and periodic solids. TURBOMOLE uses Gaussian basis sets and has been designed with robust and fast quantum-chemical applications in mind, ranging from homogeneous and heterogeneous catalysis to inorganic and organic chemistry and various types of spectroscopy, light–matter interactions, and biochemistry. This Perspective briefly surveys TURBOMOLE’s functionality and highlights recent developments that have taken place between 2020 and 2023, comprising new electronic structure methods for molecules and solids, previously unavailable molecular properties, embedding, and molecular dynamics approaches. Select features under development are reviewed to illustrate the continuous growth of the program suite, including nuclear electronic orbital methods, Hartree–Fock-based adiabatic connection models, simplified time-dependent density functional theory, relativistic effects and magnetic properties, and multiscale modeling of optical properties.


INTRODUCTION
TURBOMOLE is a collaborative, multinational software development project aiming to provide highly efficient and stable computational tools for quantum-chemical simulations of molecules, clusters, periodic systems, and solutions.The software suite is optimized for widely available, inexpensive, and resource-efficient hardware, such as multicore workstations and medium-size compute clusters.TURBOMOLE specializes in electronic structure methods with an outstanding accuracy− cost ratio, such as density functional theory including the random phase approximation (RPA), GW-Bethe−Salpeter equation (BSE) methods, second-order Møller−Plesset (MP2) theory, and coupled-cluster (CC) methods.The code is based on Gaussian basis sets and has been pivotal for the development of many fast and low-scaling algorithms in the past three decades, such as integral-direct methods, the resolution-of-the-identity (RI) approximation, fast multipole methods, imaginary frequency integration, Laplace transform, and pair natural orbital methods.
The development of TURBOMOLE was started in the late 1980s by Reinhart Ahlrichs and his group.Integral-direct algorithms and non-Abelian point group symmetry were among the first distinctive capabilities of TURBOMOLE, which initially focused on Hartree−Fock (HF) methods, with subsequent extensions for second-order MP2 perturbation theory 1−7 and time-dependent HF (TDHF) response properties. 8A major milestone was the relatively early adoption of density functional theory (DFT) using newly designed quadrature algorithms 9 and that of time-dependent DFT (TDDFT) shortly afterward. 10With extensions to meta-GGA, 11 RPA, and other fifth-rung methods, 12−16 as well as current density 17 and local hybrid functionals, 18,19 critical performance improvements, 20−24 and a plethora of available analytical derivatives of ground-and excited-state energies, 25−39 TURBOMOLE has become one of the leading allpurpose molecular (TD)DFT codes.The development and implementation of the RI approximation for Coulomb (RI-J 40,41 ) and exchange contributions, (RI-K 42,43 ), as well as its generalization to post-HF theories such as MP2 and CC2 44−46 and its multipole-accelerated version for extended systems (MARI-J), were other critical innovations. 47−68 More recent additions include explicitly correlated wave function methods up to CCSD(T) and BCCD(T), 68,69 efficient pair natural orbital (PNO) approaches, 70−73 solvation models and embedding, 74−83 two-component relativistic methods, 84−86 GW-BSE type methods 52,87,88 real-time (RT) TDDFT, 89 and nonadiabatic molecular dynamics. 90o ensure continuity and coordinate the development, maintenance, and distribution independent of individual developers or groups, TURBOMOLE GmbH, a limited liability company located in Karlsruhe, Germany, was founded in 2007.TURBOMOLE GmbH has adopted an irrevocable bylaw preventing the distribution of dividends to ensure that all profits are reinvested into the project.TURBOMOLE GmbH distributes fee-based end-user licenses itself and through partners, as well as free developer licenses and access to the source code based on project proposals. 91ere we focus on recent developments and provide illustrative applications to chemistry and materials science.For an overview of existing features, as well as development, licensing, and distribution, the reader is referred to refs 92 and 93 and the TURBOMOLE Web site. 91

BRIEF FEATURE OVERVIEW
The program suite consists of a series of modules with a broad range of methods from universal force field to fast semiempirical methods, state-of-the-art DFT and MP2, and coupled-cluster and post-HF methods for ground and excited states.For convenience, the use of modules is facilitated by various tools such as the scripts woelfling, raman, vcd, and genetic.pyfor reaction path optimization, 103 vibrational Raman spectra, 29 vibrational circular dichroism spectra (VCD), 104 and genetic algorithms, respectively. 105Moreover, the graphical user interface TmoleX is of great help for running calculations and visualizing results. 106most all time-consuming parts are parallelized for multicore systems or clusters using OpenMP 107 for sharedmemory parallelization (SMP) 46,62,95−97 and the messagepassing interface 108 (MPI) for parallelization across multiple nodes, 98−102 as outlined in Table 1.The older Fork-SMP 94 is available as a fallback for some modules.Starting with the latest release (V7.7),support for graphics processing units (GPUs) has become available. 109he list of parallelized HF/DFT modules includes those for molecular self-consistent field (SCF) energy (dscf and ridft) and gradient calculations (grad and rdgrad), response properties such as vibrational frequencies (aoforce), NMR/EPR spectra (mpshift), excited-state properties (escf, egrad), and electron transport properties (evib).For these modules, the OpenMP version is recommended for most calculations due to its cost−benefit ratio in terms of computer hardware.Accordingly, post-Kohn− Sham methods (rirpa) and calculations with periodic boundary conditions at the HF/DFT level (riper), as well as the CCSD and CCSD(T) program ccsdf12, are only parallelized with OpenMP.In contrast, the implementation of approximate CC methods in ricc2 and pnoccsd widely supports the OpenMP and MPI standard and combinations thereof for large-scale calculations on multiple nodes.

Local Hybrid Functionals for Strong Correlation and Range-Separated Local Hybrid Functionals.
Local hybrid functionals (LHs) 110,111 with a position-dependent exact-exchange (EXX) admixture governed by a local mixing function (LMF) have been part of TURBOMOLE since release V7.2 and have since been extended in various ways over the past years.Using seminumerical integration techniques, such LHs have been implemented in an efficient way in the code, 18 with functionalities that exceed by far those available in any other quantum-chemistry package that contains LHs.Beyond ground-state SCF 112 and nuclear gradients, 33 this now includes linear-response TDDFT energies 19,113 and excitedstate gradients, 35 frequency-dependent and frequency-independent polarizabilities, 113,114 NMR chemical shifts 115  Fork-SMP 94 and the OpenMP version 46,62,95−97 are restricted to calculations on a single node.MPI 5,98−101 and OpenMP/MPI hybrid 102 implementations allow for the use of multiple nodes.The availability of first-and second-order derivatives as well as excitation energies is also indicated.−128 Here, we focus on two recent extensions.Essentially, the aims of these works have been to conserve the established advantages of LHs and improve other aspects.We start with the fundamental goal to escape the often invoked "zero-sum game" 129,130 between reduced self-interaction errors and delocalization errors or "fractional-charge errors" (FCEs) on one side and minimizing static-correlation errors or "fractionalspin errors" (FSEs) on the other side. 131The enhanced EXX admixture usually helps minimize FCEs, and LHs have been shown to achieve this goal while retaining some of the important left−right correlation in bonding regions. 110On the other hand, larger EXX admixtures usually are detrimental in cases with large FSEs, such as dissociating or stretched bonds or many transition-metal systems with appreciable static correlation.Standard LHs so far have not been a way out of this dilemma, at least not to a larger extent.Relevant real-space approaches to reduced FSEs are Becke's B13 functional 132 and a modified approach by Kong and Proynov (KP16/B13), 133 which have both been implemented self-consistently into a local developer's version of TURBOMOLE. 134ircumventing the numerical difficulties and poor SCF convergence in many cases of the B13 and KP16/B13 functionals, the idea of a local strong-correlation factor has recently been transferred to the LH framework.Initial attempts were still based on relatively simple first-generation LHs but did already show that FSEs can be reduced when multiplying the LH term for nondynamical correlation (NDC) by a somewhat modified KP16/B13-type q AC factor. 134Most recently, the more advanced scLH22t functional has been constructed. 135It is based on the more recent and overall more accurate LH20t functional. 120Using a damping factor for smaller NDC contributions, an almost complete decoupling between the underlying LH20t and the added q AC factor could be obtained.That is, the optimized parameters of LH20t, as well as its excellent performance for weakly correlated situations (e.g., for GMTKN55 main-group energetics), are retained, but FSEs and the related spin-restricted dissociation curves of covalent bonds are dramatically improved. 135otably, the q AC factor forms part of a new LMF (Figure 1).
The second extension of LHs has been the implementation and construction of range-separated local hybrids (RSLHs), combining the ideas of local hybridization in real space and range-separated hybrids in interelectronic distance space. 136,137hat is, instead of full-range semilocal exchange, short-range exchange is mixed in.This corresponds to the use of a longrange EXX admixture governed by a suitable range-separation parameter ω.Based on earlier work for locally range-separated hybrids, 138 RSLHs have been implemented for ground-state SCF (modules dscf/ridft) and gradients (modules grad/rdgrad), as well as for linear-response TDDFT (module escf). 139Using this implementation, the ωLH22t RSLH has been constructed and optimized.It retains most of the good performance of LH20t for main-group and transitionmetal energetics, as well as for core, Rydberg, and triplet valence excitations.At the same time, however, it decisively corrects errors for excitations with appreciable charge-transfer character and improves on the potential-energy curves of three-electron cations known to be affected strongly by selfinteraction errors. 139Most recently, ωLH22t has been demonstrated to provide unprecedented accuracy in quasiparticle energies for applications in molecular electronics and organic photovoltaics, without the usual system-dependent tuning of range-separated hybrids. 140

Inclusion of the Current Density in DFT.
The kinetic energy density τ(r) is a commonly used ingredient in many functionals to detect iso-orbital regions or describe the inhomogeneity of the electron density ρ(r).For its extension τ(r,t), used in the time-dependent Kohn−Sham formalism (TDKS), it has been shown 17 that this quantity is not invariant under a gauge transformation in the external potential.S u b s t i t u t i o n o f τ b y i t s g e n e r a l i z a t i o n 1 4 1 , where j P is the paramagnetic current density, restores gauge invariance.This leads to additional terms in the magnetic orbital rotation Hessian in linear-response TDDFT calculations, accounting for the response of j P .While the original implementation of these terms in TURBOMOLE dates back to 2012 (release V6.4), we highlight four recent important updates.
First, an incorrect prefactor 2 has been removed from the original implementation with V7. 6. 119,144 TDDFT calculations employing τ-dependent functionals performed with previous versions erroneously overcorrected for the effects of the current density response and should be reassessed, although average changes are on the order of 0.03 eV. 144econd, a recent investigation 145 reveals that the effect restoring gauge invariance has on the final excitation energies can be significantly larger than previously assumed depending on the functional and type of excitation.In one particular investigation of d−d excitations in nickel(II) complexes, restoring gauge invariance shifts the excitation energies with the M06-2X functional by more than 0.4 eV closer to the experimental reference values, as shown in Figure 2. 145 A broader analysis reveals that the importance of imposing gauge invariance can be linked to the derivative of the exchange energy integrand with respect to τ. 145 Moreover, n → π* excitations are significantly more affected by restoring gauge invariance than most π → π* excitations with the exception of π → π ⊥ * excitations, where the dominantly contributing molecular orbitals (MOs) are perpendicular (⊥) to each other. 145These findings suggest that a reassessment of previously reported TDDFT results obtained with τ-dependent functionals is warranted, particularly for cases that are potentially more sensitive due to the choice of the functional, the type of excitation, or both.Gauge invariance is restored by default with τ-dependent functionals at moderate (nonhybrid functionals) or negligible (hybrid functionals) additional computational cost. 17Recently, excited-state gradients and quadratic response properties for τ-dependent meta-generalized gradient approximations (mGGAs) have been implemented, enabling gauge invariant computations of excited-state equilibrium structures, relaxed dipole moments, (dynamic) hyperpolarizabilities, and two-photon absorption cross sections. 147These developments will be available in a future TURBOMOLE release.
Third, the inclusion of j P is required for gauge invariance of magnetic properties and implemented for magnetizabilities, 119 NMR coupling constants, 119 and NMR shifts of closedshell 119,148 and open-shell systems, 118 as well as EPR hyperfine coupling constants 149 and g-tensors. 118Recent findings indeed hint at the inclusion of the current density response also being crucial for NMR and EPR properties. 118,119,148Especially for open-shell systems, neglecting the current density leads to large deviations, 118 as shown in Figure 3.
Finally, it was shown that the inclusion of the paramagnetic current is also crucial in relativistic two-component calculations. 151Contrary to all previous cases, the modifications outlined above must already be taken into account for the ground state in the presence of spin−orbit coupling (SOC).SOC gives rise to an internal magnetic field, inducing a paramagnetic current that already has a nonvanishing contribution at the energy-level. 151Accordingly, if properties such as light−matter interactions are targeted in the presence of SOC, the interplay of the ground-and excited-state paramagnetic currents must be taken into account.This gives rise to highly nonlinear j P -dependent terms, which have a profound impact on many properties. 151,152Current density functional theory (CDFT) for τ-based functionals is available for ground-state energies (module ridft) and gradients (module rdgrad) and linear response properties (modules escf and mpshift). 151Given the profound impact of j P , we therefore strongly recommend the exclusive use of the current-dependent forms.

Methods for Finite Magnetic Fields.
−155 For these applications, the magnetic field is usually treated perturbatively, as it is orders of magnitude smaller than the electronic interactions responsible for the formation of a chemical bond.
−167 Consequently, spectra obtained from such interstellar objects can only be interpreted using quantum-chemical calcula-   150 We list the mean signed percentwise deviation (MSPD), the mean absolute percentwise deviation (MAPD), and the root-mean-square deviation (RMSD).Reprinted from ref 118 under a CC BY license.Copyright 2022 the Authors.tions 168−172 In such extreme conditions, entirely new types of chemical bonding, spin-phase transitions, and other exotic phenomena have been shown to occur. 162,173,174The effects of arbitrarily weak or strong magnetic fields on atoms or molecules may be computed using the finite magnetic field approach, which is implemented for Hartree−Fock, 175 CDFT, 152,176 GW/BSE, 177,178 RPA, 179 and CC2. 178Through the calculation of electronic ground states, molecular gradients, and electronic excitations, a wide variety of applications for molecular spectroscopy in magnetic fields are now accessible.
Due to the efficiency of our implementation, systems containing dozens of atoms can be routinely computed in explicit magnetic fields. 152,175,176,178To highlight the capabilities of our approach, we calculated the MCD spectrum of ZnDiNTAP, 176,180 a zinc tetraazaporphyrin with two fused naphthalene units, Figure 4a, using CAM-B3LYP and a mixed def2-TZVP (Zn)/def2-SVP (all other atoms) basis set.Similarly to the experiment, 180 Figure 4b, a magnetic field of 5 T was applied.The resulting spectrum is shown in Figure 4c.
Minor differences can all be attributed to solvation effects and vibronic coupling, which were neglected in our calculation. 176urthermore, the MCD spectrum of ZnDiNTAP in an explicit magnetic field of 1000 T is shown in Figure 4c.While some of the bands, particularly in the fingerprint region, are not affected by nonlinear effects induced via such a strong field, the two Q bands are significantly shifted (650 → 693 nm and 536 → 514 nm).
Other applications of the finite magnetic field ansatz include higher-order properties such as magnetizabilites and hypermagnetizabilities through the use of numerical derivatives.
Moreover, molecules that are naturally prone to exhibiting "magnetic" effects, such as aromatic compounds (see also section 3.7), may show a nonlinear response even to weak magnetic fields. 152This is not captured by a perturbative approach but can be routinely investigated using our implementation.

EPR Properties and Single-Molecule Magnets.
Over the past decade, TURBOMOLE has pioneered the in silico study of f-element chemistry.To provide accurate descriptions of compounds containing these heavy elements, robust DFT routines are leveraged to deliver a balanced treatment of dynamic correlation, static correlation, solvation, and relativistic effects.More recently, these developments have enabled the discovery of new species with novel chemistry, subsequently necessitating new and improved computational methods capable of describing them.
For example, a series of recently discovered Ln-based singlemolecule magnets (SMMs), [La(OAr*) 3 ] − , [Lu(NR 2 ) 3 ] − , and [Lu(OAr*) 3 ] − (OAr* = 2,6-Ad 2 -4-t-Bu-C 6 H 2 O, Ad = adamantyl, t-Bu= tert-butyl, R = SiMe 3 with Me = methyl), were found to exhibit exotic EPR properties, with [Lu-(OAr*) 3 ] − producing hyperfine coupling (HFC) constants of unforeseen magnitude and furthermore demonstrating extended magnetic coherence facilitated by a hyperfine clock transition.A primary investigation of these species with nonrelativistic HFC operators attributed the large hyperfine coupling constants to a sizable Fermi contact contribution from the highest occupied MO (HOMO) of each system.However, the predictions of hyperfine coupling constants and g-tensor values themselves produced errors of roughly one  order of magnitude, strongly suggesting the need for more rigorous methods. 181uch improved predictions of EPR spectra are possible with relativistic exact two-component (X2C) theory, 38,39 including spin−orbit effects up to the noncollinear two-component (2c) DFT framework. 38,39,117For meta-generalized gradient approximations and local hybrid functionals, this also includes the paramagnetic current density in the ground state (cf.section 3.2). 151This rigorous formulation can be truncated to the scalar-relativistic limit 63 or a perturbative ansatz 149 to study the individual contributions of each term to the EPR parameters.Gauge origin invariance of the g-tensor calculations is ensured by the gauge including atomic orbitals, 36,39,149 which are crucial for systems with a spatially distributed spin density.These methods were implemented, and their performance is shown here for the aforementioned SMM [Lu(OAr*) 3 ] − shown in Figure 5.The all-electron relativistic methods lead to good agreement with the experiment, as shown in Table 2.The HFC constant of [Lu(OAr*) 3 ] − is dominated by the scalar-relativistic contribution due to the localization of the spin-density in the 6s/d HOMO producing a large Fermi contact interaction.The importance of the paramagnetic spin−orbit contribution increases with the number of unpaired electrons and the scalar formulation, as well as the spin−orbit perturbation theory (SOPT) break down for systems such as [TbPc 2 ] − with six unpaired electrons, Pc = bis(phthalocyaninato). 38,39,117,151he self-consistent 2c methods are thus pivotal.
With the next release version, the EPR Euler transformations for the HFC, g-tensor, and electric-field gradient as well as the nuclear quadrupole interaction tensor will further become available for users. 39

NMR Coupling Constants Across the Periodic Table of Elements.
NMR spectroscopy is key to the analysis and structure determination not only for organic compounds but also for inorganic systems consisting of heavy elements.The NMR spin−spin coupling constant describes the splitting of the signals or peaks in the NMR spectra and is a characteristic property driven by the chemical environment.Formally, the coupling tensor is obtained as the mixed derivative of the energy with respect to the corresponding nuclear magnetic moments, which are introduced via the principle of minimal coupling.NMR couplings are available within a nonrelativistic 60 scalar X2C 183 and the spin−orbit X2C framework. 86All functional classes up to local hybrids are supported and include the current density for gauge invariance. 61,119,151or systems made up of light elements, the nonrelativistic treatment is sufficient.Here, the coupling constant is generally partitioned into the Fermi-contact (FC), spin−dipole (SD), paramagnetic spin−orbit (PSO), and diamagnetic spin−orbit (DSO) contributions.The FC, SD, and PSO terms necessitate the solution of the response equations, whereas the DSO term is computed with the ground-state density.Typically, the FC contribution is the leading term, and accurate coupling constants require large basis sets.Thus, a nuclear selection scheme and locally dense basis sets are often applied to largescale calculations.
Systems containing heavy elements such as Sn, Pb, Pd, and Pt require the inclusion of relativistic effects, 61,86,183,185−187 i.e., methods based on the Dirac equation are introduced.For such methods, the FC, SD, and PSO terms are coupled due to spin−orbit interaction, and they come with drastically increased computational demands.Nevertheless, when using a local X2C ansatz, 84−86 large-scale calculations are possible, as illustrated in Figure 6 for the Karplus relationship of Me 3 Sn− CH 2 −CHR−SnMe 3 , where R is different substituents (Me = CH 3 ).The relativistic DFT approach reproduces the experimental findings with fairly good agreement.Improvements are possible with the correlation kernel augmented BSE (cBSE) based on the Green's function GW ansatz.Here, the DFT response equations are replaced with their BSE counter parts. 61,113o demonstrate the efficiency, the calculation of the Sn−P coupling constants of [({SIDipp}P) 2 Sn] (SIDipp = 1,3bis(2,6-di-isopropylphenyl)-imidazolidin-2-ylidene) with 137 atoms 187 takes about 44 min (PBE) and 55 h (PBE0) using 12 OpenMP threads of an Intel Xeon Gold 6212U CPU (2.40 GHz). 86nd a scalar-relativistic treatment 36,104 are available for the orbital contribution, including the response of the current density. 118,119,148The paramagnetic NMR (pNMR) shielding tensor σ I tot for a nucleus I reads with S denoting the spin, μ e denoting the Bohr magneton, γ I denoting the gyromagnetic ratio of nucleus I, k B denoting the Boltzmann constant, and T denoting the temperature.Here, the orbital contribution σ orb is the straightforward open-shell generalization of the closed-shell limit. 63,118Additionally, a temperature-dependent contribution arises, which includes the HFC tensor A I of nucleus I and the g-tensor g already discussed in section 3. 4. Both the HFC and the g-tensor depend on spin−orbit coupling.For the calculation of 1 H/ 13 C pNMR spectra of large molecules, a perturbative treatment of spin− orbit coupling is preferred over the 2c ansaẗze due to lower computational costs. 149The viability of the pertubative ansatz in the X2C framework is demonstrated for two negatively charged Ru(III) complexes in Figure 7, which depicts the good agreement between calculated results and the experimentally measured 189 pNMR 1 H and 13 C shifts of the two compounds.
The calculation of properties depending on the density in the vicinity of the nuclei requires additional tight basis functions.Thus, the pcJ, 190 pcS, 191 and pcH 192 basis sets are recommended for nonrelativistic calculations.For relativistic calculations, the x2c-s basis sets were developed. 193,194he efficiency of the pNMR implementation is similar to that of the closed-shell case, as shown in Figure 8. Coulomb integrals can be calculated with the RI-J/MARI-J approximations for the Coulomb contribution 62,63 and the seminumerical scheme for exchange integrals. 18,109.7.Ring Currents of Heavy-Metal Clusters.Aromatic compounds, such as benzene, show characteristic signals at 7 ppm in the 1 H NMR spectra.This shift is a consequence of the cyclic electron delocalization associated with the π-orbitals, which deshield the nuclei due to an induced ring current. 195his magnetically induced current density may be calculated indirectly with the nucleus-independent chemical shift 196 (NICS) or directly using TURBOMOLE's interface to the GIMIC program, which was reworked for release V7.−204 The occurrence of ring currents and the concept of aromaticity are not restricted to cyclic conjugated hydrocarbons and related organic compounds.All-metal systems may also sustain a ring current, and these systems are therefore classified as all-metal aromatic compounds.For instance, the endohedral [Th@Bi 12 ] 4− cluster features a nonlocalizable πorbital around the {Bi 12 } torus, which leads to a ring current. 205Figure 9 shows a streamlined representation of this ring current, whose strength amounts to about 25 nA/T . 119,205spite featuring only two delocalized π-electrons, almost the same ring current strength as that in porphine is induced.
Here, the thorium atom in the torus center is only needed for stability and the synthesis, but it does not take part in the ring current.Furthermore, the open-shell variant [U@Bi 12 ] 3− shows a strong ring current, 63 as the same π-orbital is occupied by two electrons. 205,206Moreover, prismatic {Bi 6 }-based clusters such as [(CpRu) 3 Bi 6 ] − show ring currents of more than 25 nA/T. 207Therefore, all-metal clusters help to push the frontiers of aromaticity.

Characterization of Novel Electronic Configurations in f-Block Elements by DFT and RPA Methods.
Electronic structure calculations of large f-element complexes contain a great deal of complexity due to the competition between metal oxidation states, d-and f-shell occupations, spin   coupling, and relativistic effects. 208Computational studies in these compounds have to scan a range of formal electronic occupations and apply a range of techniques to ensure convergence to the desired electronic configuration, including Fermi smearing with suitable parameters, 209,210 and a combination of damping, level shifting, 211 and direct inversion in the iterative subspace (DIIS) extrapolation. 212Moreover, the stability of the ground-state reference is a concern in calculations of molecular properties, for example, electronic absorption spectra. 213FT results helped characterize the structures and properties of nontraditional Ln 2+ complexes possessing the 4f n 5d 1 configuration in the [Ln(C 5 H 4 SiMe 3 ) 3 ] − (Ln = Ce−Nd, Gd− Er), 214−217 [Ln{N(SiMe 3 ) 2 } 3 ] − (Ln = La, Gd), 218,219 and [Ln(Cp iPr 5 ) 2 ] (Ln = Tb, Dy) series. 220The preference for the 4f n 5d 1 configuration relative to the traditional 4f n+1 occupation of Ln 2+ results from the stabilization of the Ln 5d z 2 orbital by the trigonal ligand environment or extremely bulky ligands.Nontraditional Ln 2+ complexes show a characteristic intense absorption band in the visible range due to excitations from the occupied Ln 5d orbital.The prediction of UV−vis spectra of low-valent lanthanide complexes, in particular those with a nontraditional configuration, is improved by including diffuse augmentation in lanthanide basis sets. 221DFT and RPA methods were employed to examine the strong ferromagnetic coupling between the Ln 3+ centers in [(C 5 Me 5 ) 2 Ln(μ-S) 2 Mo(μ-S) 2 Ln(C 5 Me 5 ) 2 ] − (Ln = Y, Gd, Tb, Dy) and the Mo→Ln electronic excitations in the near-infrared spectral region. 222Excited-state studies of [Ln(C 5 Me 5 ) 2 (C 5 Me 4 H)] and [Ln(C 5 Me 5 ) 2 (η 3 -C 3 H 4 )] complexes (Ln = Y, Dy, Lu) using TDDFT elucidated their unexpected photochemical activation, which was used to reduce dinitrogen and sulfur and to polymerize isoprene. 223,224omputational studies of neutral actinide complexes [An-(Cp iPr 5 ) 2 ] (An = Th, U, Pu, Am, Bk, No, Lr) (pentaisopropylcyclopentadienyl = Cp iPr 5 ) using DFT predicted ground states with a linear ligand arrangement of S 10 symmetry and significant An 6d orbital occupation for An = Th, U, Lr. 225 The calculations were carried out with the TPSS exchange− correlation functional, 226 Stuttgart−Cologne scalar-relativistic small-core effective core potentials (ECPs), 227 and the corresponding basis sets. 228,229Mixed 5f/6d occupation was predicted in the corresponding Pu complex, while the An = Am, Bk, No complexes were found to have 5f n+1 configurations.The Pu and Am complexes showed a slight deviation from the perfectly symmetric structure, with Cp−M−Cp bending angles of 11°and 12°, respectively.The simulated absorption spectra showed intense peaks in the UV−vis range due to the metal− ligand charge transfer excitations from the An 6d orbital shown in Figure 10.Comparisons with the previously experimentally known Ln analogs (Ln = Dy, Tb) 220 suggested that the synthesis of the predicted actinocene complexes was thermodynamically feasible.The computational predictions received experimental confirmation for An = U while the results were still under review.Layfield and co-workers reported the synthesis of the linear S 10 -symmetric "secondgeneration" uranocene [U(Cp iPr 5 ) 2 ]. 230 The U−Cp centroid distance was determined from crystallographic studies as 2.504 Å, in good agreement with the computational result of 2.483 Å.The measured UV−vis spectra showed broad and intense absorption, as predicted by TDDFT calculations.
3.9.Damped Response Theory.Assessing light−matter interactions in extended or exotic systems is an active field of development in TURBOMOLE. 53,119,152,176,178,231,232While root-by-root linear-response methods have been crucial for many applications, they are not suitable for spectrally dense systems or core excitation due to the high number of excited states (roots).Damped response theory, 233,234 or the equivalent complex polarization propagator approach, 235 provides a convenient framework to formulate resonance convergent response functions, circumventing these problems.It provides a convenient route to directly compute a variety of absorptive and dispersive effects in both UV−vis and X-ray frequency regions, which is particularly advantageous for large systems, and in frequency regions with high densities of states, as it does not require to solve eigenvalue equations for all contributing states and individual transitions matrix elements between them. 236−238,241−254 TURBOMOLE is, to the best of our knowledge, the only program to date that offers damped linear response functions at the RI-CC2 level of theory 232,253 (available since release V7.7).
3.9.1. Damped Response for Multiscale Modeling.DFTbased damped response implementations in TURBOMOLE cover IR spectroscopy, VCD, and Raman spectroscopy, as well as absorption and electronic circular dichroism (ECD) in the visible and ultraviolet spectral range.Furthermore, this approach has recently been extended to the modern framework of the GW-BSE method, being especially useful for core excitations. 113,255t a given complex frequency ω ex = ω + iγ, where ω and γ are the real and imaginary parts of the external field, respectively, the coupled perturbed equation 254,256 is first solved.A and B are defined as ε p marks orbital or quasiparticle energies, v pq,rs is a Coulomb integral, f pq,rs is the exchange-correlation kernel, if present, and K pq,rs is a exchange integral.The precise kinds of f pq,rs and K pq,rs depend on the chosen method. 254p ν describes the external perturbation, e.g., an electric or magnetic field. 254,256The polarizability can then be calculated as Damped polarizabilities and the related magnetizabilities have further emerged as the quantum-chemical cornerstone of the transition matrix (T-matrix) based approach for multiscale modeling of light−matter interactions. 254,257,258This way, the functionalization of molecular structures within optical devices is possible.In that regard, optical cavities filled with molecular materials or metasurfaces of cylinders consisting of molecular materials can be designed.These devices exhibit tailored optical properties for a variety of applications, such as enhancing the circular dichroism of a chiral molecule.Through the multiscale approach depicted in Figure 11, the properties of a molecular unit and the macroscopic sample can be distinguished and combined to achieve a specific effect.Combining TURBOMOLE and, for example, the multilayered periodic general Mie method (mpGMM), 259 simulations and predictions of the light−matter interactions of layer-structured materials ranging from a few to hundreds of nanometers are now routinely possible. 254,260To target molecular materials of arbitrary shape, the T-matrix approach was furthermore coupled with state-of-the-art homogenization techniques. 261ombining classical electrodynamics and quantum mechanics has proven to be a worthwhile approach in the field of light− matter interactions, and TURBOMOLE will remain at the forefront of these developments.
3.9.2. Damped Response Theory for One-Photon Absorption and CD Spectra with RI-CC2.For RI-CC2, the equations for the damped response 237,262 of the cluster amplitudes are recast in a form that only involves effective matrices in the space of single excitations, avoiding the storage of parameters for the double excitation space.
1 A 1 Subscripts Re and Im represent real and imaginary components, respectively.The effective matrices in the equations above are where A SS , A SD , and A DS are, respectively, the singles−singles, singles-doubles, and doubles-singles blocks of the CC2 Jacobi matrix.
frequency-shifted orbital energy differences, and ω and γ are again the real and imaginary parts of the frequency of the external field, respectively.The effective right-hand sides are where ξ S x and ξ D x are the single and doubles parts of the righthand sides in the nonpartitioned form, respectively. 232,253The partitioned formulation that avoids the need to store double excitation amplitudes and four-index integrals allows applications to system sizes otherwise not accessible at the CC2 level.
As an illustrative application, we computed the UV−vis absorption spectrum and electronic circular dichroism spectrum of a donor−acceptor cyclophane 263 shown in Figure 12.The absorption spectrum was obtained from calculations of the imaginary damped dipole polarizability and the ECD spectrum from the imaginary part of the optical rotation tensor in the velocity gauge.The asymmetric form 253 of the damped linear response function was used in the calculations (the symmetric form is also available).
Cyclophanes are well-studied examples of strained aromatic organic compounds (hydrocarbons) that exhibit planar chirality.The UV−vis spectrum is a typical example where the traditional state-wise approach converges only slowly with the number of states, as seen by the large difference between the results for 14 and 59 states that have been included for comparison.The problem does not appear in the damped response calculations.
Ongoing work is concerned with extending the implementation to induced and nonlinear spectra like MCD and RIXS. 253.10.Nonadiabatic Molecular Dynamics Simulation for Spectroscopic Observables.Many photophysical and photochemical processes involve multiple electronic excited states coupled by radiative and nonradiative transitions.Efficient simulations of these processes by nonadiabatic molecular dynamics (NAMD) have recently become possible with transition dipole moments 30 and nonadiabatic couplings between excited states 32,264 computed within quadratic response TDDFT.The TURBOMOLE implementation of the multistate fewest-switches surface hopping (FSSH) algorithm enables simulations of molecular systems with 50− 100 atoms and simulation times of >10 ps. 265,266In addition, time-resolved fluorescence (TRF) and transient absorption (TA) 34 spectra can be simulated for comparison with experimental results.
The simultaneous treatment of multiple electronic excited states enables the examination of Kasha's rule. 267According to Kasha's rule, singlet excited states energetically located above S 1 undergo ultrafast decay to the S 1 state and thus are not directly involved in fluorescence or photoactivated reactions.However, exceptions to Kasha's rule are well-known in molecules in the gas phase, such as azulene and pyrene. 268n our recent study, we used the NAMD implementation in TURBOMOLE to investigate the dynamics of several polycyclic aromatic hydrocarbons including pyrene, azulene, and bicyclo[6.2.0]decapentaene (BCDP, an isomer of azulene) at the PBE0/def2-SVP 269,270 level.Azulene was found to exhibit non-Kasha behavior due to emission from the S 2 state, in agreement with experiment. 271BCDP obeys Kasha's rule and emits only from S 1 .−275 Multistate NAMD simulations describe the non-Kasha behavior as a combination of S 1 → S 2 transitions and the change in the diabatic character of the S 1 state from dark (L b ) to bright (L a ) at points of near degeneracy between the S 1 and S 2 states.The high-energy shoulder in the fluorescence spectrum of pyrene can be understood as Stick spectra are the results from state-wise calculations, and the red and orange lines are obtained by convoluting these computed stick spectra with a Lorentzian broadening with a half-width-at-halfmaximum of ≈1000 cm −1 including, respectively, the lowest 14 and 59 states.originating from excited states with diabatic bright (L a ) character.The S 2 lifetime in pyrene was computed by an exponential fit of the decay of the state population as 63 fs, in agreement with the experiment value of 85 fs in methanol. 276he S 1 lifetime of azulene was computed to be 2.2 ps in comparison to the experimental result of 1.4 ps in cyclohexane. 277The computed lifetime of the S 1 state of BCDP was found to be 0.8 ps.
The NAMD trajectories were also used to obtain the TA spectrum of pyrene as an ensemble average of the Gaussianbroadened excited-state absorption and emission spectra (Figure 13).The experimental TA spectrum of pyrene shows an intense band at 580 nm and the growth of a steady state signal at around 450 nm corresponding to S 2 → S n transient absorption that decays rapidly and a S 1 → S n transition (n > 4), respectively. 276,278NAMD simulations predict S 2 → S 4 and S 1 → S 4 bands at 1500 and 800 nm, respectively.The time evolution of the S 1 and S 2 states is in good agreement with experiment, while the absorption maxima (λ(S 2 → S 4 ) = 1500 nm and λ(S 1 → S 4 ) = 800 nm) are red-shifted due to truncation of the electronic excitation space.
3.11.Generating Function Methods for Vibrationally Resolved Electronic Spectroscopy.The vibrational structure of electronic spectra gives detailed information about molecular structure 279 and excited-state phenomena such as internal conversion and intersystem crossing. 280,281−287 However, this approach fails for rigid molecules with high vibrational frequencies. 288In this case, special care must be taken to include the quantum nature of nuclear vibrations.−292 The radless module 288,293,294 makes use of the generating function method to compute vibrationally resolved absorption and emission spectra, as well as photoelectron ionization spectra.Spectra can be computed within the global harmonic approximation, which only requires equilibrium geometries for initial and final structures as well as vibrational spectra of both structures.The method accounts for the full Duschinsky rotation, 295 taking into account differences in initial and final state structures and vibrational modes.Due to its efficiency, the method can be applied to large molecules, such as polyaromatic hydrocarbons. 123,124,293,296,297An extension of the module further allows the computation of emission and absorption spectra arising from singly occupied vibrationally excited initial states, allowing the simulation of single vibronic level (SVL) fluorescence 293 and vibrationally promoted electronic resonance (VIPER) spectra. 298n addition, the newly implemented semiclassical thawed Gaussian approximation (TGA) 291,299 offers an efficient way of evaluating the time-correlation function without resorting to the global harmonic approximation.The relation between vibronic spectroscopy and semiclassical dynamics stems from the wavepacket autocorrelation picture of the dipole time correlation function, first popularized by Heller. 290In TGA, an initial Gaussian wave function is evolved under an effective local harmonic potential constructed at each step about the center of the wavepacket.As a result, its Gaussian form is conserved; the center of the wavepacket follows a classical trajectory, while its width is adjusted according to the instantaneous Hessian of the potential energy surface (PES).Whereas in the original ab initio TGA 300−303 the Hessian of the potential energy is updated over time, in the single-Hessian version, 304−306 implemented in TURBOMOLE as part of the radless module, 294 the Hessian is kept constant throughout the dynamics.Therefore, the overall additional cost compared to the conventional harmonic approximation is that of a single ab initio trajectory in the final electronic state, which is simulated using the frog module.Since the trajectory experiences the true anharmonic PES, the method can account for anharmonicity at least approximately, although it cannot describe more subtle quantum dynamics, such as wavepacket splitting or tunneling.The TGA approach has proven to be especially useful in systems with a large displacement between the ground-and excited-state minima and in systems with a double-well-shaped PES along a low-frequency mode in the final electronic state.In such molecules, the harmonic approach typically fails because the global harmonic PES constructed around one of the wells is not adequate.Moreover, in contrast to the global harmonic methods, the TGA results often depend weakly on the choice of the Hessian, as illustrated in Figure 14.Overall, the implementation in TURBOMOLE combines these advantages of TGA with accurate and efficient excited-state electronic structures, such as ADC(2) and CC2 methods.
3.12.Molecular Properties from Self-Consistent GKS-spRPA.The generalized Kohn−Sham semicanonical projected random phase approximation (GKS-spRPA) provides a route for obtaining one-particle energies at the RPA level of theory. 16hese one-particle energies provide accurate estimates of ionization potentials (IPs) and electron affinities due to a correct description of orbital correlation and relaxation effects.Its computational cost was reduced from N ( ) 6 to N ( ) . Experimental absorption spectrum of 1,2,4,5-tetrafluorobenzene (black) compared with the spectra simulated using the (adiabatic Hessian, AH) global harmonic method ("Harmonic", blue dotted line) and single-Hessian TGA with either the adiabatic Hessian (AH, red dashed line, evaluated in the excited electronic state at its corresponding optimized geometry) or the initial Hessian (IH, green dash-dotted line, ground-state Hessian computed at the ground-state optimized geometry).All electronic structure calculations, including geometry optimization, energies, and forces for the ab initio dynamics and Hessians, were performed at the CC2/def2-TZVP level.Adapted with permission from ref 294.Copyright 2022 American Chemical Society.
using well-known analytic continuation (AC) techniques. 307−309 3.12.1.Applications to Nonvalence Anionic States and Xray Emission Spectroscopy.The AC version of GKS-spRPA retains a high accuracy across energy scales from valence to core-ionization energies.The versatility of the AC GKS-spRPA was shown by its application to problems involving very weakly bound anionic states and very strongly bound core-hole ionization energies.
Nonvalence anionic states of molecules are weakly bound states where the excess electron occupies a diffuse orbital.The excess electron is bound by a combination of long-range electrostatic and correlation effects.When electrostatics are sufficient to bind an excess electron, the anionic states are referred to as nonvalence electrostatic-bound (NVEB), and when correlation effects are necessary the states are referred to as nonvalence correlation-bound (NVCB).Using a model water tetramer cluster, the GKS-spRPA approach was shown to provide electron affinities (EAs) from the lowest unoccupied molecular orbital (LUMO) energies that were within 10 meV of those provided by the EOM-CCSD(T)a* approach, see Figure 15.The high accuracy of the GKS-spRPA approach is due to the correct description of long-range correlation effects.The AC approach introduces errors of less than 5 meV for nonvalence states.
The AC GKS-spRPA method can also be used for simulating nonresonant X-ray emission (XE) spectra using just the information from the one-particle eigenspectrum. 310XE energies, ΔE, computed by taking the difference between core and valence IPs and the oscillator strengths, f osc., are evaluated within a frozen orbital approximation where ϕ c and ϕ v are the core orbital and valence GKS-spRPA orbitals involved in the XE process, respectively.The AC GKS-spRPA approach was used in conjunction with the scalarrelativistic (SR) X2C approach and uncontracted basis sets to estimate highly accurate XE spectra for molecules containing second-and third-period elements, for example, see Figure 16.
Using uncontracted basis sets, the XE energies were found to have MAEs of 0.7 eV for both second and third period-based XE.The X2C-based AC GKS-spRPA approach thus enables the simulation of nonresonant X-ray emission within a simple one-particle picture while avoiding the use of empirical shifts or core-hole reference states.The latter is an appealing aspect, since issues related to variational instability, which are present in core−hole reference based methods, are avoided.

Orbital Ordering in Quinacridone.
The accuracy of GKS-spRPA 16 for ionization potentials (IPs) compared to nonselfconsistent RPA 23 is due to partial satisfaction of functional self-consistency, which requires that the Kohn− Sham (KS) density equal the interacting density defined as the functional derivative of the ground-state energy with respect to the external potential. 16,317A comparison of the GKS-spRPA orbitals and orbital energies to experimental photoelectron spectra of quinacridone illustrates this point, see Figure 17.Within both the KS and the GKS approach, the HOMO energy equals the first IP and subsequent lower-lying orbital energies approximate higher principal IPs. 318,319However, this is often not true for semilocal density functional approximations.PBE 313 predicts a HOMO(−1), HOMO(−2), and HOMO(−3) ordering inconsistent with the results of optimally tuned range-separated hybrid functional (OT-RSH) calculations, yielding an accurate description of experimental photoelectron spectra and G 0 W 0 IPs. 320,321The canonical GKS-spRPA orbital ordering qualitatively agrees with the one OT-RSH one down to HOMO(−3), see Figure 17.and the SR-GKS-spRPA method. 310The computed spectrum (in blue) was obtained by Gaussian broadening of the vertical transitions (red lines) using a width parameter of 1 eV.The vertical dashed lines denote the two intense peak positions for SR-GKS-spRPA.Reprinted and adapted with permission from ref 310.Copyright 2022 American Chemical Society.

Subquadratic-Scaling PNO-CCSD(T) and PNO-CCSD(T)-F12
BCCD(T)(F12*). 68,69 When combined with a sufficiently large basis set, these methods deliver ground-state energy differences for reaction or interaction energies and potential energy surfaces to an accuracy that enables a quantative comparison with experimental results.The characteristics of the TURBOMOLE implementation are low memory and disk requirements and shared-memory parallelization, which make calculations on systems with ∼20 atoms routinely possible on modern machines.
To evaluate energies for larger molecules, the steep N ( ) scaling of the costs with system size must be overcome.Local approximations based on the short-ranged nature of electron correlation in molecules provide a route to near-linear scaling of costs with the system size.The PNO approach uses an approximate MP2 pair density to construct the set of local virtual orbitals for each pair that is best suited to describing the correlation of that pair.This, in combination with screening and local density fitting approximations, among others, reduces the scaling from N ( ) 7 to N ( ) in the asymptotic limit.For ground-state energies, PNO-MP2, PNO-CCSD, and PNO-CCSD(T) explicitly correlated variants PNO-MP2-F12, PNO-CCSD(F12*), and PNO-CCSD(T)(F12*) have been available since V7. 6. 70−73 The efficiency of the implementation is greatly improved if the PNOs are expanded in projected atomic orbitals (PAOs) rather than directly in terms of atomic orbitals (AOs).This is known as domain-based local PNO theory (DLPNO).In TURBOMOLE, the PAO domains for each pair are determined also using the approximate MP2 density rather than by analyzing the MO coefficients and are consequently much more compact.The resulting domains are called principle domains. 323Very recently, the principle domain approach has been extended to F12 theory, where principle domains and PNOs are required for every subspace in the F12 strong orthogonality projector. 324In contrast to implementations in other software packages, we do not use the simplified "A" approximation 95 for the MP2-F12 contributions.This has the benefit that the energies converging smoothly to both the canonical limit and the basis set limit can be extrapolated using PNO extrapolation techniques. 325o illustrate the characteristics of the TURBOMOLE implementation and what is now possible, we report timings of PNO-CCSD(T) and PNO-CCSD(T)(F12*) calculations on a sequence of alkane chains and rock salt crystal fragments in Figure 18.The default tight (10 −7 ) PNO truncation threshold was used in all cases.The def2-TZVPP basis and def2-TZVP were used for C n H 2n+2 and Na n Cl n , respectively, and calculations were run on a 48 core Intel processor with 200 Gb of memory.For the linear systems, the observed scaling is subquadratic, with C 128 H 256 taking 15 h to complete.For the globular systems, the observed scaling is subcubic, with Na 50 Cl 50 taking 45 h to complete.The F12 calculations are 2− 3× more costly than non-F12 calculations but provide energies

Journal of Chemical Theory and Computation pubs.acs.org/JCTC
Perspective close to the basis set limit without requiring basis sets with a large number of AOs per atom.

Real-Time TDDFT for Molecules.
While linearresponse TDDFT can be used to study excitations under weak fields, a better understanding of nonlinear excited state dynamics at the femto-and sub-femtosecond time scales requires using RT-TDDFT.
RT-TDDFT for molecular systems has been implemented 89 within the riper module [55][56][57][58]327,328 using the Magnus propagator 329 and the predictor-corrector scheme 330 for time integration.The implementation utilizes density fitting and continuous fast multipole method (DF-CFMM) techniques 56 to speed up the KS matrix evaluation and scales almost quadratically with system size. Previously the implementation was used to simulate absorption spectra over a wide range of frequencies by perturbing the molecules with weak electric fields. 89 ecently, the code has been extended to simulate high harmonic generation (HHG) spectra under intense laser pulses and was utilized to complement and analyze the pulse-induced electron dynamics in the organic semiconductor molecules tetraphenylporphyrin (TPP) and zinc-tetraphenylporphyrin (ZnTPP). 326Figure 19 shows the good agreement of the experimental absorption and HHG spectra with the calculated ones.The difference in the first few harmonics is due to the fact that the experimental spectra contain contributions from the quartz substrate. Ovall, the HHG and absorption spectra of both the porphyrins are quite similar due to the very similar electronic structure.Our simulations combined with experiments showed that π → π* excitation plays a major role in the harmonic generation process in porphyrins.It was also discovered that resonant excitation leads to an early onset of nonperturbative behavior for the fifth harmonic, and similar effects are expected in Brunel harmonic generation with other organic materials.326 3. 15. Deelopments of the DFT-Based Embedding Implementations.3.15.1. Froen Density Embedding Implementation.Environment effects on molecular systems beyond static Coulomb potentials can be treated in TURBOMOLE with the conductor-like screening model (COSMO), polarizable embedding (PE), or frozen density embedding (FDE).Developments on COSMO and PE in TURBOMOLE have been presented in the previous review.92 In contrast to the latter approaches, FDE is an atomistic QM/ QM embedding model, which compared to QM/MM schemes uses a purely quantum-mechanical description of the total system and does not require any system-specific parametrizations prior to a calculation.
The implementation of FDE in TURBOMOLE was until recently restricted to just two subsystems, limiting its applicability.It has now been extended to handle arbitrarily many subsystems and to an FDE variant free of intermediates evaluated on the supermolecular orbital basis.For the embedding potentials, the newly implemented approach uses only electron densities and electrostatic potentials of the subsystems, which are computed on an integration grid generated for the total system.It is available for HF and DFT within the dscf and ridft programs.In combination with the ricc2 program, it can be used to compute groundstate energies with MP2 and CC2 within the perturbation to the energy (PTE) or post-SCF reaction field coupling approaches.Excitation energies can be calculated within the frozen solvent approximations either with or without the kernel contributions.An important feature of the new implementation is the possibility to include pseudopotentials to improve the embedding potential.This allows for a more accurate description of the Pauli repulsion, 83 which is particularly important for electronically excited states to prevent unphysical charge spill-out.Besides all of the new features in FDE, it should also be mentioned that the current implementation is still limited to closed-shell subsystems.
To demonstrate its capabilities, we computed the first eight excitation energies of acetone at the CC2/aug-cc-pVDZ level embedded in 237 water molecules within the frozen solvent approximation.The active subsystem also includes two water molecules beside acetone and was treated with the HF method during the freeze and thaw (FaT) cycles, while the remaining  326 water molecules were treated as separate subsystems using DFT with the PBE exchange correlation functional 313 and a gridsize of 3. 9,331 For the embedding potential, the PBE exchange correlation and the LC94 kinetic energy functional 332 with a gridsize of 3 were used.A comparison of the spectrum of solvated acetone with that of isolated acetone is shown in Figure 20.
Recently, another update scheme of the subsystem densities was implemented, where the densities of all subsystems are updated simultaneously at the end of a FaT cycle instead of successive updates immediately after each subsystem calculation.Thus, the updated densities of each subsystem only contribute to the embedding potentials in the next FaT cycle, offering the advantage that interim results no longer depend on the enumeration of the subsystems; this allows the subsystem calculations to be performed simultaneously to improve efficiency in parallel runs.In Table 3, the wall times for both update schemes are compared for an reduced example of the above system containing only 78 water molecules.The timings show that this kind of parallelization leads to a decrease of the wall time in the case of many small subsystems where it is inefficient to parallelize single subsystem calculations over many CPU cores.
Future developments on the implementation of the FDE method aim to increase the applicability to systems containing large subsystems by combining it with PNO-based methods 333 and to reduce further the computation time of the Coulomb contribution to the embedding potential.

FDE and Projection-Based
Embedding for Molecules and Solids.The FDE implementation described above is suitable for molecular and weakly interacting subsystems, as it employs embedding potentials based on approximate kinetic energy density functionals (KEDFs). 334cently, a DFT-based embedding scheme that treats both molecular and periodic systems on equal footing has been implemented within the riper module. 82This implementation supports both FDE and projection-based embedding 335 (PbE) via a level-shift projection operator 336 (LSPO).PbE along with FaT cycles can be used to perform exact DFT-in-DFT embedding for molecules and solids and reproduce the exact total DFT energies, even for strongly interacting subsystems.Similar to the implementation in section 3. 15.1, the embedding scheme is also coupled with correlated wave function (CW) methods and additionally with RT-TDDFT, enabling CW-in-DFT and RT-TDDFT-in-DFT calculations, respectively, on a cluster embedded in a molecular/periodic environment.However, here the CW-in-DFT calculations are performed only by adding the DFT-based embedding potential as a static term to the HF core potential of the active subsystem and obtaining the converged HF reference orbitals for post-SCF calculations.RT-TDDFT-in-DFT, on the other hand, does support updating a portion of the embedding potential during the time evolution in select cases.
As an illustrative application of the work, Figure 21 and Table 4 show that the solvatochromic shift in the first excitation energy, calculated using the CC2-in-DFT method (molecule-in-periodic) for an acetone molecule (active) solvated with 113 water molecules (periodic environment) in a cubic box with 3D periodicity, is in remarkable agreement with the shifts calculated using standard CC2 46 for acetone +  a T(X) represents the computational wall time obtained using X CPU cores.n FaT gives the number of FaT iterations where it is necessary to reach convergence.at only a fraction of the computational cost.While beneficial, it is also important to acknowledge the limitations of the implementation as well, such as being restricted to only two closed-shell subsystems and using only LDA/GGA functionals for the embedding potential.

Periodic Hartree−Fock Exchange.
In ref 328, we have presented a robust implementation of the periodic Hartree−Fock exchange in TURBOMOLE's riper module. 56,57Without precautions, exchange matrix elements may be divergent, arising from an artificial periodicity of the density matrix.This periodicity of the density matrix is introduced in practical calculations by the discretization of wavevectors.The finite k-mesh determines in turn the size of the Born−von Kaŕmań supercell.We have demonstrated 328 that a minimum image convention 338 removes the divergence for discrete kmeshes.While calculations with periodic HF exchange may be unstable for small supercells, stable SCF calculations and convergence of total energies are typically achieved for sufficiently large sizes of the supercells.The size of the supercell or k-mesh that is required for a reliable energy depends on the locality of the density matrix and hence both the electronic structure of the studied material and the chosen basis set.For selected insulators and semiconductors, we have demonstrated that HF total energies converge exponentially with the number of k-points, 328 see Figure 22.
Through our implementation of periodic exchange, conventional Hartree−Fock calculations can be carried out with TURBOMOLE for periodic systems of any dimension.In addition, DFT calculations with hybrid functionals can now be performed routinely for semiconductors and insulators, and we showed successful applications of PBE0 313,344 and HSE06 345 hybrid and range-separated hybrid functionals. 328As the next important step, analytical gradients shall be added for structure optimization.Furthermore, the Hartree−Fock exchange infrastructure that is available may be used in the development of new electronic structure methods for periodic systems that require exchange.

Nuclear Electronic Orbital Method.
Proton-coupled electron transfer (PCET) reactions are an important class of reactions that cannot be adequately described within the Born−Oppenheimer approximation. 346A remedy to this problem is the nuclear electronic orbital (NEO) methods, which treat not only electrons but also the protons of selected hydrogen atoms quantum mechanically. 347,348This is of particular importance for reactions that include proton transfer, such as, for example, acid−base reactions.The corresponding effects become especially important when the proton transfer is coupled to the electronic structure, such as in photoacids and photobases.In an initial proof-of-principle implementation, nuclear electronic orbitals were made available for the Hartree−Fock method (NEO-HF) and second-order Møller−Plesset perturbation theory (NEO-MP2) in a development version of TURBOMOLE.Furthermore, for the NEO-HF method, analytical gradients have since been implemented to allow for structure and basis set optimization.
Figure 23a shows the nuclear orbitals of the trans-Zundel isomer H 9 O 4 + . 349The nuclear orbital energies, as calculated by the NEO-HF method, give an estimate of the binding energies of the protons.While the four outer ones have energies from −450 to −436 kcal/mol, the central one is the least stable with an energy of −417 kcal/mol.If the outer two water molecules are removed, the energies change to between −400.5 and −399.6 kcal/mol for the four outer protons and −371 kcal/ mol for the central one.This hints at the ionic cluster being stabilized by the outer water molecules.
The protons of the neutral glycine molecule in Figure 23b have orbital energies between −633 and −546 kcal/mol, with the least stable one being the one at the carboxyl group.The zwitterionic glycine structure in Figure 23c has been optimized as a positively charged system with six classical protons.A NEO-HF calculation with five protons resulted in proton orbital energies between −663 and −558 kcal/mol.This demonstrates that protonation and deprotonation at certain sites can be elegantly investigated by the ab initio occupation of nuclear orbitals at the respective sites instead of placing classical nuclei according to intuition.
Implementing these methods in TURBOMOLE allows for the use of highly efficient schemes that already exist for purely electronic methods.The NEO methods can also use various existing programs to analyze the results.is one of the most used approaches for wavefunction-based correlation energies, as well as being used in the double-hybrid (DH) DFT. 353,354Nevertheless, the MP2 method shows several limitations, e.g., it overestimates the correlation energy in large systems 355 and diverges for systems with a vanishing gap. 356−358 Another more recent path in this direction is the Møller− Plesset or Hartree−Fock adiabatic connection (HFAC) method. 359In the HFAC approach, the correlation energy E c is given as a nonlinear function of E MP2 , E x HF (the HF exchange energy), and two semilocal functionals of the HF density (W c = W c [ρ HF ] and W c ′ = W c ′[ρ HF ]).The latter are derived from the strong-correlation regime, 360,361 and we have The HFAC method de facto includes an infinite-order resummation of the MP correlation series thanks to the interpolation with the strong-correlation limit, as in the more conventional AC based on DFT. 362The nonlinear function F can be approximated by modeling the HFAC curve at various coupling strengths 359 using known exact asymptotic conditions. 361,363Consequently, F satisfies two important limits: where G is a nonlinear function whose form depends on the choice of F. For well-behaving approximations of F, the condition in eq 11b yields a finite energy whenever E MP2 → − ∞, thus removing one main limitation of MP2 and DH functionals for systems with a vanishing gap.The condition (eq 11a) is an exact condition, 359 which is violated in all the regularized MP2 methods. 356,357Thus, the HFAC method allows us to overcome the main drawbacks of the MP2 approach within a well-defined theoretical framework at the small extra cost of a post-HF semilocal DFT calculation.Some working approximations of F have been proposed (e.g., ISI, 362,364 RevISI, 365 and MPACF1 363 ), and they have been implemented in TURBOMOLE together with the currently available DFT approximations for W c and W c ′. 360,366 Note that eq 10 is not size-consistent for systems composed of different species of fragments (as F is a nonlinear function).However, a size-consistent correction (SCC) 367 can be readily computed with TURBOMOLE at no additional costs, allowing the calculation of dissociation curves.
Two examples of applications where MP2 and DH functionals fail whereas the HFAC implementation can be readily used are displayed in Figure 24.
Figure 24a shows the dissociation curve for the coronene dimer, a prototype for a wide class of problems that are hardly tractable with high-level correlated wave function methods but are poorly described by MP2 because of the lack of high-order correlation contributions. 355In constrast, HFACM methods are very close to the available reference data at the equilibrium geometry.The MPACF1 functional has been tuned on dispersion complexes and on average exhibits an error 2−5× smaller than MP2. 3639][130][131]357,364 All methods work well close to the equilibrium geometry. However for a larger separation where the energy-gap closes, MP2 and DH functionals rapidly diverge.HFAC methods, however, remain well behaved, yielding a finite interaction energy, see eq 11b.The exact result is not reproduced, as the available HFAC functionals are approximated and do not take into account the recent theory developments.361 Further development and testing are thus required, and the HFAC implementation in TURBOMOLE represents an efficient platform to this end.

4.3.
Approximate TDDFT Approaches.Time-dependent density functional theory is still the most-used approach for the calculation of excitation energies of molecular and extended systems, thanks to its favorable accuracy/computational-cost ratio.Nevertheless, the computational cost of first-principles TDDFT calculations limits routine calculations to systems with a few hundreds of atoms, depending on the choice of the exchange-correlation (XC) functional.−377 Another efficient, though approximated, path is to perform a semiempirical tight-binding linear-response (TBLR) approximation, 378−380 using first-principles KS orbitals and eigenvalues.TBLR methods speed up the TDDFT calculation by about two orders of magnitude.TBLR is accompanied by a loss in accuracy of about 0.1−0.2eV, which is comparable to the overall TDDFT accuracy. 381,382−380 Instead, in the TDDFT-as method, the latter approximation is not employed and, moreover, the calculation of the semilocal XC kernel contribution on the grid is not required, 383 as it can be modeled/approximated by the same exponent α of the s-type Gaussian auxiliary basis  368,369 and the equilibrium distance is R 0 = 3.458 Å.All calculations were performed with the hfacm script.The coronene dimer results are based on a complete basis set extrapolation from cc-pVQZ results.H 2 calculations employed the aug-cc-pVQZ basis set.The strong-correlation functional used in the HFAC results is the hPC model. 366unction.However, the exponent α needs to be optimized for each atom type separately.
Here, we shortly report on results for two computationally expensive cases using the PBE 313 functional.First, a 120 atom silver nanoparticle, Ag 120 , with T d symmetry calculated using the def-SVP 41 basis set.Second, a fullerene with C 1 symmetry calculated using the def2-TZVP 270 basis set.Both systems contain only a single type of atom.We optimized α by minimizing the root-mean square (RMS) averaged excitation energy error E avg , avoiding state flipping, 383 and considering 400 t 2 (600 a) excited states for the silver nanoparticle (fullerene).The optimization yields α Ag = 0.036 with E avg = 5 meV and α C = 0.18 with E avg = 12 meV for Ag 120 and C 60 , respectively, as reported in Table 5.The amazing accuracy is also retained when the maximum error on all excitation energies E max , including many optically dark states in C 60 , is considered.This is demonstrated in Table 5 and Figure 25, where we report the absorption spectra of the two systems considered.The TDDFT-as absorption spectra in Figure 25, reported on a log scale to also highlight states with low oscillator strengths, can be hardly distinguished from the reference TDDFT results in a wide energy range for both systems investigated.Compared to TBLR approaches, 383 the accuracy of TDDFT-as is therefore increased by an order of magnitude, showing that it is a highly competitive approach.
The computational speed up of the TDDFT-as method is shown in Table 5.The XC part is completely neglected, and the RI-Coulomb part is reduced by a factor of 50−70 due to the strongly reduced dimension of the auxiliary basis set.
More recently, TDDFT-as was extended to hybrid functionals and defined across the periodic table, forming the TDDFT-ris model. 384Using the PBE0 hybrid functional 344 on a test set comprising small-to-large organic molecules, TDDFT-ris has an average error of E avg = 60 meV 384 compared to an average error of E avg = 240 meV for sTDDFT. 378,379Thus, the TDDFT-as and TDDFT-ris methods, both fully available in the next release, are efficient and accurate approximations of standard TDDFT, providing a significantly less empirical alternative to TBLR approaches.Thanks to the flexible and efficient implementation, accurate simulations of the absorption spectra of large nanoparticles and organic molecules are available at a fraction of the computational cost of standard TDDFT.Examples for such codes are the multilayered periodic general Mie method code (mpGMM), JCMsuite, and COMSOL Multiphysics, which describe light−matter interactions of complex optical devices made form novel materials.In the foreseeable future, the current developments and implementations will be converted into a fully fledged workflow for optical material simulations.This will make the ab initio-based Tmatrix approach available to a broader scientific community interested in a bottom-up approach of simulating complex artificial molecular materials and photonic devices.Future work will be dedicated to nonlinear optical properties, as this topic is currently seeing increasing interest in the scientific community.The change in the molecular dipole moment Δμ i (polarization) upon exposure to the oscillating external electric field E i at the excitation frequencies is often expressed as a power series of the incident field E.
In eq 12, α ij denotes the polarizability, β ijk denotes the first hyperpolarizability, γ ijkl denotes the second hyperpolarizability and so on.Currently, optical multiscale studies are limited to linear response, taking into account α ij . 254,261To take into account nonlinear effects, i.e., β ijk and/or γ ijkl , the additionally arising quadratic (and/or cubic) response terms of eq 12 need to be taken into account.While TURBOMOLE already allows calculation of the first hyperpolarizabilities β ijk for real frequencies, ongoing work is dedicated to expanding this toward general complex frequencies.Ultimately, this will allow not only studies of the nonlinear light−matter interactions on In the TDDFT(-as) calculations, the 1s core orbital and the 4s4p orbitals were kept frozen in C 60 and Ag 120 , respectively.the individual molecular level but also the construction of "hyper-T-matrices".The latter can be used to investigate for example second-harmonic generation (SHG) efficiencies, macroscopic second-order susceptibilities, and two-photon absorption of photonic devices made from molecular materials.

Relativistic Effects and Magnetic Properties of Periodic Systems.
As is evident from the many sections in this Review focusing on molecules, TURBOMOLE was initially developed to study finite molecular systems.However, the code infrastructure was extended to support calculations with periodic boundary conditions almost 15 years ago, 55 and developments for molecular systems can be transferred to the periodic code.Recently, a two-component DFT framework was implemented for ground-state calculations, supporting energies and various plotting options for bands and the electron density on a grid. 385The reason to use such a framework for periodic systems is twofold.First, it allows for the inclusion of spin−orbit coupling in a variational ansatz for relativistic effects.Second, the 2c formalism is necessary to study magnetic properties and arbitrary spin alignments, i.e., ferromagnetic, antiferromagnetic, and noncollinear spin textures.
A pilot application to the band gaps of AgI is shown in Table 6.Here, relativistic effects substantially affect the band energies.That is, the nonrelativistic approach shows a large deviation toward the four-component (4c) ansatz directly based on the Dirac equation.In contrast, the scalar 1c approach yields substantial improvement.The spin−orbit 2c approach, employing ECPs, further improves upon these results and is in good agreement with the four-component reference.This shows that the 2c framework serves as an excellent approximation with drastically reduced computational demands.Here, a Kramers-restricted (time-reversal symmetric) approach is available for closed-shell systems, whereas a Kramers-unrestricted approach (breakdown of timereversal symmetry) is used for open-shell systems.For the latter, a noncollinear treatment of the electron spin is applied.
Current endeavors cover the extension of the 2c formalism to energy gradients 56 and the stress tensor, 59 as well as the inclusion of the current density, 151 see section 3.2.This will allow us to perform structure optimizations and sophisticated studies with τ-dependent functionals, as the current density is of crucial importance for materials such as Weyl semimetals 392 and magnetic Hopfions. 393Another strong point of the 1c implementation for periodic systems is the availability of HF exchange, which can be applied at a reasonable cost as localized basis functions are used. 328This allows for the use of generally all available global and range-separated hybrid functionals in a stable and convergent framework.Extension of this feature to the 2c framework will allow for a more precise description of band gaps and other properties.For example, for magnetic properties it was shown that the amount of HF exchange incorporated is crucial. 119Therefore, a robust implementation of 2c HF exchange in the periodic framework will be useful in determining related quantities also for materials in the solid state and nanostructures.

OUTLOOK
The quality of a code strongly correlates with the health and functioning of its developers' community. 394TURBOMOLE developers are organized in small units pursuing their own scientific agendas, as illustrated by this Review.While this is a typical and, to a degree, necessary modus operandi for large scientific coding projects, the need to secure original authorship and demonstrate scientific independence often conflicts with sharing plans and code, taking collective responsibility, and avoiding "technical debt".As a result, TURBOMOLE has historically not been particularly easy to use, contribute to, or interface with other codes.TURBO-MOLE GmbH was founded precisely to address these issues and has provided a framework to advance common goals and improve code quality.Nevertheless, incentives to collaborate and adopt sustainable coding practices remain few and far between.The future of the TURBOMOLE project will vitally depend on whether the conditions set by the environment, i.e., academic institutions, funding agencies, reviewers, the developers, and not least the users, foster a thriving and collaborative community, which incentivizes continued investment in the code base.He co-organized the collaboration for this paper and wrote and revised parts of the manuscript.

Notes
The

Figure 1 .
Figure 1.Spin-restricted dissociation curve of P 2 with the original (LH20t) and sc-corrected t-LMF (scLH22t) along the bond axis for different distances on the dissociation curve (right).See ref 135 for related graphics for other molecules.

Figure 2 .
Figure 2. Effect of restoring gauge invariance (GINV) for various τdependent functionals on the average deviation of vertical TDDFT excitation energies from the experimental absorption maxima for five Ni(II) complexes. 146Reprinted from ref 145 with permission.Copyright 2022 AIP Publishing.

Figure 4 .
Figure 4. MCD spectrum of ZnDiNTAP.(a) Molecular structure.(b) Experimental spectrum.Reprinted with permission from ref 180.Copyright 2007 RSC Publishing.(c) Spectra as calculated in finite magnetic fields of 5 and 1000 T. Adapted with permission from ref 176.Copyright 2022 the Authors.

3 . 6 .
Paramagnetic NMR Shieldings and Shifts.NMR spectroscopy is also an important technique for the characterization of open-shell chemical compounds.NMR shielding tensors and chemical shifts describe the positions of the peaks in NMR spectra.In the closed-shell case, only the temperatureindependent orbital contribution is relevant to the calculation

Figure 8 .
Figure 8.Total wall times for SCF and NMR calculations with respect to the number of glucose units (6-31G* basis).Wall times were measured on a single thread of an Intel Xeon Gold 6212U CPU (2.40 GHz).Reprinted with permission from ref 63.Copyright 2021 the Authors.

Figure 9 .
Figure 9. Aromatic ring current of [Th@Bi 12 ] 4− .The color scheme (red to blue) indicates strong to weak currents.Data from ref 205

Figure 11 .
Figure 11.Scheme for T-matrix-based multiscale modeling of light−matter interactions using damped response polarizabilities as outlined in ref 254.Reprinted with permission from ref 254 under a CC BY-NC-ND license.Copyright 2022 the Authors.

Figure 12 .
Figure 12.ECD (upper panel) and absorption spectra (lower panel) of (4S p )-4,7-dicyano-12,13,15,16-tetramethyl[2.2]-paracyclophane computed at the CC2/aug-cc-pVDZ level.The MP2-optimized structure and the experimental CD spectrum have been taken from ref 263.Results from damped response theory are plotted as blue squares; the blue line is a cubic spline fit to these computed points.Stick spectra are the results from state-wise calculations, and the red and orange lines are obtained by convoluting these computed stick spectra with a Lorentzian broadening with a half-width-at-halfmaximum of ≈1000 cm −1 including, respectively, the lowest 14 and 59 states.

Figure 13 .
Figure 13.Computed (left) and experimental (right) TA spectra of pyrene in the gas phase.Band origins of S 1 and S 2 were shifted by −0.42 and −0.02 eV, respectively.Experimental spectra reprinted from ref 278.Copyright 1999 Hindawi Publishing Corporation.Distributed under a CC-BY license.
Methods.Since 2010, TURBOMOLE has supported accurate wave function methods for computing the ground-state energies of molecular configurations, and up to CCSD(T) and BCCD(T) are available, together with their explicitly correlated counterparts CCSD(T)(F12*) and

Figure 15 .
Figure 15.EAs of a water tetramer cluster for different dimer separations (R) computed using GKS-spRPA, ΔHF, ΔCCSD(T), and EOM-CCSD(T)a* methods.The inset shows the arrangement of the tetramer with a 70% isosurface of the LUMO at R = 4.047 Å. aug-cc-pVDZ basis sets were used for O and H atoms, and a 7s7p set of basis functions located at the center of the cluster was used.The shaded (unshaded) area of the plot corresponds to NVCB (NVEB) anionic states.Reprinted with permission from ref 309.Copyright 2021 American Chemical Society.

Figure 16 .
Figure16.XE spectra of S 8 from experiments311,312 and the SR-GKS-spRPA method.310The computed spectrum (in blue) was obtained by Gaussian broadening of the vertical transitions (red lines) using a width parameter of 1 eV.The vertical dashed lines denote the two intense peak positions for SR-GKS-spRPA.Reprinted and adapted with permission from ref310.Copyright 2022 American Chemical Society.

Figure 17 .
Figure 17.Plots of the five highest occupied orbitals at isovalues of ±0.02 shown for (a) PBE and (b) GKS-spRPA (with the PBE potential).The structure of quinacridone was optimized with PBE-D3 313−315 and cc-pVTZ basis sets.316Reproduced from ref 317 with the permission of AIP Publishing.

Figure 19 .
Figure 19.(a) Experimental linear absorption spectra and (b) HHG spectra of 100 nm thick TPP and ZnTPP samples.(c) Calculated linear absorption spectra and (d) HHG spectra of TPP and ZnTPP molecules.326

Figure 22 .
Figure 22.Self-consistent total energy differences |E(N k 1/3 ) − E(∞)| per primitive cell for the PBE and HF methods.LiH is calculated in the rocksalt structure with a lattice constant of 4.084 Å 339 and is described with the basis set from ref 340.For Si, we use the diamond structure with a lattice constant of 5.430 Å 341 and the pob-TZVP 342 and def-SVP 343 basis sets.E(∞) is approximated by the energy obtained with a 31 × 31 × 31 k-mesh, since |E(25) − E(31)| < 4 × 10 −10 a.u.for LiH and Si.Reprinted with permission from ref 328.Copyright 2018 American Chemical Society.

Figure 23 . 4 4. 2 .
Figure 23.Nuclear orbitals as isosurfaces of the total densities with a cutoff of 0.0001 calculated with NEO-HF in the dscf program using the def2-TZVP basis for electrons and DZSPDN 348 for protons.White, classical proton; pink, center for basis set of a quantum proton; red, oxygen; and blue, nitrogen.(a) trans-Zundel isomer of H 9 O 4

Figure 24 .
Figure24.Dissociation energy curves for (a) a coronene dimer and (b) H 2 , as computed with different methods.For the coronene dimer, reference data are taken from the literature,368,369 and the equilibrium distance is R 0 = 3.458 Å.All calculations were performed with the hfacm script.The coronene dimer results are based on a complete basis set extrapolation from cc-pVQZ results.H 2 calculations employed the aug-cc-pVQZ basis set.The strong-correlation functional used in the HFAC results is the hPC model.366

4 . 4 .
Multiscale Modeling Extensions for the Nonlinear Optical Response of Molecular Materials.Further ongoing work concerning multiscale modeling of optical molecular materials is envisaged in the direction of integrating the automatic construction of T-matrices, as used in the multiscale approach in section 3.9.1.These T-matrices or effective material parameters derived thereof are often used in scattering codes and finite-element method-based programs.

Figure 25 .
Figure 25.TDDFT and TDDFT-as absorption spectra (in log scale) for (a) Ag 120 and (b) fullerene determined using the PBE XC functional and a Lorentzian broadening of 20 meV.

7 .
He wrote and revised parts of the manuscript.R. Grotjahn performed investigations on the importance of the current-density correction in TDDFT and implemented the excited-state gradients and quadratic response properties for current-dependent MGGAs.Furthermore, R.Grotjahn contributed to the development of the ωLH22t functional.M. Kaupp supervised the work of S. Furst and A. Wodynśki as well as part of the work of R. Grotjahn and supported work on local hybrid functionals, including range-separated local hybrids, strong-correlation corrections, and aspects pertaining to current dependence.M. Kehry was involved in the design and programming steps of the underlying (2c) damped response modules for GW-BSE and 2c TDDFT.He wrote and revised parts of the manuscript.J. H. Andersen and D. A. Fedotov contributed to the development of the dampedresponse RI-CC2 functionalities.J. H. Andersen and C. Haẗtig carried out the illustrative damped RI-CC2 calculation here reported.C. Haẗtig and S. Coriani conceived, conceptualized, and supervised the damped-response RI-CC2 project and wrote the text in section 3.9.2.M. Krstićand B. Zerulla devised the multiscale modeling approach and designed and programmed the interface to external Maxwell-solver based programs.They wrote and revised parts of the manuscript.F. Mack contributed to the implementation of the NMR coupling constants (nonrelativistic and 2c) at the DFT and GW-BSE levels.S. Majumdar contributed to the research discussed in section 3.10 and was involved in conceptualization, execution, and writing.G. S. Phun contributed to writing section 3.10.A. Rajabi and D. Rappoport were involved in both the conceptualization and execution of the research discussed in section 3.8.S. M. Parker developed and supervised the TDDFT-ris method.F. Pauly and his group contributed to method development on periodic Hartree−Fock exchange and relativistic effects in periodic systems and the corresponding sections of this review.He wrote and revised parts of the manuscript.A. Pausch supported Y. J. Franzke and C. Holzer in the implementation of the two-component CDFT framework.A. Pausch and C. Holzer developed the methods for finite magnetic fields.He wrote and revised parts of the manuscript.E. Perlt performed implementations of the NEO methods and contributed to writing of section 4.1.T. Schrader performed NEO calculations and wrote section 4.1.M. Sharma, under the supervision of M. Sierka, contributed to the RT-TDDFT code and extended it to perform HHG simulations and also wrote the corresponding section 3.14.Additionally, M. Sharma implemented the molecular and periodic DFT-based embedding coupled with RT-TDDFT and wave function methods and wrote section 3.15.2.B. Samal and B. D. Nguyen contributed to writing section 3.12.V. K. Voora developed the AC version of the GKS-spRPA method and contributed to the writing of section 3.12.A. Wodynśki carried out the implementation, optimization, and evaluation of strong-correlation-corrected local hybrid functionals.J. M. Yu contributed to the conception and execution of the research detailed in section 3.4.R. Treß developed the frozen density embedding implementation and carried out the calculation reported in section 3.15.1, and C. Haẗtig supervised the project.F. Furche conceived, oversaw, and participated in the development of the material in sections 3.8 and 3.10 and parts of the material in sections 3.2, 3.12, 3.4, and 4.1.He also initiated and organized the collaboration for this paper and wrote and revised parts of the manuscript.M. Sierka supervised the work of M. Sharma described in sections 3.14 and 3.15.2,as well as conceived, oversaw, and participated in the development of the RT-TDDFT code (section 3.14).He coorganized the collaboration for this paper and wrote and revised parts of the manuscript.F. Weigend supervised the work of F. Bruder, S. Gillhuber, and F. Mack and supported the EPR, pNMR shift, NMR coupling, and ring current studies.

Table 1 .
and Available Parallelizations of Various Modules Shown in Version 7.7 a

Table 3 .
Wall Times for an Acetone Molecule in a Water Cluster for the Two Different Update Schemes a

Table 4 .
Solvatochromic Shift ΔE of Acetone (Ac) S 1 Excitation in Water Calculated via CC2 on Ac + (H 2 O) n a T denotes the wall time of the CC2 solvatochromic shift calculations.Journal

Table 5 .
Optimized Exponent of the s-Type Auxiliary Basis Function, RMS-Averaged (E avg ) and Maximum (E max ) Errors on Excitation Energies, and Computational Cost (on a Single-Core Intel Xeon Gold 6132) for the First Davidson Step for the RI-J and XC Part for the Two Systems Considered a

Corresponding AuthorsTable 6 .
386rgy Band Gaps (in eV) of Three-Dimensional AgI (Lattice Constant 6.169 Å, Rocksalt Structure386) Obtained for Various k-Points with the PBE Functional a for this paper and wrote and revised parts of the manuscript.C.Holzer contributed to the conception and implementation of the current density framework for onecomponent and two-component ansaẗze (finite magnetic field, SCF, NMR, EPR, and TDDFT).Additionally, C. Holzer implemented the initial graphics processing (GPU) unit support in the modules listed in Table1.C. Holzer further contributed to the conception and implementation of the (2c) GW-BSE method in general and the related solvers.Furthermore, the 2c damped response and static 2c CPKS framework were designed and implemented together with Y. J. Franzke and M.Kehry.Additionally, the implementation and interfacing for the T-matrix based multiscale modeling of optical properties was codesigned by C.Holzer.He coorganized the collaboration for this paper and wrote and revised parts of the manuscript.Y. J. Franzke and C. Holzer have rewritten the local hybrid code from 2020−2022, i.e., increased efficiency, the just in time (JIT) framework, corrected memory handling, and so on.They implemented the general two-component version (Kramers-unrestricted LMFs and calibration function, XC kernel, and CDFT).T. Begusǐćand E. Tapavicza implemented the TGA method for vibronic spectra.F. Bruder added spin−orbit perturbation theory for EPR and pNMR using the X2C routines developed above.He wrote and revised parts of the manuscript.F. Della Sala implemented the HFACM approach and developed the TDDFT-as/TDDFT-ris methods.E. Fabiano developed and implemented the HFACM approach.S. Furst carried out the implementation, optimization, and evaluation of rangeseparated local hybrid functionals, including the new ωLH22t functional.S. Gillhuber implemented the nonrelativistic EPR and pNMR extension.Y. J. Franzke added the scalar X2C part.S. Gillhuber contributed to the research of section 3. ration authors declare the following competing financial interest(s): Principal Investigator Filipp Furche has an equity interest in TURBOMOLE GmbH.The terms of this arrangement have been reviewed and approved by the University of California, Irvine, in accordance with its conflict of interest policies.Christof Hattig and David P. Tew have an equity interest in TURBOMOLE GmbH.Marek Sierka and Florian Weigend have an equity interest in TURBOMOLE GmbH and serve as its chief executive officers.■ACKNOWLEDGMENTS All past and present developers' contributions to the TURBOMOLE project are gratefully acknowledged.A list of TURBOMOLE contributors is available on the TURBO-MOLE website. 91.J. Franzke was supported by fellowships from Fonds der Chemischen Industrie (FCI, German Chemical Industry Fonds), Deutscher Akademischer Austauschdienst (DAAD, German Academic Exchange Service), and TURBOMOLE GmbH.C. Holzer and M. Krstićgratefully acknowledge funding by Volkswagen Stiftung.T. Begusǐćand E. Tapavicza acknowledge scientific support from J. Vanícěk in the development and application of the TGA method.F.Della Sala acknowledges the financial support from ICSC−Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union− NextGenerationEU−PNRR.D. A. Fedotov and S. Coriani acknowledge support from the European Unions Horizon 2020 research and innovation program under the Marie Skłodowska-Curie European Training Network COSINE (grant agreement no.765739).J. H. Andersen and S. Coriani acknowledge financial support from the Independent Research Fund Denmark-DFF-FNU RP2 (grant no.7014-00258B).S. Gillhuber is supported by a fellowship from Fonds der Chemischen Industrie (FCI no.110160).R. Grotjahn acknowledges support via a Walter-Benjamin postdoctoral fellowship funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), no.501114520.C.Hättig acknowledges support by the Deutsche Forschungsgemeneinschaft (DFG) via project Ha 2588/10-1.The Kaupp group has been supported by the Deutsche Forschungsgemeinschaft (DFG) via projects KA1187/14-1 and KA1187/14-2.M. Kehry acknowledges financial support by the DFG through the Transregional CRC 88 "Cooperative Effects in Homo-and Heterometallic Complexes" (project C1).F.Mack acknowledges support from TURBOMOLE GmbH and from the DFG through the CRC 1176 (Project Q5).The material in section 3.8 and parts of the material in section 3.12 is based upon work supported by the US National Science Foundation under CHE-2102568.The material in section 3.2 and in section 3.10 is based upon work supported by the US Department of Energy, Office of Basic Energy Sciences, under award number DE-SC0018352.A. Pausch was supported by a fellowship from Fonds der Chemischen Industrie and Studienstiftung des deutschen Volkes (German Academic Scholarship Foundation).E. Perlt and T. Schrader acknowledge support from the Carl Zeiss Foundation within the CZS Breakthroughs Program.M. Sierka and M. Sharma gratefully acknowledge financial support from Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the CRC 1375 NOA, project A4, and from the Carl Zeiss Foundation within the CZS Breakthroughs Program.B. Samal and V. K. Voora were supported by the Department of Atomic Energy, Government of India, under project no.RTI2001.F. Weigend acknowledges support from the DFG through the Collaborative Research Centre (CRC) 1573 (Project Q).J.M. Yu acknowledges support from the US National Science Foundation under Grant DGE-1839285.B. Zerulla acknowledges support by the KIT through the "Virtual Materials Design" (VIRTMAT) project.