The Jahn-Teller effect in mixed aqueous solution: the solvation of Cu²⁺ in 18.6% aqueous ammonia obtained from ab initio quantum mechanical charge field molecular dynamics

Simulation protocol

The QMCF-MD simulation was performed in a cubic simulation cell with a side length of 31.75 Å consisting of one Cu²⁺ ion, 815 ammonia and 184 water molecules. This setup corresponds to 18.6% aqueous ammonia at the experimental density of 0.927 g cm⁻³ [54]. The simulation was carried out in a canonical (NVT) ensemble under periodic boundary conditions and the minimum image convention was applied to all interactions. A second order Adams–Bashforth predictor-corrector algorithm was used to integrate the equation of motions with a time step of 0.5 fs. The temperature was set to 298.15 K corresponding to the experimental temperature employed to study 18.6% aqueous ammonia solutions [54]. To maintain constant temperature conditions along the simulation the Berendsen algorithm with a bath relaxation time of 0.1 ps was employed [55]. The Barker-Watts reaction field approach was used to account for the cutoff of long-range electrostatic interactions in the simulation [56].

The radius of the QM region was set to 6.5 Å, and a smoothing region with a thickness of 0.2 Å was applied to ensure continuous transitions of molecules between the QM and MM zones [57]. The calculation method of this simulation was Hartree-Fock, the respective basis sets applied to Cu²⁺ and all solvent molecules (H₂O, NH₃) were LANL2DZ-ECP and Dunning DZP, respectively. The flexible BJH-CF2 water model [58] was employed to describe all water interactions, enabling the description of explicit hydrogen motion. A flexible ammonia potential model [59] was used to describe the respective intra- and inter-molecular potential contributions. The interaction between NH₃ and H₂O in the MM region is described using a previously published pair potential [60]. The starting configuration for the QMCF MD simulation was obtained via classical MD employing published pair potentials for the Cu²⁺–O and Cu²⁺–N interactions [61]. First, the system was equilibrated for several hundred picoseconds using classical MD until a stationary temperature control was achieved after including heating and cooling conditions. After confirming stationary total energy and temperature, the system was re-equilibrated using QMCF MD for 4 ps to reach rapid equilibration. Next, data sampling was carried out for 20 ps. All quantum mechanical calculations in the QMCF MD simulation were executed using TURBOMOLE 5.9 [62], [63]. The VMD package was employed to visualize the trajectory data [64].

The NBO analysis was performed using Gaussian 09 [65], applying the respective molecular coordinates extracted from the QMCF MD simulation trajectory. The specific frame with the most stable geometry was employed as starting configuration for the NBO analysis, which was subsequently energy minimized at the HF level of theory using the LANL2DZ ECP basis set for Cu(II) and Dunning DZP for the solvent molecules. Finally, the NBO analysis was carried out, and the associated Wiberg bond index and stabilization energies were calculated.

Analysis

The structural properties of the Cu²⁺-solvation are analyzed via radial and angular distribution functions (RDFs, ADFs), as well as the respective coordination number distributions (CNDs). In addition, to investigate the exchange dynamics of water and ammonia ligands in the second solvation shell, the MRTs are calculated by counting the number of successful ligand exchange via the “direct method” [38]:

t_sim represents the simulation time, N_av is the average coordination number of ligands in the respective shell, and Nex0.5  is number of exchanges with a minimum ligand displacement time of 0.5 ps.

The frequencies of the octahedral vibrational modes are calculated via Fourier transform of the respective velocity autocorrelation functions (VACFs) C(t) as defined in Eq. 2.

with N and N_t representing the number of particles and time origins t_I, respectively, and v_J corresponds to a particular velocity components of particle J. By projecting the velocities onto the fundamental modes of the octahedral coordination complex, a decomposition of the vibrational characteristics according to the fundamental modes being A_1g, E_g, T_1u and T_2g is achieved [66].

A second-order perturbation theory analysis of the Fock matrix is performed to calculate the donor-acceptor interactions in the NBO basis [67]. The interactions are based on the loss of occupancy from the localized NBOs in the idealized Lewis structure into empty non-Lewis orbitals. For each donor NBO (i) and acceptor NBO (j), the stabilization energy associated with the delocalization i→j is estimated as [67].

with q_i being the donor orbital occupancy, ε_i and ε_j are the associated orbital eigenenergies, and F(i, j) is the respective off-diagonal NBO Fock matrix element.

Structural properties

The radial distribution functions (RDFs) of the Cu²⁺-ligand interaction (i.e. Cu–N,O) and the respective individual Cu–O and Cu–N contributions along with the respective running integration are depicted in Fig. 1. To improve the clarity, the Cu-ligand RDF is shown at an offset of +0.5. The details of the structural data are summarized in Table 1.

Fig. 1:

Radial distribution functions (RDFs) of Cu²⁺ in 18.6% aqueous ammonia. Solid lines show the RDFs, and the dashed lines indicate the associated integration numbers.

The Cu–N,O RDF shows a narrow first shell peak in the range of 1.97–3.17 Å (peak maximum at 2.19 Å) composed of six ligands followed by broad second shell peak ranging from 3.17 to 5.45 Å (peak maximum at 4.55 Å) being on average populated by 13.5 ligand molecules. Decomposition of the total RDF into the contributions arising from the Cu–N and Cu–O interaction shows that the first shell is exclusively occupied by ammonia molecules. In contrast, the second shell is on average populated by 8.8 H₂O molecules, whereas the average number of NH₃ ligands is only about 4.7. The Cu-O and Cu-N pair distribution also show, that the ammonia content tends to decrease upon increasing Cu-ligand distance, while at the same time the pair density of H₂O molecules is increasing.

Close inspection of the first shell peak demonstrates a pronounced tailing towards larger distances, which has been attributed to the influence of the Jahn-Teller (JT) distortion of the solvation complex [68], [69]. In order to obtain detailed insight into the impact of the JT effect and the associated ligand dynamics, the pair distances observed for the three trans-NN pairs present in the first solvation shell have been plotted against the simulation time (see Fig. 2a–d). It can be seen that the three NN-distances undergo structural changes between stationary distances being at approximately 4.3 and 5.0 Å at the picosecond time scale. Whenever the distance between ligands associated to an elongated NN-pair is shortened, this process is compensated by elongating the distance of another NN-pair. Thus, on average two NN-pairs are found at distances close to 4.3 Å, whereas the third pair adopts an elongated distance. The observed ligand separations of 4.3–5.0 Å correspond to stable ion-N distances amounting to approx. 2.15 and 2.50 Å, respectively. However, it can be seen that in many instances the instantaneous N–N distances are at an intermediate state in the range of 4.5–4.7 Å (corresponding to Cu–N distances of 2.25–2.35 Å) arising from the conformational changes of the solvation complex occurring on the picosecond scale. This may explain why the accurate determination of the axial and equatorial distance is particularly challenging. Decomposition of the first shell peak of the Cu–N distribution (see Fig. 2e) demonstrates that each of the three N-N pairs indeed show highly similar contributions to the overall Cu–N RDF.

Fig. 2:

(a–c) Time evolution of the instantaneous N–N distance of the three trans-pairs in the first solvation shell. The dashed lines mark the two stable N–N distances of 4.3 and 5.0 Å arising due to the Jahn-Teller distortion, corresponding to average Cu–N distances of 2.15 and 2.5 Å, respectively. (d) Snapshot depicting the N-atoms of the three trans-pairs in red, blue and green. Second shell NH₃ and H₂O molecules are highlighted in light green and light blue, respectively. (e) Decomposition of the first shell RDF peak into the contributions arising from the respective trans-pairs.

A comparison between the structural data obtained from this simulation study and data reported in the literature is presented in Table 2. The structural properties are consistent with the experimental study based on EXAFS [7] reporting the solvation of Cu²⁺ in aqueous ammonia as an octahedral structure with the presence of four equatorial and two elongated axial ligands. The remaining deviation in the ion-nitrogen distances can be attributed to the challenging nature of the dynamic Jahn-Teller effect on the one hand and the possibly influence of perchlorate as counter ions at the comparably high concentrations employed in the experimental studies [7]. The distance of the Cu²⁺–NH₃ interaction in aqueous ammonia is shorter than that observed in the aqueous case [6], [15], which was shown to be longer than in pure liquid ammonia solutions [7], [12]. These findings highlight the impact of the presence of water on the interaction between Cu²⁺ and NH₃ ligands in aqueous ammonia, which is in line with the hard/soft Lewis acids and bases (HSAB) classification [70] of water being a hard ligand compared to NH₃, which is typically attributed as being soft [71], [72].

In contrast to the zero RDF intensity separating the first and second shell peak in the Cu–N RDF, which is associated to the absence of first shell ligand exchange, the broad second shell peaks in the Cu–N and Cu–O RDFs imply that a higher solvent mobility linked to the exchange of ligands with the surrounding bulk phase is observed. Since the coordination numbers obtained from the integration of the radial distribution functions correspond to the ensemble average over the whole simulation trajectory, the respective variation has been investigated using coordination number distribution (CND) plots shown in Fig. 3. As expected, the coordination number in the first solvation shell is found to be exclusively six.

Fig. 3:

Coordination number distributions (CNDs) of the (a) first and (b, c) second solvation shell of Cu²⁺ in 18.6% aqueous ammonia obtained from the QMCF MD simulation.

In the second solvation shell, a broad distribution of coordination numbers of NH₃ ligands ranging from 2 to 8 with the highest probability of 33.0% being observed for coordination number 5. This broad distribution results due to the large distances between the ion and ligands that cause a high ligand mobility and rapid ligand exchange between the second shell and the bulk phase. The associated fluctuating time series of second shell coordination numbers shown in Fig. 4 point towards the occurrence of a large number of ligand exchange events. Despite the rapid fluctuation of the individual time series of water and ammonia, the respective sum displays a notably smaller variation over the simulation trajectory.

Fig. 4:

Time series of the coordination number of NH₃ and H₂O ligands in the second solvation shell along with the respective sum.

In order to obtain further details on the structural properties in the first solvation shell, the ligand-ion-ligand angular distribution function (ADF) was evaluated. Because only ammonia ligands are present within the first solvation shell, only the N–Cu²⁺–N angle had to be calculated (Fig. 5). Two peaks ranging from 69.5°–114.5° and 155.5°–179.5° with maxima at 90.5° and 174.5° are visible, which correspond to the octahedral geometry of the solvation complex. Thus, despite the dynamical character of the Jahn-Teller effect leading to strong variations of the instantaneous Cu²⁺–N distances, the average angle distribution analysed over the whole simulation trajectory corresponds to that of an ideal octahedral solvation complex.

Fig. 5:

Angular distribution function (ADF) of the N–Cu²⁺–N angle in the first solvation shell with the visualization of the [Cu(NH₃)₆]²⁺ solvation complex.

Ligand dynamics

The ligand dynamics in the solvation shells is evaluated employing the time evolution of all ligand visiting the first and second solvation layers of Cu²⁺ (see Fig. 6). This dynamical analysis confirms that only N-atoms are present in the first shell and no ligand exchange between first and second solvation shell is observed. In contrast, ligands in the second shell show a high degree of mobility and frequently undergo exchange reactions with the bulk phase. The mean residence time (MRT) of second shell ligands was determined via direct monitoring of ligand exchange events [38]. The resulting MRT values were determined as 2.2 ps and 2.6 ps for ammonia and water ligands, respectively. Thus, although ammonia forms a highly stable first solvation shell, the second shell MRT value is shorter than that of water, which can be attributed to the weak hydrogen bond properties of the NH₃ molecules compared to H₂O.

Fig. 6:

Time evolution of (a) the Cu²⁺–O and (b) the Cu²⁺–N distance of all ligand molecules registered in the first or second solvation layer. The black dashed line represent the boundary between first and second solvation shell, the blue dashed lines indicate the boundary between the second shell and the bulk phase.

The vibrational power spectra of the octahedral A_1g, E_g, T_1u and T_2g-modes of the [Cu(NH₃)₆]²⁺ complex obtained via Fourier transform of the projected velocity auto-correlation functions are depicted in Fig. 7. The A_1g-mode corresponding to the symmetric stretch vibration enables a direct analysis of the associated ion-ligand force constants. Two peaks at 200 and 385 cm⁻¹ are visible in the spectrum of the A_1g-mode, with corresponding average force constants of 27 N m⁻¹ and 100 N m⁻¹, representing the different bond strength of the equatorial and axial ion-ligand interactions present in the Jahn-Teller distorted complex. This finding agrees well with the force constants for Cu²⁺ in aqueous solution [15] reported as 63 and 121 N m⁻¹.

Fig. 7:

The vibrational power spectra of the octahedral A_1g, E_g, T_1u and T_2g-modes of the [Cu(NH₃)₆]²⁺ complex obtained via Fourier transform of the projected velocity auto-correlation functions (VACFs).

Natural bond orbital (NBO)

The application of the NBO analysis was employed as an additional characterization to evaluate the properties of the ion-ligand interaction and the associated stabilization energy. Recently, this analysis has been successfully employed to characterize the ion-ammonia interaction in studies of K⁺ and Rb⁺ in pure liquid ammonia [36], [41] as well as Cu⁺ in aqueous ammonia solution [73].

In this work, the NBO analysis is applied to the solvated complex [Cu(NH₃)₆]²⁺. The Cu²⁺–NH₃ binding is dominated by orbital interactions between the unoccupied valence non-bonding (LP*) of Cu²⁺ and (i) the lone-pair orbitals (LP) of the N-atoms as well as (ii) the bonding N–H orbitals (BD) of ammonia. The associated Wiberg bond indices (WBIs) and the stabilization energy of specific donor-acceptor NBOs are listed in Table 3. The values for the WBIs in the range from 0.15 to 0.23 imply that although the ion-ligand interaction is dominated by electrostatic contributions, a residual covalent contribution appears to be present in the Cu²⁺–NH₃ binding.

The nature of the Jahn-Teller distorted solvation complex is visible in the Wiberg bond indices as well. While the four equatorial ammonia ligands show consistent bond indices of 0.23, the associated values for the axial NH₃ pair only amount to 0.15 and 0.17, respectively. The same trend is visible in the respective stabilization energy estimated for individual orbital interactions. While values in the range from 230 to 260 kJ mol⁻¹ are observed for the four equatorial ligands, significantly lower stabilization energies of 128 and 159 kJ mol⁻¹ are obtained for the NH₃ molecules in axial position. As expected, the main contribution to the donor-acceptor interaction stems from the LP of the N-atoms acting as donor NBO, however, the contributions from the N–H bonding orbitals still are significant. This analysis demonstrates that the analysis of natural bond orbitals can also be employed to study complex chemical interactions in a Jahn-Teller distorted coordination complex.