Developing orbital-dependent corrections for the non-additive kinetic energy in subsystem density functional theory

Subsystem density functional theory (sDFT)^1–3 is based on the commonly used Kohn–Sham density functional theory (KS-DFT) and adopts the idea of partitioning the total molecular system into subsystems based on the electron density. This approach can provide a very favorable computational scaling, allowing one to compute large molecular systems composed of up to a few thousand atoms.⁴ However, due to the non-additive nature of the exchange–correlation (XC) and kinetic energies, the density partitioning gives rise to new terms in the sDFT energy expression. As long as the XC energy is given by a pure functional of the density, the corresponding non-additive XC contribution is trivial to compute. Unfortunately, the orbital-dependent non-additive kinetic energy expression in a monomer basis is unknown and requires an additional approximation to be made.⁵ 

Several different strategies were developed over the past decades to account for this energy contribution or to avoid the problem altogether. Among those are decomposable approximations based on the use of explicit density-dependent kinetic energy functionals (e.g., see Ref. 5), the projection-based embedding theory^6–8 enforcing external orthogonality between subsystem orbitals and ensuring that the non-additive kinetic energy is equal to zero, and the potential reconstruction technique^9–12 that is employed to obtain accurate embedding potentials. Unfortunately, the use of explicit kinetic energy functionals is associated with limitations such as the inability to describe strongly interacting molecules and to cut through covalent bonds, whereas projection-based embedding and potential reconstruction techniques often lead to a large increase in the computational cost of sDFT. Therefore, the problem of accurately approximating non-additive kinetic energy contributions in a cost-efficient way persists. For more information on this topic, we refer to the recent review of sDFT in Ref. 13.

As opposed to projection-based embedding, which enforces external orthogonality between subsystems, an approximate strategy employing Slater determinants composed of non-orthogonal Kohn–Sham-like molecular orbitals (MOs) for the direct evaluation of expectation values of quantum operators was demonstrated within the frozen-density embedding diabatization (FDE-diab) technique.^14,15 This approach was originally developed by Pavanello et al. and received increasing attention over the past years (see Refs. 16–21). Although not being formally exact, it was successfully used for electron- and hole-transfer simulations,^22,23 computations of spin densities,^16–18 and isotropic components of hyperfine coupling constants.²⁰ To the best of our knowledge and somewhat surprisingly, this strategy has never been tested in the context of sDFT for computations of kinetic energy contributions. Being inspired by the previous success of FDE-diab, we make the first important step toward filling this gap by developing orbital-dependent approximations for the non-additive kinetic energy, which are based on non-orthogonal Kohn–Sham-like MOs. In this regard, our main priority does not lie in the construction of a formally exact theory, as opposed to many existing approaches, but rather in creating an alternative route to inexpensive, practical, and fully self-consistent sDFT computations applicable to large molecular systems.

This work is organized as follows. The formal theory behind the sDFT method as well as derivations of orbital-dependent approximations for the non-additive kinetic energy are outlined in Sec. II. Computational details for numerical tests conducted are given in Sec. III. Subsequently, assessments of approximations made and proof-of-principle computations employing the new approximations are presented in Sec. IV. Conclusions to the results and associated discussions are given in Sec. V.

In the following, we briefly outline the theory behind sDFT in Sec. II A. For more details on this topic, the interested reader is referred to reviews in Refs. 5 and 13. Subsequently, new approximations for the non-additive kinetic energy are described in Sec. II B.

Although the approximate non-additive kinetic energy $Tsnad[ρA,ρB]$ is the main error source in sDFT computations (e.g., see Ref. 43), the overall performance of the method also depends on the chosen XC functional. In practical computations, the so-called conjoint functionals,^44,45 a pair of XC and kinetic energy functionals sharing the same form of enhancement factor, are often applied. In this regard and in the context of this work, it was unclear whether the development of new approximations for the non-additive kinetic energy alone would lead to inconsistencies in evaluating the embedding potential. Therefore, similar approximations were derived for the non-additive XC contributions in Secs. S3 and S4.2 in the supplementary material and tested in Sec. IV B.

All computations presented in this work were carried out in a locally modified version of the Serenity program.^46–48 The geometry of the T-shaped Be⁺⋯H₂ electrostatic complex was taken from Ref. 49 and used without further structure optimization. Molecular clusters of small solvent molecules such as water⋯water (H₂O⋯H₂O), water⋯methanol (H₂O⋯CH₃OH), water⋯acetone [H₂O⋯(CH₃)₂O], and methanol⋯methanol (CH₃OH⋯CH₃OH) were optimized with KS-DFT using the PW91 XC functional^50,51 and the valence triple-ζ polarization def2-TZVP basis set.^52,53 The resulting molecular structures are shown in Fig. 1. Subsequently, the sets of displaced structures were created by varying the intermolecular O⋯H and Be⁺⋯H₂ bond distances while keeping other degrees of freedom fixed. No further structure optimization was performed on the resulting geometries.

The generated molecular structures were used in subsequent KS-DFT and sDFT single-point calculations. To that end, the same def2-TZVP basis set^52,53 was employed for all molecular clusters except for Be⁺⋯H₂, which was computed using a smaller 3-21G basis,⁵⁴ unless stated otherwise. The XC PW91 functional^50,51 was consistently applied in both KS-DFT and sDFT computations. We defined and tested several different sDFT computational protocols varying in (i) the choice of the electron density (ρ_Φ or ρ) used for the evaluation of non-additive XC contributions, (ii) the non-additive kinetic energy approximation employed, and (iii) the truncation level M of the Neumann expansion. Thus, standard sDFT computations using PW91^50,51 and PW91k⁵⁵ functionals to account for non-additive contributions are referred to as sDFT/PW91k in the following. sDFT/EA[M] denotes computational protocols, which employ orbital-dependent embedding approximations of up to the Mth order of the Neumann series for the non-additive kinetic energy. If additionally the electron density ρ_Φ obtained from the Mth order truncation of the series from Eq. (31) of the supplementary material is applied to evaluate non-additive XC contributions in conjunction with the PW91 XC functional, the notation sDFT/EA[M, M] is used. In all these cases, fully self-consistent computations were carried out and freeze-and-thaw cycles²⁵ were performed until full convergence of electron densities, i.e., until the sum of absolute element-wise differences between density matrices from subsequent cycles was below the default convergence threshold of 1.0 × 10⁻⁶ a.u.

In what follows, we first analyze the convergence of the Neumann series for a number of molecular clusters in Sec. IV A. To that end, overlap matrices generated with standard sDFT computations employing density-dependent approximations for the non-additive kinetic energy are used. Then, in Sec. IV B, the performance of newly proposed approximations is analyzed in detail for the test case of the Be⁺⋯H₂ electrostatic complex. Finally, in Sec. IV C, a semi-empirical approach to calculating interaction energies is proposed and demonstrated.

For larger molecular systems, spectral radii R(A) can be computed numerically. Such sDFT/PW91k computations for the molecular clusters H₂O⋯H₂O, H₂O⋯CH₃OH, CH₃OH⋯CH₃OH, and H₂O⋯(CH₃)₂O at different intermolecular displacements are presented in Fig. 2. As one can see, the spectral radii R(A) are well below 1 for all molecular clusters at all investigated displacements, thus signifying the convergence of the Neumann series. Furthermore, it can be seen that the spectral radii R(A) tend to zero for larger displacements. This is due to inter-subsystem overlap integrals tending to zero and, hence, A becoming the zero-matrix 0.

Calculations of the truncation error Δ₂ from Eq. (28) for a set of molecular complexes at different intermolecular displacements relative to the equilibrium structure are demonstrated in Fig. 3. As can be seen, very similar results are obtained for all complexes. At intermolecular separations larger than about 2 Å, the overlap matrix S and its inverse S⁻¹ become identity matrices I. This situation is accurately described by the Neumann series truncated at the zero order M = 0, since the corresponding expansion term (I −S)⁰ is also equal to I, whereas all higher-order terms M > 0 yield zero matrices 0. As a result, the error Δ₂ is equal to zero. At shorter displacements, orbital overlaps grow and S⁻¹ start deviating from I. Therefore, the zero-order expansion term is no longer sufficient for describing S⁻¹. Using the Neumann series truncated at the first order M = 1, the error Δ₂ can be kept around zero for intermolecular displacements as short as 0.5 Å. However, higher expansion terms are required for even shorter distances. We find the second expansion order M = 2 sufficient for our applications as it yields very small errors at the equilibrium distance and is less computationally demanding than the third-order expanded series.

For the initial numerical testing of our new approach, we employ a small electrostatic complex of the beryllium cation Be⁺ with a hydrogen molecule H₂. This complex was previously investigated theoretically (e.g., see Refs. 49 and 67–70) and experimentally.⁷¹ It is known that Be⁺ and H₂ are bound by weak (D_e ≈ 0.4 eV) electrostatic and induction interactions resulting in a T-shaped molecular geometry. The ground and first excited electronic states are well-separated from each other, allowing us to apply single-reference electronic-structure methods. These make it a convenient example to test the performance of sDFT-based approaches. In fact, this compound was previously employed for numerical tests in Ref. 16.

Before presenting the results generated with the new computational scheme, we point out two potential issues. First, by introducing new correction terms dependent on the non-orthogonal MOs, we, in principle, incorporate a new electron density $ρΦ(r⃗)$ in sDFT as seen from Eq. (31) of the supplementary material. This new density is not equal to $ρ(r⃗)$⁠, which is given in Eq. (1) and is being optimized/relaxed within sDFT to find the minimum of energy. This creates an inconsistency in the overall approach and means that the new scheme can no longer be considered formally exact as opposed to the standard KS-DFT and sDFT. The second issue stems from the fact that $ρΦ(r⃗)$ does not necessarily integrate to the number of electrons N in the total molecular system when being represented as a truncated Neumann series. This is easy to see from Eq. (31) of the supplementary material, where the zero-order term $ρ(0)(r⃗)$ is equal to the electron density $ρ(r⃗)=ρB(r⃗)+ρA(r⃗)$ and, therefore, by definition, integrates to the number of electrons N. As a consequence, integration of $ρΦ(r⃗)$ results in N only if the integral of the sum of higher-order expansion terms is equal to zero. This requirement is satisfied for the infinitely large series, but does not necessarily hold for all possible truncated expansions. To further analyze this aspect, we computed integrals of $ρΦ(r⃗)$ with sDFT/EA[M] for different truncation orders M as seen in Fig. S2 in the supplementary material. Our results demonstrate that the first-order truncated density expansion, M = 1, does not correctly reproduce the number of electrons in Be⁺⋯H₂ for intermolecular displacements shorter than about 1.5 Å. The deviation reaches about −1 electron for the equilibrium distance (i.e., at the displacement of 0.0 Å). This also shows that the density correction ρ⁽¹⁾ could be negative. However, already for M = 2, the correct number of electrons N = 5 is obtained for all displacements. Furthermore, no negative density areas were found when analyzing the second-order expanded ρ_Φ. Higher-order terms were found to have negligible contributions to the number of electrons.

As the next step, we analyze the performance of new approximations by computing non-additive energy contributions and potential energy curves. To that end, KS-DFT and sDFT/PW91k approaches are used as reference. The results are shown in Fig. 4. As one can see from Fig. 4 (top left), the non-additive energy contributions computed with sDFT/PW91k have opposite signs. The non-additive kinetic energy $Tsnad$ is always positive and, to a large extent, cancels the negative $EXCnad$ contributions to the interaction energy. In fact, it was proven that the non-additive kinetic energy $Tsnad$ computed as a functional of electron density is always non-negative.⁷² The corresponding sDFT/PW91k potential energy curve is qualitatively correct, but strongly underestimates the interaction strength when compared to KS-DFT as seen from Fig. 4 (bottom left). On the contrary, both sDFT/EA[1] non-additive energies are negative and result in qualitatively incorrect and quickly descending potential curves, see Fig. 4 (top left and bottom left). This is probably a consequence of the first-order expansion term of the Neumann series from Eq. (15) featuring a negative sign. This assumption is supported by sDFT/EA[1,1] computations, demonstrated in Fig. 4 (top right and bottom right), which lead to both non-additive contributions having opposite signs to those from sDFT/PW91k. In this case, there is partial cancellation between non-additive contributions. However, the corresponding sDFT/EA[1,1] potential energy curve is overly repulsive and still qualitatively incorrect. Note that the sign of $Ts(1)$ depends on the truncation level M since our non-additive corrections are employed self-consistently. It is negative in sDFT/EA[1] computations but becomes positive when higher-order correction terms are included, i.e., for sDFT/EA[M] with M > 1. This fact is demonstrated in Sec. S7 in the supplementary material. All sDFT/EA[2] and sDFT/EA[2,2] non-additive energy contributions show a much better agreement with sDFT/PW91k results, as shown in Fig. 4 (top left and top right, respectively), and reproduce signs correctly. However, the deviations are still too large to correctly reproduce the shape of the associated potential energy curves. In addition, sDFT/EA[2] and sDFT/EA[2,2] have convergence issues in self-consistent field procedures at intermolecular displacements shorter than 0.0 Å. As can be seen from Fig. 4 (bottom left and bottom right), both curves have a qualitatively incorrect non-bonding character. It should also be noted that potential energy curves and non-additive XC energy contributions $EXCnad$ are very similar for the sDFT/EA[2] and sDFT/EA[2,2] approaches. Therefore, the use of the electron density ρ_Φ for evaluations of the non-additive XC contribution does not significantly change the outcome of computations when second- or higher-order expansion terms are employed. In addition, we assessed the performance of the sDFT/EA[3] computational protocol. The obtained interaction energies were found to be very similar to those from sDFT/EA[2] with the largest deviation of about 0.003 eV. Therefore, we conclude that the Neumann series is sufficiently well converged at the second order of truncation and the use of even higher-order terms is not likely to lead to considerable improvements. Finally, we analyzed the basis set and XC functional dependencies of the sDFT/EA[M] and sDFT/EA[M, M] computational schemes. To that end, similar computations of interaction energies for Be⁺⋯H₂ were performed using the double-ζ def2-SVP and triple-ζ def2-TZVP^52,53 basis sets, and pairs of XC and kinetic energy functionals such as (i) LDA and TF,^73,74 and (ii) BP86^75,76 and LLP91K.⁴⁴ The results of this analysis are shown in Sec. S8 of the supplementary material. A strong dependency of standard sDFT computations on XC and kinetic energy functionals was reported before in the literature (e.g., see Ref. 77) and was, therefore, expected to be observed for sDFT/EA[M] and sDFT/EA[M, M] as well. However, qualitatively same and unsatisfactory results, to those presented in Fig. 4, were obtained with sDFT/EA[M] and sDFT/EA[M, M] showing only rather minor dependencies on the basis set and XC functional applied.

As demonstrated in Sec. IV B, the use of our approximations constructed for the non-additive kinetic energy led to qualitatively incorrect results. Since it was also shown that the Neumann series converges sufficiently well already at the second order, errors in the interaction energy are probably due to other assumptions made. First, as outlined above in Sec. II B, we assumed that the non-interactive kinetic energy of the total molecular system can be computed according to Eq. (10). Second, and as was pointed out previously, we incorporated an inconsistency by introducing a new density $ρΦ(r⃗)$⁠. Analyzing these assumptions in more detail is not a trivial task, which clearly goes beyond the scope of this work. Instead, we adopt a more pragmatic approach and show how qualitatively correct results could still be obtained by introducing purely empirical parameters in orbital-dependent expressions for the non-additive kinetic energy. This decision is motivated by sDFT/EA[M] results, which show very little dependence on the basis set and XC functional used. Therefore, a set of parameters found for one molecular system, and a specific basis set and functional might be transferable to other cases. To analyze this hypothesis, several parametric forms based on Eq. (18) were tested. Among those are the non-additive kinetic energies $Tsnad$ being represented as α(T⁽¹⁾ + T⁽²⁾), αT⁽¹⁾ + α²T⁽²⁾, and αT⁽¹⁾ + βT⁽²⁾. Note that computations were still performed self-consistently and parameters α and β were used as scaling factors for the associated energy- and Fock-matrix contributions. The best results were obtained with the latter fit expression, setting α to −1.0 and finding β by minimizing deviations between sDFT/EA[2] and KS-DFT interaction energy curves of Be⁺⋯H₂. To that end, the PW91 XC functional and 3-21G basis set were employed. The results of this minimization procedure (with β = 0.17) are shown in Fig. 5. It can be argued that the use of a positive α parameter is a more natural choice. However, since $Ts(1)$ and $Ts(2)$ obtained from sDFT/EA[2] are positive, as was mentioned before in Sec. IV B, such a fit does not result in a bound state.

As can be seen from Fig. 5 (left), the fitted sDFT/EA[2] interaction energy curve (denoted as sDFT/EA[2]/fit) shows qualitatively correct behavior and outperforms sDFT/PW91k in reproducing the well-depth D_e value. The corresponding equilibrium distance is shorter than that from KS-DFT but agrees well with that from the original coupled cluster computation from Ref. 49. This improvement in the performance of sDFT/EA[2]/fit is, of course, not surprising since the fitting was performed on the very same molecular cluster. However, deviations from KS-DFT results are still considerable. At the displacement of 0.0 Å, the difference between sDFT/EA[2]/fit and KS-DFT is about 0.1 eV. It is also interesting to note that the integrated sDFT/EA[2]/fit density error, computed according to Eq. (24), is smaller than that from standard sDFT/PW91k at short intermolecular displacements as shown in Fig. 5 (right).

Subsequently, the parameters found for Be⁺⋯H₂ were used without re-optimization for the H₂O⋯H₂O, H₂O⋯CH₃OH, H₂O⋯(CH₃)₂O, and CH₃OH⋯CH₃OH molecular complexes. To that end, the larger def2-TZVP basis set was employed. The results of these computations are shown in Fig. 6. As one can see, all sDFT/EA[2]/fit interaction energy curves show qualitatively correct behavior. Furthermore, issues with converging the self-consistent field procedure are no longer observed. In all four cases, the use of sDFT/EA[2]/fit results in slightly larger equidistant intermolecular distances than those from KS-DFT. Furthermore, the value of the well-depth D_e is underestimated for H₂O⋯H₂O and CH₃OH⋯CH₃OH by about 0.015 eV and overestimated for H₂O⋯CH₃OH and H₂O⋯(CH₃)₂O by about 0.015 and 0.024 eV, respectively. Computing root-mean-square deviations of sDFT/PW91k and sDFT/EA[2]/fit interaction energies from reference KS-DFT results (on 26 equidistantly separated grid points for displacements from −0.5 to 2.0 Å), we obtain errors below about 0.04 eV in all cases. sDFT/PW91k outperforms sDFT/EA[2]/fit for three compounds, namely, H₂O⋯H₂O, H₂O⋯(CH₃)₂O, and CH₃OH⋯CH₃OH, by only about 0.01 eV, whereas sDFT/EA[2]/fit shows a higher accuracy than sDFT/PW91k by 0.006 eV in the case of H₂O⋯CH₃OH. Furthermore, we computed sDFT/PW91k and sDFT/EA[2]/fit integrated density errors relative to KS-DFT. The results of this analysis are presented in Sec. S9 of the supplementary material. As one can see, sDFT/EA[2]/fit outperforms sDFT/PW91k in all cases except for the H₂O⋯H₂O complex. Therefore, we conclude that the proposed semi-empirical approach is robust and has the potential to be transferable between different molecular systems and basis sets. However, a thorough benchmark study is required to find optimal parameters applicable to a broad range of molecular systems and interactions. In this regard, the extended molecular test sets S22x5⁷⁸ and S66x8⁷⁹ are especially attractive. Furthermore, other parametric forms could be investigated. However, this goes beyond the scope of this work and will be conducted elsewhere.

In this work, we presented an alternative route to constructing inexpensive approximations for the non-additive kinetic energy contribution in sDFT. The use of Slater determinants composed of non-orthogonal Kohn–Sham-like MOs for computing the kinetic energy of the total molecular system as an expectation value and the Neumann expansion of the inverse MO overlap matrix constitute the core of this methodology. By deriving the first few terms of the Neumann-expanded kinetic energy expression and taking the corresponding functional derivatives, we constructed a series of orbital-dependent approximations to the non-additive kinetic energy, which can be directly and self-consistently incorporated in sDFT. We also pointed out and discussed the differences and similarities of the obtained expressions with those from the projection-based embedding theory, which employs the Huzinaga operator.⁶ For testing purposes, similar approximations for the non-additive XC energy contributions were derived as well. However, it should be noted that current derivations were carried out for the case of two subsystems and would require introducing additional approximations to be formulated in the general case of multiple subsystems.

Subsequently, we studied the behavior of the Neumann series in detail and discussed the necessary and sufficient conditions for its convergence based on the eigenvalue analysis. For larger molecular systems, an alternative inexpensive technique for performing this analysis was proposed. Furthermore, we demonstrated that the Neumann expansion converges sufficiently well already at the second-order truncation level for molecular systems investigated in this work, namely, water⋯water, water⋯methanol, water⋯acetone, and methanol⋯methanol clusters, and for a large range of intermolecular displacements. The inclusion of higher-order terms affected the results slightly and was found to be important only for very short intermolecular displacements and strongly interacting molecular systems. Therefore, we conclude that the Neumann series is an efficient and robust tool for approximating the inverse MO overlap matrix and avoiding expensive matrix inversion operations.

The derived approximations were applied for computations of potential energy curves of the Be⁺⋯H₂ electrostatic complex and compared against standard KS-DFT and sDFT approaches, which employed explicit functionals of density. Although corrections to the non-additive kinetic and XC energies expanded to the second-order showed an agreement with the corresponding energy contributions from sDFT, the resulting potential energy curves were qualitatively incorrect. The inclusion of higher-order correction terms, as well as the use of different XC functionals and basis sets, did not lead to improved results. In fact, very little dependence of our computational approach on the choice of the XC functional and basis set was observed.

This led us to the idea of introducing empirical parameters into the derived expressions and optimizing them such that deviations to KS-DFT potential energy curves for Be⁺⋯H₂ are minimized. As expected, the use of these new semi-empirical approximations resulted in improved accuracy and better agreement with the KS-DFT reference for the Be⁺⋯H₂ complex. Most importantly, we demonstrated that the very same parameters can be employed for calculations of other molecular clusters while using a larger basis set still resulting in quantitatively correct interaction energies. Thus, for water⋯water, water⋯methanol, water⋯acetone, and methanol⋯methanol complexes, the average deviations from KS-DFT energies were about 0.04 eV. For comparison, the standard sDFT employing the decomposable PW91k⁵⁵ approximation led to comparable accuracy and was even slightly outperformed by our semi-empirical approach in the case of the water⋯methanol complex. Furthermore, our approach demonstrated smaller integrated density errors than sDFT for all complexes except for water⋯water. Therefore, based on these proof-of-principle computations, we conclude that the obtained semi-empirical approximations have the potential to be transferable between different molecular systems.

In conclusion, this work is an important step toward developing novel orbital-dependent approximations for the non-additive kinetic energy in sDFT. Although the current computational protocol requires the use of empirical parameters to correctly reproduce potential energy curves, it also shows a very good agreement with reference KS-DFT results for a set of molecular complexes, weak dependency on the basis set and XC functionals employed, and a high potential of optimized parameters to be applicable to other types of chemical compounds without re-optimization. Furthermore, the use of these semi-empirical approximations comes with a cubic computational cost with respect to the number of atomic orbitals in the active subsystem. Therefore, this approach is only slightly more expensive than sDFT with explicit kinetic energy functionals and could be applied to large molecular systems. However, finding a suitable set of empirical parameters applicable to a broad range of molecular systems and interaction types could be a challenging task and requires more thorough benchmarking, which will be conducted elsewhere.

See supplementary material for (S1) detailed derivations of non-additive kinetic energy contributions, (S2) derivations of corresponding functional derivatives, (S3) new approximations for the non-additive exchange–correlation contributions, (S4) expressions in atomic-orbital representation and associated computational cost, (S5) additional theory aspects on matrix norms, (S6) results of density integration, (S7) consideration of energy contributions to the non-additive kinetic energy, (S8) analysis of basis set dependence, and (S9) integrated density errors.

Funding was provided by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) (Project No. 545861628). L.S.E. acknowledges the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes) for its support. The authors thank Dr. Jan P. Götze for helpful discussions. They also thank the HPC Service of FUB-IT, Freie Universität Berlin, for computing time.

The authors have no conflicts to disclose.

Larissa Sophie Eitelhuber: Conceptualization (supporting); Data curation (supporting); Formal analysis (equal); Investigation (equal); Methodology (supporting); Software (supporting); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Denis G. Artiukhin: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.