The ¹²⁵Te Chemical Shift of Diphenyl Ditelluride: Chasing Conformers over a Flat Energy Surface

Abstract

The interest in diphenyl ditelluride (Ph₂Te₂) is related to its strict analogy to diphenyl diselenide (Ph₂Se₂), whose capacity to reduce organic peroxides is largely exploited in catalysis and green chemistry. Since the latter is also a promising candidate as an antioxidant drug and mimic of the ubiquitous enzyme glutathione peroxidase (GPx), the use of organotellurides in medicinal chemistry is gaining importance, despite the fact that tellurium has no recognized biological role and its toxicity must be cautiously pondered. Both Ph₂Se₂ and Ph₂Te₂ exhibit significant conformational freedom due to the softness of the inter-chalcogen and carbon–chalcogen bonds, preventing the existence of a unique structure in solution. Therefore, the accurate calculation of the NMR chemical shifts of these flexible molecules is not trivial. In this study, a detailed structural analysis of Ph₂Te₂ is carried out using a computational approach combining classical molecular dynamics and relativistic density functional theory methods. The goal is to establish how structural changes affect the electronic structure of diphenyl ditelluride, particularly the ¹²⁵Te chemical shift.

1. Introduction

Diphenyl ditelluride (Ph₂Te₂) was first synthesized in 1894 [1,2,3]. Since then, the study on organotellurides has been rapidly increasing. In contrast to the historical notion that compounds of tellurium are difficult to handle and very sensitive to air and light, besides having a very pungent odor, many of these substances proved stable to air and light and were odorless [4,5]. Like organoselenides, organotellurides are currently widely used in organic synthesis to obtain useful intermediates, e.g., via transmetallation or cross-coupling reactions, and in the interconversion of functional groups [6,7,8]. In particular, Ph₂Te₂ proved to be efficient in catalytic oxidations of olefins [9]. In this reaction, the active catalyst is phenyl tellurol which forms upon Te–Te bond cleavage and can be oxidized to phenyltellurenic and phenyltellurinic acid. In fact, the chalcogen–chalcogen bond is soft and polarizable and can be easily broken [10], much more easily after initial oxidation, as reported for diphenyl diselenide (Ph₂Se₂) [11,12].

By analogy to Ph₂Se₂, Ph₂Te₂ has also been considered for its antioxidant potential, mimicking the enzymatic activity of glutathione peroxidase (GPx) [13,14]. Rocha and co-workers have reported evidence of thiol peroxidase activity of diphenyl ditelluride in mice [15,16]. The GPx-like mechanism of peroxide degradation was assessed in vitro with different thiols as reducing agents. Differently from what was found in the case of analogous diselenides—for which the active species is regarded to be a selenol/selenolate—in presence of tellurium the reaction seems to be associated with reversible changes in the chalcogen’s oxidation state [15]. Mechanistic differences in the oxidation step of GPx enzymatic cycle [17,18] have been recently elucidated in silico also by modeling the catalytic pocket of Se-GPx and of its Te-mutant [19].

Interesting studies about the toxicity of Ph₂Te₂ are available in the literature [20,21,22,23,24]. However, the toxicological data on Ph₂Te₂ and organotellurium compounds, in general, are not exhaustive and conflicting results were obtained, the majority of which lean towards the view that organotellurides have a higher toxicity compared to their selenium analogs [25].

Ph₂Te₂ is a chiroptical molecule, due to its peculiar structural features (Scheme 1). The dihedral angle Ψ has an equilibrium value of −89°, close to the value of H₂O₂; conversely, the absolute values of the dihedral angles Φ_1,2 may be both close to 0° or both close to 90°, defining the orientation of the phenyl rings in two minimum energy geometries, respectively. In the crystallographic structure [26], the values of the Φ_1,2 dihedral angles are 96° and 90°. This arrangement is found to be lower in energy, as revealed by accurate calculations at different levels of theory [10]. The conformation with Φ_1,2 ≈ 0 is another minimum, although lying very close in energy to the former one, as we will see below. Importantly, in solution the phenyl rings are almost free to rotate. The fast conformational motions of the phenyl rings were observed experimentally [27] as well as quantified through density functional theory (DFT) calculations [10,28]. The barrier in gas phase is about 5.5 kcal mol⁻¹ for the rotation about Ψ and 0.9 kcal mol⁻¹ for the rotation about Φ_1,2 at B3LYP/cc-pVTZ(-PP) level of theory [25].

Nowadays, computational protocols based on DFT are widely used for the investigation of structural [29] and spectroscopic properties of organic and organometallic molecules [30]. In particular, the prediction of NMR chemical shifts has been the subject of several studies, especially concerning ¹H and ¹³C chemical shifts of organic molecules [31,32,33], also for the case of natural substances including heavy atoms [34]. In the latter case, relativistic effects must be accounted for, because of the well-known HALA (heavy atom on light atom) effect [35,36]. Relativistic effects are even more important when the sought after NMR properties are those of the heavy atom itself [37], as is the case of ¹²⁵Te-NMR properties in organotellurides.

¹²⁵Te-NMR parameters have been the subject of computational works [38,39,40]. For example, Rusakova et al. recently investigated several organotellurides and organoditellurides using a full four-component DFT relativistic approach. They observed that a very good agreement between calculated and experimental data could be obtained after inclusion, besides relativistic effects, of solvent and vibrational contributions [41]. However, diphenyl ditelluride was not included in the set of studied compounds.

In a recent paper [42], Elder and Vargas-Baca reported an experimental investigation of the aggregation behavior in solution of some derivatives of diphenyl ditelluride. They used ¹²⁵Te-NMR spectroscopy, supported by relativistic DFT calculations. Based on the concentration dependence of the δ (¹²⁵Te), the authors proposed that aggregation occurs in low polarity solvents, e.g., hexane, toluene, and dichloromethane. However, the observed variations were of the order of ≈ 10 ppm/mol L⁻¹, that is varying the concentration from 10⁻² M to 10⁻¹ M, would roughly change the ¹²⁵Te chemical shift of the investigated compounds by at most 1 ppm. The effect of the temperature, in the investigated range, was comparable, of the order of 0.2 ppm/K. Nevertheless, the authors stated in their Conclusions that “the influence of concentration on the measured δ¹²⁵Te, however, is dwarfed by the effects of conformational changes or the dielectric properties of the medium” [42].

Therefore, in this work, we will analyze the issues related to the calculation of the chemical shift of ¹²⁵Te of diphenyl ditelluride; emphasis will be given to the dependence of δ(¹²⁵Te) on conformational effects and the difficulties related with an accurate modeling of the conformational energy landscape, this one being an essential requisite for a reliable prediction of NMR properties. By using purely quantum as well as combined classic and quantum chemistry methods, we discuss the factors affecting the results, that is the geometry and the electronic effects. Particularly, the analysis of the dependence of the paramagnetic component of the shielding constant on selected relevant geometrical parameters such as the Te–Te and Te–H bond lengths in a simpler model compound, provides general insight for the calculation of δ(¹²⁵Te) in organotellurides.

2. Results

Recent computational studies on Ph₂Te₂ and its selenium analog [10,28] have revealed the presence of two stable minima. In fact, we have located two conformers on the diphenyl ditelluride potential energy surface (PES), displaying very similar values for the C–Te–Te–C dihedral (Ψ, −86° and −89°, respectively) but differing values of the C–C–Te–Te dihedrals (Φ_1,2 = 100° and 10°, respectively). We refer to these two structures as the ‘open’ (Figure 1, left, with Φ_1,2 dihedrals of 100°) and the ‘closed’ (Figure 1, right, with Φ_1,2 dihedrals of 10°) conformer, due to the arrangement of the phenyl rings that resemble open and closed doors, respectively. They lie very close in energy as, at the ZORA-OPBE/TZ2P level of theory, the closed conformer is found only 0.5 kcal mol⁻¹ above the open one, which is the most stable, in agreement with the crystallographic structure [26]. Both fully optimized structures were used for the calculation of the ¹²⁵Te isotropic shielding constants (σ) (If not differently stated, a single value for the ¹²⁵Te shielding constant/chemical shift is provided for Ph₂Te₂, which is the average of the almost identical values computed for each Te nucleus), which resulted in 2095.4 and 2390.2 ppm for the open and closed configurations, respectively, calculated in chloroform at the COSMO-ZORA-OPBE/QZ4Pae level of theory on the ZORA-OPBE/TZ2P optimized geometries (this level of theory will be referred to as COSMO-ZORA-OPBE/QZ4Pae//ZORA-OPBE/TZ2P); they correspond to δ(¹²⁵Te) equal to 574.8 and 280.0 ppm, respectively, using Me₂Te as a reference. These computed values reveal that the phenyl orientation strongly affects the ¹²⁵Te chemical shift, as the difference amounts to 294.8 ppm. For a comparison with the experimental value of the δ(¹²⁵Te) of Ph₂Te₂, which is 420.8 ppm [43], the Boltzmann-weighted mean with an open to closed ratio of 70/30 was computed, i.e., 486.4 ppm. The inclusion of spin orbit coupling [41,44] leads δ(¹²⁵Te) of 476.8 and 184.5 ppm in the open and closed conformer, respectively (level of theory: COSMO-SO-ZORA-OPBE/QZ4Pae//ZORA-OPBE/TZ2P); the Boltzmann averaged δ(¹²⁵Te) is 389.1 ppm. In the case of the open and closed conformers, for which the values of the C–Te–Te–C dihedral angle are similar, spin orbit effects can be regarded as a constant correction (see Figure S2). However, the spin orbit effect is significant, i.e., about 100 ppm, and it is due to HAHA (heavy atom on heavy atom) effect analogous to the HALA effect discussed in the Introduction: since the molecule of interest has two Te atoms and the reference, Me₂Te only one, such effect on the shielding constant cannot be cancelled in the calculation of the chemical shift. Nonetheless, the most significant difficulty in the calculation of the δ(¹²⁵Te) of Ph₂Te₂ is intrinsic in the nature of the molecular PES, being related to the small energy difference between the conformers. We note, in fact, that even an error as small as 0.1 kcal mol⁻¹ in the calculation of the relative energy, affects the population ratio at room temperature by a factor of about 1.2. It is worth mentioning that calculations of the δ(¹²⁵Te)—with a full four component relativistic method, also taking into account solvation effects—were performed for a series of mono tellurides in a recent paper [41]. The authors report a mean absolute error of 23 ppm for a series of six compounds. Taking into account the presence of two Te nuclei and the conformational issues, the inclusion of scalar relativistic effects for our system is sufficient and is not the most important error source.

With the aim of understanding the relation between the structure and the value of δ(¹²⁵Te), a classical molecular dynamics (MD) simulation was employed to obtain a set of structures spanning the conformational space of Ph₂Te₂. The generalized AMBER force field (GAFF) parameters required to completely define our molecule were generated using a computational protocol recently devised by some of us (see Materials and Methods section for further details). We could not use the GAFF parameters of (Ph₂Te₂) reported in [28] because they were obtained using quantum mechanical calculations performed at a different level of theory, i.e., B3LYP/cc-pVTZ(-PP). As demonstrated in this work, different functionals lead to slightly different equilibrium bond lengths and angles, which have a strong impact on the calculations of the shielding constant of the ¹²⁵Te nucleus. Thus, it is of fundamental importance that the force field parameters are derived consistently.

A 500 ns MD simulation on a single Ph₂Te₂ molecule, solvated with chloroform in a cubic box with a surrounding buffer of 12 Å, was carried out and 1500 structures were extracted at intervals of 300 ps. A structural analysis of the collected snapshots showed values of the Ψ dihedral angle clustered around 90° or −90° during the whole simulation (Figure 2), which are very close to the value obtained from the single molecule quantum mechanical optimization (−86° and −89° for the open and closed conformer, respectively). In contrast, the values of Φ_1,2 dihedrals are more spread out, encompassing the whole 360° range. This is in agreement with the fact that Ψ and Φ have very different rotational barriers, estimated to be about 9 and 0.5 kcal mol⁻¹ in gas phase at ZORA-OPBE/TZ2P [10]. During the MD simulation, the Te–Te bond distance fluctuates around a mean value of 2.67 Å, which is close to the average of the corresponding bond lengths in the open (2.69 Å) and the closed conformer (2.66 Å) optimized at the ZORA-OPBE/TZ2P level, with a standard deviation of 0.07 Å.

For each sampled snapshot, the ¹²⁵Te isotropic shielding constants (σ) of both Te nuclei of the isolated Ph₂Te₂ molecule were calculated at the COSMO-ZORA-OPBE/QZ4Pae level of theory and averaged. The mean value of the calculated δ(¹²⁵Te) was 501.9 ppm with a standard deviation of 269.5 ppm, clearly denoting that the δ(¹²⁵Te) is highly sensitive to the geometrical structure. δ(¹²⁵Te) values fall between a minimum of −225.5 ppm and a maximum of 1503.4 ppm (see Figure S1), an interval covering a good part of the known ¹²⁵Te-NMR chemical shift range [40]. Interestingly, we have analyzed the ¹²⁵Te shielding constants, also extracting a decreasing number of snapshots and found that the result does not change when using 1000 snapshots and the deviation remains smaller than 10 ppm up to an average based on 250 snapshots. Since, as noted above, the spin-orbit correction is non-negligible, we also calculated the mean value of the shielding constant using the spin-orbit relativistic Hamiltonian: the average δ(¹²⁵Te) based on the previous 250 MD snapshots at the COSMO-SO-ZORA-OPBE/QZ4Pae level of theory is 403.8 ppm. This result is in very good agreement with the experimental value, being underestimated by only 17 ppm.

A benchmark study was conducted to assess the effect of the functional and of the basis set on the geometry of the conformers and thus on δ(¹²⁵Te) values. In a previous computational study by some of us [10], the focus was centered on the deviation of the computed Ph₂Te₂ geometry from the crystallographic structure, to assess the best method for geometry optimization. Therefore, the computations were carried out starting from the crystallographic conformer—i.e., the open one—which was identified correctly as a minimum using 25 different functionals. The conclusion was drawn that the choice of the functional is not so crucial to predict a fairly good minimum energy structure of Ph₂Te₂. Conversely, in the same study, differences among the methods arose when computing reaction energies. For homolytic Te–Te bond dissociation, only few functionals provide good agreement with experimental data—including OPBE, OLYP, and B3LYP—and among these ones only OLYP provides the correct dissociation energy trend when comparison is made between Ph₂Se₂ and Ph₂Te₂.

For our analysis, the structures of both conformers are needed to obtain δ(¹²⁵Te) values close to the experimental one. Therefore, a systematic study on both conformers was carried out (i) to locate them on the PES testing different functionals and basis sets; and (ii) to evaluate the structural and energy discrepancies among them when changing the computational method. The structures of the closed and open conformers were fully optimized using 23 different functionals. The relative single point energies and the significant geometrical parameters, i.e., the Te-Te bond length (R) and the values of the Ψ/Φ_1,2 dihedral angles are reported in

Table 1. Relative energies of Ph₂Te₂ conformers (kcal mol⁻¹), Te–Te bond lengths (R, Å) and amplitudes of the Ψ and Φ_1,2 dihedral angles (°) calculated with 23 different functionals combined with the TZ2P basis set; scalar ZORA effects were included in the calculations. Open conformer is taken as reference in the ΔE calculation.

Functional	R	Ψ	Φ1,2	R	Ψ	Φ1,2	ΔE
LDA	2.69	−79	81	2.65	−82	17	−2.1
BLYP	2.78	−87	97	2.74	−93	17	0.7
BP86	2.74	−86	98	2.70	−90	11	-a
OLYP	2.73	−87	98	2.70	−91	15	1.0
PBE	2.73	−85	93	2.70	−89	15	−0.4
OPBE	2.69	−87	100	2.66	−89	10	0.5
RPBE	2.75	−86	95	2.71	−90	15	0.5
revPBE	2.74	−86	97	2.71	−90	14	0.4
HTBS	2.71	−86	97	2.67	−89	10	0.1
PW91	2.73	−86	89	2.69	−88	13	−0.3
mPBE	2.74	−85	94	2.70	−89	15	−0.1
mPW	2.74	−86	98	2.70	−89	11	- a
TPSS	2.73	−85	95	2.69	−90	16	0.0
M06-L	2.76	−60	96	2.70	−83	19	−1.5
B3LYP	2.74	−90	93	2.71	−97	27	0.9
PBE0	2.70	−88	90	2.67	−93	21	0.1
OPBE0	2.67	−89	96	2.64	−92	15	0.6
mPW1PW	2.71	−90	91	2.68	−94	23	-a
M06	2.73	−79	94	2.70	−93	26	-a
TPSSH	2.71	−86	94	2.68	−91	19	−0.3
BLYP-D3(BJ)	2.77	−69	89	2.70	−90	25	0.1
BP86-D3(BJ)	2.74	−60	95	2.69	−82	17	−0.1
PBE-D3(BJ)	2.73	−77	87	2.69	−84	17	−0.6

^a The difference is smaller than 0.1 kcal mol⁻¹.

.

Importantly, both conformers are found on the PES employing all the selected functionals and the energy differences between the open and closed structures are within a 1 kcal mol⁻¹, except when LDA and M06-L are used, for which the open-closed energy difference amounts to −2.1 and −1.5 kcal mol⁻¹, respectively. The Te–Te bond length (experimental value: 2.71 Å [26]) results systematically longer in the open conformer with all the tested functionals. However, the differences between these bond lengths in the two conformers is of hundredths of Å and changes slightly when changing functional, as it goes from a minimum of 2.67 Å in the case of the open structure and of 2.64 Å in the closed conformer when using OPBE0 to a maximum of 2.78 and 2.74 Å for the open and closed conformers when using BLYP. A peculiar behavior is found for the dihedral angle Ψ when looking at the results obtained with the functionals empirically corrected for dispersion. In these cases, Ψ is found to be appreciably smaller in the open conformer, resulting in structures in which the phenyl rings are almost stacked. Nevertheless, this geometrical difference does not significantly affect the energies of the conformers which are found to be very close to each other.

Additionally, full geometry optimizations of both conformers were carried out with the ZORA-OPBE functional combined with basis sets of different size to check whether the basis set choice may affect the final geometry. This test has confirmed once again the flatness of the PES of Ph₂Te₂, as with the smallest basis sets (i.e., SZ and DZP of single and double-ζ quality, respectively) the closed conformer was never located, but all attempts ended up in the open conformation. Increasing the number of the functions (i.e., employing the TZ2P and QZ4P basis sets) resulted in a more accurate description of the PES of Ph₂Te₂, which correctly returned both optimized structures although with minimal relative energy differences.

These results indicate that the main issue to consider in order to calculate a δ(¹²⁵Te) close to the experimental value is rooted in the accurate description of the molecular geometry and of the PES of Ph₂Te₂ and all functionals perform quite similarly. While we cannot act on the nature of the PES, it is desirable to understand how the calculated δ(¹²⁵Te) is affected by changes in the relevant geometrical parameters, which are related not only to the different conformation, but more in general to even small differences on the optimized geometries provided by the different methods. In order to tackle this aspect, a decomposition analysis of the gas phase shielding constant was carried out (see the Materials and Methods section). The most important results come from the paramagnetic part of the shielding constant (σ_p) as it is the one that displays the greatest variation upon geometrical modification of the structure. Moreover, for simplicity and in order to have a much broader view on the importance of the electronic environment of ¹²⁵Te, the shielding constant analysis was first performed using a smaller model compound—i.e., PhTeH—focusing on the lengthening of the Te–H bond and on the rotation of the dihedral angle H–Te–C–C.

We found that the total ¹²⁵Te shielding constant in PhTeH (2646.7 ppm) is higher than those of both conformers of Ph₂Te_2., indicating that in the monochalcogenide the Te nucleus is more shielded. According to the scheme proposed by Schreckenbach et al. [45,46], we decomposed the isotropic shielding constant into its diamagnetic (σ_d) and paramagnetic contributions (σ_p), which amount to 5385.6 and −2738.9 ppm, respectively. σ_p can be expressed as the sum of two contributions, one arising from the coupling of pairs of occupied orbitals (σ_p^occ-occ) and the other from the coupling of an occupied and an empty orbital (σ_p^occ-virt) [46]. The interesting information comes from σ_p^occ-virt, since it is possible to discern which pairs of filled and empty molecular orbitals possess a strong magnetic coupling and thus contribute the most to the isotropic shielding tensor, that is to the chemical shift. This analysis showed that the sum of all the possible σ_p^occ-virt terms is slowly converging to the total σ_p. The decomposition of σ_p shows that the main term comes from the magnetic coupling of the HOMO–LUMO pair (about 32 % of the total σ_p). The orbital analysis reveals that the HOMO is an antibonding orbital composed mainly of the Te p_y orbital (Figure 3A, left), whereas the LUMO is also an antibonding orbital located in the xz plane resulting from a combination of the Te p_z and p_x orbitals, the p_x orbital on the carbon directly bonded to Te and the s orbital of the hydrogen atom linked to Te (Figure 3A, right). This composition is favorable, since after application of the magnetic coupling operator, the HOMO orbital is rotated along the z- or x- axis and can overlap substantially with the LUMO (see Materials and Methods).

Upon stretching the Te–H bond, σ(¹²⁵Te) decreases when going from the equilibrium value (1.66 Å) to 2.66 Å, i.e., the value changes from 2647.7 ppm to −2455.6 ppm (Figure 4A). This last value, when decomposed, results in a σ_d of 5389.6 ppm and a σ_p of −7843.2 ppm. Comparing these values to those of the equilibrium structure, it is clear that the paramagnetic part is responsible for the decrease of the σ(¹²⁵Te). The decomposition of σ_p confirmed that the stretching of the Te-H bond enhances those terms arising from the occupied orbitals and the LUMO, as the energy of this anti-bonding orbital decreases with increasing Te–H distance (Figure 4A). Focusing on the largest contribution to σ_p—i.e., the one arising from the HOMO–LUMO pair (Figure 3B), the term σ_p^occ-virt changes from −891.4 ppm in the relaxed structure to −4998.8 ppm in the deformed geometry, while the energy gap between the two molecular orbitals shrinks from 3.49 to 1.10 eV.

The profile of the σ(¹²⁵Te) as a function of the H–Te–C–C dihedral shows the presence of three minima at 0°, 92.6°, and 180° values of the dihedral angle (Figure 4B). The two structures at 0° and 180° are equivalent by symmetry and therefore have almost identical σ(¹²⁵Te)—i.e., 2646.4 and 2645.6 ppm—whereas the relative minimum at 92.6° has a ¹²⁵Te shielding constant value of 2673.7 ppm. The analysis of the σ_p^occ-virt showed that when the dihedral angle H–Te–C–C is 0° or 180°, the empty molecular orbital mostly involved in the largest terms is the LUMO, which has the correct symmetry to allow favorable overlap, after rotation, with HOMO and HOMO-3 orbitals. As the dihedral value increases, the overlap decreases, resulting in a smaller—i.e., less negative—σ_p. The relative minimum at 92.6° arises thanks to the presence of unfilled orbitals that have the correct symmetry to overlap with the HOMO after rotation around the Te–C bond, namely LUMO + 2 and LUMO + 3 (see Supplementary Materials Table S3). However, these unfilled orbitals lie at a higher energy than LUMO and the increased energy gap implies an overall less negative σ_p, leading to a higher σ.

For Ph₂Te₂, we decomposed the ¹²⁵Te σ of both conformers into its diamagnetic and paramagnetic contributions. The results confirmed that σ_d changes only slightly for the two conformers—i.e., 5387.2 ppm for the closed conformer and 5388.5 ppm for the open one—and is almost identical to that computed for PhTeH, i.e., 5385.6 ppm. On the other hand, σ_p is appreciably different in the two conformers—i.e., −3495.1 ppm for the open one and −2990.7 ppm for the closed one. The decomposition of σ_p showed that the sum of all the σ_p^occ-virt is slowly converging as for PhTeH. In addition, the multitude of terms makes it difficult to outline distinctive features that differentiate the chemical shifts. Therefore, we focused on the most relevant contributions. Particularly, we first addressed the σ_p^occ-virt ascribed to the HOMO–LUMO pair, which has the smaller energy gap (3.05 eV in the closed conformer and 2.72 eV in the open conformer) and should therefore provide a substantial contribution to σ_p. We found that in both conformers these terms are similar and represent approximatively 10% of the total σ_p, i.e., −327.9 ppm (first Te nucleus) and −324.9 ppm (second Te nucleus) in the closed conformer, and −300.7 ppm (identical in each Te nucleus) in the open one. The HOMOs of both molecules display strong and either p_x or p_y atomic orbital character on the Te nuclei (Figure 5A). Conversely, the LUMOs have net Te p_z atomic character. Either B_y or B_x field components will rotate the p_x or p_y to p_z, thus leading to a non-null overlap integral with the LUMO. The component of the applied magnetic field is important for the evaluation of the integral in Equation (2) (see Materials and Methods): due to the mutually orthogonal nature of the p orbitals on the Te nuclei in the HOMO, the use of B_y or B_x determines a strong magnetic coupling only on one of the two chalcogens.

The largest contribution to σ_p in the closed conformer comes from the HOMO-6–LUMO orbital pair (−994.7 ppm (first Te nucleus) and −992.4 ppm (second Te nucleus)), whereas for the open structure the orbital pair HOMO-1–LUMO + 6 gives the most substantial contribution (−707.4 ppm (first Te nucleus) and −707.2 ppm (second Te nucleus)). These orbital pairs are still significantly localized on Te nuclei with p atomic character (see Figure 5), but the origin of the strong magnetic coupling, despite their similar energy gaps, is different. In the closed conformer, the chalcogen nuclei have mixed p_z and p_x/p_y character in the HOMO-6, these latter atomic contributions suitable for overlap with the Te p_z based LUMO. Conversely, in the open conformer, both HOMO-1 and LUMO + 6 in the open conformer have p_x (first Te nucleus), p_y (second Te nucleus) and p_y (first Te nucleus), p_x (second Te nucleus) character, so that magnetic coupling occurs simultaneously at both Te nuclei.

The deformation of the Te–Te bond in Ph₂Te₂ has a strong effect on the σ(¹²⁵Te). From the values of 2095.4 and 2390.2 ppm calculated for the two minima, we computed a value of 2515.5 ppm when the bond is compressed to 2.60 Å and of 1329.2 ppm when the bond length is increased from the equilibrium geometries to 2.90 Å. The σ(¹²⁵Te) profile along the stretching of the Te–Te bond shows an almost linear decrease with the bond length (Figure 6A). Importantly, starting from the closed conformer, as the bond lengthens, a rotation of the phenyl groups occurs, resulting in the switch of the structure from closed to open after 2.66 Å. This exactly corresponds to the sudden drop of the shielding constant in Figure 6A.

The decomposition of σ(¹²⁵Te) computed for the structure with the largest Te–Te distance shows that the value of σ_p changes from −3433.1 ppm in the minimum energy structure of the open conformer to −4054.5 ppm. Inspection of the most important contributions of σ_p^occ-virt in the deformed geometry shows that the pairs of orbitals that play an important part in determining the value of the ¹²⁵Te paramagnetic shielding constant involve the LUMO as the empty orbital (Table S4). The LUMO of Ph₂Te₂ is an anti-bonding σ MO mainly composed by the out-of-phase overlap of two p_z orbitals on the two chalcogen atoms. Upon stretching of the Te–Te bond, it remains almost unchanged (Figure 7), but decreases in energy leading to a progressively smaller HOMO-LUMO gap (Figure 6A). Thus, since the σ_p^occ-virt are weighted by a factor that is inversely proportional to the energy gap between the involved MOs, the stabilization of the LUMO implies an increase of the magnetic coupling of all pairs of MOs including the LUMO as the empty molecular orbital. This leads to a larger (negative) σ_p and consequently to a higher ¹²⁵Te chemical shift.

Finally, we investigated the effects of a rotation around the Te-Te bond on the ¹²⁵Te shielding constant. Results showed a decrease of σ(¹²⁵Te) from an initial value of 2169.3 ppm (Ψ = 0°) to a relative minimum of 2061.7 ppm recorded at Ψ = 76°, followed by an increase up to 2190.8 ppm (Figure 6B). Right after the minimum point, a non-monotonous increase can be observed, which is caused by the slightly different values of Φ₁ and Φ₂ at each point during the completely relaxed torsional scan. Nevertheless, the profile of σ(¹²⁵Te) follows the trend of the Te–Te distance, since when the chalcogen atoms get closer or further, the shielding constant decreases or increases, respectively (Figure 6B). This behavior follows from Schreckenbach analysis since the elongation of the inter-chalcogen bond leads to a stabilization of the LUMO and thus to a larger (negative) σ_p and consequently to a higher ¹²⁵Te chemical shift, as explained in the previous paragraph.

The HOMO–LUMO pair provides a very small σ_p^occ-virt contribution in the structures with Ψ equal to 0° or 180° due to their shape (Figure 8) that results in almost no HOMO–LUMO overlap after application of the magnetic moment operator due to the cancellation of the positive and negative orbital phase. Conversely, the largest contributions to σ_p arise from the HOMO-2–LUMO and HOMO-1–LUMO pairs when Ψ = 0° with values of −716.8 and −698.1 ppm respectively, and from the pair HOMO-1–LUMO when Ψ = 180°, with a value of −801.7 ppm. Although the energy gap is larger than for HOMO–LUMO, the magnetic coupling is strong due to the composition of these orbital pairs in terms of Te p atomic orbitals. Comparing the two extreme structures obtained in the torsional scan with the fully optimized open conformer of Ph₂Te₂, it is clear how the intermediate value of Ψ found in the minimum energy structure, with the two phenyl rings arranged in skewed fashion, favors not only the strongest orbital interactions, but also enhances the magnetic couplings, resulting in the structure with the lowest isotropic shielding constant. These analyses have been carried out also with the inclusion of spin-orbit effects: the two limit structures, Ψ = 0° and Ψ = 180°, have been chosen because of the biggest difference between relativistic scalar and spin-orbit ¹²⁵Te shielding constants. The larger, and therefore the main, σ_p contributions are localized as the same shown above using the scalar ZORA approximation (Tables S5 and S6). This allowed us to entirely perform the analysis within a scalar approximation without losing important magnetic coupling contributions and obtaining reliable data that are qualitatively comparable with those obtained with the extended spin-orbit Hamiltonian.

3. Materials and Methods

All geometry optimizations were carried out with the Amsterdam density functional (ADF) program (version 2017, SCM, Amsterdam, The Netherlands) developed by Baerends and co-workers [47,48,49]. In agreement with a recent benchmark study [10], the OPBE functional [50,51,52,53] was chosen in combination with the TZ2P basis set, which is a large and uncontracted set of Slater-type orbitals (STOs) of triple-ζ quality, augmented with two sets of polarizations functions for each atom: 2p and 3d for hydrogen, 3d and 4f in the case of carbon, and 5d and 4f for tellurium. The frozen core approximation was employed: up to 1s for carbon and up to 4p for tellurium. Relativistic effects were included in the calculations using the scalar relativistic zeroth-order regular approximation (ZORA) [54]. This level of theory is denoted ZORA-OPBE/TZ2P. The scans of selected geometrical parameters were carried out at the ZORA-OPBE/TZ2P level of theory with the (linear) transit protocol implemented in ADF.

For the computation of the δ(¹²⁵Te), the calculations of the shielding tensors were carried out with the NMR module of the ADF program (version 2017, SCM, Amsterdam, The Netherlands) [45], using the OPBE functional in combination with the QZ4P all-electron basis set [55] in chloroform (ε = 4.8) employing the COSMO solvation model [56,57]. This level of theory, which gave very good results for the δ(⁷⁷Se) calculations on organic mono and diselenides [11], is here denoted COSMO-ZORA-OPBE/QZ4Pae. As a reference, in order to convert the isotropic shielding constant to δ(¹²⁵Te), using the relationship δ = σ_ref − σ, (CH₃)₂Te was employed [43] with a calculated σ_ref of 2670.3 ppm at the COSMO-ZORA-OPBE/QZ4Pae level of theory. Spin-orbit effects have been taken into account in selected cases and are denoted in this paper as COSMO-SO-ZORA-OPBE/QZ4Pae.

The analysis of σ was carried out in the gas phase at the ZORA-OPBE/QZ4Pae level of theory. The shielding constant can be decomposed as the sum of the diamagnetic (σ_d) and the paramagnetic (σ_p) contributions [46]. The spin-orbit term was neglected because it was not explicitly considered in these calculations. The diamagnetic term depends mainly on the core–shell electrons and is therefore almost constant for a given nucleus in different compounds and cancels out when calculating chemical shifts. The paramagnetic term can be broken down into the sum of two contributions

$${\sigma }_p={\sigma }_p^{\mathrm{occ}-\mathrm{occ}}+{\sigma }_p^{\mathrm{occ}-\mathrm{virt}},$$

(1)

in which σ_p^occ-occ is a sum over all pairs of occupied orbitals and σ_p^occ-virt is a sum over all pairs of orbitals consisting of an occupied orbital and virtual orbital; each term is weighted by the inverse of the energy gap and is proportional to the magnetic coupling between the involved molecular orbitals

$$\langle {\mathsf{\varphi }}_{\mathrm{virt}}|\widehat{B}|{\mathsf{\varphi }}_{\mathrm{occ}}\rangle ,$$

(2)

in which $\widehat{B}$ is the magnetic moment operator (Scheme 2). σ_p^occ-virt is the largest contribution to σ_p and is mostly responsible for the variation in chemical shift.

To run MD simulations, a novel set of GAFF (generalized AMBER force field) parameters for Ph₂Te₂ was obtained. A recent computational protocol devised by some of us was employed for this purpose [28]. The two energetically relevant conformers of Ph₂Te₂, fully optimized using OPBE [50,51,52,53] combined with the TZ2P basis set for all the atoms with the ADF program, were used as starting structures for all the scans (see below). Charges were calculated using R.E.D. Server Development [60] which is a web-based interface that employs PyRED program (version SEP-2015, Université de Picardie-Jules Verne, Amiens, France and Sanford Burnham Prebys Medical Discovery Institute, La Jolla, USA) [61] and R.E.D. tools (version SEP-2015, Université de Picardie-Jules Verne, Amiens, France and Sanford Burnham Prebys Medical Discovery Institute, La Jolla CA, USA) [62]. Calculations of the charges were carried out, as recommended [63], at the Hartree-Fock level of theory using the 6-31G(d) basis set for all the atoms except tellurium for which the Stuttgart Dresden effective core potential basis set [64,65] was employed; the Gaussian09 package (rev D.01, Gaussian, Inc.: Wallingford CT, USA) was used for these calculations [66]. To obtain highly reproducible charges, two different orientations were taken into account for each conformer, using the rigid-body re-orientation algorithm implemented in R.E.D. Tools. In order to obtain the final charge values for Ph₂Te₂, a simultaneous restrained electrostatic potential (RESP) [67] fit of four structures was performed.

Constrained scans (i.e., calculations where only one relevant parameter is changed and the rest of the molecule is kept frozen) were carried out at ZORA-OPBE/TZ2P to obtain the force constants of the harmonic contributions (i.e., a_i and b_i terms in Equation S1). Stretching and bending parameters were thus computed for Te–C and Te–Te bonds and for Te–Te–C and Te–C–C angles, respectively. Relaxed scans of the most stable conformer were performed for the dihedrals C–Te–Te–C and C–C–Te–Te. Moreover, although generalized parameters for the dihedral Te–C–C–H (in which the chalcogen atom is in a terminal position) were already present in GAFF, we decided to improve the data by refitting the coefficients. The energies of the resulting ensemble of 72 structures were simultaneously fit with Paramfit [68], a tool included in AmberTools15 software [69]. These new parameters, which are included in the Supplementary Materials, were used in MD simulations of a single molecule of Ph₂Te₂ using the Assisted Model Building with Energy Refinement (AMBER) suite of programs (version 2015, University of California, San Francisco, San Francisco, CA, USA) [69]. The molecule was initially solvated in chloroform in a cubic box with a surrounding buffer of 12.0 Å. First, the potential energy was minimized by keeping the diphenyl ditelluride atoms fixed and relaxing the solvent; afterwards, a full minimization was done. The temperature was first increased using the Langevin thermostat up to 100 K for 20 ps and, subsequently, the Berendsen barostat was added for 1 ns to reach the production temperature (310 K) and pressure (1 bar). A pre-equilibration (200 ns) and production dynamics (500 ns) were conducted by switching to the Monte Carlo barostat. The cutoff distance for the nonbonded interactions was set to 10.0 Å; the shake constraints were applied to all hydrogen atoms, and a time step of 2 fs was used.

4. Conclusions

In this work, the effects of structural changes on the ¹²⁵Te chemical shift of Ph₂Te₂ were investigated using a combined quantum mechanical and classical computational methodology. Relativistic density functional theory-based ¹²⁵Te-NMR calculations were performed using fully optimized conformers of Ph₂Te₂ as well as on structures extracted from a classical MD trajectory. At the scalar relativistic level, the agreement with the experimental result is satisfactory in both instances, but a non-negligible discrepancy can be observed. However, the inclusion of spin-orbit corrections leads to a much better agreement already when considering only the two energy-minimized open and closed conformers; the calculated result approaches the experimental values within just 17 ppm when the average geometry obtained from the dynamics is considered. Nonetheless, the MD simulation highlighted the high flexibility of Ph₂Te₂ molecule, affording a very wide span of calculated chemical shifts, from −225.5 to 1503.4 ppm as a function of the phenyl rings orientation and Te–Te bond lengths (see Figure S1 in Supplementary Materials). Again, the time-averaged ¹²⁵Te chemical shift is, therefore, strongly dependent on the very accurate modeling of the whole potential energy surface, which suffers from the intrinsic approximations of classical force fields. Thus, the take-home message here is that conformational effects in diphenyl ditelluride (and closely related derivatives), are dominant over other issues like solvent, temperature, and aggregation, in the calculations of NMR properties. Carefully taking into account the conformational landscape and the spin-orbit correction leads to a very satisfactory agreement.

Finally, a detailed decomposition of the calculated isotropic shielding constant was performed on the stable minima as well as on geometrically ad hoc deformed Ph₂Te₂ molecules to gain insight on the origin of the strong variation of ¹²⁵Te chemical shift. Calculations on a simpler model monochalcogenide—i.e., PhTeH—were also included to derive more general conclusions about the relationship between the geometry and the shielding constant in organotellurides.

The analyses conducted on PhTeH and Ph₂Te₂ show that, in general, a stretching of the bond in which Te is involved implies a lowering of the LUMO energy, a larger (negative) σ_p, and consequently a higher ¹²⁵Te chemical shift. Conversely, the rotation about a bond in which Te is involved—i.e., about the Te–Te bond (in Ph₂Te₂ ) or the C–Te bond (in PhTeH)—leads to larger (negative) σ_p, whenever a favorable magnetic coupling of occupied/empty molecular orbital pairs is established; since in general these interactions involve high energy unoccupied molecular orbitals, the variation of ¹²⁵Te σ (and δ) is less significant, due to the large energy gaps of these pairs.