Low-cost calculation and analysis of 2D IR spectra of model diiron trinitrosyl complexes in the NO stretch region with vibrational perturbation theory

Abstract

Two-dimensional infrared spectroscopy offers unique capabilities for probing vibrational coupling in complex metal–ligand systems. In this paper, we combine two-dimensional infrared spectroscopy with vibrational perturbation theory to investigate vibrational coupling in a diiron trinitrosyl complex across three stable redox states. Although these systems are challenging for electronic structure methods, we demonstrate that key features of experimental 2D IR spectra can be accurately reproduced using reduced-dimensional anharmonic calculations with a small harmonic frequency scaling. Analysis reveals that N–O stretching modes maintain high locality across all redox states, with coupling patterns that directly reflect variations in Fe–N bond strength. Using curvilinear coordinate analysis, we demonstrate these differences result from systematic changes in cubic anharmonic force constants rather than mode delocalization. Our results establish N–O stretches as sensitive probes of metal–ligand bonding strength, expanding the toolkit for studying biologically relevant nitrosyl complexes.

---

1. Introduction

Iron nitrosyl complexes play important roles in biological systems, particularly in the storage and transport of nitric oxide (NO).^1,2 The affinity of NO for sulfur in biosulfur molecules, such as glutathione or cysteine, and for iron as in iron–sulfur clusters leads to a panoply of biologically significant structures, several of which contain all three constituents: Fe, S, and NO.³ Many of these structures appear to be a factor in intricate processes of NO delivery to targets of biological or physiological significance, ultimately landing on significant sites such as on heme iron in blood pressure–regulating processes or on iron in a form of nitrile hydratase, (N₂S₂)Fe(NO) (where N₂S₂ is derived from the backbone of a cysteine–serine–cysteine tripeptide unit), a natural catalyst of industrial significance for the hydration of nitriles.^4,5 As the free neutral NO radical is highly reactive, carrier molecules are required to transport NO from its biosynthetic preparations; in fact, the Fe(NO)₂ unit is touted as “a working form of nitric oxide.”⁶ M–NO systems are less straightforward to analyze spectroscopically than other strong infrared absorbers, such as M–CO or M–CN, in that the NO ligand itself may exist in three redox levels: NO⁺, NO˙, or NO⁻.^7–9 Frontier molecular orbitals of the M–NO unit are thus dependent on the preferred redox level of the metal to which it is bound, as well as geometry. These orbitals are also highly dependent on the ligands within the first and even second coordination sphere, L_xM(NO). Unlike the metal carbonyls and metal cyanides, there are few homoleptic M(NO)_n derivatives with n > 2.

The focus of this paper is on further understanding the chemical physics pertaining to the oxidation states in a [FeFe]-hydrogenase model system, [(NO)Fe(N₂S₂)Fe(NO)₂]^+/0/− (N₂S₂ = N,N′-bis(2-mercaptoethyl)-1,4-diazamethylethane), shown in Fig. 1, which we refer to as [Fe₂(NO)₃] for the rest of the paper.^7,9 The [Fe₂(NO)₃] model system displays two reversible redox events in its cyclic voltammogram and can form two redox congeners upon chemical reduction, forming cationic, neutral, and anionic species while maintaining the two-sulfur-bridged Fe₂(NO)₃ scaffold.^7–9 This and related complexes have previously been studied using linear infrared (IR) spectroscopy.^7–9 These previous studies were able to determine through harmonic frequency calculations and ¹⁵N-labeling that the central vibrational peak arises from the mononitrosyl.

Fig. 1 Experimental X-ray diffraction (XRD) structures of the iron nitrosyl complex of interest with noted Fe–Fe distances and mononitrosyl Fe–N–O bond angles, with indications for states of charge and electron spin. The {Fe(NO)}^D indicates D number of electrons in the frontier orbitals of the system according to the Enemark–Feltham notation.¹⁰

While linear IR spectroscopy has provided initial insights into the vibrational states of these complexes, two-dimensional infrared (2D IR) spectroscopy can provide a direct report of the mixing and coupling that exists between the vibrational chromophores of a system from the measured cross-peaks.^11–15 2D IR has been used extensively to study the vibrational couplings and dynamics of metal carbonyls as they have served as important systems in the technique's development.^15–23 Metal carbonyl systems' prevalence to 2D IR spectroscopy stems from their strong transition dipole moments, narrow lineshapes, and absorbance in an uncongested portion of the infrared region.²⁰ This long history shows that 2D IR is a powerful tool for studying the structure and dynamics of transition metal complexes, but limited work has been done for other ligand stretches like nitrosyls.^24–27 Metal nitrosyl systems share common traits with metal carbonyls in their infrared absorption spectrum, posing an interesting system for nonlinear spectroscopic techniques to quantify and making 2D IR a promising tool for analyzing these Fe₂(NO)₃ [FeFe]-hydrogenase model systems.

There are multiple available methods for modeling 2D IR spectra from the appropriate third-order nonlinear response functions that vary in utility based on the system being studied. For systems with large numbers of weakly coupled vibrational modes and significant configurational changes on the timescale of the experiment, such as liquid water and proteins, Numerical Integration of the Schrödinger equation (NISE) theory is often used with time-dependent Hamiltonians from spectroscopic maps to capture system dynamics and solvatochromism.^12,28–31 For systems with smaller mode spaces and in the higher coupling regime, Redfield theory methods applied to a static geometry provide a well-suited approach for calculating 2D IR line shapes, but they require knowledge of the system's underlying energy levels. Multiple approaches exist for obtaining the necessary anharmonic vibrational frequencies and transition dipole moments, including direct experimental fits,^16,18 bilinearly coupled Morse oscillators combining experimental fits with anharmonicities from DFT,²² fourth-order Taylor series expansions of the potential energy,¹⁷ and the localized-mode vibrational self-consistent method with the localized-mode vibrational configuration interaction method.¹⁵ Here, we have chosen to use second-order vibrational perturbation theory (VPT2), which allows us to start at the harmonic normal mode picture and consider the anharmonicity to be a weak perturbation with the potential of the system expressed as a Taylor series. Corrections to the harmonic frequencies are expressed in terms of the cubic and quartic force constants, which are often calculated by finite difference on the Hessian matrix. Obtaining the necessary partial quartic force field for VPT2 requires only 2N + 1 Hessian calculations, where N is the number of modes for which anharmonic frequencies and transition moments are required.³² This relatively low computational cost makes it well suited to systems without already parametrized potential energy surfaces, and previous studies have shown its successful application to modeling 2D IR spectra of metal carbonyls.^19,20,23

In this study, we calculate and interpret the vibrational coupling between iron nitrosyl ligands in the N–O stretch region for the redox series of a diiron trinitrosyl complex with comparisons of the experimental and theoretical 2D IR spectra. We utilize a reduced-dimensional treatment of vibrational perturbation theory providing a low cost route to calculate the eigenstates of the anharmonic vibrational Hamiltonian. We continue our analysis in curvilinear internal coordinates and oblique coordinates, where we decompose the original normal mode coordinates into specified internal coordinates representative of the system. Finally, we quantify our comparisons with the experiment to relate systematic changes in the electronic structure of the system between spin and redox state that still preserve the overall geometry (i.e., common bond distances, angles, and dihedrals) with the vibrational transitions shown in the experiment. The complex of interest exists in varying electron spin and redox states, providing a useful series to probe with these techniques. Our analysis reveals that the strength of Fe–N interactions plays a crucial role in determining the 2D IR cross-peak patterns, with different electronic structure methods capturing these interactions with varying degrees of accuracy. These findings provide new insights into how metal–ligand bonding affects vibrational mode coupling in biologically relevant iron nitrosyl complexes.

2. Methods

2.1. Experimental methods

The complex was synthesized using previously described methods,⁷ with full details in Section S1A. The 2D IR spectra of the complex in THF were collected using a custom-built spectrometer described in prior work,³³ with full details in Section S1B.

2.2. Computational details

To model the 2D IR spectra of these complexes, we adapted Redfield theory methods for calculation of the lineshape to ultimately interpret anharmonic treatments as outlined throughout prior 2D IR studies.^11,13,16,18

All electronic structure calculations were performed using the Gaussian 16 suite at two levels of density functional theory.³⁴ First, the TPSSh functional^35,36 was used with the 6-311++G(d,p) basis set^37–39 on the non-iron atoms and the Watchers–Hay basis set^40–42 on the iron, under the 6-311++G(d,p) designation within Gaussian 16 (referred to as the TPSSh method for the rest of the text). Second, the BP86 functional^43,44 was used with the SDD pseudopotential⁴⁵ on the iron and the 6-311++G(d,p) basis set^37–39 for the carbon, hydrogen, oxygen, nitrogen and sulfur atoms (referred to as the BP86 method). The choices for the levels of theory were based on prior studies of this and related complexes.^7–9,46 Those readers interested in the electronic structure of this complex are directed towards previous work by Ghosh et al.⁷ Additional information on geometric differences between the experimental XRD structures and our calculated geometries can be found in Table S13.

2.3. Anharmonic vibrational calculations

In prior theoretical 2D IR studies, VPT2 has successfully been used to evaluate the necessary anharmonic corrections to the frequencies and transition moments.^19,20,23 Using VPT2 has the benefits of being computationally inexpensive and widely available in electronic structure packages. However, most implementations of VPT2 do not easily enable modification of the underlying Hamiltonian or access to transition moments for transitions from states other than the ground state, both of which are immediately available when constructing a local mode model. To overcome these limitations, we use a recently developed implementation of vibrational perturbation theory based on the automatic reduction of corrections to polynomial forms.⁴⁷ This allows us to both introduce physically motivated scaling factors to adjust the parameters obtained directly from electronic structure calculations and to directly evaluate all necessary transition moments.

As seen in Fig. S3, neither level of theory provides good agreement with the linear IR spectra. However, anharmonic effects are less sensitive to the level of electronic structure theory than the harmonic effects are.^48,49 We introduce a scaling factor to account for these discrepancies at the harmonic level, providing an effective harmonic frequency that causes the anharmonic fundamental frequencies to agree with the center frequencies in the experimental linear FTIR spectrum. To calculate this effective harmonic frequency, the anharmonicity returned from our VPT2 calculations was added to the FTIR center frequency. These effective harmonic frequencies were used as the starting point for scaled VPT2 calculations to calculate eigenstates of the anharmonic vibrational Hamiltonian. The scaling factors and relevant frequencies are given in Section S6 and Tables S7 (effective harmonic frequencies) and S8 (frequency scalings). Other than the effective harmonic frequencies, no other corrections were introduced to the vibrational calculations. In our VPT2 calculations, we consider the N–O stretching modes, the highest vibrational frequency, and low-frequency rocking motions. The choice of mode space was due to the presence of degeneracies within the unscaled frequency calculations.

2.4. 2D IR spectra calculation from reduced dimensional anharmonic calculations

For calculating the 2D IR spectra from the VPT2 results, the relevant perturbatively calculated transition frequencies and dipoles for the N–O stretch fundamentals, overtones, and combinations were extracted. The third-order molecular response was evaluated using the appropriate response functions to generate the rephasing and nonrephasing pathways. The reader is directed to the Section S2 for the response functions.

The spectra were simulated with a population time (t₂) of 10 ps to match experiment and a dephasing time (T₂) of 1 ps, which approximately aligns with the experimentally determined dephasing dynamics for other metal–nitrosyl complexes.^24–26 The spectra were obtained by taking the Fourier transform of the sum of the third-order nonlinear response functions found in Section S2. The formalism presented in this work can be used to generate the rephasing and the nonrephasing 2D IR spectra for any number of coupled oscillators.⁵⁰ The simulated 2D IR spectra in this work were calculated with perpendicular 〈XXZZ〉 polarization conditions to match the experiment. We define the double-sided Feynman diagrams and the analytical expressions for the response functions with 〈XXZZ〉 polarization in Section S2.

2.5. Curvilinear internal coordinate representation

Utilizing an internal coordinate representation is advantageous in exploring the physical origins of cross-peak intensity. Specifically, we identified nine internal coordinates that were most important for exploring coordinate mixing in the N–O stretch modes.

In Fig. 2, each of the individual stretches and bends is defined as [N–O]_1–3 for the local N–O stretches, [Fe–N–O]_1–3 for the Fe–N–O bends, and [Fe–N]_1–3 for the local Fe–N stretches, with the same order preserved for each set of three redox congeners. The ordering is consistent with the dominant contribution to the normal mode.

Fig. 2 Structural labels for internal coordinate and oblique coordinate decomposition.

To simplify the comparison to a local mode model, we added an oblique scaling transformation, first introduced by Zúñiga⁵¹ and generalized by one of the authors,⁵² that preserves internal coordinate identity while making the internal coordinates equivalent to normal modes up to a pure rotation. This has the added benefit of making the internal coordinate Wilson and matrices⁵³ equivalent, allowing the internal coordinate Hessian to act as an effective Hamiltonian. By using this coordinate choice, specific terms in the effective Hamiltonian may be correlated with physically interpretable couplings that drive the structure of the 2D IR spectrum. We also made comparisons between the internal coordinate Hessians to understand the coupling in the potential term.

3. Results and discussion

3.1. Analysis of normal vibrational mode locality

The normal modes found with both levels of theory have a relatively high degree of locality that was assigned in past work to specified portions of the complex.⁷ Here, we computed the percentages that each specific stretch and bend contributes to the normal mode. In Section S3, we report the decompositions in Tables S3 (cation), S4 (neutral), and S5 (anion); we briefly summarize the results in Table 1.

Table 1 Displacement contribution percentages to normal modes from dominant local N–O stretch and corresponding Fe–N stretch, where i denotes the specific numbering from Fig. 2 equivalent to the normal mode number in increasing frequency. Complete percentages can be found in Tables S3–S5

Electronic Structure	BP86 Method			TPSSh Method
Displacement	Normal Mode 1	Normal Mode 2	Normal Mode 3	Normal Mode 1	Normal Mode 2	Normal Mode 3
Cation [N–O]i	83.5	86.2	82.7	82.8	86.0	80.3
Neutral [N–O]i	64.0	68.3	81.5	80.2	82.6	84.1
Anion [N–O]i	67.1	73.0	76.9	84.3	87.0	84.2
Cation [Fe–N]i	11.5	11.4	11.7	10.2	9.5	9.4
Neutral [Fe–N]i	12.4	7.5	13.8	14.1	8.3	11.4
Anion [Fe–N]i	14.9	10.9	14.5	12.6	10.1	11.4

---

Due to the locality of the system, we mostly discuss the modes with local-mode-like terminology (e.g., [N–O]_i with normal mode i). For most normal mode vibrations, the second largest contributor to each normal mode is the Fe–N stretch adjacent to each N–O stretch, with the exception of longer-range delocalization between [N–O]₁ and [N–O]₂. The BP86 method predicts that the movement of the [N–O]₂ and the [N–O]₁ each contribute on a similar magnitude to the Fe–N stretching in normal mode 1 and normal mode 2, respectively, for select redox states. In one case, the additional N–O accounts for nearly 20% of the normal mode (see Tables in Section S3). Despite the locality of the modes, the use of a nonlinear spectroscopic technique provides a probe into the vibrational couplings between the nitrosyl stretches across redox states. The degree of N–O mode locality shows a systematic trend across redox states: the cation maintains the highest locality (>80% for all modes with both methods), while the neutral species exhibits the most delocalization, particularly with the BP86 method where cross-mode mixing reaches nearly 20%.

3.2. 2D IR spectra with unscaled harmonic frequencies

We utilized a recently developed implementation of vibrational perturbation theory to evaluate the relevant transition frequencies and transition dipole moments.⁴⁷ The 2D IR spectra, calculated from anharmonic frequencies and transition dipole moments evaluated at the minimum energy structures obtained with both TPSSh and BP86 methods, are reported in Fig. S4. In the theoretical spectra, there is a wider frequency range (Δω) between the mode 1 and mode 3 diagonal peaks when compared to the experiment. In Fig. S4B, for the anionic and the neutral species, there is significant overlap between signals associated with the first two vibrational modes, which is consistent with the linear IR signals in Fig. S3. For the neutral species, both the TPSSh method and BP86 method calculations produced spectra where the range of the frequencies is Δω = 230 cm⁻¹, while the experimental 2D IR has a much smaller overall distribution of frequencies, Δω = 110 cm⁻¹. There is significant overlap in the theoretical spectra such that each of the peaks cannot be distinguished. While the calculated frequencies of the fundamental transitions for the neutral and the anionic species differ appreciably relative to experiment, the overall structure of the 2D IR spectrum is in reasonable agreement for the cationic species, with the exception of a few cross-peaks.

3.3. Impact of normal mode frequency scaling

To address the disagreement in fundamental frequencies between calculation and experiment, we have followed a procedure laid out in prior literature.²² In this approach, we introduce a correction to the fundamental frequencies to agree with experiment, while leaving all anharmonicities and transition moments unaffected. As much of the fundamental physics and couplings in the system are encoded in cross-peak splittings, and as these splittings are only weakly perturbed by these corrections, this has the effect of allowing an improved visual comparison of the most important spectral features while avoiding overfitting of the spectrum. Similar effects have been seen previously in other anharmonic treatments.¹⁹

Our extension to this approach is to include these corrections directly in the harmonic frequencies, which are then used to evaluate the anharmonicities and transition moments. Specifically, we first evaluated the anharmonicity using the unscaled frequencies and added this to the center FTIR frequencies to obtain an effective harmonic frequency for use in subsequent calculations. As the frequency shifts are on the order of tens of wavenumbers out of frequencies on the order of 1600 cm⁻¹, the anharmonic corrections to the frequencies do not change significantly, and the final anharmonic frequencies agree very well with the experimental FTIR values as seen in Fig. 3. This results in a different scaling factor for each vibrational frequency between 0.96 and 1.10, shown in Table S8. It should be noted that this is predominantly for visual comparison to the experiment, allowing for greater insight, and none of the subsequent analysis changes qualitatively if these scaling factors are not used. We note additional information on the differences in cross-peak splittings between the scaled normal mode calculations and unscaled normal mode calculations in Tables S9–S11.

Fig. 3 FTIR spectra for the cationic, neutral, and anionic species with the unscaled (BP86), scaled (BP86), and experimental spectrum.⁷

3.4. 2D IR with scaled harmonic frequencies

In Fig. 4, we report the 2D IR spectra calculated with scaled harmonic frequencies. Compared to the unscaled results (Fig. S4), frequency scaling reduces the diagonal distribution for the anionic and neutral species, better matching the experimental spectral range, as expected. The cationic species shows improved frequency distribution but still lacks the experimentally observed cross-peaks between modes 1 and 2 (cross-peaks 1 and 3). The neutral species exhibits cross-peaks above and below the diagonal, including the prominent experimental cross-peak in the bottom right. The anionic species shows good diagonal agreement, with the BP86 method nearly matching experiment except for the bleach cross-peak between modes 2 and 3 (cross-peak 4). While there are unique scaling factors for each normal mode frequency, the choice of scaling factor does not affect the ability to interpret important features. The cross-peak splittings in Fig. 4 and Fig. S4 are very close, in most cases differing within half of a wavenumber, with the exception of cross-peaks 2 and 5 in the neutral species for the TPSSh method. Applying the scaling factors to each independent harmonic frequency does not affect the ability to interpret the spectra and is solely a simple mechanism for the reader to visually identify the splittings and crosspeaks presence or absence.

Fig. 4 A. Experimental 2D IR spectra, B. Theoretical 2D IR spectra at BP86/6-311++G(d,p) (C, H, O, N, S) and BP86/SDD (Fe) calculated with scaled harmonic frequencies, C. Theoretical 2D IR spectra at TPSSh/6-311++G(d,p) calculated with scaled harmonic frequencies. The numbering on the experimental 2D IR spectrum for the cation denotes the cross-peak numbering convention for analysis.

Even with the scaling factors, there are differences between the experimental 2D IR spectra (Fig. 4A) and the theoretical 2D IR spectra calculated with both the BP86 (Fig. 4B) and TPSSh (Fig. 4C) methods. From visual comparison, there is near excellent agreement for some of the cross-peaks and poor agreement for some of the others. In general, even when agreement on some of the cross-peaks present is poor, the general structure is preserved. Specifically, cross-peaks are not present for cross-peaks 1 and 3 for the cationic species with either method, as well as cross-peaks 1 and 3 for the anionic TPSSh spectrum. Ground state bleach cross-peaks are not present for cross-peaks 5 and 6 for the neutral BP86 method spectrum, cross-peaks 1, 3, and 4 for the neutral TPSSh spectrum, cross-peak 4 for the anionic BP86 spectrum, and cross-peaks 4 and 6 for the anionic TPSSh spectrum. To both explain the presence of cross-peaks, as well as the absence, we turn to key items that define the relevant vibrational transition frequencies: the potential derivatives and the magnitudes of the transition moments.

3.5. Analysis of cross-peak frequency differences from cubic potential derivatives and transition moment magnitudes

The calculated 2D IR spectra (Fig. 4) are qualitatively different between the BP86 and TPSSh methods, even after our experimentally informed scaling factors. This can be most directly seen through the presence or absence of cross-peaks. For example, for the anionic species, the BP86 method produces the cross-peak shown in the experiment between modes 1 and 2, while the TPSSh method fails to produce this cross-peak. To understand this difference, we turn to our VPT2 analysis. In the anionic species, where the three N–O stretches have similar frequencies, the cross-peak splitting is given by the difference between the frequency of the transition from the fundamental |010〉 to the combination state |110〉 and the frequency of the transition from ground |000〉 to the fundamental, |100〉. This splitting is given exactly in terms of the anharmonic matrix, with

However, for three stretches with similar frequencies, the dominant contribution to the splitting is given by

where f_i,j,k is the third derivative of the potential involving modes i, j, and k, f_i,j,k,l is the corresponding fourth derivative, and [small omega, Greek, macron] is the average stretch frequency.

For stretches that are well-approximated as Morse potentials, the fourth derivative contribution is expected to be approximately half the magnitude of the cubic contribution (with opposite sign).⁵⁴ Therefore, as we analyze the origins of cross-peak splittings, we focus entirely on the cubic derivatives. Moreover, as f_1,2,3 is much smaller than f_1,2,2 or f_1,1,2 and the denominator is much larger, we focus on f_1,1,2 and f_1,2,2. When generalized to modes i and j, the off-diagonal terms, f_i,i,j and f_i,j,j, have large contributions to cross-peak splittings. The cross-peak splitting in the three-mode system depends on the frequency differences between overtone and combination bands, which are determined by the cubic force constants (see Section S8 for detailed mathematical derivations and additional cross-peak analysis).

In Fig. 5, we report the cubic derivatives of the potential expressed in dimensionless normal mode coordinates. When examining the force constants in Fig. 5, the f_1,1,2 force constants for the anionic system at each method stand out: −209.289 cm⁻¹ with the BP86 method, and −37.791 cm⁻¹ for the TPSSh method. The f_1,2,2 force constants follow a similar pattern: −131.949 cm⁻¹ for the BP86 method, and −48.616 cm⁻¹ for the TPSSh method. These differences in the cubic force constants are factors for the frequency differences that give the splitting between the cross-peaks, specifically in cross-peaks 1 and 3 in Fig. 4.

Fig. 5 Magnitudes of cubic force constants for the three N–O normal mode coordinate stretches for the cationic, neutral, and anionic spin-redox states. Subscripts i, j, k indicate the normal modes in increasing order with the most dominant contribution for normal mode 1 being from [N–O]₁.

Having established the role of cubic force constants, we now examine the frequency differences and transition dipole magnitudes for specific cross-peaks between the 2D IR spectra in Fig. 4. Here, we focus on identifying signals with varying cross-peak intensity between each of the two electronic structure methods and map back to the transition frequency differences and transition dipole moment magnitudes (Section S9) with the cubic derivatives in Fig. 5.

For cross-peak 3 in Fig. 4 for the anionic species, this cross-peak is present with the BP86 method but not with the TPSSh method. Examining the frequency differences, there is roughly a 5.73 cm⁻¹ cross-anharmonicity for the BP86 method and a 0.35 cm⁻¹ cross-anharmonicity for the TPSSh method. When combined with the differences in the cubic derivatives noted above, these control the splitting in the 2D IR spectra. There is similar reasoning for cross-peak 1. For the BP86 method, there is roughly a 5.73 cm⁻¹ cross-anharmonicity while there is roughly a 0.35 cm⁻¹ cross-anharmonicity for the TPSSh. For cross-peak 3 for the neutral species, the cross-anharmonicity is 11.085 cm⁻¹ for the BP86 spectrum and 3.167 cm⁻¹ for the TPSSh spectrum. An additional item of note is the absence of cross-peak 3 for the cationic species in both the TPSSh method and BP86 method. The f_1,1,2 and f_1,2,2 force constants do not differ greatly between the two methods: −34.837 cm⁻¹, −6.835 cm⁻¹, −49.161 cm⁻¹, −49.339 cm⁻¹, respectively.

While the cross-peak splitting is directly correlated with the cubic derivatives (as discussed above), the presence or absence of cross-peaks in the theoretical 2D IR spectrum can also be attributed to the magnitude of the relevant transition dipole moments. For cross-peak 3 for the neutral species, there is an 11.085 cm⁻¹ splitting for the BP86 method, and a 3.17 cm⁻¹ splitting for the TPSSh method. Although there is the small splitting for the cross-peak in the TPSSh spectrum, there are differences in the transition moment magnitudes, 0.246 (ground state bleach) and 0.405 (excited state absorption). The stronger transition moment for the excited state absorption drives the more pronounced blue absorption peak. The subtitles in differences in the transition moment magnitudes can have strong effects in the presence of smaller cross-peak splitting.

3.6. Analysis of curvilinear internal coordinate Hessians

To gain further insights on the differences between the electronic structure methods and spin–redox states, we finish by returning to the vibrational couplings at the harmonic level, with the use of curvilinear internal coordinates.⁵⁵ Normal mode coordinates can be delocalized throughout the molecule. A curvilinear internal coordinate representation allows for the decomposition of the Cartesian normal modes into various stretches and bends in our user-defined internal coordinate system allowing us to separate any large-amplitude motion from the smaller vibrational motions of interest.^55,56 Thus, we can examine how the N–O stretches ([N–O]_1–3), Fe–N–O bends ([Fe–N–O]_1–3), Fe–N stretches ([Fe–N]_1–3), change as spin–redox state changes, without additional couplings to other components of the system (see Fig. 2 for labeling scheme). In Section S10, we report the upper 9 × 9 block of the curvilinear coordinate Hessian matrices. The entire Hessian can be found along with the other code and data for this paper in the Zenodo repository.

To better visualize the role each N–O stretch plays in the signal collected in the experiment, we report below the curvilinear internal coordinate bilinear couplings from the Hessian between each independent N–O stretch with additional N–O stretches, Fe–N–O bends, and Fe–N stretches for the BP86 method, as well as the difference between the BP86 and TPSSh methods to better evaluate the differences present.

In Fig. 6, we compare the couplings of the mononitrosyl unit N–O bond vector and the two N–O bond vectors in the dinitrosyl unit to each of the other key curvilinear motions. Here, we see that [N–O]₂ (mononitrosyl) couples more strongly to the [Fe–N–O]₂ bend in the mononitrosyl than [N–O]₁ or [N–O]₃ couple to their respective Fe–N–O bends. This is expected because this N–O stretch and Fe–N–O bend are on the mononitrosyl portion of this coordination complex. Additional differences are the couplings between each local N–O and Fe–N stretch. The [N–O]₁/[Fe–N]₁, and [N–O]₃/[Fe–N]₃ coupling terms (boldened matrix elements in Section S10) in the Hessians reflect the potential couplings between the local N–O stretch and local Fe–N stretch in each of the NO ligands in the dinitrosyl portion of the complex in curvilinear internal coordinates. The [N–O]₁/[Fe–N]₁ couplings show large variations between spin–redox state and electronic structure method (TPSSh followed by BP86 in each case): 9.34 cm⁻¹ and 82.11 cm⁻¹ for the cationic species, 92.31 cm⁻¹ and 71.01 cm⁻¹ for the neutral species, and 1.25 cm⁻¹ and 64.40 cm⁻¹ for the anionic species. These differences are evident in the bottom portions of the figure. The large differences also occur in the [N–O]₃/[Fe–N]₃ coupling terms (again, TPSSh followed by BP86): 16.34 cm⁻¹ and 86.43 cm⁻¹ for the cationic species, 22.23 cm⁻¹ and 78.02 cm⁻¹ for the neutral species, and −6.02 cm⁻¹ and 68.99 cm⁻¹ for the anionic species. The more significant mixing between these important coordinates could be indicative of the larger differences in the force constants in the normal mode coordinates in Fig. 5.

Fig. 6 Curvilinear internal coordinate bilinear couplings at the quadratic level for both electronic structure methods obtained from the internal coordinate Hessian.

3.7. Analysis of scaled internal coordinate hamiltonians with the oblique transform

To simplify the analysis of curvilinear internal coordinate contributions to the normal modes, we also scaled the set of internal coordinates via an oblique transformation, described above in Section 2.5. As previously noted, the Wilson and matrices⁵³ are equivalent after this scaling transformation and act as an effective harmonic Hamiltonian expressed in internal coordinates. We report the oblique coordinate model Hamiltonians with the two electronic structure methods for the cationic, neutral, and anionic species in Section S11.

In Fig. 7, we show the N–O and Fe–N stretching frequencies and the N–O/Fe–N bilinear coupling from the oblique coordinate Hamiltonian matrices.⁵² These effective Hamiltonians provide insight into the origin of the differences in mode mixing demonstrated in Table 1, which drive the differences in the cubic and quartic force constants responsible for the presence or absence of spectral features using the BP86 and TPSSh methods. As an example, for the anion species treated with the BP86 method, the lowest frequency normal mode is composed of approximately 80% motion of [N–O]₁ (see Fig. 2) and 15% motion of [N–O]₂, and the intermediate frequency normal mode is comprised of approximately 80% motion of [N–O]₂ and 12% motion of [N–O]₁. This leads directly to the appearance of cross-peaks 1 and 3 in Fig. 4, while for the TPSSh method, there is no appreciable mixing and hence no cross-peak.

Fig. 7 Oblique coordinate model Hamiltonian matrix elements for the N–O stretch, Fe–N stretch, and coupling between the N–O and Fe–N stretches for the BP86 level of theory and the difference between the BP86 and TPSSh matrix elements.

The effective Hamiltonians provide us a method for determining what physical features drive these important mode mixings. By either replacing the Fe–N and N–O local frequencies or off-diagonal couplings from the BP86 method with the corresponding terms from the TPSSh method and rediagonalizing, we can determine which terms are most important. When doing this, we find that replacing the off-diagonal elements has only a small effect, while replacing the local frequencies allows one to transform from the BP86 mixing percentages to the TPSSh percentages. This results from the fact that the BP86 N–O stretch frequencies are significantly lower than those at the TPSSh level for the anion system, while the Fe–N frequencies are significantly higher. When considering that the off-diagonal couplings vary by much smaller amounts, this provides a potential mechanism for the observed mixing. In particular, as the Fe–N interaction increases in strength, the corresponding N–O interaction decreases in strength, which leads to higher frequency Fe–N stretches and lower frequency N–O stretches. Therefore, the presence or absence of cross-peaks is a potential reporter of Fe–N interaction strengths.

3.8. Comparison of spectral features across redox states

The 2D IR cross-peak patterns reveal systematic changes across the three spin–redox states that reflect both the electronic structure variations and geometric changes as the system is oxidized or reduced. The cross-peak intensities increase substantially from the anionic to cationic species, with the theoretical spectra showing overall excellent agreement for the anionic and neutral species, but missing important cross-peaks between modes 1 and 2 and modes 2 and 3 for the cationic species (labeled cross-peaks 1 and 6).

To better illustrate how differences in the anharmonic couplings affect cross-peak splitting, we conducted an interpolation of the cubic and quartic derivatives between the values obtained at the BP86 and TPSSh levels (Fig. 8). Here, we calculated the anharmonic frequencies in the absence of degeneracies by interpolating the cubic and quartic derivatives starting from those obtained from the BP86 method to the TPSSh method. From the calculated frequencies of the specified vibrationally excited states, we calculated the cross-peak splitting. This interpolation used the scaled normal mode frequencies from the BP86 method, see Section S6.

Fig. 8 Theoretical cross-peak splitting as a function of interpolated cubic and quartic derivatives from the BP86 method to the TPSSh method.

We report cross-peak splittings as functions of an interpolation between cubic and quartic derivatives. The frequencies reflect nearly accurate cross-peak splittings, compared to the splittings in Fig. 4 with all of the pure BP86 information. As we move from the BP86 method higher-order potential derivatives to the TPSSh method, the cross-peak splitting reflects the theory cross-peak splitting with the exception of the cross-peak splitting between neutral cross-peaks 2 and 5. We attribute this to possible degeneracies created in our frequency scaling. The interpolation scheme allows us to envision what possible cross-peak splittings can be seen at some level of theory that could exist between the two methods. When comparing to the experiment, there are differences between the theoretical cross-anharmonicities by a factor of approximately 4 or 5. This is indicative of further delocalization between each vibrational mode in experiment.

To provide further assistance in understanding these discrepancies, we examined correlations between the experimental bond lengths and the computed vibrational frequencies that utilized the scaled internal coordinate Hamiltonians (Section S11), which are presented in Section S12. The isolated mononitrosyl Fe–N distance has an excellent inverse correlation with the local Fe–N stretch frequencies. As the Fe–N distance increases from 1.666 Å in the cationic species to 1.702 Å in the anionic species, the corresponding frequencies decrease systematically from 740 cm⁻¹ to 677 cm⁻¹. This inverse relationship is consistent with Badger's rule, and similar trends have been documented in the literature.^57–60

However, when examining the cationic species, significant discrepancies emerge for the N–O stretches. As an example, the experimental structures show nearly equivalent N–O distances (1.162 Å for the mononitrosyl unit, and 1.165 Å for both in the dinitrosyl unit), but the BP86 calculations predict three distinct distances: 1.170 Å, 1.180 Å, and 1.165 Å. Bond length changes of 0.001 Å can alter vibrational frequencies by at least 10 cm⁻¹, so these structural errors directly impact predicted cross-peak patterns. Our mode-specific scaling approach attempted to mitigate these effects, but such an approach cannot work perfectly for nearly degenerate states, unless all couplings are otherwise correct. When we artificially symmetrize the N–O frequencies in the scaled Hamiltonian, we find that the mixing between [N–O]₁ and [N–O]₂ increases to 8%, which would be sufficient to generate the missing cross-peaks in Fig. 4.

4. Conclusion

We have shown that it is possible to compute 2D IR spectra for a series of complex diiron trinitrosyl complexes via a low-cost, experimentally-informed approach that applies minimal scaling to correct the fundamental frequencies based on linear IR, then evaluates the anharmonic response using second-order vibrational perturbation theory with the corrected frequencies. The use of VPT2 provides a direct connection between anharmonic force constants and the presence or absence of cross-peaks in the 2D IR spectrum. From this, it is possible to infer specific possible physical mechanisms for the cross-peak appearance by using an effective reduced-dimensional Hamiltonian approach in curvilinear local mode coordinates. This approach can explain the magnitude of the anharmonic force constants. Specifically, we note that the presence of cross-peaks in the anionic and neutral redox states is a potential reporter of the strength of the Fe–N interaction, while for the cationic species the levels of electronic structure utilized in this work (based on prior literature results for IR spectroscopy) are unable to recover any appreciable cross-peak intensity in key regions.

The vibrational analysis can be connected to a detailed discussion of the structural differences between the experimental crystal structures and the minimum energy geometries at the different levels of theory. In addition, we provide a complete breakdown of the relevant bilinear couplings between the most important curvilinear internal coordinates of the system for describing the molecular vibrations and their connection to the anharmonic force constants for the relevant normal modes.

Beyond a demonstration of a simple, effective approach to the simulation of multidimensional spectra, these results highlight the ability of 2D IR spectra in the N–O stretching region to provide insights into underlying physical processes, complementing well-established 2D IR studies of C–O and C–N stretches. The connection between spectral features and specific Fe–N interactions has implications towards investigating different metal–ligand interactions in nitrosyl-containing systems where these bonds play crucial roles in biological function, as well as other coordination complexes. The combination of minimally adjusted simulation and detailed analysis of effective Hamiltonian models is extensible to any similar future studies. Finally, the compact analytic form for determining cross-peak splittings between modes of similar frequency introduced in this work may be inverted to provide potential insight into anharmonic couplings directly from experiment.

Conflicts of interest

There are no conflicts to declare.

Data availability

All relevant data and electronic supplementary information (SI) can be found at the linked Zenodo repository: https://doi.org/10.5281/zenodo.17065451. Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp03578c.

Acknowledgements

H. A. M. acknowledges support from the Department of Defense (DoD) through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program. C. R. B. acknowledges the Robert A. Welch Foundation (F-1891) and the Department of Energy (DE-SC0023221). A. F. M. and D. P. T. acknowledge the Robert A. Welch Foundation (A-2049). M. A. B. and D. P. T. acknowledge support from Research Corporation for Science Advancement, through the Cottrell Scholar Program (CS-CSA-2023-054). M. Q., S. S., and M. Y. D gratefully acknowledge the Robert A. Welch Foundation (A-092A). This work was partially supported by the Reconfigurable Materials Inspired by Nonlinear Neuron Dynamics (REMIND) Energy Frontier Research Center, funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DE-SC0023353 to M.Y.D for the synthesis and characterization of complexes conducted in the Department of Chemistry at Texas A&M University. Additional laboratory and facility support was provided by the National Science Foundation under Grant No. CHE-2102159 to M. Y. D. The Welch Foundation supported laboratory supplies under Grant No. A-0924 to M. Y. D. The authors also thank Dr. Nattamai Bhuvanesh for help in solving the X-ray crystal structures. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.

References

1. S. M. Andrabi, N. S. Sharma, A. Karan, S. M. S. Shahriar, B. Cordon, B. Ma and J. Xie, Adv. Sci., 2023, 10, 2303259.
2. D. C. Pectol, S. Khan, R. B. Chupik, M. Elsabahy, K. L. Wooley, M. Y. Darensbourg and S.-M. Lim, Mol. Pharm., 2019, 16, 3178–3187.
3. D. C. Pectol, S. Khan, M. Elsabahy, K. L. Wooley, S.-M. Lim and M. Y. Darensbourg, Inorg. Chem., 2020, 59, 16998–17008.
4. S. Sanfui, M. Jana, N. Small, D. J. Darensbourg and M. Y. Darensbourg, J. Biol. Inorg. Chem., 2025, 30, 161–168.
5. H. Yamada and M. Kobayashi, Biosci., Biotechnol., Biochem., 1996, 60, 1391–1400.
6. A. F. Vanin, Biophysics, 2024, 69, 937–954.
7. P. Ghosh, S. Ding, M. Quiroz, N. Bhuvanesh, C.-H. Hsieh, P. M. Palacios, B. S. Pierce, M. Y. Darensbourg and M. B. Hall, Chem. – Eur. J., 2018, 24, 16003–16008.
8. C.-H. Hsieh, S. Ding, Ö. F. Erdem, D. J. Crouthers, T. Liu, C. C. L. McCrory, W. Lubitz, C. V. Popescu, J. H. Reibenspies, M. B. Hall and M. Y. Darensbourg, Nat. Commun., 2014, 5, 3684.
9. S. Sanfui, M. Quiroz, J. Li, Y. Ha, F. Yang, J. Guo, N. Bhuvanesh, B. S. Pierce, P. B. Balbuena, P. A. Lindahl, M. B. Hall and M. Y. Darensbourg, Adv. Sci., 2025, e13976.
10. J. Enemark and R. Feltham, Coord. Chem. Rev., 1974, 13, 339–406.
11. S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford University Press, New York, 1995.
12. C. R. Baiz, B. Basiak, J. Bredenbeck, M. Cho, J.-H. Choi, S. A. Corcelli, A. G. Dijkstra, C.-J. Feng, S. Garrett-Roe, N.-H. Ge, M. W. D. Hanson-Heine, J. D. Hirst, T. L. C. Jansen, K. Kwac, K. J. Kubarych, C. H. Londergan, H. Maekawa, M. Reppert, S. Saito, S. Roy, J. L. Skinner, G. Stock, J. E. Straub, M. C. Thielges, K. Tominaga, A. Tokmakoff, H. Torii, L. Wang, L. J. Webb and M. T. Zanni, Chem. Rev., 2020, 120, 7152–7218.
13. P. Hamm and M. Zanni, Concepts and Methods of 2D Infrared Spectroscopy, Cambridge University Press, 2011.
14. M. J. Tucker, C. J. Mallon and M. Hassani, J. Phys. Chem. B, 2025, 129, 1439–1452.
15. J. Brüggemann, M. Wolter and C. R. Jacob, J. Chem. Phys., 2022, 157, 244107.
16. O. Golonzka, M. Khalil, N. Demirdöven and A. Tokmakoff, J. Chem. Phys., 2001, 115, 10814–10828.
17. A. M. Moran, J. Dreyer and S. Mukamel, J. Chem. Phys., 2003, 118, 1347–1355.
18. M. Khalil, N. Demirdöven and A. Tokmakoff, J. Chem. Phys., 2004, 121, 362–373.
19. C. R. Baiz, P. L. McRobbie, N. K. Preketes, K. J. Kubarych and E. Geva, J. Phys. Chem. A, 2009, 113, 9617–9623.
20. C. R. Baiz, P. L. McRobbie, J. M. Anna, E. Geva and K. J. Kubarych, Acc. Chem. Res., 2009, 42, 1395–1404.
21. C. R. Baiz, K. J. Kubarych and E. Geva, J. Phys. Chem. B, 2011, 115, 5322–5339.
22. C. R. Baiz, K. J. Kubarych, E. Geva and E. L. Sibert, J. Phys. Chem. A, 2011, 115, 5354–5363.
23. R. Fernández-Terán, J. Ruf and P. Hamm, Inorg. Chem., 2020, 59, 7721–7726.
24. J. F. Brookes, K. M. Slenkamp, M. S. Lynch and M. Khalil, J. Phys. Chem. A, 2013, 117, 6234–6243.
25. Y. Fujii, K. Aikawa, J. Tayama, M. Banno, K. Ohta and K. Tominaga, J. Chem. Phys., 2023, 158, 134510.
26. M. Feng, J. Zhao, P. Yu and J. Wang, J. Phys. Chem. B, 2018, 122, 9225–9235.
27. K. Adamczyk, M. Candelaresi, R. Kania, K. Robb, C. Bellota-Antón, G. M. Greetham, M. R. Pollard, M. Towrie, A. W. Parker, P. A. Hoskisson, N. P. Tucker and N. T. Hunt, Phys. Chem. Chem. Phys., 2012, 14, 7411–7419.
28. T. l C. Jansen and J. Knoester, J. Phys. Chem. B, 2006, 110, 22910–22916.
29. T. l C. Jansen, B. M. Auer, M. Yang and J. L. Skinner, J. Chem. Phys., 2010, 132, 224503.
30. C. Liang, M. Louhivuori, S. J. Marrink, T. L. C. Jansen and J. Knoester, J. Phys. Chem. Lett., 2013, 4, 448–452.
31. A. W. Smith, J. Lessing, Z. Ganim, C. S. Peng, A. Tokmakoff, S. Roy, T. L. C. Jansen and J. Knoester, J. Phys. Chem. B, 2010, 114, 10913–10924.
32. P. R. Franke, J. F. Stanton and G. E. Douberly, J. Phys. Chem. A, 2021, 125, 1301–1324.
33. S. C. Edington, A. Gonzalez, T. R. Middendorf, D. B. Halling, R. W. Aldrich and C. R. Baiz, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E3126–E3134.
34. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision B.01, 2016.
35. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
36. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
37. A. D. McLean and G. S. Chandler, J. Chem. Phys., 1980, 72, 5639–5648.
38. K. Raghavachari, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
39. T. Clark, J. Chandrasekhar, G. W. Spitznagel and P. V. R. Schleyer, J. Comput. Chem., 1983, 4, 294–301.
40. A. J. H. Wachters, J. Chem. Phys., 1970, 52, 1033–1036.
41. P. J. Hay, J. Chem. Phys., 1977, 66, 4377–4384.
42. K. Raghavachari and G. W. Trucks, J. Chem. Phys., 1989, 91, 1062–1065.
43. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
44. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822–8824.
45. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.
46. S. M. Brothers, M. Y. Darensbourg and M. B. Hall, Inorg. Chem., 2011, 50, 8532–8540.
47. M. A. Boyer and A. B. McCoy, PyVibPTn, A General Package for Vibrational Perturbation Theory, 2021,  DOI:10.5281/zenodo.5563090.
48. L. McCaslin and J. Stanton, Mol. Phys., 2013, 111, 1492–1496.
49. C. Puzzarini, M. Biczysko and V. Barone, J. Chem. Theory Comput., 2010, 6, 828–838.
50. M. Cho, Chem. Rev., 2008, 108, 1331–1418.
51. J. Zúñiga, A. Bastida and A. Requena, J. Phys. B: At., Mol. Opt. Phys., 2019, 52, 055101.
52. M. A. Boyer, I. Sibert and L. Edwin, J. Chem. Phys., 2023, 159, 234108.
53. E. Wilson, J. Decius and P. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, Dover Publications, 1980.
54. M. A. Boyer and A. B. McCoy, J. Chem. Phys., 2022, 156, 054107.
55. M. Mendolicchio, Theor. Chem. Acc., 2023, 142, 133.
56. A. B. McCoy and M. A. Boyer, J. Phys. Chem. A, 2022, 126, 7242–7249.
57. R. M. Badger, J. Chem. Phys., 1934, 2, 128–131.
58. M. A. Boyer, O. Marsalek, J. P. Heindel, T. E. Markland, A. B. McCoy and S. S. Xantheas, J. Phys. Chem. Lett., 2019, 10, 918–924.
59. M. M. Foreman, J. F. Stanton and J. M. Weber, J. Phys. Chem. A, 2023, 127, 9717–9722.
60. J. M. Finney and A. B. McCoy, J. Phys. Chem. A, 2024, 128, 868–879.

---