Introduction

The Jahn-Teller (JT) effect leads to a distortion of high-symmetry structures of polyatomic molecules with degenerate electronic states [1,2,3]. Computational work on the influence of JT and/or PJT effects on the structure and stability of JT systems in both gas phase, solution, and recently also in the solid state [4, 5] include, e.g., monitoring changes in geometry and symmetry [6,7,8], or changes in molecular properties such as natural bond orbitals (NBO)s and related NBO charges [9]. Also reactivity indices [10] including Fukui functions have been applied to analyze JT distortions [11,12,13,14,15] and dynamic processes such as bond pseudo-rotation triggered by JT distortion have been explored [16,17,18]. In astrochemistry, investigating structural and electronic changes from associated JT distortion through vibrational spectroscopy of diamondoids has also garnered increasing interest [19,20,21].

In this work we added another perspective, namely the elucidation of JT distortions based on vibrational spectroscopy. For this purpose, we applied two specific aspects of the local vibrational mode (LVM) theory developed in our group [22, 23], the characterization of normal mode (CNM) procedure [24,25,26,27] and an adiabatic connection scheme (ACS) connecting local and normal vibrational modes [22, 23, 28], exploring their potential use for (i) probing JT distortions in transition metal complexes and (ii) subsequently explaining structure and stability of transition metal ion compounds, providing a framework for a new understanding of their structural preferences connected with JT distortions.

We targeted two distinct categories displaying JT distortions, one category with moderate JT distortions accompanied with bond lengthening and one category with slight JT distortions accompanied with bond shortening, often more difficult to capture. As representative for the first category we chose the \(\mathrm {T_h}\) symmetric hexaaquachromium(III) cation, \(\mathrm {[Cr{(OH_2)}_6]^{3+}}\) (1 in Fig. 1), and its transition to \(\mathrm {[Cr{(OH_2)}_6]^{2+}}\) with \(\mathrm {D_{2h}}\) symmetry (2 in Fig. 1). As representative for the second category we chose ferrous and ferric hexacyanides, \(\mathrm {[Fe{(CN)}_6]^{4-}}\) with \(\mathrm {O_h}\) symmetry (3 in Fig. 1) and its transition to \(\mathrm {[Fe{(CN)}_6]^{3-}}\) with \(\mathrm {D_{4h}}\) symmetry (4 in Fig. 1).

Fig. 1
figure 1

Schematic representation of test examples 1–4, showcasing the point groups associated with the equilibrium geometries obtained at the PBE0/Def2-TZVP level of theory. The difference in axial bond lengths of 2 and 4 has been exaggerated for illustration purposes

Full size image

Both hexaaqua complex ions 1 and 2 are high-spin [29], with the JT distortion leading to an elongation of the axial Cr–O bonds [30, 31], see Fig. 1. Conversely, both ferrous and ferric hexacyanides 3 and 4 are low-spin [29]. Neutron diffraction studies have revealed that 3 adopts a regular octahedral configuration \(\mathrm {O_h}\) [32] whereas for 4, a d\(^5\) complex, a slight JT distortion is expected [33, 34]. As confirmed by ab initio calculations, 4 exhibits smaller Fe–ligand \(\pi \) back-donation compared to 3 [33]. X-ray diffraction analysis, supported by infrared (IR) spectroscopy and ligand field theoretical calculations, predicts in addition that the iron–ligand bonds in 4 are predominantly of \(\sigma \) type, with minimal or no \(\pi \) contribution [35]. The 3/4 transition stands for an unusual example of JT distortion, and to the best of our knowledge, an explanation of its intricacies is still pending.

Ultimately, this proof-of-concept study was aimed at the investigation of the 1/2 case as a commonly known example that exhibits a JT distortion followed by the more difficult 3/4 case, comparing common and distinct features. The CNM/ACS protocol for characteristic normal vibrational modes was employed as a novel tool, contributing to a deeper understanding of Jahn-Teller distortions by offering valuable insights into the molecular fragments responsible for symmetry disruption at the molecular level.

Computational methods

Geometry optimization and Hessian matrix calculation

All electronic structure calculations were carried out using the Gaussian 16 quantum chemistry program [36] at the density functional theory (DFT) level [37]. An ultra-fine integration grid and a tight convergence criterion were applied for the self-consistent field procedure [38]. The equilibrium geometries of the test examples 1–4, as well as their corresponding Hessian matrices and normal vibrational modes, were obtained using the PBE0 functional [39] in conjunction with the Def2-TZVP basis set [40]. PBE0 has demonstrated favorable outcomes across diverse properties of molecules and materials, encompassing spectroscopic characteristics, making it highly versatile for applications in both quantum chemistry and condensed matter physics [39]. Its performance for transition metal compounds has also been assessed [39]. The hexaaquachromium cations were calculated with unrestricted DFT (UDFT), and ferrous and ferric hexacyanides were calculated with restricted DFT (RDFT) and UDFT, respectively [41].

Local vibrational mode theory

Normal vibrational modes throughout polyatomic systems are generally delocalized [42, 43], which imposes an important limitation on determining the intrinsic bond strength through the direct use of normal mode frequencies and normal mode force constants. This is precisely where the LVM theory comes into play. LVM was originally introduced by Konkoli and Cremer [44, 45] and further developed by our group. For more details on LVM, including its theoretical foundation and broad application spectrum across chemistry and beyond, the reader is referred to two recent review articles and literature cited therein [22, 23].

Normal vibrational modes \(\textbf{d}_n\) given in internal coordinates \(q_n\) (with \(n = 1,\cdots ,N_{vib}\) and \(N_{vib}\) = (3N - 6) for non-linear N-atomic complexes and (3N - 5) for linear N-atomic complexes) and the diagonal normal mode force constant matrix \(\textbf{K}\) given in normal coordinates \(Q_n\) (being available in the output of any modern quantum chemistry software that performs a standard normal mode analysis) are transformed into their local mode counterparts leading to local mode vectors \(\textbf{a}_n\) associated with internal coordinates \(q_n\) via:

$$\begin{aligned} \textbf{a}_n = \frac{\textbf{K}^{-1}\textbf{d}_n^{\dagger }}{\textbf{d}_n\textbf{K}^{-1}\textbf{d}_n^{\dagger }}. \end{aligned}$$
(1)

By utilizing the transformation matrix \(\textbf{L}\) [42] the local mode vector \(\textbf{a}_n\) can be easily converted into Cartesian coordinates \(\textbf{x}\), resulting in \(\textbf{a}_n^x\) given by the following expression:

$$\begin{aligned} \textbf{a}_n^x = \textbf{L} \textbf{a}_n. \end{aligned}$$
(2)

The calculation of the corresponding local mode force constant \(k_n^a\) can be performed using the following expression:

$$\begin{aligned} k_n^a = \textbf{a}_n^{\dagger } \textbf{K} \textbf{a}_n. \end{aligned}$$
(3)

This enables the computation of the local mode frequency \(\omega ^a_n\):

$$\begin{aligned} (\omega ^a_n)^2 = {(4 \pi ^2 c^2)}^{-1} \frac{k_n^a}{m^a_n}, \end{aligned}$$
(4)

with \(m^a_n\) being the local mode mass.

Local mode force constants \(k^{a}\) have proven to be a reliable tool to quantify the strength of covalent chemical bonds including metal-ligand bonds and weak chemical interactions such as halogen, chalcogen, pnicogen, tetrel, or hydrogen bonding. For a summary, see Refs. [22, 23] and citations therein. We also derived a new measure for the assessment of metal-ligand bonding, the metal-ligand electronic parameter (MLEP) [46,47,48,49], applied, e.g., in recent work for the investigation of iron-ligand bonding in carboxy myoglobins and carboxy neuroglobins [50, 51] and a new assessment of non-covalent \(\pi \)–interactions in mutated aquomet-myoglobin proteins [52]. Furthermore, our recent incorporation of local mode force constants into lanthanide spectroscopy has provided a novel explanation for the inverse relationship between lanthanide-ligand strength and ligand effective polarizability, offering insights into bonding in lanthanide chemistry [25].

An adiabatic connection scheme (ACS) [28], founded on the Decius compliance matrix [53], denoted as \(\mathbf {\Gamma } = (\textbf{F}^q)^{-1}\), where \(\textbf{F}^q\) represents the Hessian matrix in internal coordinates, in conjunction with the Wilson \(\textbf{G}\) matrix, establishes a significant one-to-one correspondence between every complete and non-redundant set of local vibrational modes and normal vibrational modes:

$$\begin{aligned} ( \textbf{G}_d + \lambda \textbf{G}_{od} ) \textbf{R}_{\lambda } = ( \mathbf {\Gamma }_d + \lambda \mathbf {\Gamma }_{od} ) \textbf{R}_{\lambda }\mathbf {\Lambda }_{\lambda }, \end{aligned}$$
(5)

where \(\textbf{R} = \mathbf {\Gamma }^{-1}\textbf{D}\) with \(\textbf{D}=\textbf{B}\textbf{L}\) (\(\textbf{B}\) is the Wilson B matrix [42]). \(\textbf{G}_d\) and \(\mathbf {\Gamma }_d\), and \(\textbf{G}_{od}\) and \(\mathbf {\Gamma }_{od}\) are the diagonal and off-diagonal parts of Decius and Wilson matrices, respectively [42]. The local vibrational modes (\(\lambda =0\)) are adiabatically transformed into their normal mode counterparts (\(\lambda =1\)) by a shift parameter \(\lambda \) [28].

ACS serves as the fundamental basis for the characterization of normal mode (CNM) procedure [22, 23, 54, 55] which enables the decomposition of each normal vibrational mode \(\textbf{l}{_\mu }\) (in Cartesian coordinates) into local vibrational modes contributions. The CNM analysis is conducted for a complete and non-redundant set of \(N_{vib}\) local vibrational modes \(\textbf{a}_n\), offering a powerful tool for quantitatively assessing the composition of normal vibrational modes and facilitating the interpretation of vibrational spectra. The overlapping encoded by \(S_{n\mu }\) is [22, 23, 54, 55]

$$\begin{aligned} {S}_{n\mu }=\frac{\langle \textbf{a}_n^x | \textbf{F}^x | \textbf{l}_{\mu } \rangle ^2}{\langle \textbf{a}_n^x| \textbf{F}^x |\textbf{a}_n^x\rangle \langle \textbf{l}_{\mu }| \textbf{F}^x |\textbf{l}_{\mu } \rangle }. \end{aligned}$$
(6)

Thus, the local mode contribution \(C_{n\mu }\) (local mode character) of local vibrational mode \(\textbf{a}_n^x\) to normal vibrational mode \(\textbf{l}_{\mu }\) (in Cartesian coordinates) can be calculated as shown below [22, 23, 54, 55]

$$\begin{aligned} C_{n\mu } = \frac{S_{n\mu }}{\sum _{m=1}^{N_{vib}}S_{m\mu }}. \end{aligned}$$
(7)

Our CNM analysis has facilitated a new way of analyzing and interpreting vibrational spectra stretching from smaller molecular systems to large biomolecules [22, 25,26,27, 54,55,56,57,58], being described with hybrid quantum chemistry and molecular mechanics (QM/MM) [59,60,61,62] thanks to our novel LModeAGen protocol, which automatically determines a complete non-redundant set of local mode parameters [26] and which has been coupled for QM/MM systems with our generalized subsystem vibrational analysis (GSVA), extracting intrinsic fragmental vibrations of any fragment/subsystem from the whole system via the evaluation of the corresponding effective Hessian matrix [63, 64].

Fig. 2
figure 2

CNM plots capturing metal–ligand (a, b), and O–H (c, d), stretching vibrations for \(\mathrm {[Cr{(OH_2)}_6]^{3+} (T_h)}\) and \(\mathrm {[Cr{(OH_2)}_6]^{2+} (D_{2h})}\) at the PBE0/Def2-TZVP level of theory

Full size image
Fig. 3
figure 3

CNM plots capturing metal–ligand (a, b), and C–N (c, d), stretching vibrations for \(\mathrm {[Fe{(CN)}_6]^{4-} (O_h)}\) and \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\) at the PBE0/Def2-TZVP level of theory

Full size image

For the test examples 1–4, the generation of non-redundant and complete sets of local vibrational modes was achieved using the LModeAGen protocol [26]. Subsequently, the local mode analysis was conducted using the standalone LModeA package [65].

Electron density analysis

The electron density analysis was performed with the AIMAll program [66]. Applying Bader’s quantum theory of atoms in molecules (QTAIM) [67,68,69], we complemented LVM data with bond properties calculated at the bond critical point \(\textbf{r}_b\) of the electron density \(\rho (\textbf{r})\) between the two atoms forming a chemical bond or weak chemical interaction.

The Cremer-Kraka criterion [70, 71] was applied to assess the covalent character of the chemical bond and/or weak chemical interaction which relies on the local energy density \(H(\textbf{r})\) given by the equation:

$$\begin{aligned} H(\textbf{r}) = G(\textbf{r}) + V(\textbf{r}) \,, \end{aligned}$$
(8)

where \(G(\textbf{r})\) is the kinetic energy density (positive, destabilizing), and \(V(\textbf{r})\) is the potential energy density (negative, stabilizing). If the value of \(H(\textbf{r}_b)\) is negative the bond/interaction is predominately covalent. Conversely, if it is positive the bond/interaction is primarily electrostatic [70, 71].

Results and discussion

Figures 2 and 3 provide an overview of the CNM results, displaying the local mode contributions within the characteristic frequency ranges associated with metal–ligand and ligand stretching vibrations for each of the test examples 1–4 from Table 1. In the hexaaquachromium complexes, there is an overall increase in the localization of the characteristic normal vibrational modes corresponding to Cr–O and O–H stretching vibrations, as reflected in Fig. 2(a)–(d). In contrast, ferrous and ferric hexacyanides, depicted in Fig. 3(a)–(d), show an opposite trend for Fe–C stretching vibrations. Overall, the normal vibrational modes remain delocalized in both scenarios, as depicted in Figs. 2 and 3. Such delocalization is a characteristic feature observed in systems with significant symmetry, providing valuable insights into the vibrational behavior and the functional group interactions in such compounds. Considering the significance of a reliable bond strength descriptor for understanding coordination sites, the observed delocalization of the characteristic normal vibrational modes analyzed through CNM provides compelling evidence for the utilization of local mode properties to assess the intrinsic strength of metal–ligand bonding in JT systems.

Table 1 presents a comprehensive collection of the properties associated with the local mode parameters of the same metal–ligand and ligand stretching vibrations featured on the CNM plots. The table comprises bond length (l), local mode force constant (\(k^a\)), local mode frequency (\(\omega ^a\)), electron density [\(\rho (\textbf{r}_b)\)], energy density [\(H(\textbf{r}_b)\)], and delocalization index (DI) for the complex ions (test examples 1–4) and the isolated ligands. See Supplementary Information (SI) for Cartesian coordinates of the optimized geometries and the QTAIM atomic charges.

These values were obtained through calculations performed at the PBE0/Def2-TZVP level of theory. On the one hand, the expected equivalence of the properties associated with the local mode parameters of the same metal–ligand and ligand stretching vibrations in both reference complex ions, \(\mathrm {[Cr{(OH_2)}_6]^{3+} (T_h)}\) and \(\mathrm {[Fe{(CN)}_6]^{4-} (O_h)}\), is observed and depicted in Table 1. On the other hand, the breaking of symmetry leads to differentiation between the properties describing the two axial and four equatorial metal–ligand local mode parameters, as observed in the complex ions where JT distortions occur, \(\mathrm {[Cr{(OH_2)}_6]^{2+} (D_{2h})}\) and \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\).

Table 1 shows that the axial local mode parameters, Cr1O4 and Cr1O6, provide insights into the octahedral elongated geometry of \(\mathrm {[Cr{(OH_2)}_6]^{2+} (D_{2h})}\). JT distortion in this compound results in the observed substantial Cr1O4 and Cr1O6 bond length elongation followed by a general reduction in the local force constant, local frequency, electron density, and delocalization index of metal–ligand bond parameters, with the Cr1O4 and Cr1O6 parameters displaying the lowest values. Interestingly, the Cremer-Kraka criterion indicates an increase in axial covalent character and a striking shift to electrostatic character for the equatorial bonds after distortion. On the contrary, due to JT distortion, the local force constant, local frequency, electron density, and delocalization index of water ligands experience a considerable overall increase, with a slightly enhanced covalent character for all the O–H bonds. The results for the isolated form of H\(_2\)O are shown for comparison with the bonded forms, where an increase in bond length and a decrease in the absolute value of all the other computed properties are observed in both cases.

Table 1 Selected local mode parameters (param.) followed by the corresponding values of bond length (l), local mode force constant (\(k^a\)), local mode frequency (\(\omega ^a\)), electron density [\(\rho (\textbf{r}_b)\)], energy density [\(H(\textbf{r}_b)\)], and delocalization index (DI), for the complex ions and the isolated ligands, calculated at the PBE0/Def2-TZVP level of theory
Full size table
Fig. 4
figure 4

ACS plots relating the frequencies of a complete and non-redundant set of local modes (left) with those of normal modes (right) for \(\mathrm {[Cr{(OH_2)}_6]^{3+} (T_h)}\) (a), and \(\mathrm {[Cr{(OH_2)}_6]^{2+} (D_{2h})}\) (b). The pairwise breakdown of the corresponding separate metal–ligand frequency range is provided on the bottom (c, d). Results obtained at the PBE0/Def2-TZVP level of theory

Full size image

In the case of \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\), Table 1 also shows that Fe1C4 and Fe1C6 local stretching modes serve as indicators of the octahedral compressed geometry. Despite the subtle compression, there is a significant overall increase in the local force constant, local frequency, and electron density for both metal–ligand and ligand \(\mathrm {CN^-}\) bond parameters, with Fe1C4 and Fe1C6 parameters showing slightly smaller values compared to the equatorial local mode parameters. Additionally, the energy density values demonstrate that both metal–ligand and ligand bonds become notably more covalent upon distortion, with no shift in bond nature. It has to be noted that the bond length is not always a qualified bond strength descriptor. Numerous cases have been reported in the literature illustrating that a shorter bond is not always a stronger bond [72,73,74,75,76,77,78]. This also applies to the metal–ligand bonds in \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\), despite the delocalization index indicating otherwise.

Comparing the isolated and bonded forms of the \(\mathrm {CN^-}\) ligand reveals that in \(\mathrm {[Fe{(CN)}_6]^{4-} (O_h)}\) a slight increase in bond length and a decrease in the local force constant, local frequency, electron density, the absolute value of energy density, and delocalization index occur as expected due to the \(\pi \) back-donation phenomenon discussed in the literature [33] (see Table 1 for details). When compared to \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\), both the electron and energy density change slightly, especially for the axial ligands, with a slight decrease in bond length. Still, there is an increase in both the local force constant and local frequency upon coordination with Fe(III). The weaker back-donation effect found in \(\mathrm {[Fe{(CN)}_6]^{3-}}\) can account for this behavior, as it implies a reduced transfer of electron density from the metal center to the ligands, affecting the CN bond characteristics. The difference in vibrational properties between \(\mathrm {[Fe{(CN)}_6]^{4-} (O_h)}\) and \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\) is in accordance with the literature [33, 79].

Fig. 5
figure 5

ACS plots relating the frequencies of a complete and non-redundant set of local modes (left) with those of normal modes (right) for \(\mathrm {[Fe{(CN)}_6]^{4-} (O_h)}\) (a), and \(\mathrm {[Fe{(CN)}_6]^{3-} (D_{4h})}\) (b). The pairwise breakdown of the corresponding separate metal–ligand frequency range is provided on the bottom (c, d). Results obtained at the PBE0/Def2-TZVP level of theory

Full size image

Furthermore, Figs. 4 and 5 present the ACS plots for test examples 1–4. In these plots, the parameter \(\lambda \) distinguishes between local modes (\(\lambda =0\)) and normal modes (\(\lambda =1\)). Both \(\mathrm {[Cr{(OH_2)}_6]^{3+}}\) and \(\mathrm {[Cr{(OH_2)}_6]^{2+}}\) exhibit 51 normal vibrational modes, which are adiabatically connected to 6 groups of local vibrational modes, as shown in Fig. 4. Similarly, \(\mathrm {[Fe{(CN)}_6]^{4-}}\) and \(\mathrm {[Fe{(CN)}_6]^{3-}}\) possess 33 normal vibrational modes, which are adiabatically connected to 5 groups of local vibrational modes, as depicted in Fig. 5.

From Fig. 4(a) and (c), it is evident that the \(\mathrm {[Cr{(OH_2)}_6]^{3+}}\) complex ion with the \(\mathrm {T_h}\) point group symmetry exhibits equivalence among the \(\mathrm {OH_2}\) ligands, leading to the expected result. Consequently, several degenerate local vibrational modes can be observed, which can be categorized as follows: six degenerate Cr–O bond stretching modes and twelve degenerate O–H bond stretching modes, as highlighted in Table 1, in addition to six degenerate local O–Cr–O modes, twelve degenerate local H–O–Cr modes, three degenerate local O–Cr–O–O modes, and twelve degenerate local H–O–Cr–O modes. It can be observed from Table 1 and Fig. 4(c) and (d) that the expected sixfold degeneracy of the Cr–O bond local vibrational modes in \(\mathrm {[Cr{(OH_2)}_6]^{3+}}\) is broken, forming three groups of twofold degeneracy in \(\mathrm {[Cr{(OH_2)}_6]^{2+}}\). Figure 4(b) and (d) clearly illustrate the symmetry breaking due to the JT elongation captured by ACS. The axial local vibrational modes, denoted by Cr1O6 and Cr1O4, occur as the lower-frequency twofold degenerate group, whereas the four equatorial local vibrational modes split into two very closely separated groups due to the equilibrium geometry being a \(\mathrm {D_{2h}}\) minimum. Of particular importance is the local vibrational mode evidence of the strength of the JT distortion, seen by the pronounced frequency separation between the axial and equatorial groups of local vibrational modes in Fig. 4(d). This finding agrees with the marked difference in bond length tabulated in Table 1.

It is worth mentioning that crystallographic data and EPR studies have indicated that \(\mathrm {[Cr{(OH_2)}_6]^{2+}}\) preferably exhibits an octahedral elongated structure [30, 31]. Furthermore, it is well-established [6, 29, 31] that octahedral JT elongated structures can be understood in terms of the d-electron stabilization energy, where a greater electrostatic repulsion is observed between the metal ion and the two ligands along the z-axis, leading to their stabilizing elongation. It is also remarkable that the local mode force constants (\(k^a\)) are smaller for Cr1O4 and Cr1O6 bonds, indicating that they are weaker than the Cr–O equatorial bonds. This trend in \(k^a\) can be attributed to the elongation of the axial Cr–O bonds, which serves to decrease the repulsion between the d-orbitals along the z-axis and the ligands. By elongating the axial bonds, the d-electron distribution around the chromium ion is modified, leading to reduced energy levels and increased stability. This effect is in line with the concept of JT elongation, where the distortion of the octahedral geometry results in lower energy states for the system. Moreover, the observed decrease in \(\rho (\textbf{r}_b)\) upon JT distortion further supports the weakening of the axial bonds. Additionally, the negative value of \(H(\textbf{r}_b)\) for Cr–O axial bonds reveals that the electron density at the bond critical point is still significant, indicating a considerable bonding interaction despite the weakening. Overall, the findings reported in this study align with the known principles of JT elongation in octahedral complexes and provide valuable insights into the structural and electronic properties of \(\mathrm {[Cr{(OH_2)}_6]^{2+}}\).

Figure 5(a) and (c) provide clear evidence that the \(\mathrm {[Fe{(CN)}_6]^{4-}}\) complex ion, exhibiting \(\mathrm {O_h}\) point group symmetry, demonstrates equivalence among the \(\textrm{CN}\) ligands, as expected. Consequently, several degenerate local vibrational modes can be observed, classified as follows: six degenerate Fe–C bond stretching modes and six degenerate C–N bond stretching modes, as shown in Table 1, in addition to six degenerate local C–Fe–C modes, twelve degenerate local N–C–Fe modes, and three degenerate local C–Fe–C–C modes. From Table 1 and Fig. 5(c) and (d), it can be observed that the sixfold degeneracy of the Fe–C bond local vibrational modes in \(\mathrm {[Fe{(CN)}_6]^{4-}}\) is broken in \(\mathrm {[Fe{(CN)}_6]^{3-}}\), resulting in two groups of twofold and fourfold degeneracy. The symmetry breaking due to the JT compression captured by ACS is clearly illustrated in Fig. 5(b) and (d).

The axial local vibrational modes, denoted by Fe1C6 and Fe1C4, occur as the lower-frequency twofold degenerate group, while the four equatorial local vibrational modes occur as the fourfold degenerate group, in line with the equilibrium geometry of \(\mathrm {D_{4h}}\) symmetry. The local vibrational mode evidence of the weakness of the JT distortion in \(\mathrm {[Fe{(CN)}_6]^{3-}}\), as mentioned in the literature [33], can be observed by the small frequency separation between the axial and equatorial groups of local vibrational modes in Fig. 5(d). This finding is in agreement with the minor difference in bond length tabulated in Table 1.

As an important development, the innovative method of probing JT distortions through ACS and CNM, as demonstrated in this study, is poised to undergo further expansion to encompass the vibrational analysis of coordination compounds with diverse sizes and complexities. This advancement opens up a wide range of potential applications that can be tailored to investigate different metal ions, ligands, chemical environments, and coordination spheres, thereby enriching our understanding of coordination chemistry. Work is in progress to extend the protocol described herein to post-Hartree-Fock methods.

Conclusions

We developed in this work a new way for assessing JT distortions via a combination of two important features of LVM, namely the characterization of normal mode (CNM) procedure and an adiabatic connection scheme (ACS) connecting local and normal vibrational modes. The viability of this new CNM/ACS protocol was successfully demonstrated through the examination of classical coordination compounds at the PBE0/Def2-TZVP level of theory. Our CNM/ACS results provide new insights into the changes in symmetry resulting from JT distortions, which sets the stage for applying this approach to a diverse array of systems with varying size and complexity. More specifically, the following results emerged from our study, unravelling (i) distinct pairwise patterns of ACS and CNM variations across different oxidation states, and (ii) subtle differences between systems 1/2, representing the category with moderate JT distortions accompanied with bond lengthening, and 3/4, representing the category with slight JT distortions accompanied with bond shortening:

  • CNM results illustrated that the characteristic normal vibrational modes exhibit overall delocalization for both scenarios 1/2 and 3/4, confirming what we generally observe for systems with high symmetry. There is a slight increase in localization in the transition from 1 to 2 for both Cr–O and O–H stretching vibrations. In contrast, there is a reverse trend for Fe–C stretching vibrations in the transition from 3 to 4.

  • The ACS analysis clearly reflected the symmetry breaking via the differentiation between axial and equatorial groups following redox changes for the pairs 1/2 and 3/4 as a consequence of the JT distortions. The ACS profile of 4 exhibits a slight frequency separation between the axial and equatorial groups, whereas in the case of 2, a more pronounced distinction is observed. This is in accordance with the fact that 2 undergoes a moderate JT distortion compared to 4, which experiences a minor JT distortion. These findings suggest the use of ACS patterns as sensitive indicators of the extent of JT distortions, based on the frequency separation between axial and equatorial groups displayed in the ACS profile.

  • The Cremer-Kraka criterion revealed an enhancement in axial covalent nature and a noticeable shift toward electrostatic character for the equatorial bonds following distortion from 1 to 2. Moving from 3 to 4, the energy density data indicate a significant increase in covalence for both metal–ligand and ligand bonds upon distortion, while their bond nature remains consistent. Comparison with the isolated CN\(^-\) molecule confirms the distinct \(\pi \) back-donation disparity between 3 and 4, as documented in the literature [32,33,34,35].

In summary, the ACS profile emerges as a powerful tool to investigate and confirm the subtle and less-explored JT distortion in 4 and similar systems, often more difficult to capture, and in this way providing an answer to the current debate on the existence of a JT distortion in this system. As demonstrated in this proof-of-concept study our new CNM/ACS protocol opens up new possibilities for future JT studies in coordination chemistry. Its adaptability to different metal ions, ligands, chemical environments, and coordination spheres can be utilized to gain a better understanding of coordination compounds. The versatility of this innovative approach allows for a broad range of applications, opening the avenue for further CNM/ACS JT explorations stretching into organic chemistry [80], biochemistry [81] and beyond [82, 83]. Work is in progress to extend the protocol described herein to post-Hartree-Fock methods.