TURBOMOLE: Today and Tomorrow

1TURBOMOLE 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.

2The 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− and time-dependent HF (TDHF) response properties. A major milestone was the relatively early adoption of density functional theory (DFT) using newly designed quadrature algorithms and that of time-dependent DFT (TDDFT) shortly afterward. With extensions to meta-GGA, RPA, and other fifth-rung methods,− as well as current density and local hybrid functionals,, critical performance improvements,− and a plethora of available analytical derivatives of ground- and excited-state energies,− TURBOMOLE has become one of the leading all-purpose molecular (TD)DFT codes. The development and implementation of the RI approximation for Coulomb (RI-J,) and exchange contributions, (RI-K,), as well as its generalization to post-HF theories such as MP2 and CC2− and its multipole-accelerated version for extended systems (MARI-J), were other critical innovations. The RI methods are still cornerstones in many modern implementations and outstanding features of the program suite.,,− More recent additions include explicitly correlated wave function methods up to CCSD(T) and BCCD(T),, efficient pair natural orbital (PNO) approaches,− solvation models and embedding,− two-component relativistic methods,− GW-BSE type methods,, real-time (RT) TDDFT, and nonadiabatic molecular dynamics.

3To 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.

4Here 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 and and the TURBOMOLE Web site.

5The 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.py for reaction path optimization, vibrational Raman spectra, vibrational circular dichroism spectra (VCD), and genetic algorithms, respectively. Moreover, the graphical user interface TmoleX is of great help for running calculations and visualizing results.

6Almost all time-consuming parts are parallelized for multicore systems or clusters using OpenMP for shared-memory parallelization (SMP),,− and the message-passing interface (MPI) for parallelization across multiple nodes,− as outlined in Table 1. The older Fork-SMP is available as a fallback for some modules. Starting with the latest release (V7.7), support for graphics processing units (GPUs) has become available.

7The 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.

8Local hybrid functionals (LHs), 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, with functionalities that exceed by far those available in any other quantum-chemistry package that contains LHs. Beyond ground-state SCF and nuclear gradients, this now includes linear-response TDDFT energies, and excited-state gradients, frequency-dependent and frequency-independent polarizabilities,, NMR chemical shifts and spin–spin coupling constants,, EPR hyperfine couplings and g-tensors,, the related NMR shieldings of paramagnetic systems, magnetizabilities, and quasiparticle states based on GW. From the variety of possible applications and evaluations, many of which have been touched upon in the 2020 overview of TURBOMOLE and in a 2019 comprehensive review of local hybrids, we highlight in particular the outstanding performance of LHs for mixed-valence systems− and for phosphorescence spectra.− We also point to further recent LH publications and to reviews by other authors.,,,−

9Here, 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”, between reduced self-interaction errors and delocalization errors or “fractional-charge errors” (FCEs) on one side and minimizing static-correlation errors or “fractional-spin errors” (FSEs) on the other side. The 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. On 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 and a modified approach by Kong and Proynov (KP16/B13), which have both been implemented self-consistently into a local developer’s version of TURBOMOLE.

10Circumventing 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 qAC factor. Most recently, the more advanced scLH22t functional has been constructed. It is based on the more recent and overall more accurate LH20t functional. Using a damping factor for smaller NDC contributions, an almost complete decoupling between the underlying LH20t and the added qAC 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. Notably, the qAC factor forms part of a new LMF (Figure 1).

11The 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., That is, instead of full-range semilocal exchange, short-range exchange is mixed in. This corresponds to the use of a long-range EXX admixture governed by a suitable range-separation parameter ω. Based on earlier work for locally range-separated hybrids, 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). Using this implementation, the ωLH22t RSLH has been constructed and optimized. It retains most of the good performance of LH20t for main-group and transition-metal 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 self-interaction errors. Most 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.

12The 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 that this quantity is not invariant under a gauge transformation in the external potential. Substitution of τ by its generalization−, where jP 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 jP. While the original implementation of these terms in TURBOMOLE dates back to 2012 (release V6.4), we highlight four recent important updates.

13First, an incorrect prefactor 2 has been removed from the original implementation with V7.6., 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.

14Second, a recent investigation 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.

15A broader analysis reveals that the importance of imposing gauge invariance can be linked to the derivative of the exchange energy integrand with respect to τ. 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. These 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. Recently, 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. These developments will be available in a future TURBOMOLE release.

16Third, the inclusion of jP is required for gauge invariance of magnetic properties and implemented for magnetizabilities, NMR coupling constants, and NMR shifts of closed-shell, and open-shell systems, as well as EPR hyperfine coupling constants and g-tensors. Recent findings indeed hint at the inclusion of the current density response also being crucial for NMR and EPR properties.,, Especially for open-shell systems, neglecting the current density leads to large deviations, as shown in Figure 3.

17Finally, it was shown that the inclusion of the paramagnetic current is also crucial in relativistic two-component calculations. Contrary 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. Accordingly, 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 jP-dependent terms, which have a profound impact on many properties., Current 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). Given the profound impact of jP, we therefore strongly recommend the exclusive use of the current-dependent forms.

18Quantum-chemical calculations are routinely carried out for various types of molecular properties in magnetic fields, including NMR, EPR, and magnetic circular dichroism (MCD) spectroscopy (see also sections 3.4, 3.5, 3.6, and 3.9.2 of this Review).− 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.

19Other applications necessitate the use of a more general approach, particularly if the magnetic field becomes strong enough to compete with the electronic interactions within a molecule (>1000 T).− Such conditions may be found in the vicinity of interstellar objects like magnetic white dwarfs and cannot be reproduced in a laboratory.− Consequently, spectra obtained from such interstellar objects can only be interpreted using quantum-chemical calculations− In such extreme conditions, entirely new types of chemical bonding, spin-phase transitions, and other exotic phenomena have been shown to occur.,, The 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, CDFT,,GW/BSE,, RPA, and CC2. Through 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.

20Due to the efficiency of our implementation, systems containing dozens of atoms can be routinely computed in explicit magnetic fields.,,, To highlight the capabilities of our approach, we calculated the MCD spectrum of ZnDiNTAP,, 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,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. Furthermore, 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).

21Other 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. This is not captured by a perturbative approach but can be routinely investigated using our implementation.

22Over 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.

23For example, a series of recently discovered Ln-based single-molecule magnets (SMMs), [La(OAr*)3]−, [Lu(NR2)3]−, and [Lu(OAr*)3]− (OAr* = 2,6-Ad2-4-t-Bu-C6H2O, Ad = adamantyl, t-Bu= tert-butyl, R = SiMe3 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.

24Such improved predictions of EPR spectra are possible with relativistic exact two-component (X2C) theory,, including spin–orbit effects up to the noncollinear two-component (2c) DFT framework.,, For meta-generalized gradient approximations and local hybrid functionals, this also includes the paramagnetic current density in the ground state (cf. section 3.2). This rigorous formulation can be truncated to the scalar-relativistic limit or a perturbative ansatz 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,,, 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 [TbPc2]− with six unpaired electrons, Pc = bis(phthalocyaninato).,,, The self-consistent 2c methods are thus pivotal.

25With 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.

26NMR 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 scalar X2C and the spin–orbit X2C framework. All functional classes up to local hybrids are supported and include the current density for gauge invariance.,,

27For 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 large-scale calculations.

28Systems containing heavy elements such as Sn, Pb, Pd, and Pt require the inclusion of relativistic effects,,,,−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,− large-scale calculations are possible, as illustrated in Figure 6 for the Karplus relationship of Me3Sn–CH2–CHR–SnMe3, where R is different substituents (Me = CH3). 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.,

29To demonstrate the efficiency, the calculation of the Sn–P coupling constants of [({SIDipp}P)2Sn] (SIDipp = 1,3-bis(2,6-di-isopropylphenyl)-imidazolidin-2-ylidene) with 137 atoms takes about 44 min (PBE) and 55 h (PBE0) using 12 OpenMP threads of an Intel Xeon Gold 6212U CPU (2.40 GHz). Notably, using a single NVIDIA A100 GPU, the PBE0 timing can be reduced to 3.5 h.

30NMR 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 temperature-independent orbital contribution is relevant to the calculation of the shielding tensor. Both a nonrelativistic treatment,, and a scalar-relativistic treatment, are available for the orbital contribution, including the response of the current density.,, The paramagnetic NMR (pNMR) shielding tensor σItot for a nucleus I readswith S denoting the spin, μe denoting the Bohr magneton, γI denoting the gyromagnetic ratio of nucleus I, kB 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., Additionally, a temperature-dependent contribution arises, which includes the HFC tensor AI 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 1H/13C pNMR spectra of large molecules, a perturbative treatment of spin–orbit coupling is preferred over the 2c ansätze due to lower computational costs. The 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 pNMR 1H and 13C shifts of the two compounds.

31The calculation of properties depending on the density in the vicinity of the nuclei requires additional tight basis functions. Thus, the pcJ, pcS, and pcH basis sets are recommended for nonrelativistic calculations. For relativistic calculations, the x2c-s basis sets were developed.,

32The 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, and the seminumerical scheme for exchange integrals.,

33Aromatic compounds, such as benzene, show characteristic signals at 7 ppm in the 1H 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. This magnetically induced current density may be calculated indirectly with the nucleus-independent chemical shift (NICS) or directly using TURBOMOLE’s interface to the GIMIC program, which was reworked for release V7.7 and now supports open-shell calculations.− The latter approach is more flexible and also applicable to complicated multicyclic systems.−

34The 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@Bi12]4– cluster features a nonlocalizable π-orbital around the {Bi12} torus, which leads to a ring current.Figure 9 shows a streamlined representation of this ring current, whose strength amounts to about 25 nA/T.,

35Despite 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@Bi12]3– shows a strong ring current, as the same π-orbital is occupied by two electrons., Moreover, prismatic {Bi6}-based clusters such as [(CpRu)3Bi6]− show ring currents of more than 25 nA/T. Therefore, all-metal clusters help to push the frontiers of aromaticity.

36Electronic 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. Computational 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,, and a combination of damping, level shifting, and direct inversion in the iterative subspace (DIIS) extrapolation. Moreover, the stability of the ground-state reference is a concern in calculations of molecular properties, for example, electronic absorption spectra.

37DFT results helped characterize the structures and properties of nontraditional Ln2+ complexes possessing the 4fn5d1 configuration in the [Ln(C5H4SiMe3)3]− (Ln = Ce–Nd, Gd–Er),− [Ln{N(SiMe3)2}3]− (Ln = La, Gd),, and [Ln(CpiPr5)2] (Ln = Tb, Dy) series. The preference for the 4fn5d1 configuration relative to the traditional 4fn+1 occupation of Ln2+ results from the stabilization of the Ln 5dz2 orbital by the trigonal ligand environment or extremely bulky ligands. Nontraditional Ln2+ 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. DFT and RPA methods were employed to examine the strong ferromagnetic coupling between the Ln3+ centers in [(C5Me5)2Ln(μ-S)2Mo(μ-S)2Ln(C5Me5)2]− (Ln = Y, Gd, Tb, Dy) and the Mo→Ln electronic excitations in the near-infrared spectral region. Excited-state studies of [Ln(C5Me5)2(C5Me4H)] and [Ln(C5Me5)2(η3-C3H4)] complexes (Ln = Y, Dy, Lu) using TDDFT elucidated their unexpected photochemical activation, which was used to reduce dinitrogen and sulfur and to polymerize isoprene.,

38Computational studies of neutral actinide complexes [An(CpiPr5)2] (An = Th, U, Pu, Am, Bk, No, Lr) (pentaisopropylcyclopentadienyl = CpiPr5) using DFT predicted ground states with a linear ligand arrangement of S10 symmetry and significant An 6d orbital occupation for An = Th, U, Lr. The calculations were carried out with the TPSS exchange–correlation functional, Stuttgart–Cologne scalar-relativistic small-core effective core potentials (ECPs), and the corresponding basis sets., Mixed 5f/6d occupation was predicted in the corresponding Pu complex, while the An = Am, Bk, No complexes were found to have 5fn+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) 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 S10-symmetric “second-generation” uranocene [U(CpiPr5)2]. 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.

39Assessing light–matter interactions in extended or exotic systems is an active field of development in TURBOMOLE.,,,,,, While 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,, or the equivalent complex polarization propagator approach, 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. It also allows to compute nonlinear transition properties in the vicinity of additional resonances, e.g., the resonant inelastic X-ray scattering (RIXS) transition amplitudes, as well as properties at imaginary frequencies, like the C6 dispersion coefficient.,−

40A number of implementations of damped response theory have been presented in the past two decades for time-dependent density functional theory and a variety of wavefunction methods.,,,−,− 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, (available since release V7.7).

41DFT-based 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.,

42At a given complex frequency ωex = ω + iγ, where ω and γ are the real and imaginary parts of the external field, respectively, the coupled perturbed equation,is first solved. A and B are defined asεp marks orbital or quasiparticle energies, vpq,rs is a Coulomb integral, fpq,rs is the exchange-correlation kernel, if present, and Kpq,rs is a exchange integral. The precise kinds of fpq,rs and Kpq,rs depend on the chosen method.pν describes the external perturbation, e.g., an electric or magnetic field., The polarizability can then be calculated as

43Damped 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.,, This 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), simulations and predictions of the light–matter interactions of layer-structured materials ranging from a few to hundreds of nanometers are now routinely possible., To target molecular materials of arbitrary shape, the T-matrix approach was furthermore coupled with state-of-the-art homogenization techniques. Combining 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.

44For RI-CC2, the equations for the damped response, 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.Subscripts Re and Im represent real and imaginary components, respectively. The effective matrices in the equations above arewhere ASS, ASD, and ADS are, respectively, the singles–singles, singles-doubles, and doubles-singles blocks of the CC2 Jacobi matrix. Δijab = εa – εi + εb – εj – ω are 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 arewhere ξSx and ξDx are the single and doubles parts of the right-hand sides in the nonpartitioned form, respectively., The 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.

45As an illustrative application, we computed the UV–vis absorption spectrum and electronic circular dichroism spectrum of a donor–acceptor cyclophane 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 of the damped linear response function was used in the calculations (the symmetric form is also available).

46Cyclophanes 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.

47Ongoing work is concerned with extending the implementation to induced and nonlinear spectra like MCD and RIXS.

48Many 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 and nonadiabatic couplings between excited states, 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., In addition, time-resolved fluorescence (TRF) and transient absorption (TA) spectra can be simulated for comparison with experimental results.

49The simultaneous treatment of multiple electronic excited states enables the examination of Kasha’s rule. According to Kasha’s rule, singlet excited states energetically located above S1 undergo ultrafast decay to the S1 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.

50In 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, level. Azulene was found to exhibit non-Kasha behavior due to emission from the S2 state, in agreement with experiment. BCDP obeys Kasha’s rule and emits only from S1. Previous studies assigned the high-energy shoulder in the pyrene fluorescence spectrum to non-Kasha emission as a result of the reverse S1 → S2 internal conversion.,− Multistate NAMD simulations describe the non-Kasha behavior as a combination of S1 → S2 transitions and the change in the diabatic character of the S1 state from dark (Lb) to bright (La) at points of near degeneracy between the S1 and S2 states. The high-energy shoulder in the fluorescence spectrum of pyrene can be understood as originating from excited states with diabatic bright (La) character. The S2 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. The S1 lifetime of azulene was computed to be 2.2 ps in comparison to the experimental result of 1.4 ps in cyclohexane. The computed lifetime of the S1 state of BCDP was found to be 0.8 ps.

51The NAMD trajectories were also used to obtain the TA spectrum of pyrene as an ensemble average of the Gaussian-broadened 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 S2 → Sn transient absorption that decays rapidly and a S1 → Sn transition (n > 4), respectively., NAMD simulations predict S2 → S4 and S1 → S4 bands at 1500 and 800 nm, respectively. The time evolution of the S1 and S2 states is in good agreement with experiment, while the absorption maxima (λ(S2 → S4) = 1500 nm and λ(S1 → S4) = 800 nm) are red-shifted due to truncation of the electronic excitation space.

52The vibrational structure of electronic spectra gives detailed information about molecular structure and excited-state phenomena such as internal conversion and intersystem crossing., In the case of conformationally flexible molecules, molecular dynamics sampling often allows the shape of electronic spectra to be accounted for.− However, this approach fails for rigid molecules with high vibrational frequencies. In this case, special care must be taken to include the quantum nature of nuclear vibrations. Efficient theoretical simulation of molecular vibrationally resolved electronic (i.e., vibronic) spectra has been advanced by switching from the time-independent approach, which requires a tedious sum-overstates evaluation of Franck–Condon overlaps, to the time-dependent formalism (generating function approach), where the spectrum is given by the Fourier transform of an appropriate time-correlation function.,− The radless module,, 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, 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.,,,, An 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 and vibrationally promoted electronic resonance (VIPER) spectra.

53In addition, the newly implemented semiclassical thawed Gaussian approximation (TGA), 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. In 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− the Hessian of the potential energy is updated over time, in the single-Hessian version,− implemented in TURBOMOLE as part of the radless module, 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.

54The generalized Kohn–Sham semicanonical projected random phase approximation (GKS-spRPA) provides a route for obtaining one-particle energies at the RPA level of theory. These 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 to using well-known analytic continuation (AC) techniques.−

55The 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.

56Nonvalence 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.

57The 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. XE energies, ΔE, computed by taking the difference between core and valence IPs and the oscillator strengths, fosc., are evaluated within a frozen orbital approximationwhere ϕ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 scalar-relativistic (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.

58The accuracy of GKS-spRPA for ionization potentials (IPs) compared to nonselfconsistent RPA 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., A 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., However, this is often not true for semilocal density functional approximations. PBE 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 G0W0 IPs., The canonical GKS-spRPA orbital ordering qualitatively agrees with the one OT-RSH one down to HOMO(−3), see Figure 17. Furthermore, the negative GKS-spRPA HOMO energy of 7.07 eV is close to the experimental IP of 7.23 eV, whereas the negative PBE HOMO energy of 4.92 eV is significantly too small.

59Since 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 BCCD(T)(F12*)., 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.

60To evaluate energies for larger molecules, the steep 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 to 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.−

61The 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. Very 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. In contrast to implementations in other software packages, we do not use the simplified “A” approximation 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.

62To 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 CnH2n+2 and NanCln, 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 C128H256 taking 15 h to complete. For the globular systems, the observed scaling is subcubic, with Na50Cl50 taking 45 h to complete. The F12 calculations are 2–3× more costly than non-F12 calculations but provide energies close to the basis set limit without requiring basis sets with a large number of AOs per atom.

63While linear-response 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.

64RT-TDDFT for molecular systems has been implemented within the riper module−,, using the Magnus propagator and the predictor-corrector scheme for time integration. The implementation utilizes density fitting and continuous fast multipole method (DF-CFMM) techniques 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. Recently, 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).Figure 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. Overall, 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.

65Environment 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. 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.

66The 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 ground-state 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, 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.

67To 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 water molecules were treated as separate subsystems using DFT with the PBE exchange correlation functional and a gridsize of 3., For the embedding potential, the PBE exchange correlation and the LC94 kinetic energy functional 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.

68Recently, 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.

69Future 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 and to reduce further the computation time of the Coulomb contribution to the embedding potential.

70The 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). Recently, a DFT-based embedding scheme that treats both molecular and periodic systems on equal footing has been implemented within the riper module. This implementation supports both FDE and projection-based embedding (PbE) via a level-shift projection operator (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.

71As 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 for acetone + (H2O)20, acetone + (H2O)35, and acetone + (H2O)48 clusters 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.

72In ref, we have presented a robust implementation of the periodic Hartree–Fock exchange in TURBOMOLE’s riper module., Without 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 Kármán supercell. We have demonstrated that a minimum image convention removes the divergence for discrete k-meshes. 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, see Figure 22.

73Through 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, and HSE06 hybrid and range-separated hybrid functionals. As 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.

74Proton-coupled electron transfer (PCET) reactions are an important class of reactions that cannot be adequately described within the Born–Oppenheimer approximation. A 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., This 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.

75Figure 23a shows the nuclear orbitals of the trans-Zundel isomer H9O4+. The 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.

76The 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.

77Implementing 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. In future developments, NEO–DFT will be implemented for molecular and periodic systems.

78MP2 is one of the most used approaches for wavefunction-based correlation energies, as well as being used in the double-hybrid (DH) DFT., Nevertheless, the MP2 method shows several limitations, e.g., it overestimates the correlation energy in large systems and diverges for systems with a vanishing gap. For this reason, several regularized and/or scaled MP2 methods have been developed.−

79Another more recent path in this direction is the Møller–Plesset or Hartree–Fock adiabatic connection (HFAC) method. In the HFAC approach, the correlation energy Ec is given as a nonlinear function of EMP2, ExHF (the HF exchange energy), and two semilocal functionals of the HF density (Wc = Wc[ρHF] and Wc′ = Wc′[ρHF]). The latter are derived from the strong-correlation regime,, and we haveThe 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. The nonlinear function F can be approximated by modeling the HFAC curve at various coupling strengths using known exact asymptotic conditions., Consequently, 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 EMP2 → – ∞, thus removing one main limitation of MP2 and DH functionals for systems with a vanishing gap. The condition (eq 11a) is an exact condition, which is violated in all the regularized MP2 methods., Thus, 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.

80Some working approximations of F have been proposed (e.g., ISI,, RevISI, and MPACF1), and they have been implemented in TURBOMOLE together with the currently available DFT approximations for Wc and Wc′., 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) can be readily computed with TURBOMOLE at no additional costs, allowing the calculation of dissociation curves.

81Two examples of applications where MP2 and DH functionals fail whereas the HFAC implementation can be readily used are displayed in Figure 24.

82Figure 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. In 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.

83Figure 24b shows the H2 dissociation in a restricted formalism, a prototype of a strongly correlated system.−,, 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. Further development and testing are thus required, and the HFAC implementation in TURBOMOLE represents an efficient platform to this end.

84Time-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. Different methods and algorithms have been developed in order to increase the TDDFT efficiency,− most prominently the RI-J technique.− Another efficient, though approximated, path is to perform a semiempirical tight-binding linear-response (TBLR) approximation,− 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.2 eV, which is comparable to the overall TDDFT accuracy., More recently, it has been shown that, for semilocal XC functionals, the TBLR approaches can be considered as an approximation of the RI-TDDFT scheme with only one s-type Gaussian basis function per atom in the RI auxiliary basis set (TDDFT-as) and with the three index RI integrals replaced with a Löwdin approximation.− 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, as it can be modeled/approximated by the same exponent α of the s-type Gaussian auxiliary basis function. However, the exponent α needs to be optimized for each atom type separately.

85Here, we shortly report on results for two computationally expensive cases using the PBE functional. First, a 120 atom silver nanoparticle, Ag120, with Td symmetry calculated using the def-SVP basis set. Second, a fullerene with C1 symmetry calculated using the def2-TZVP basis set. Both systems contain only a single type of atom. We optimized α by minimizing the root-mean square (RMS) averaged excitation energy error Eavg, avoiding state flipping, and considering 400 t2 (600 a) excited states for the silver nanoparticle (fullerene). The optimization yields αAg = 0.036 with Eavg = 5 meV and αC = 0.18 with Eavg = 12 meV for Ag120 and C60, respectively, as reported in Table 5. The amazing accuracy is also retained when the maximum error on all excitation energies Emax, including many optically dark states in C60, 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, the accuracy of TDDFT-as is therefore increased by an order of magnitude, showing that it is a highly competitive approach.

86The 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.

87More recently, TDDFT-as was extended to hybrid functionals and defined across the periodic table, forming the TDDFT-ris model. Using the PBE0 hybrid functional on a test set comprising small-to-large organic molecules, TDDFT-ris has an average error of Eavg = 60 meV compared to an average error of Eavg = 240 meV for sTDDFT., Thus, 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.

88Further 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. 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 T-matrix 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 Ei at the excitation frequencies is often expressed as a power series of the incident field E.

89In 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., To 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 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.

90As 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, 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. The 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.

91A 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 time-reversal symmetry) is used for open-shell systems. For the latter, a noncollinear treatment of the electron spin is applied.

92Current endeavors cover the extension of the 2c formalism to energy gradients and the stress tensor, as well as the inclusion of the current density, 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 and magnetic Hopfions. Another 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. This 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. Therefore, 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.

93The quality of a code strongly correlates with the health and functioning of its developers’ community. TURBOMOLE 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. TURBOMOLE 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.