Mechanism of reaction between titanocene pentasulfide and sulfenyl chloride: The effect of continuum solvation on the energy surface.

An investigation of the kinetics of the reaction between titanocene pentasulfide and sulfenyl chloride[1] leading to the formation of the S₇ allotrope of sulfur was accompanied by supporting DFT calculations which led to the conclusion that of five possible mechanisms for the reaction, the most probable corresponded to a variant of the concerted σ-bond metathesis (Scheme 1, Mechanism IV, R = Cl). Here we take a closer look at the DFT predictions from the point of view of the impact of continuum solvation on the calculated mechanism.

Scheme 1.  Possible reaction mechanisms.

The original study used the new r²scan-3c/Def2-mTZVPP composite DFT functional[2] as implemented in the ORCA 6 program code,[3] for which the (gas phase) mechanistic reported pathway was obtained as shown in scheme 2 below. Here, we re-label the succeeding steps as TS3 and TS4 (see Scheme 3 below, rather than TS2 from scheme 2) to form the final product P, for reasons that will become apparent below.

Scheme 2.  Reaction  scheme and energy profile.[1]

Methodology

In this study, we used an older functional, MN15L[4] for R=Cl (scheme 1), which we have found very effective for the transition elements[5] and which – like r²scan-3c – was designed to be accurate for multi-reference and single-reference systems and for noncovalent interactions, This functional, unlike r²scan-3c, is implemented in the Gaussian 16 program code and hence has the advantage of allowing computed intrinsic reaction coordinate (IRC) data to be usefully visualised using e.g. Gaussview.^‡ Transition state geometries were initially obtained from the supporting information given in the original article [1] and re-optimised in the Gaussian program using MN15L/Def2-TZVPP. The thermochemical energies shown in Table A1 were all obtained using GoodVibes[6] (see manual) with an entropic quasi-harmonic treatment frequency cut-off value of 2.0 wavenumbers[7] and an enthalpic quasi-harmonic treatment frequency cut-off value of 50.0 wavenumbers.[8] In the table, harmonic values are indicated as e.g. hG and quasi-harmonic values as qh-G; the difference between these two is relative small. If desired, other harmonic cut-off values can be obtained, as well as at other molar concentrations, using log files obtained from the repository DOIs indicated below, via the command line:

python3 -m goodvibes -q --fs 2 --fh 50 -c 0.0409 logfilename.

Results

Part 1: Gas phase model

The intrinsic reaction coordinate (IRC) deriving from transition state TS1 computed using Gaussian 16, and without inclusion of a solvent continuum model (a gas phase model) is shown below (Table 1, Figure 1).[9] It leads directly to the “half-way” product Int2, with no intervening intermediate such as that shown in Scheme 2 (there labelled Int1[1]). So here, the TS1 IRC conflates the originally reported TS1 and TSint[1] as shown in scheme 2, with the conflation point occuring at an IRC value of ~9. This point can also be seen below as a prominent “hidden intermediate” in the gradient norm plot at the same IRC value. A gradient norm at this point of not quite zero is what makes it a “hidden” rather than a “real” intermediate. The conflation point also ~corresponds to a minimum in the dipole moment plot. Here,[9] the stepsize between points in the IRC calculation was selected as 20 (in units of 0.01 Bohr) and the Hessian was recalculated every 5 steps. The equivalent parameters for the IRCs as noted – but not visualised – in the original article[1] were not stated; it is entirely possible that differences in either these parameters, or the algorithms used to compute the IRC or indeed the use of a different functional could account for this slight difference in behaviour. Slight, because TSint is shown as a very shallow intermediate (Scheme 2[1]).

Figure 1. Computed geometry of TS1 in the gas phase at the MN15L/Def2-TZVPP level.

Table 1. IRC for TS1 in a gas phase MN15L/Def2-TZVPP model.
Total Energy
Gradient norm
Dipole moment

Part 2: Continuum solvent model (scheme 3).

FAIR data for all the calculations conducted using a dichloromethane continuum solvent are summarised here, [10] with individual calculations indicated as a reference to a FAIR repository dataset  (Table A1). The computed geometry of TS1 changes from the initial values of a) the new S₇-S₈ bond 2.195 → 2.356Å, b) S₈-Cl 2.918 → 2.616Å and c) S₇-Ti 2.572 → 2.562Å (Table 2, atom numbering shown in  Figure 1).

Scheme 3.  Revised reaction scheme showing the formation of the ion-pair Int0 rather than Int1 computed using geometries optimised with continuum solvation effects included.

When same potential energy surface is computed with a continuum solvent model (Table 3), the “hidden intermediate” present for the gas phase model in the original IRC profile of TS1 (Table 1) now becomes a real ion-pair intermediate (labelled here Int0 in  scheme 3 to distinguish it from Int1 in scheme 2) with a discrete chloride anion. Thus Int0 occurs at an IRC value of ~3.5 in the plots below, although it has a very small exit barrier via a transition state, here labelled TS2 (different from the TS2 labelled in scheme 2 above). The plots below  (Table 3) are the result of concatenating two separate IRC plots for TS1 and TS2.

Table 3. IRCs for TS1 + TS2 in a continuum solvent model

The gas phase (left) and continuum solvent (right) IRCs are summarised below (Figure 2) to enable a visual comparison of the two potential energy surfaces.

Figure 2. A comparison of the computed IRC energy profiles in the gas phase (left) and continuum dichloromethane solvent (right).

This indicates a change from Mechanism IV under gas phase conditions (scheme 2) to one closer to mechanism II with continuum solvent; an S_N2 like displacement of chloride ion at sulfur, followed in a second step by Cl…Ti bond formation and Ti…S cleavage, Scheme 3. This can also be summarised by the following plots of bond lengths in Figure 3.

Figure 3. Selected bond length variation for the concatenated IRC profile for TS1 + TS2 in a continuum solvent. See Figure 1 for atom numberings.

The overall results can be summarised in Table 4 indicting that both the original r²scan-3c/Def2-mTZVPP and the MN15L functional used here are both in reasonable agreement with the experimental results obtained from kinetic studies. Also noteworthy is that for substituents such as e.g. 2b, the enthalpy of activation may actually be negative, with the resulting positive value for the free energy of activation being a consequence of a very negative entropy of activation.

Conclusions

The overall conclusion is simple; when the possibility of ion-pair formation on a reaction potential energy surface is present, it matters how the geometries of all the species involved are obtained. A gas phase geometry optimised model is likely to disfavour the formation of such an ion-pair, whereas a continuum solvent geometry optimised model is more likely to promote such species. This effect can be seen especially in Figure 2, where the two potentials are shown side by side. Such promotion of ion-pairs is likely to reduce any concerted behaviour of the reaction into a two-stage stepwise process. Thus the concerted nature of mechanism IV (Scheme 1) morphs into more stepwise behaviour (Mechanism II, Scheme 1 modified as in Scheme 3). Overall however, while the resulting calculated energetic barriers, although reduced by inclusion of solvation, are unlikely provide definitive evidence of which mechanism actually prevails, the greater change in dipole moment in the stepwise behaviour is more consistent with the experimentally observed effects of solvent identity on the rate.

Appendix

All energies for the computed species, obtained with optimisation with a dichloromethane continuum solvent model. Values for the default concentration of 0.0409M are shown for completeness and to enable facile reconciliation with those included in the published FAIR datasets.

^a0.0409M ^b1.0M ^cInt0 is ion pair comprising S⁺ and Cl^– See Table 3 and scheme 3.

Footnotes

^‡Currently at least, Gaussian is better supported by associated visualisation programs then ORCA.

This post has DOI: 10.59350/1a48f-rj714

References

1. P.H. Helou de Oliveira, P.J. Boaler, G. Hua, N.M. West, R.T. Hembre, J.M. Penney, M.H. Al-Afyouni, J.D. Woollins, A. García-Domínguez, and G.C. Lloyd-Jones, “Kinetics of sulfur-transfer from titanocene (poly)sulfides to sulfenyl chlorides: rapid metal-assisted concerted substitution”, Chemical Science, vol. 15, pp. 11875-11883, 2024.
2. S. Grimme, A. Hansen, S. Ehlert, and J. Mewes, “r2SCAN-3c: A “Swiss army knife” composite electronic-structure method”, The Journal of Chemical Physics, vol. 154, 2021.
3. F. Neese, F. Wennmohs, U. Becker, and C. Riplinger, “The ORCA quantum chemistry program package”, The Journal of Chemical Physics, vol. 152, 2020.
4. H.S. Yu, X. He, and D.G. Truhlar, “MN15-L: A New Local Exchange-Correlation Functional for Kohn–Sham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids”, Journal of Chemical Theory and Computation, vol. 12, pp. 1280-1293, 2016.
5. E.M. Richards, L. Casarrubios, J.M. D’Oyley, H.S. Rzepa, A.J.P. White, K. Goldberg, F.W. Goldberg, J.A. Bull, and S. Díez‐González, “Bidentate NHC‐Containing Ligands for Copper Catalysed Synthesis of Functionalised Diaryl Ethers”, Advanced Synthesis & Catalysis, vol. 367, 2024.
6. I. Funes-Ardoiz, and R.S. Paton, “GoodVibes: version 2.0.3”, 2018.
7. S. Grimme, “Supramolecular Binding Thermodynamics by Dispersion‐Corrected Density Functional Theory”, Chemistry – A European Journal, vol. 18, pp. 9955-9964, 2012.
8. Y. Li, J. Gomes, S. Mallikarjun Sharada, A.T. Bell, and M. Head-Gordon, “Improved Force-Field Parameters for QM/MM Simulations of the Energies of Adsorption for Molecules in Zeolites and a Free Rotor Correction to the Rigid Rotor Harmonic Oscillator Model for Adsorption Enthalpies”, The Journal of Physical Chemistry C, vol. 119, pp. 1840-1850, 2015.
9. H. Rzepa, “TS1 Cp2TiS5 + S2Cl2, gas phase, Def2-QZVPP on Ti, G = -4943.798931 IRC”, 2025.
10. H. Rzepa, “Mechanism of reaction between titanocene pentasulfide and sulfenyl chloride”, 2025.
11. H. Rzepa, “CYPTIS5 Def2-QZVPP on Ti, TZVPP on rest G = -3227.179639 (Def2-TZVPP = -3227.141839)”, 2025.
12. H. Rzepa, “S2Cl2, G = -1716.667668”, 2025.
13. H. Rzepa, “TS1 for Sn2 by S on S-Cl Def2-QZVPP on Ti, G = -4943.825172 (TS2 -4943.827519)”, 2025.
14. H. Rzepa, “Ion pair Int SCRF(CPCM), G = -4943.831670”, 2025.
15. H. Rzepa, “TS2 Def2-QZVPP on Ti, TZVPP on rest G = -4943.827519”, 2025.
16. H. Rzepa, “TS1, IRC => Product G = -4943.872859”, 2025.
17. H. Rzepa, “C2pTiS5. TS3, second SN2 with Def2-QZVPP on Ti, G = -4943.858887”, 2025.
18. H. Rzepa, “C2pTiS5. TS3, second SN2 with Def2-QZVPP on Ti ==> Int3 G = -4943.859300 ==> Int3 G =0.106 oz -4943.864446”, 2025.
19. H. Rzepa, “Cp2TsS5 + S2Cl2, TS4, G = -4943.862461”, 2025.
20. H. Rzepa, “C2pTiS5. Final Product, G = -4943.902398 S7 -2787.152347”, 2025.
21. H. Rzepa, “C2pTiS5. Final Product, G = -4943.902398 Cp2TiCl2 G =-2156.754402 + -2787.152347 = -4,943.906749”, 2025.
22. H. Rzepa, “2a, Def2-QZVPP on Ti, G = -2946.7928059”, 2025.
23. H. Rzepa, “TS2 for Sn2 by S on S-Cl Def2-QZVPP on Ti, G = -4943.825172 (TS1 -4943.827519) => 2a G = -4663.439222”, 2025.
24. H. Rzepa, “TS1, Def2-QZVPP on Ti, TZVPP on rest G = -4943.827519 => 2a G = -4663.446325”, 2025.
25. H. Rzepa, “2b CYPTIS4C(NMe2)2 G = -3135.870870”, 2025.
26. H. Rzepa, “TS2 for Sn2 by S on S-Cl Def2-QZVPP on Ti, => 2b G = -4852.523215”, 2025.
27. H. Rzepa, “TS1, Def2-QZVPP on Ti, TZVPP on rest G = -4943.827519 => 2b G = -4852.528126”, 2025.
28. H. Rzepa, “6m. CYPTIS4C(CN)2 G = -3052.669074”, 2025.
29. H. Rzepa, “TS2 for Sn2 by S on S-Cl Def2-QZVPP on Ti, => 6m G = -4769.302718”, 2025.
30. H. Rzepa, “TS1, Def2-QZVPP on Ti, TZVPP on rest => 6m”, 2025.