Florian
Bogdain
a,
Sebastian
Mai
b,
Leticia
González
bc and
Oliver
Kühn
*a
aInstitute of Physics, University of Rostock, Albert-Einstein-Str. 23-24, 18059 Rostock, Germany. E-mail: oliver.kuehn@uni-rostock.de
bInstitute of Theoretical Chemistry, Faculty of Chemistry, University of Vienna, Währinger Straße 17, 1090 Vienna, Austria
cVienna Research Platform on Accelerating Photoreaction Discovery, University of Vienna, Währinger Straße 17, 1090 Vienna, Austria
First published on 2nd July 2025
A protocol for generating potential energy surfaces and performing photoinduced nonadiabatic multidimensional wave packet propagation is presented. The workflow starts with the parameterization of a linear vibronic coupling (LVC) Hamiltonian using the Green's function – Bethe–Salpeter equation (BSE@GW) approach. In a second step, the LVC model is used as input for multi-layer multi-configurational time-dependent Hartree (ML-MCTDH) wave packet propagation. To facilitate automated ML tree generation, a spectral clustering algorithm is applied based on a correlation matrix obtained from nuclear coordinate expectation values of a full-dimensional time-dependent Hartree (TDH) simulation. The performance of the protocol is tested on the photoinduced spin-vibronic dynamics of a transition metal complex, [Fe(cpmp)]2+. For this example, it is shown that BSE@GW provides a more robust description of the character of the transitions contributing to the absorption spectrum compared to TD-DFT. Furthermore, the LVC parameterization is tested against explicit calculations of potential energy curves to find the validity of the linear approximation over a wide range of normal mode elongation. Finally, the flexibility of spectral clustering is used to generate different ML trees, resulting in very different numerical efficiencies for ML-MCTDH propagation. In terms of electronic structure and dimensionality, [Fe(cpmp)]2+ is a challenging example, suggesting that the new protocol should be applicable to a wide range of systems.
Electronic excitations in transition metal complexes can be classified as being of metal/ligand centered (MC/LC) or charge transfer type, e.g. metal-to-ligand charge transfer (MLCT) or vice versa (LMCT). Stepping from rather successful 4d/5d to 3d transition metals leads to a drastically altered electronic structure as far as the interplay between these types of transitions is concerned. This often comes along with a much reduced MLCT lifetime, which is disadvantageous, e.g. for catalytic applications.7 This is a challenge for ligand design and in particular systems with carbene ligands showed MLCT lifetimes in the nanosecond regime have been reported.8
The challenge to theory provided by transition metal compounds with typical organic ligands is twofold. First, the electronic structure is notorious for the many complications arising from static and dynamic correlations. Second, dynamics simulations for typical systems will have to cope with hundreds of nuclear degrees of freedom (DOFs). In terms of electronic structure theory, density functional theory (DFT) with tailored exchange–correlation (XC) functionals has proven a low-cost yet reasonably accurate alternative to expensive wave function calculations.9 In particular non-empirically tuned range-separated functionals10,11 perform rather well, although often the actual agreement, e.g. with absorption spectra is hard to predict clear beforehand.12 A particular challenge is provided by the interplay between MC and MLCT transition as shown in ref. 13. Absorption spectra at the Franck–Condon geometry are often the target when it comes to the assessment of electronic structure methods. However, as also demonstrated in ref. 13 using trajectory surface hopping (TSH)14 different parameterizations of a range-separated functional can lead to rather different photodynamics even though the absorption spectra are similar.
A viable alternative to common time-dependent density functional theory (TD-DFT) calculations is the Green's function, Bethe–Salpeter equation (BSE@GW) approach.15–18 It has demonstrated notably strong performance for solids19–21 and interfaces,22 but also delivered accurate results for single-point calculations across a wide range of molecular systems.23–26 Recently, we have provided a systematic investigation of electronic excitations of first-row transition metal compounds using BSE@GW.27 It was shown that BSE@GW often outperforms TD-DFT in terms of the prediction of absorption spectra. Even more important, however, is the observation that the mixing of MLCT and MC states can be accounted for with results being robust against variations of the XC functional.
In the context of dynamics simulations, a rather successful approximation to the PES is the linear vibronic coupling (LVC) model,1,28,29 which builds on the availability of gradients on the PES. While being standard for TD-DFT, computations of gradients at the BSE@GW level of theory remain scarce and are restricted to the numerical determination.30 Furthermore, to our knowledge there exists no straightforward and easy-to-use implementation of numerical gradients at the BSE@GW level of theory. This work aims at filling this gap, yielding a LVC model that paves the way to perform TSH31 or full quantum mechanical dynamics simulations in the form of MCTDH.32–34
The multi-configuration time-dependent Hartree (MCTDH) method has proven a versatile approach to high-dimensional quantum dynamics. Introduced in 1990 by Meyer and coworkers35,36 there have been numerous applications, e.g., to nonadiabatic electron-vibrational dynamics,37 including transition metal compounds,1,38,39 excitation energy transfer in pigment-protein complexes,40,41 but also to vibrational dynamics42,43 (for an overview, see also ref. 44). It should be noted that for efficiency MCTDH requires to have at hand a Hamiltonian in the form of a sum of products. Therefore, the LVC model discussed above is ideally suited for the combination with MCTDH dynamics. A boost concerning the dimensionality came with the multi-layer extension (ML-MCTDH) introduced by Wang and Thoss45 and later generalized in ref. 46 and 47. While in MCTDH the wave packets are expanded into Hartree products of single-particle functions (SPFs), in ML-MCTDH the idea is to combine strongly correlated DOFs into multidimensional particles, which in turn are expanded in MCTDH form. This is repeated until one- or low-dimensional particles are obtained. The resulting wave packet has a tree-like (ML-tree) structure that needs to be converged with respect to the number of SPFs assigned to the branches. This results in the challenge to find cost-efficient ML-trees, since the numerical performance in terms of CPU time can vary vastly for different trees.48 While there are simple guidelines for combining modes (e.g. similar frequencies), standardized routines for the systematic generation of reliable ML-tree structures remain scarce. Recently, Mendive-Tapia et al. proposed to use hierarchical clustering for generating ML-trees, focusing on vibrational dynamics.49 In this work we will present a method based on spectral clustering,50 yielding ML-trees for nonadiabatic electron-vibrational dynamics.
The computational workflow presented in this paper starts with the determination of a full-dimensional LVC coupling Hamiltonian using BSE@GW (Section 2.2) that serves as an input for the generation of ML-trees using spectral clustering (Section 2.3). The new protocol will be applied to the [Fe(cpmp)]2+ complex (cpmp = 6,2′-′carboxypyridyl-2,2′-methylamine-pyridyl-pyridin) shown in Fig. 1. The design of this pseudo-octahedral complex follows two strategies for obtaining long-lived MLCT states, i.e. forming a close to octahedral coordination sphere such as to maximize the ligand field splitting and destabilizing the MC states by electron-rich ligands while simultaneously stabilizing the MLCT states by π-acceptor groups.51
![]() | ||
Fig. 1 Chemical structure and geometry of the [Fe(cpmp)]2+ complex investigated in this work (for synthesis and experimental characterization, see ref. 51). |
Previous research has studied the excited state dynamics of this system using both classical TSH13 and ML-MCTDH52 dynamics simulations. In both cases, TD-DFT was the underlying level of theory at which the PES was calculated. The range-separated LC-BLYP(α, ω) functional was used, where the parameters α and ω were tuned non-empirically. The results presented in Section 4 primarily serve to illustrate the new protocol rather than introduce novel insights into the photophysics of this particular complex. However, this paper includes a valuable comparison of the absorption spectrum assignments for TD-DFT and BSE@GW, along with a detailed investigation of the validity of the LVC model.
![]() | (1) |
The matrix elements are given by55
![]() | (2) |
![]() | (3) |
![]() | (4) |
![]() | (5) |
![]() | (6) |
![]() | (7) |
![]() | (8) |
![]() | (9) |
![]() | (10) |
This LVC Hamiltonian may also be extended by an additional term, accounting for the interaction with an external electric field, (t), in dipole approximation (coordinate independent diabatic dipole vector matrix elements
mn)
![]() | (11) |
E(t) = E0![]() | (12) |
![]() | (13) |
![]() | (14) |
![]() | (15) |
![]() | (16) |
The extension to more layers is straightforward. In this case, each layer introduces a new set of logical coordinates, each derived from a subset of the coordinates combined in the preceding layer. The resulting structure of the wave packets can be visualized using so-called ML-trees.46 Equations of motion for SPFs and coefficient vectors can be obtained by applying the Dirac–Frenkel time-dependent variational principle.47 Convergence can be controlled by inspecting the so-called natural orbital populations related to the SPFs. Below we will define thresholds for the largest population of the least occupied natural orbital, pNO.
Finally, we note that the limit where the wavepacket ψm(Q,t) is expressed as a simple Hartree product for all physical coordinates is called the time-dependent Hartree (TDH) approach.
![]() | (17) |
Next, we define the adjacency matrix W of an associated weighted graph with matrix elements
Wξξ′ = |Cξξ′| − δξξ′. | (18) |
In this work, we relied on a spectral clustering algorithm,50,68–71 in contrast to ref. 49 which used hierarchical clustering. How the two methods compare in the context of ML-MCTDH efficiency will require further investigation. We prefer spectral clustering as it provides a means to directly control the size of the clusters, which will be important for the discussion below. In addition, it circumvents the problem of predefining a metric that dictates the structure of the obtained dendrogram in hierarchical clustering. In short, spectral clustering proceeds as follows. First, the symmetrized Laplacian is computed, for which the lowest k eigenvalues and eigenvectors are determined. Afterward, each row of the matrix, whose columns are the eigenvectors, is treated as a data point. Data points are grouped into k clusters using k-means. Given a clustering at the level of the full correlation matrix, the procedure is repeated until subclusters of either size one or two are obtained.
The final result is a ML-tree for which the number of SPFs at each vertex has to be obtained. This is done by iteratively propagating the equations of motion while increasing the number of SPFs until a predefined threshold for pNO has been reached. For the specific examples below this was done in the time interval [0:1000 fs].
All quantum chemical calculations were performed using TURBOMOLE version 7.8.1.72,73 In particular, the BSE@GW and TD-DFT calculations were done with the ESCF module74,75 using the TDA. Self-consistent field calculations were performed in the RIDFT76 module, making use of the resolution of identity (RI) approximation. SOC matrix elements were calculated with TURBOMOLE's PROPER module using the effective one-electron operator option, which includes a screened nuclear charge.64,77 The range-separated functional LC-BLYP was used with parameters taken from ref. 13. The equilibrium geometry was obtained through DFT geometry optimization using the LC-BLYP (α = 0.2, ω = 0.08) functional, which resulted in a C2 symmetric structure, although symmetry was not enforced and has not explicitly been used in subsequent calculations.
All calculations were performed using a def2-TZVP basis set. The auxiliary basis sets (jbas and cbas) were always of the same quality as the regular ones. Convergence criteria were set to $scfconv = 7 and $rpaconv = 5. The m5 grid option was always chosen. UV-Vis spectra were computed by broadening vertical transition energies by a Lorentzian function with a full-width-half-maxima of 0.25 eV. The length gauge was always chosen as the height of the absorption peak. The analysis of the transition density matrices was performed with the THEODORE program package version 3.2.78,79
All calculations regarding the LVC model were carried out in a modified version of the SHARC 3.0 program package.14,80 To do so, a new interface between SHARC and TURBOMOLE's ESCF/PROPER module was developed, which follows the established SHARC routines. It enables one to carry out a full LVC parametrization32 of the PES using both the BSE@GW and TD-DFT levels of theory. The computation of the overlap matrix was also performed in SHARC's WFOVERLAP module.66 Our LVC model included the lowest 10 singlet (including the singlet ground state) and 14 triplet states. This covers the range of excitation energies up to about 2.5 eV, thus enabling a comparison with the transient absorption experiments conducted to study the MLCT states in this range and reported in ref. 13. We opted for a def2-TZVP basis set as we observed that the BSE@GW calculations demonstrate a high sensitivity to basis set quality (see Fig. S1 in the ESI†). Both the RI and TDA were employed and the first 49 core molecular orbitals were discarded to calculate the wave function overlaps. A wave function truncation threshold of 0.998 was chosen. The same convergence criteria and grid settings for the BSE@GW calculations were chosen as those used for the spectra calculations. For this model, we used a displacement factor of 0.05. To calculate the GW quasiparticle energies, contour deformation techniques (CD)81–83 were used with an orbital range of [157,177] for which the correlation part of the self-energy was calculated. The eigenvalue-only self-consistent GW scheme (evGW) was used to converge quasiparticle energies. The range-separated functional LC-BLYP (α = 0.2, ω = 0.08) was taken as a starting point.
ML-MCTDH computations were performed with the Heidelberg MCTDH package Version 8.6.5.35,36,84,85 For the primitive basis in eqn (16) a ten-point harmonic oscillator discrete variable representation was used. The equations of motion were solved employing variable mean fields, together with a 5th-order Runge–Kutta solver. The ML-trees were generated using spectral clustering as implemented in ref. 86.
All computations were performed on an Intel Xeon Gold 6326 CPU@2.90 GHz, which served as the benchmark system for comparing the computational demand of the different calculation setups.
Second, we will perform ML-MCTDH quantum dynamics simulations using the obtained LVC Hamiltonian. Here, the emphasis will be put on the performance of the spectral clustering algorithm for ML-tree generation.
![]() | ||
Fig. 3 Absorption spectra, according to TD-DFT and BSE@GW using the LC-BLYP functional for different α and ω values. The gray-shaded area refers to the experimental data from ref. 51. |
In contrast, the BSE@evGW framework provides an improvement of the absorption spectrum when comparing with the TD-DFT results for (α = 0.2, ω = 0.08), albeit it does not perform as well as the TD-DFT (α = 0, ω = 0.14) protocol. At first glance this appears disappointing, particularly recalling the superior performance for similar Fe-compounds reported in ref. 27. In that work, it was also shown that there is only a small dependence of the BSE@GW spectra on the XC-functional, implying for the present case that along this line there is perhaps little room for improvement of the result shown in Fig. 3. In fact, the only variable to impact the resulting spectra appears to be the amount of exact exchange used as seen in the ESI,† Fig. S12 and S13.
However, the correct reproduction of the absorption spectrum is only one part of the picture. In addition, the CT characters must also be described properly, to obtain an accurate characterization of the underlying electronic transitions. To scrutinize this point, an analysis of the TDM was performed for the singlet and triplet excitations contributing to the absorption spectrum in Fig. 3. Here the singlet ground state geometry was assumed. We will compare the order of the states with CASPT2 results on the lowest excitations, which gave for singlets/triplets that the MLCT transitions is energetically below/above the MC one.13 Note that in ref. 13 the active space was designed with a focus on the lowest singlet and triplet transitions only. That is, there is no CASPT2 grade information on the absorption spectrum and its assignment.
The results in Fig. 4 show that LC-BLYP (α = 0, ω = 0.14) provides a reasonable description of the absorption spectrum, even though it incorrectly predicts the ordering of the excited triplet states, placing the MLCT state below the MC state, which does not align with previous CASPT2 computations.13 In fact, an increase of the short-range exact exchange contribution (α) is required to obtain the correct ordering, which simultaneously worsens the alignment with the experiments. Further, it bears the risk of overcompensation up to the point where the MC states are placed below the MLCT in the singlet absorption spectrum. We obtained the correct ordering of the lowest MLCT/MC transitions for the parameters (α = 0.1, ω = 0.11).
![]() | ||
Fig. 4 Analysis of the singlet and triplet TDM for the different tuning parameters and methods used in Fig. 3. The geometry corresponds to a vertical excitation from the singlet ground state. Note that the pronounced mixing a MC and MLCT states evidences the deviation from octahedral symmetry of this complex. Numerical values for transition energies are given in Tables S1 and S2 of the ESI.† |
In comparison, BSE@evGW yields the correct ordering in both the lowest singlet and triplet states without requiring the separate tuning of parameters, aligning well with observations for similar systems.27 This advantage of the BSE@GW framework over TD-DFT motivates its usage as a basis to construct an LVC model and perform subsequent dynamics simulations.
A histogram of these obtained values is given in Fig. 5. It shows that most of the coupling constants are small, i.e. below about 0.07 eV for κ and λ and below about 0.025 eV for SOC. While this is comparable to the TD-DFT results reported in ref. 52, the majority of κ values are a bit larger (below about 0.12 eV) in TD-DFT. Above the mentioned values, some scattered coupling parameters do not have exact counterparts in TD-DFT. Notice that the TD-DFT results had been obtained using LC-BLYP (α = 0.2, ω = 0.08), which according to our discussion in the previous section (cf.Fig. 4) shows a different composition of transitions as compared with BSE@evGW. This will have, of course, an influence on the coupling parameters, in accordance with the present findings. For instance, the MC character will be reflected in the shift of the PES. In light of the discussion of Fig. 5, which concluded on a more robust description of the character of transitions by BSE@evGW, the κ and λ values should be considered more reliable than the published TD-DFT ones.52
It is interesting to ask whether the LVC model itself is applicable. The general quality of the κ and λ parameters can be estimated by comparing the resulting adiabatic PES (obtained by diagonalization of the LVC Hamiltonian) with single-point calculations for a set of geometries, displaced along particular normal modes. In Fig. 6 respective PES cuts along the mode with ξ = 25 having a frequency of 186 cm−1 are shown. It is a symmetric breathing mode meaning that a positive displacement corresponds to an elongation of the Fe–N bond lengths (results for other modes can be found in the ESI,† Fig. S2).
Overall, one observes that the calculated LVC PES, as obtained with BSE@evGW (Fig. 6(a, b)) as well as TD-DFT (Fig. 6(c, d)), decently align with the single-point calculations, especially around the equilibrium geometry, where the PES show harmonic characteristics (for a more anharmonic case, see Fig. S2(a) in ESI†). In particular, we observe a good agreement for the lowest triplet states, which are often difficult to describe properly. Besides the fairly good performance of the LVC approximation as such, the PES of BSE@evGW and TD-DFT show distinct differences. For instance, the lowest excited singlet state potential curve is considerably more displaced (large κ) in TD-DFT as compared with BSE@evGW. This is not surprising given the fact that according to Fig. 4, TD-DFT predicts a dominant MC (anti-bonding) character. The displacement of the potential curve for the lowest triplet state, on the other hand, is rather similar in BSE@evGW and TD-DFT, in accord with its comparable character in Fig. 4. Finally, we notice that BSE@evGW predicts a higher density of states in this energy range as compared with TD-DFT.
As has been highlighted for the case of TD-DFT and different LC-BLYP tuning parameters in ref. 13 already, differences in coupling parameters can lead to rather different dynamics, even though the absorption spectra are similar. We will not further elaborate on a comparison of BSE@evGW and TD-DFT in terms of the quantum dynamics. Instead, the focus of the next section will be on general aspects of the performance of the ML-tree generation using spectral clustering.
Based on this procedure, we generated five different trees for detailed investigation. All trees are shown in the ESI,† Fig. S3–S8. The first four were binary trees, where W was decomposed into two clusters in each layer, but with a varying threshold for entries. In this instance, we have chosen the four thresholds ΔW 0 (I), 0.1 (II), 0.2 (III), and 0.3 (IV). Entries of W that fell below these values were set to 0 accordingly, meaning that the two modes are considered uncorrelated. In principle, a small threshold will give a larger set of isolated modes, leading the clustering algorithm to generate a deeper ML-tree (cf. ESI,† Fig. S3–S6). This occurs because the likelihood increases that, at a specific node, the clustering algorithm groups a small set of normal modes into one branch, while passing the majority of normal modes further down to the next layer. This way a small threshold in the adjacency matrix may be seen as a tool to deepen the ML-tree by choice. However, increasing the threshold could lead to a situation where most entries of the adjacency matrix are located around the diagonal. This may aid the spectral clustering algorithm in identifying highly correlated pairs and thus reducing the depth of the ML-tree (cf. ESI,† Fig. S6). For the present models (I) to (IV) we have 10, 11, 18, and 17 layers, respectively. In general, having deeper ML-trees could be desirable as related research has shown that this may lead to a lower computational effort.87 However, at the same time, this could increase the number of SPFs, which may counterbalance the positive effect on efficiency. The remaining tree (V) was chosen such that the full adjacency matrix was divided into clusters of different sizes. The first layer was decomposed into 3 clusters, followed by two for all remaining layers. Of course, the present choice of ML-trees is just a subset of in principle possible trees, picked by intuition and serving the purpose of comparison.
Below we will discuss the performance and convergence of the different ML-trees in terms of the total (singlet + triplet) MLCT population. The connection to the overall dynamics is provided in Fig. 8. These results were obtained with tree (IV), which gave the best performance as will be detailed below. In Fig. 8 one observes the expected behavior for this type of system. Initially, the laser pulse excited the 1MLCT states which rapidly convert into 3MLCT states due to intersystem crossing. Internal conversion subsequently leads to a population relaxation in the triplet manifold towards the low-lying 3MC states. Depopulation of MLCT states in the excitation window takes place within a few hundred femtoseconds, which is in accord with the experimental observation by Moll et al.51 In passing we note that the states with the highest energy in the given interval are only marginally populated, thus giving support for the restriction to 10 singlet and 14 triplet states for the given excitation conditions.
![]() | ||
Fig. 8 Energy-resolved ML-MCTDH population dynamics of diabatic spin-free states for the 63-dimensional LVC model and for tree (IV) with pNO converged to 0.2%. |
In what follows, the convergence behaviour of the five ML-trees is investigated by gradually decreasing the threshold for pNO. In Fig. 9(a) the total MLCT population for tree (IV) is shown for pNO ranging from 5% to 0.2%. Only minor differences are observed between the two curves with the tightest convergence criteria. We therefore take the curve with the tightest convergence criteria as a reference to compare them against the remaining models. It is interesting to note that only the reference calculation reproduces dynamics in accordance with what would be expected for this type of system. In contrast, higher convergence thresholds lead to oscillations as well as long-lived populations in the higher-lying excited triplet states (for an example, see ESI,† Fig. S1).
![]() | ||
Fig. 9 Convergence behavior for MLCT diabatic populations for models (a) (IV), (b) (I), and (c) (V). The graphs are labeled according to the convergence threshold set pNO × 1000. The tightest convergence criteria for model (IV) of 0.2% will be set as a reference. Notice that in all curves there is a slight initial rise, which is discussed in ESI,† Section S7. |
The required convergence threshold and the speed at which convergence is achieved depend strongly on the structure of the ML-tree. This can be seen by comparing the different panels of Fig. 9. In all cases, the not sufficiently converged results show a slower decay of the MLCT population, i.e. in particular a trapping of population in high-lying excited triplet states. Increasing the convergence threshold will lead to the dynamics trending towards the reference results, although with varying speeds (to make sure that this is not a feature of the particular model, it has been checked for a 42-dimensional model, see ESI,† Fig. S11).
It appears that even the tightest convergence thresholds used remain insufficient for model (I), as rather large changes are observed. The same behavior can be observed for systems (II) and (III), which all appear to converge slower than (IV) (cf. ESI,† Fig. S9). Thus we can conclude that deepening of the ML-tree structure by setting a threshold for the adjacency matrix may indeed improve convergence with respect to pNO, but at the time being it is not clear whether there is an optimal choice. Model (V) (having a comparatively flat tree) in Fig. 9(c) demonstrated by far the highest computational demand, combined with the slowest convergence behavior. In passing, we note that an ML-tree with 4 clusters at the top, followed by 3 in the second layer and 2 for all remaining ones was also tested. However, the convergence was too slow and the computational demand too high to obtain results in a time comparable to all other models.
The CPU time requirement of all ML-trees is shown in Fig. 10. Here, we notice that all binary ML-trees result in a comparable computational demand at a given threshold pNO. This observation should be highlighted in particular, as the ML-trees of the first four models have a completely different structure and, depending on the depth of the structure, a significantly higher amount of SPFs (see ESI,† Fig. S3–S8). We conclude that a binary ML-tree structure should consequently be preferred over other related structures for this system to achieve a reliable low computational cost, as all results trend toward the converged reference calculation with comparable computational cost.
![]() | ||
Fig. 10 CPU time of the different ML-tree structures for different convergence thresholds. All results are normalized to the model (IV) simulation for a threshold of pNO = 0.5%. |
The new protocol was applied to the photoinduced dynamics of an Fe(II) complex. Here it should be noted that BSE@evGW gave a reasonable absorption spectrum if compared with experiment and TD-DFT using an optimized LC-BLYP functional. More importantly and in contrast to TD-DFT, however, BSE@evGW provided a reasonable assignment of the lowest singlet and triplet transitions. Since the characters of transitions determine vibronic couplings, the BSE@evGW derived LVC parameters can be considered superior if compared to TD-DFT. ML-MCTDH dynamics after pulsed excitation yielded results for the MLCT populations in accord with experiment. The performance of different ML-trees was rather different in terms of convergence with respect to the number of SPFs and associated computational demands. Although there is not yet a strict method for tree optimization, for the present case it has been shown that spectral clustering with binary trees and using a threshold for the adjacency matrix leads to efficient ML-MCTDH setups.
The investigated transition metal complex can be considered a challenging case in terms of electronic structure and dimensionality of spin-vibronic dynamics. This suggests that the presented protocol should be applicable to a broad range of molecular systems.
Footnote |
† Electronic supplementary information (ESI) available: Further information on the choice of the basis set, the validity of LVC model, the structure of the ML-trees, the adjacency matrix of model (IV), the population dynamics for a larger convergence threshold, a 42-dimensional model, and some discussion of accuracy of state couplings. See DOI: https://doi.org/10.1039/d5cp01208b |
This journal is © the Owner Societies 2025 |