Blog

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[cite]10.1039/d4sc02737j[/cite] 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[cite]10.1063/5.0040021[/cite] as implemented in the ORCA 6 program code,[cite]10.1063/5.0004608[/cite] 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.[cite]10.1039/d4sc02737j[/cite]

Methodology

In this study, we used an older functional, MN15L[cite]10.1021/acs.jctc.5b01082[/cite] for R=Cl (scheme 1), which we have found very effective for the transition elements[cite]10.1002/adsc.202400909[/cite] 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 [cite]10.1039/d4sc02737j[/cite] and re-optimised in the Gaussian program using MN15L/Def2-TZVPP. The thermochemical energies shown in Table A1 were all obtained using GoodVibes[cite]10.5281/zenodo.1435820[/cite] (see manual) with an entropic quasi-harmonic treatment frequency cut-off value of 2.0 wavenumbers[cite]10.1002/chem.201200497[/cite] and an enthalpic quasi-harmonic treatment frequency cut-off value of 50.0 wavenumbers.[cite]10.1021/jp509921r[/cite] 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).[cite]10.14469/hpc/15219[/cite] It leads directly to the “half-way” product Int2, with no intervening intermediate such as that shown in Scheme 2 (there labelled Int1[cite]10.1039/d4sc02737j[/cite]). So here, the TS1 IRC conflates the originally reported TS1 and TSint[cite]10.1039/d4sc02737j[/cite] 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,[cite]10.14469/hpc/15219[/cite] 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[cite]10.1039/d4sc02737j[/cite] 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[cite]10.1039/d4sc02737j[/cite]).

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, [cite]10.14469/hpc/15204[/cite] 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.

Table 2. Calculated geometries for TS1 – Ts4 with solvation
#	Geometry	#	Geometry
TS1		Int0
TS2		Int2
TS3		TS4

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.

Table 4. Solvent DCM model: Mechanism II through to Ion-Pair intermediate Int0 rather than the original Int1.
Source	ΔH	ΔS	ΔG
1a: Article	6.0	-45.7	19.6
1a: This work	3.4	-31.5	12.8
1a: Expt/1M	0.0	-49.7	14.8
2a: Article	1.7	-45.7	15.3
2a: This work	0.8	-39.5	12.6
2a: Expt/1M	-2.3	-48.2	12.6
2b: Article			~10.8
2b: This work	-1.4	-34.1	8.80
6m: Article			~25.0
6m: This work	8.95	-39.2	20.6

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.

Table A1. Calculations at the MN15L/Def2-TZVPP; Def2-QZVPP on Ti/CPCM=Dichloromethane level with geometry optimisation.
1a (scheme 2)
Species	h-H	T.h-S	h-G	qh-H	T.qh-S	qh-G
Z=S	-3227.1212^a -3227.1212^b	0.05845 0.05543	-3227.179639 [cite]10.14469/hpc/15150[/cite] -3227.176621	-3227.12140 -3227.12140	0.05845 0.05543	-3227.17985 -3227.17683
S2Cl2	-1716.63131^a -1716.63135^b	0.03637 0.03335	-1716.667668 [cite]10.14469/hpc/15196[/cite] -1716.664658	-1716.63135 -1716.63135	0.03637 0.03335	-1716.66771 -1716.66470
Sum	-4943.75254^a -4943.75254^b	0.09481 0.08878	-4943.847307 -4943.841279	-4943.75275 -4943.75275	0.09481 0.08878	-4943.84756 -4943.84152
TS1	-4943.74534^a -4943.74534^b	0.07984 0.07682	-4943.825172 [cite]10.14469/hpc/15179[/cite] -4943.822153	-4943.74734 -4943.74734	0.07984 0.07682	-4943.82717 -4943.82415
ΔTS1	4.52	-31.53	13.88	3.39	-31.52	12.79
ΔTS1^b	4.52	-25.17	12.00	3.39	-25.17	10.90
Int0^c	-4943.75128	0.08039	-4943.831670 [cite]10.14469/hpc/15205[/cite]	-4943.75298	0.08039	-4943.83337
ΔInt0	0.76	-30.36	9.81	-0.15	-30.36	8.91
TS2	-4943.74951	0.07801	-4943.827519 [cite]10.14469/hpc/15154[/cite]	-4943.75100	0.07801	-4943.82901
ΔTS2	1.88	-35.33	12.42	1.10	-35.36	11.64
Int2	-4943.79255	0.08031	-4943.872859 [cite]10.14469/hpc/15206[/cite]	-4943.79432	0.08031	-4943.87463
ΔInt2	-25.13	-30.53	-16.03	-26.09	-30.53	-16.99 (-11.0 lit)
TS3	-4943.77886	0.08003	-4943.858887 [cite]10.14469/hpc/15248[/cite]	-4943.78065	0.08003	-4943.86068
ΔTS3	-16.54	-30.12	-7.27	-4943.78065		-8.23
Int3	-4943.78456	0.079887	-4943.864446 [cite]10.14469/hpc/15254[/cite]	-4943.78619	0.079887	-4943.86607
ΔInt3	-20.12	-31.42	-10.76	-20.98	-31.42	-11.62
TS4	-4943.78439	0.07807	-4943.862461 [cite]10.14469/hpc/15251[/cite]	-4943.78599	0.07807	-4943.86406
ΔTS4	-20.01	-35.23	-9.51	-20.86	-35.23	-10.36
S7 as product	-2787.10568	0.04667	-2787.152347 [cite]10.14469/hpc/15578[/cite]	-2787.10608	0.04667	-2787.15275
Cp2TiCl2	-2156.70450	0.04990	-2156.754401 [cite]10.14469/hpc/15579[/cite]	-2156.70456	0.04990	-2156.75447
Product P	-4943.81017	0.09657	-4,943.906749	-4943.81064	0.09657	-4,943.90722
ΔProduct P	-36.19	3.70	-37.30	-36.33	3.70	-37.43 (-36.9 lit)
2a:
Z=CMe2	-2946.72598^a -2946.72598^b	0.06683 0.063808	-2946.7928059 [cite]10.14469/hpc/15193[/cite]-2946.789786	-2946.72674 -2946.72674	0.06683 0.063806	-2946.79356 -2946.790542
S2Cl2	-1716.63131^a -1716.63135^b	0.03637 0.03335	-1716.667668 [cite]10.14469/hpc/15196[/cite] -1716.664658	-1716.63135 -1716.63135	0.03637 0.03335	-1716.66771 -1716.66470
Sum	-4663.35729 -4663.35729	0.10319 0.097158	-4663.4604739 -4,663.454444	-4663.35808 -4663.35808	0.10319 0.097156	-4663.46127 -4663.455242
TS1	-4663.35477^a -4663.35477b	0.08445 0.081429	-4663.439222 [cite]10.14469/hpc/15195[/cite]-4663.436203	-4663.35677 -4663.35677	0.08445 0.081429	-4663.44122 -4663.438201
ΔTS1^a	1.58	-39.45	13.34	0.82	-39.46	12.59
ΔTS1^b	1.58	-33.10	11.45		-33.10	10.69
TS2	-4663.36425	0.08207	-4663.446325 [cite]10.14469/hpc/15194[/cite]	-4663.36574	0.08207	-4663.44782
ΔTS2	-4.37	-44.44	8.88	-4.80	-44.44	8.44
2b:
Z=C(NMe2)2	-3135.79512 -3135.79512	0.07575 0.072730	-3135.870870 [cite]10.14469/hpc/15197[/cite]-3135.867852	-3135.79611 -3135.79611	0.07575 0.072730	-3135.87186 -3135.86884
S2Cl2	-1716.63131^a -1716.63135^b	0.03637 0.03335	-1716.667668 [cite]10.14469/hpc/15196[/cite] -1716.664658	-1716.63135 -1716.63135	0.03637 0.03335	-1716.66771 -1716.66470
Sum	-4852.4264 -4852.4264	0.11212 0.10608	-4852.538538 -4852.53251	-4852.42745 -4852.42745	0.11212 0.10608	-4852.5396 -4852.53354
TS1	-4852.42719 -4852.42719	0.09602 0.09300	-4852.523215 [cite]10.14469/hpc/15199[/cite] -4852.520196	-4852.42964 -4852.42964	0.09591 0.09289	-4852.52555 -4852.522536
ΔTS1^a	-0.49	-33.86	9.62	-1.37	-34.11	8.80
ΔTS1^b	-0.49	-27.53	7.27	-1.37	-27.76	6.91
TS2	-4852.43669	0.09143	-4852.528126 [cite]10.14469/hpc/15198[/cite]	-4852.43823	0.09143	-4852.52966
ΔTS2	-6.44	-43.54	6.53	-6.76	-43.54	6.22
6m:
Z=C(CN)2	-3052.59951 -3052.59951	0.06956 0.06654	-3052.669074 [cite]10.14469/hpc/15200[/cite] -3052.666055	-3052.60050 -3052.60050	0.06956 0.06654	-3052.67007 -3052.66705
S2Cl2	-1716.63131^a -1716.63135^b	0.03637 0.03335	-1716.667668 [cite]10.14469/hpc/15196[/cite] -1716.664658	-1716.63135 -1716.63135	0.03637 0.03335	-1716.66771 -1716.66470
Sum	-4769.23082 -4769.23082	0.10593 0.09989	-4769.336742 -4769.330713	-4769.23185 -4769.23185	0.10593 0.09989	-4769.3378 -4769.33175
TS1	-4769.21540^a -4769.21540^b	0.08731 0.08423	-4769.30271 [cite]10.14469/hpc/15203[/cite] -4769.29970	-4769.21759 -4769.21759	0.08731 0.08423	-4769.30490 -4769.30188
ΔTS1^a	9.68	-39.18	21.36	8.95	-39.18	20.63
ΔTS^b	9.68	-32.83	19.46	8.95	-32.83	18.74
TS2	-4769.22451	0.08582	-4769.310328 [cite]10.14469/hpc/15201[/cite]	-4769.22626	0.08582	-4769.31208
ΔTS2	3.96	-42.32	16.58	3.51	-42.32	16.13

^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

December 16, 2025

Reinvestigating the reported transition state structure of a concerted triple H-tunneling mechanism.

Substituting a deuterium isotope (²H) for a normal protium hydrogen isotope can slow the rate of a chemical reaction if this atom is involved in the reaction mode. The magnitude of the effect, referred to as a kinetic isotope effect or KIE is normally 2-7, but higher values of 20 or even more^♥ are sometimes observed due to a phenomenon known as proton tunnelling. So a recent report[cite]10.1021/acscentsci.5c00943[/cite] of a ¹H/²H of ~2440 for the following palladium catalysed reaction caught my eye:

When the protium in the solvent methanol and the hydrogen gas were replaced by deuterium, the rate of the reaction slowed by ~2400. This immediately begs one question: what was the % of deuterium incorporated into the ²H₂ and MeOD? It would have to be >99.994% to eliminate any contribution from the presumably faster reacting ¹H isotope, and this level of deuteriation is some ask! Leaving this issue aside, the authors then carried out some DFT modelling to come up with a proposed mechanism (below), which they refer to as a concerted triple hydrogen transfer reaction (the curly arrows by the way are mine; arrows are shown in the graphical abstract for this article but they are likely not curly electron arrows but simply schematic). The large value of the KIE was then attributed in part to a novel form of triple hydrogen tunnelling.

My second reality check was to search the crystal structure database for instances of the proposed catalyst containing a Pd-H substructure. Nine examples of compounds with such Pd-H bonds emerge, but none have the H-Pd(OR)₃ motif shown above, which is likely to be a transient catalytic species rather than a stable isolable one. This species (FOTBAR)[cite]10.1021/ic00306a034[/cite] with one OPh and two P ligands on the Pd is the closest match; although the trans relationship of the Pd-H and Pd-O bonds might preclude it functioning as a catalyst according to the mechanism above.

My next check related to the DFT procedure used, which was reported as B3LYP with apparently a 6-311+G(d,p) basis set, but no dispersion correction added. We had previously observed[cite]10.1002/adsc.202400909[/cite] that functionals such as B3LYP are not particularly well suited for transition metal modelling, preferring a newer variety such as MN15L.[cite]10.1021/acs.jctc.5b01082[/cite]

Finally, we also recollected our experience in modelling KIE effects using relatively modest basis sets such as 6-31G(d,p) and 6-311+G(d,p)[cite]10.1039/D3DD00246B[/cite] where we showed that the calculated KIE were inaccurate. Basis sets of eg Def-TZVPP or better were found to be essential. So here I test this hypothesis for a small selection of functional and basis sets as an initial exploration.

The calculations are published here (Table below).[cite]10.14469/hpc/15569[/cite] Row 1 shows the values given in the article,[cite]10.1021/acscentsci.5c00943[/cite] and for which a free energy of activation of 27.0 kcal/mol was indicated. Attempting to replicate this here, the main article declares that a 6-31G* basis set was adopted for the H, C, and O atoms… and the LANL2DZ basis set was adopted for Pd atoms. The supporting information records this instead as 6-311+G(d,p) for these atoms (Table S5. — B3LYP/6-311+G(d,p)/Lanl2dz level) which was used here. The basis used for Ti was not noted in either article or ESI; here it was set to 6-31G(d,p).[cite]10.1021/acs.jcim.9b00725[/cite] Using the Gaussian 16 program, my calculation gave the results shown in row 2, with geometry optimisation starting from the coordinates given in the ESI, giving a final RMS force of 0.000008 au – this value is not available for comparison with the original article, nor is the final total energy of the system. The imaginary transition state mode is 940 cm^-1 compared to the reported value of 1306 cm^-1, a not insignificant difference and which may arise from the reported basis set uncertainty. The bond lengths also differ somewhat, but the angle subtended at the Pd-H-C system is more or less linear. The newly computed free energy of activation is significantly lower. Re-modelling, but now including the effects of a methanol solvent also induces some significant changes in the geometry, but only a small change in the imaginary mode to 979 cm^-1 (entry 3).

Changing the functional from B3LYP to MN15L (entry 4) significantly reduces the imaginary mode value and here the effect of improving the basis set quality (entry 5) is large, reducing the imaginary transition state mode to 497 cm^-1 Entry 6 shows the values for the r²scan-3c functional discussed in an earlier post,[cite]10.59350/bc8j8-dtj11[/cite] revealing a transition state mode similar to the others. The free energy barriers range from 27.0 kcal/mol quoted in the article[cite]10.1021/acscentsci.5c00943[/cite] down to 18.3 (entry 3) with the r²scan-3c functional being rather higher. Given that this reaction proceeds at temperatures of 253 – 298K, one might expect a barrier closer to the lower end of this range rather than the reported computed value of 27.0 kcal/mol and in this regard, the value for the r²scan-3c functional seems quite reasonable.

The transition state mode vibrational vectors are quite similar (entries 3 and 6 shown respectively below), indicating that the PD-H-C and adjacent O-H-O contributions are quite similar, whilst the final third transfer has a smaller contribution. This shows that the three transfers are not exactly synchronous, and hence any tunnelling contributions for the three transfers are unlikely to be the same.

row	Method	r_Pd-H	r_H-C	r_O₁-H, Å	r_H-O₂	r_H-O₂	r_H-O₃	α Pd-H-C, °	ν_i	ΔG
1	B3LYP/6-311+G(d,p)/Lanl2dz gas phase^†	1.694	1.297	1.269	1.158	1.153	1.273	166.1	1306	27.0
2	B3LYP/6-311+G(d,p)/Lanl2dz gas phase^‡[cite]10.14469/hpc/15570[/cite],[cite]10.14469/hpc/15572[/cite]	1.677	1.303	1.258	1.151	1.141	1.278	177.1	940	18.6
3	B3LYP/6-311+G(d,p)/Lanl2dz/ SCRF=methanol[cite]10.14469/hpc/15574[/cite],[cite]10.14469/hpc/15575[/cite]	1.736	1.241	1.295	1.132	1.181	1.229	177.2	979	18.3
4	MN15L/6-311+G(d,p)/Lanl2dz/ SCRF=methanol[cite]10.14469/hpc/15589[/cite],[cite]10.14469/hpc/15575[/cite]	1.703	1.283	1.360	1.096	1.082	1.383	143.8	497	25.8
5	MN15L/Def2-TZVPP/ SCRF=methanol	1.633	1.363	1.306	1.120	1.073	1.405	151.3	935	28.0
6	r²scan-3c/Def2-mTZVPP/ SCRF=methanol	1.639	1.360	1.212	1.199	1.129	1.290	136.7	985	21.3

Before a transition state model can be used to infer the KIE for isotopic substitution, it has to be tested against e.g. crystal structures and variation in more accurate basis sets and density functionals. The geometry of the transition state should also be optimised to high accuracy. Whether the KIE reported (~2440) would survive modelling at these more accurate levels remains to be seen. Or indeed whether such an exceptionally high value is directly related to the synchrony of the three hydrogen transfer shown above (“triple hydrogen tunneling”).

^♥The largest value I know of that has been claimed for a KIE is the phenomenal value of ~10¹⁶[cite]https://doi.org/10.1016/j.proeng.2017.03.024[/cite] ^†[cite]10.1021/acscentsci.5c00943[/cite] SI Table S7 etc.

November 21, 2025

Short B-H…H-O Interactions in crystal structures – a short DFT Exploration using B3LYP+D4 and r2scan-3c

In the previous post,[cite]10.59350/x5k75-t2m40[/cite] I was commenting that the transition state for borohydride reduction of a ketone contained some close contacts between the hydrogen of the borohydride and the hydrogen of water. A systematic search of the CSD reveals a modest number of such contacts have been observed in crystal structures (Table). Since it is always good to have a reality check for constructed transition states, here I take a look at some of compounds showing the closest H…H contacts in the experimental database of structures.

The DFT procedures I used to calculate the geometries of the examples tabled below[cite]10.14469/hpc/15566[/cite] were

the classical B3LYP/Def2-TZVPP method, but enhanced with the D4 dispersion correction – the latter developed as a successor to the often used D3+BJ predecessor.
The DFT method named r²scan-3c, a composite described by its developers as the “Swiss army knife of DFT methods” and “r²SCAN-3c Works Well On Everything”.[cite]10.26434/chemrxiv.13333520.v2[/cite] rather grandly quotesThe specific features are described as “The unaltered r2SCAN functional is combined with a tailor-made triple-ζ Gaussian atomic orbital basis set as well as with refitted D4 and geometrical counter-poise corrections for London-dispersion and basis set superposition error“.[cite]10.1063/5.0040021[/cite] The method scales efficiently to several hundred atoms. Results for both types of DFT calculation are collected here.[cite]10.14469/hpc/15566[/cite]

Table. Calculated and observed BH…HO interactions
Molecule	r_H…H using B3LYP+D4/ Def2-TZVPP	r_H…H using r2SCAN-3c/ Def2-mTZVPP	Crystal structure, r_H…H, Å	angle, BH…H, °	angle, OH…H, °
JATMUN	1.531	1.544	1.519 (1.513) [cite]10.1016/j.jorganchem.2005.02.037[/cite],[cite]10.5517/cc8gcz1[/cite]	95.1	158.6
OLEVIL	1.612	1.654	1.806 (1.67) [cite]10.1016/j.tetlet.2010.12.106[/cite],[cite]10.5517/ccvkd70[/cite]	107.7	171.4
OLEVEH	1.623	1.666	1.857 (1.757) [cite]10.1016/j.tetlet.2010.12.106[/cite],[cite]10.5517/ccvkd6z[/cite]	105.1	176.6
OWUKID	1.685	1.770	1.871 (1.778) [cite]10.1002/cplu.202000427[/cite],[cite]10.5517/ccdc.csd.cc252bp9[/cite]	131.7	157.8
SASVAS	1.757	1.872	1.901 (1.807) [cite]10.1002/chem.200902915[/cite],[cite]10.5517/cct86qz[/cite]	94.3	155.0
BOTFOJ	1.741	1.810	1.976 (1.856) [cite]10.5517/ccdc.csd.cc2k93ww[/cite]	129.6	171.0
MOPXOG	1.888	1.966	1.990 (1.872) [cite]10.1021/om5005989[/cite],[cite]10.5517/cc136fzq[/cite]	95.4	159.5
FOLREF	1.881	1.966	1.998 (1.908) [cite]10.1021/ol500970x[/cite],[cite]10.5517/cc12tj5l[/cite]	110.2	160.8

OLEVIL@r2SCAN-3c
OLEVIL @B3LYP+D4

In general, the B3LYP+D4 method predicts slightly shorter H…H contacts than does r²SCAN-3c. Comparison with experiment is tricky, since OH and BH distances obtained directly from a classical crystal structure refinement tend to emerge as ~0.1A too short (see eg [cite]10.59350/5dy8w-0zs92[/cite] for more detailed discussion). A simple correction for these values is shown in parentheses in the table above. However, given that the angle subtended at the hydrogen atom varies enormously, this correction may too simplistic. Better would be if the original crystallographic data could be re-refined using the non-spherical atom model model described in [cite]10.59350/5dy8w-0zs92[/cite] A provisional conclusion without such treatment might be that r²SCAN-3c is somewhat more accurately predicting the H…H distances.

OLEVIL is interesting because it contains two OH groups, with only one interacting with a proximate BH group. Calculating (r2scan-3c) ν_OH gives values of 3477 cm^-1as perturbed by the close BH vs a BH unperturbed value of 3801 cm^-1 or Δν 324 cm^-1. The corresponding values using B3LYP+D4 are 3700 vs 3803, Δν 103 cm^-1. This large difference in perturbation predicted by these two DFT methods could be easily tested experimentally. Unfortunately, the experimental information reported for this compound[cite]10.1016/j.tetlet.2010.12.106[/cite] does not contain the OH stretching values, which might have been a good test of the claim that “r²SCAN-3c Works Well On Everything”.

This selection of compounds illustrates how one aspect of a transition state can be given a reality check by comparing a key interaction with experimentally determined crystal structures.

October 27, 2025

The mechanism of borohydride reductions. Part 2: 4-t-butyl-cyclohexanone – Dispersion induced stereochemistry.

Part one of this topic was posted more than ten years ago.[cite]10.59350/aqrgh-jw887[/cite] I clearly forgot about it, so belatedly, here is part 2 – dealing with the stereochemistry of the reduction of tert-butyl-cyclohexanone by borohydride in water. The known stereochemistry is nicely summarised in this article, along with an extensive history of possible explanations of the reasons for the stereochemical preference.[cite]10.1080/10610270701268815[/cite] Put simply, the hydride nucleophile attacks the carbonyl from an axial rather than equation direction with a ratio of 10:1 (ΔΔG 1.37 kcal/mol). So does the model I previously proposed[cite]10.59350/aqrgh-jw887[/cite] support this and give any indication of why the stereochemistry is axial?

The calculated transition states are shown below (click on image to get interactive 3D model). Note also the unusual calculated B-H…H-O hydrogen bonded distances of ~1.8Å. A search of the CSD (crystal structure database) reveals surprisingly few such examples, but one interesting one is as short as ~1.73Å[cite]10.1021/ja982959g[/cite]

DFT calculations were conducted using Gaussian 16 at the B3LYP/Def2-TZVPP/SCRF=water level[cite]10.14469/hpc/15559[/cite] for transition state location and energies (including D3+BJ dispersion) and using ORCA 6.1[cite]10.1002/wcms.70019[/cite] for calculating D4 dispersion energies.[cite]10.1039/D0CP00502A[/cite]

Lithium Borohydride

Sodium borohydride

As well as computing the free energy difference ΔΔG^‡ between the two transition states, I also looked at the total dispersion energy contributions to these energies at the third generation D3+BJ and the fourth generation D4 levels, shown below as the difference between the axial and equatorial transition states.

Counter ion ΔD3 dispersion ΔD4 dispersion ΔG^‡_axial ΔG^‡_equatorial ΔΔG_h^†

Lithium borohydride 1.15 0.932 12.13 12.93 0.80

Sodium borohydride 1.22 0.745 13.31 14.16 0.85

^†Harmonic energies.

The calculations reveal that using either Li or K as the counterion to borohydride, the axial transition state is lower in energy than the equatorial, by around 0.85 to 0.80 kcal/mol in total free energy. These values correspond to an axial/equatorial ratio of around 3.9-4.2:1, a little bit lower than observed experimental value of 10:1.[cite]10.1080/10610270701268815[/cite] The D3+BJ and D4 dispersion energies follow the same trend, with the suggestion that the Li ion is slightly more stereoselective than the Na ion, perhaps due to the more compact nature of the transition state.

So we might conclude that the stereochemical preference for axial hydride delivery to tert-butyl-cyclohexanone could be explained entirely by the differing dispersion energy contributions of the two transition states. This in one way is a deeply unsatisfactory explanation, since the dispersion energy difference is the total sum of many individual dispersion contributions and is largely unpredictable until calculated. Chemists like simple rules and this is not apparently amenable to such a simple rule. Indeed perhaps the rule should simply be to always compare the computed dispersion energies of two (or more) isomeric transition states before seeking other explanations!

A source of error in the calculated free energies could be the treatment of the vibrational modes as harmonic. A quasi-harmonic treatment is available[cite]10.12688/f1000research.22758.1[/cite] which should give a more reliable estimate of the relative free energy. Using this approach (for settings see [cite]10.14469/hpc/15565[/cite]) ones gets the values below. Although the discrimination is slightly reduced, the overall prediction of axial hydride attack is still confirmed.

system> hG qhG

LiBH4 ax -884.325907 -884.328233 (ΔΔG -0.62)

LiBH4 eq -884.324625 -884.327246

NaBH4 ax -1039.096789 -1039.100009 (ΔΔG -0.78)

NaBH4 eq -1039.095445 -1039.098771

October 21, 2025
Alternative reactions of the N≡N “triple bond” in a nitric oxide dimer: forming the trimer N₃O₃.

In the previous post[cite]10.59350/rzepa.29626[/cite] I mooted the possibility that a high energy form of the dimer of nitric oxide 1 might nonetheless be able to be detected using suitable traps (such as hydrogenation or cycloaddition). However, an interesting alternative is that this species could be trapped by nitric oxide itself. According to [cite]10.1016/0022-1902(65)80196-8[/cite] in an article entitled “Decomposition of nitric oxide at elevated pressures” the rate of this termolecular reaction 3NO → N₂O + NO₂ are said to obey third order kinetics. One plausible mechanism for this process is shown below.

So the question now arises as to whether this process would always intervene to prevent any trapping of species 1. Time for some calculated energies.[cite]10.14469/hpc/15549[/cite] The reaction occurs in two steps, TS1, forming Int, which then gives the product via TS2. The combined IRC for these processes is shown below,^† and the total energy barrier from 1 + NO for this process is small (~7 kcal/mol).

The NN bond length along the course of these two steps is shown below,^† indicting that it starts and ends with an ~N≡N bond, as implied in the scheme above.

The dipole moment response is shown below.^†

Animation for the IRC for TS1 is shown below:

The intermediate, N₃O₃ has a spin density covering all six atoms (click on image below to see interactive 3D model) so is difficult to represent by e.g. a single valence bond structure as shown above.

The IRC for TS2 is animated below:

The high value for ΔG^‡ ~39.5 kcal/mol arises in large part because of the loss of entropy in reducing three molecules to one at the transition states – and can be compared with e.g. the value of ΔG^‡ 31.0 kcal/mol previously quoted for the bimolecular dimerisation of two molecules of NO.[cite]10.59350/rzepa.29429[/cite]

The story is not quite over yet, since other conformations for this reaction are still being explored. Currently it appears there is only a small free energy window between forming species 1 and it reacting further in the manner 3NO → N₂O + NO₂ as shown above, so it seems unlikely at this stage that 1 could be easily detected.

^† The slight discontinuity visible at ~IRC +1 is due to a conformational change between the end point of the IRC for TS1 and the start point for TS2.

DOI: 10.59350/rzepa.29665

September 3, 2025
Hydrogenating the even more mysterious N≡N triple bond in a nitric oxide dimer.

Previously[cite]10.59350/rzepa.29429[/cite] I looked at some of the properties of the mysterious dimer of nitric oxide 1 – not the known weak dimer but a higher energy form with a “triple” N≡N bond. This valence bond isomer of the weak dimer was some 24 kcal/mol higher in free energy than the two nitric oxide molecules it would be formed from. An energy decomposition analysis (NEDA) of 1 revealed an interaction energy[cite]10.14469/hpc/15468[/cite] of +4.5 kcal/mol for the two radical fragments, compared to eg -27 kcal/mol for the equivalent analysis of the N=N double bond in nitrosobenzene dimer[cite]10.14469/hpc/15444[/cite] So here I take a look at another property of N≡N bonds via their hydrogenation energy (Scheme), mindful that the dinitrogen molecule requires forcing conditions to hydrogenate, in part because of the unfavourable entropy terms (See Wiki and also here^‡ for a calculation of ΔG₂₉₈).

Calculations at the ωB97XD/Def2-TZVPP/SCRF=water level[cite]10.14469/hpc/15516[/cite] that whilst hydrogenation of the triple bond in N₂ is strongly endo-energic, the same process for molecule 1 is exo-energic (ΔΔG -26.32 kcal/mol). The direct product is a zwitterion, but presumed rapid proton transfer to a neutral form 2 increases exo-energicity. Whilst the second hydrogenation step of N₂ is exo-energic, the equivalent second step for 1 to give 3 is now mildly endo-energic. Overall however, the thermodynamic energies of these two types of triple bond hydrogenation could not be more different.

So forming a N≡N triple bond by forcing two nitric oxide molecules to dimerise (using high pressure) in water produces a system where hydrogenation of that “difficult” N≡N bond is made very much easier thermodynamically. Time for an experiment?^♥

^‡This site reports a gas phase experimental value for ΔG -8.1 kcal/mol at 298K for this equilibrium, although the pressure is not given. The calculated value shown in the scheme above (-20.1 kcal/mol) is for 298K and 1 atm for a model using water as solvent – which might be expected to differentially solvate the product ammonia and hence promote the reaction. In the limit of low pressure (0.0001M)[cite]10.12688/f1000research.22758.1[/cite] this reduces to -13.0 kcal/mol, increases to -26.6 kcal/mol at 10M and becomes -14.3 kcal/mol at 10M/800K, illustrating how higher pressures make the reaction more exo-energic and higher temperatures less exo-energic. This was of course the problem solved in the Haber process of finding the sweet spot between pressure and temperature.

^♥Perhaps not, given the report that at high pressures, nitric oxide can become explosive.[cite]10.1016/0022-1902(65)80196-8[/cite]

August 25, 2025
The spin-offs from adding citations to blog posts.
I started adding citations to my blog posts around 2012 using the Kcite plugin.^‡ Rogue Scholar is a service that monitors registered blog sources and adds all sorts of value to the original post, including identifying such citations and creating a list of them.

I show the results for the previous blog[cite]10.59350/rzepa.29523[/cite] here.

Martin Fenner has just added some interesting new features[cite]10.53731/yjq4w-5yr32[/cite] which I thought would be useful to share with you here.
1. If you go to the Rogue Scholar archive of the post and scroll down to the References list, then click on the title of any of the references, you will get a list of all Rogue Scholar posts citing that reference: https://rogue-scholar.org/search?q=doi:10.1038/sdata.2016.18+references:10.1038/sdata.2016.18+citations:10.1038/sdata.2016.18
2. If you click on the author name in any of the entries from the previous search, you get a list of all the posts published by that person.
  https://rogue-scholar.org/search?q=orcid:0000-0002-8635-8390&sort=newest
I think this idea of adding citations to a blog post can result in a considerably enhanced discovery process – if only you could do this with journals themselves!

^‡ This is temporarily not functional due to a php update on the site. I hope to get it working again soon. Update. Thanks to Martin Fenner, the Kcite plugin is working again at version 1.7.11 and upwards.[cite]10.53731/326tr-95k32[/cite]
August 19, 2025
More on rescuing articles from a now defunct early pioneering example of an Internet journal.
Two years ago, I posted on the topic “Internet Archeology: reviving a 2001 article published in the Internet Journal of Chemistry (IJC)”.[cite]10.59350/xqerh-wam97[/cite] The IJC had been founded in 1998,[cite]10.1080/00987913.2000.10764578[/cite] in part at least to “re-invent” the scholarly journal by elevating research data to being a more integrated part of the overall article, rather than as the previously conventional addendum of SI (Supporting Information)^‡. IJC was in one sense following on from an earlier such project dating from 1995[cite]10.1080/13614579509516846[/cite] by taking it to the next level. Sadly, the pioneering IJC journal had gone off-line in 2004 and the content for around 100 articles was thought lost. It happened that I still retained the original source and associated data for one article of mine and my post[cite]10.59350/xqerh-wam97[/cite] described how I managed to get it back into more or less full working order. Now Egon Willighagen[cite]https://doi.org/10.59350/2ss5b-jpr33[/cite] has cleverly found a way of rescuing many more of these lost articles, thanks to various Web-based infrastructures:
1. From 1996 as the Internet archive (using a query such as https://web.archive.org/web/*/http://www.ijc.com/abstracts/*),
2. From 2012, WikiData (see https://www.wikidata.org/wiki/Q27211732)
3. Also from 2012, ORCID (Resarcher and collaborator) profiles. Some reserchers had the foresight (alas not me) to link their by then defunct IJC articles to their new ORCID profiles.
I link here to some examples of rescued articles as shown on Egon’s blog. I eagerly look forward to seeing what else is to come using such tools!

^‡ By around 2005, a clearer separation between the journal (the “story” or research narrative) and its associated research data was being seen as the way forward, with the data now being placed in a data repository (or Wikidata) separate from the journal, with added descriptive metadata to help make it a stand-alone object and this new entity to now be cited in the journal article (and bidirectionally the article from the data) using a persistent identifier – initially as a Handle, then as a DOI.[cite]10.1021/ci7004737[/cite] FAIR data as a concept had started to emerge from these developments, being formalised around a decade later in 2016.[cite]10.1038/sdata.2016.18[/cite]

This post has DOI:10.59350/rzepa.29523
August 19, 2025

Counter ion	ΔD3 dispersion	ΔD4 dispersion	ΔG^‡_axial	ΔG^‡_equatorial	ΔΔG_h^†
Lithium borohydride	1.15	0.932	12.13	12.93	0.80
Sodium borohydride	1.22	0.745	13.31	14.16	0.85

system>	hG	qhG
LiBH4 ax	-884.325907	-884.328233 (ΔΔG -0.62)
LiBH4 eq	-884.324625	-884.327246
NaBH4 ax	-1039.096789	-1039.100009 (ΔΔG -0.78)
NaBH4 eq	-1039.095445	-1039.098771

The even more mysterious N≡N triple bond in a nitric oxide dimer.

Previously, I pondered about the strange N=N double bond in nitrosobenzene dimer[cite]10.59350/rzepa.29383[/cite] as a follow up to commenting on the curly arrow mechanism of the dimerisation.[cite]10.59350/rzepa.28849[/cite] By the same curly arrow method, one can produce the below, showing how the simpler nitric oxide radical could potentially dimerise to a species with a N≡N triple bond!^† This involves a total of six electrons entering the N-N region, and hence raises the question of whether these all move in a single concerted/synchronous bond forming reaction, or whether they might go in (asynchronous) stages. Here are some calculations[cite]10.14469/hpc/15483[/cite]) which might shed some light on this aspect.

The structure[cite]10.1021/ja00381a052[/cite] of a nitric oxide dimer was shown in 1982 to have a very long (rather than short) N-N bond length of 2.237Å and a theoretical analysis[cite]10.1016/j.cplett.2007.11.084[/cite] showed it to be a weak complex with a very complex wavefunction showing multi-reference character.

Firstly, an IRC-based reaction path (method: uωB97XD, scrf=(cpcm,solvent=water) guess=(mix,always) def2tzvpp to allow either an open shell biradical to form and also to encourage any ion pair formation). As you can see, the (total) energy goes up to a very shallow transition state (with a tiny reverse barrier) to form a biradical with <S²> 0.628. This species, as noted existing in a very shallow energy well, has an N-N bond length of 1.725Å.

The bonding for this species is complex (analysis for a later post), but the calculated biradical spin density below shows the unpaired electrons are in the π-system (click on the image to get a 3D rotatable model).

Further contraction of the N-N length results in an IRC energy potential to a transition state with a N-N length 1.294Å across a further barrier of ~12 kcal/mol (ΔE; ΔG 13.6 kcal/mol). The overall barrier from two nitric oxide molecules is ΔG 31.0 kcal/mol with the overall thermochemistry summarised in the table. Basically, this barrier is unsurmountable at normal temperatures and the reverse barrier of ΔG 6.7 kcal/mol ensures that the N≡N triple bonded species shown above is not likely stable and will not be observed experimentally.^‡ However this product is NOT a biradical but a normal closed shell singlet molecule.[cite]10.14469/hpc/15474[/cite]

So to answer my first question, the six electrons appear to move in two stages, firstly two electrons form a weak N-N bond and then a further four electrons contract this to a triple bond. Their motion is effectively concerted, but asynchronous.

Species	ΔG	ΔH	ΔΔG	ΔΔH	T.ΔS	r_NN, Å	<S²>	DOI
2*Nitric oxide	-259.83494	-259.78839	0.0	0.0	29.2	∞	0.753	15472[cite]10.14469/hpc/15472[/cite]
Singlet biradical	-259.80716	-259.77615	17.4	7.7	19.5	1.725	0.628	15476[cite]10.14469/hpc/15476[/cite]
Triplet biradical	-259.80865	-259.77672	16.5	7.3	20.0	1.779	2.016	15475[cite]10.14469/hpc/15475[/cite]
Singlet TS	-259.78550	-259.75579	31.0 (13.6)	20.5 (12.8)	18.6	1.294	0.000	15483[cite]10.14469/hpc/15483[/cite]
Singlet N≡N dimer	-259.79614	-259.76693	24.3 (7.8)	13.5	18.3	1.114	0.000	15467[cite]10.14469/hpc/15467[/cite]

Now for a NEDA energy decomposition analysis[cite]10.1063/1.466432[/cite]

Electrical (ES+POL+SE) :  -9414.608
   Charge Transfer (CT) :  -1363.597
       Core (XC+DEF-SE) :  10782.725                      
  Total Interaction (E) :      4.520 kcal/mol.

Normally NEDA total interaction energies are -ve, but this one is positive! So the triple bond dissociation energy is not merely small, but actually negative. That is a weak triple bond and as the title implies, a very mysterious bond. In some aspects however it is conventional. Thus calculated r_NN 1.114Å and ν_NN 2604 cm^-1. However partial occupancies of NBO antibonding BD* orbitals results in a calculated Wiberg bond order of only 1.01; there is still a great deal of mystery left about this species! Probably what is fairly certain is that the closed shell single-reference wavefunction used here is not appropriate for a full explanation and more complex multi-reference procedures would have to be used to get a more complete picture of this strange non-existing^‡ little molecule. It may even be that such procedures remove the small reverse barrier noted above, thus preventing the molecule from even existing in an energy well.

^†This species does not appear to have been previously discussed or suggested, according to SciFinder/CAS.
^‡Might it exist at very high pressures in water?

To find all blog posts authored here, along with their DOIs, try https://rogue-scholar.org/search?q=orcid:0000-0002-8635-8390&sort=newest

August 18, 2025

Energy decomposition analysis of hindered alkenes: Tetra t-butylethene and others.
In the previous post,[cite]10.59350/rzepa.29383[/cite] I introduced the N=N double bond in nitrosobenzene dimer, arguing that even though it was a formal double bond, its bond dissociation energy made it nonetheless a very weak double bond! This was backed up by a technique known as energy decomposition analysis or EDA. Here I use a variant of this method known as NEDA to look at some other strained alkenes, including the famously non-existent tetra t-Butyl ethene.

The NEDA procedure gives a fragment interaction energy (decomposing it into fundamental quantum mechanically derived energies if required) with respect to a reference state for the fragments. In this case, the fragments are obtained by cutting the double bond, resulting in triplet state carbenes as the reference state. The calculations (B3LYP+GD3+BJ/Def2-TZVPP) are available here.[cite]10.14469/hpc/15463[/cite]
1. Compound 1, a relatively unstrained alkene, ΔE = -177.0 kcal/mol, R_CC 1.341Å
2. Compound 2 (PUVQUE, [cite]10.1107/S0108270198099247[/cite], [cite]10.5517/cc4cx7m[/cite]), ΔE = -164.3 kcal/mol, R_CC 1.362Å, CC torsion 16.5°
3. Compound 3 (CUBVOK, [cite]10.1016/S0040-4039(00)98504-6[/cite]) ΔE = -167.9 kcal/mol, R_CC 1.351Å, CC torsion 9.2°
4. Compound 4 (currently unknown) ΔE = -135.8 kcal/mol, R_CC 1.380Å, CC torsion 54.5°
The NEDA interaction energy is directly proportional to both the CC bond length and the C-C=C-C torsion angle. What is interesting however is the large interaction energy gap in ΔE between the two known hindered alkenes (2 and 3) and the unknown tetra-t-butyl ethene 4. It seems moving from say compound 2 by converting the two iso-propyl substituents to full t-butyl ones is just too large a change to bridge. Unless one day isolated as a very very unstable species, compound 4 seems destined not to exist!

This post has DOI: 10.59350/rzepa.29410
August 13, 2025