Triplets in the cradle: ultrafast dynamics in a cyclic disulfide

Abstract

The effect of spin–orbit coupling on the “Newton's cradle”-type photodynamics in the cyclic disulfide 1,2-dithiane (C₄H₈S₂) is investigated theoretically. We consider excitation by a 290 nm laser pulse and simulate the subsequent ultrafast nonadiabatic dynamics by propagating surface-hopping trajectories using SA(4|4)-CASSCF(6,4)-level electronic structure calculations with a modified ANO-R1 basis set. Two simulations are run: one with singlet states only, and one with both singlet and triplet states. Comparison of the two simulations suggests that the presence of triplet states depletes the singlet state population, with the net singlet and triplet populations at long times tending towards their statistical limit. Crucially, the triplet states also hinder the intramolecular thiyl radical recombination pathway via the efficient intersystem crossing between the singlet and triplet state manifolds.

---

1 Introduction

Disulfide linkages are a common chemical motif in proteins and are usually formed via covalent coupling of thiol functional groups situated on spatially adjacent cysteine amino acids.¹ As a consequence of the thermal and photolytic stability of cross-linked cysteine thiol groups, and their role as a conformational lock, disulfide linkages are important for stabilising tertiary and quaternary protein structures,^2–4 helping to maintain biological function in the native environment.

Conversely, linear alkyl disulfides are unstable with respect to heat and ultraviolet (UV) light both in solution and in the gas phase.^5–9 Cyclic disulfides, in which the disulfide bond is part of a ring, exhibit greater thermal and photolytic stability in both experimental and computational studies, although initial disulfide cleavage is still observed.^10–19 This enhanced stability has been attributed to the carbon backbone; following initial sulfur–sulfur photoinduced homolysis, the backbone tethers the two dithiyl radicals so that they may return to sufficient proximity for radical recombination to occur.

The molecule 1,2-dithiane (C₄H₈S₂) is a six-membered hydrocarbon ring in which two adjacent methylene groups have been replaced by sulfur atoms, as can be seen in Fig. 1. It serves as a convenient model for studying the photochemistry of constrained disulfide systems such as those found in proteins or other structurally-constrained environments. Its small size compared to disulfides in proteins—where polar and charged functional groups, as well as sterically bulky neighbouring functional groups, may affect the S–S bond-breaking process—renders the molecule suitable for accurate computational studies of the sub-picosecond dynamics, and allows the effect of UV light on the S–S moiety to be studied in isolation.

Fig. 1 Schematic representation of the “Newton's cradle” motion in photoexcited 1,2-dithiane, whereby the molecule repeatedly rotates about the carbon framework to bring the two sulfur atoms back into close proximity before the molecule springs open to reform the linear biradical.

The photoinduced dynamics of 1,2-dithiane upon UV excitation, and the effect of molecular geometry on the observed photodynamics of disulfides generally, have both been subject of previous experimental and computational studies.^10–20 Studies on 1,2-dithiane in particular have found that when the photodynamics is initiated on the first excited singlet electronic state S₁, ultrafast S–S homolysis occurs on a sub-100 fs timescale and results in ring-opening.^12–15 The dissociated sulfur termini remain tethered by the carbon backbone and, via torsion about the carbon backbone, encounter each other again “on the other side” of the molecule. This results in an oscillatory ring-opening/closing motion with an approximate period of 350 fs. The characteristic motion, colloquially referred to as a molecular “Newton's cradle”, is schematically illustrated in Fig. 1. Each time the two sulfur atoms return to each other, radical recombination of the S–S bond can occur on the S₀ ground electronic state. Subsequently, the molecule can spring open again or remain in the vibrationally hot cyclic ground state. The Newton's cradle dynamics is gradually damped by recombination and energy dispersion, with the latter being intramolecular vibrational energy redistribution (IVR).¹⁴

With one exception,¹³ theoretical studies have considered only singlet excited states in computational nonadiabatic dynamics studies of 1,2-dithiane.^11,14–17 However, due to the presence of the sulfur atoms, moderate spin–orbit coupling (SOC)^21,22 may be expected to play a role in the photodynamics of 1,2-dithiane. Intersystem crossing (ISC) is known to occur at short timescales in sulfur-containing molecules, allowing ISC to compete with internal conversion (IC).^23–26 Indeed, ISC has been reported even in systems where “heavy” atoms are absent.²⁷ Furthermore, the previous study that accounted for triplet states in simulations of 1,2-dithiane evidences efficient ISC from singlet states to triplet states in the short-time regime—sub-70 fs—following photoexcitation.¹³ In the present work, the effect of the triplet states on the photodynamics in 1,2-dithiane—both with respect to electronic state populations and nuclear dynamics—is considered by comparing otherwise identical nonadiabatic dynamics simulations with and without triplet states. The simulations are carried out using Tully's fewest-switches surface-hopping method,²⁸ and both simulations are propagated for 1 ps to allow both short-time and intermediate-time processes to be examined. Finally, X-ray scattering patterns for the two resulting trajectory ensembles are calculated in order to compare how any differences in nuclear dynamics between the two sets may manifest in ultrafast scattering experiments.^29–32

2 Methods

2.1 Electronic structure calculations

The electronic structure calculations are performed at the complete active space self-consistent field (CASSCF) level of theory³³ using the OpenMolcas v23.02 software package³⁴. As discussed below, we employ a (6,4) active space for the electronic structure used in the simulations. This active space, shown in Fig. 2, comprises 6 electrons distributed across 4 orbitals; the orbitals include S–S σ-bonding and σ*-antibonding orbitals (σ_S–S and ), and the two in-phase and out-of-phase combinations of lone-pair orbitals on the adjacent sulfur atoms, labelled n^u_S and n^g_S, respectively. State-averaging is performed over the four lowest energy electronic states in each spin-manifold, i.e. SA(4) for singlets only, and SA(4|4) for both singlets and triplets. For the electronic structure calculations involving triplets, the SOC Hamiltonian is calculated in the adiabatic-state basis. The spin–orbit states and associated spin–orbit coupling matrix elements (SOCMEs) are then found by diagonalising the total Hamiltonian. This is achieved with the restricted active-space self-interaction (RASSI) program within OpenMolcas, which uses atomic mean-field integrals (AMFIs).^35–37

Fig. 2 Orbitals used in the active spaces for the electronic structure calculations in 1,2-dithiane. The complete set used in the larger (10,8) active space corresponds to all orbitals shown. The four orbitals included in the smaller (6,4) active space are the subset in the grey-dashed box on the right.

In the simulations, the ANO-R1 basis set is employed.³⁸ The basis is truncated by removing the polarisation functions on all hydrogen atoms. This is justified because the hydrogen atoms are not expected to play a significant role in the ring-opening photodynamics of 1,2-dithiane upon photo-excitation at 290 nm, as supported by previous studies.^12–17 The truncated ANO-R1 basis is denoted ANO-R1(t) for brevity. The resulting electronic structure method is thus SA(4|4)-CASSCF(6,4)/ANO-R1(t), and is used henceforth unless otherwise indicated. Additional electronic structure benchmarks are provided in the SI, in which: (i) calculations using ANO-R1 and ANO-R1(t) are compared against a range of basis sets; and (ii) calculations are compared employing both the (6,4) active space and a larger (10,8) active space in which two sets of C–S σ and σ* orbitals are included in the active space—see Fig. 2.

Molecular geometries are optimised using SA(4|4)-CASSCF(6,4)/ANO-R1(t). The ground state and the first singlet excited state equilibrium structures, denoted S₀ (min) and S₁ (min), are confirmed to correspond to true minima by a frequency calculation (see SI for details). The minimum-energy conical intersection (MECI) between the S₀ and S₁ states, denoted MECI (S₀,S₁), is also located. The electronic structure calculations are benchmarked at a range of molecular geometries determined by linear interpolation in internal coordinates (LIIC) between the optimised geometries, as shown in Fig. 3. The LIIC procedure involves performing a single-point CASSCF calculation at each interpolated intermediate geometry to yield the energies of multiple electronic states. It is worth noting that the linear interpolation of molecular geometries means that all internal coordinates may change across the LIIC pathway; that is, the LIIC does not strictly correspond to a minimum energy path between the optimised geometries.

Fig. 3 Summary of static results in 1,2-dithiane. (a) Potential energy cuts for the four lowest-energy singlet states S₀–S₃ and triplet states T₁–T₄ along the LIIC pathway connecting the S₀ minimum, labelled S₀ (min), to the S₁ minimum,labelled S₁ (min), and then onwards to the minimum energy conical intersection between the S₀ and S₁ states, MECI (S₀,S₁). The excitation window used for dynamics is shown by the grey horizontal shaded region. (b) The dominant character, the corresponding CI-coefficient, and the vertical excitation energy for each electronic state at the S₀ (min) molecular geometry. (c) The optimised molecular geometries S₀ (min) (left), S₁ (min) (centre), and MECI (S₀,S₁) (right).

Electronic structure calculations with the (6,4) and (10,8) active spaces shown in Fig. 2 are compared along the LIIC pathways, noting that the differences in geometries used to create the LIIC are negligible at either level of theory (see SI). The electronic structure calculations are found to be in quantitative and qualitative agreement across the LIIC. The validity of the SA(4|4)-CASSCF(6,4)/ANO-R1(t) electronic structure calculations used in the simulations is further supported by benchmarks against other methods such as extended multi-state complete active space second-order perturbation theory, XMS-CASPT2, and multireference configuration interaction, MRCI;^39,40 these results are also included in the SI.

2.2 Dynamics

The nonadiabatic dynamics of 1,2-dithiane is simulated with the trajectory surface hopping (TSH) method, using Tully's fewest switches algorithm,²⁸ as implemented in the SHARC software package (version 2.1).^41–43

From one set of 8000 initial conditions sampled from a Wigner distribution at the S₀ (min) structure, two ensembles are prepared: one containing only the first four lowest-energy singlet states (“singlet-only”), and another containing the same set of singlet states plus the first four lowest-energy triplet states (“triplet-inclusive”).^44,45 We run 384 trajectories in each ensemble. An excitation window based on a 290 nm UV pump-pulse with 30 fs full-width-half-maximum (FWHM) pulse duration is used to excite the initial conditions for the dynamics. This window corresponds to 4.25–4.31 eV; however, a positive 0.21 eV shift is applied to this excitation window (i.e., 4.46–4.52 eV) to account for the fitting of the simulated spectrum to an experimental gas-phase absorption spectrum obtained under ambient pressure at room temperature.^46,47 The choice of these specific parameters is motivated by the anticipation that such parameters may be used in future experiments.

The UV absorption spectrum is simulated using the nuclear ensemble approach.⁴⁸ This entails summing over the absorption spectra at all initial geometries, where the absorption spectrum at each initial geometry is the convolution of the line spectrum by a Gaussian function with FWHM = 0.2 eV to represent the experimental energy resolution.

In both the singlet-only and triplet-inclusive ensembles, the trajectories are propagated for 1 ps with a nuclear timestep of 0.5 fs and an electronic timestep of 0.02 fs. Dynamics are performed in the diagonal representation to allow for favourably localised couplings between electronic states, which in principle yields more accurate dynamics compared with those performed in the diabatic and adiabatic bases.^42,49 A local diabatisation scheme is used to calculate electronic time-derivative couplings due to the greater efficiency compared with explicit calculations of nonadiabatic coupling matrix elements (NACMEs), particularly when the molecule adopts ring-opened geometries leading to dense manifolds of electronic states.⁵⁰ The speed of the quantum chemistry calculations at each timestep is further improved by only calculating the gradient of a classically unoccupied adiabatic state when the difference in energy between this state and the adiabatic state on which the dynamics is currently being propagated on is less than 0.5 eV. This is justified on account of the inverse scaling of transition probabilities with respect to the energy gap between states; population transfer is less likely between more widely spaced states. The energy-difference based decoherence correction scheme developed by Granucci and Persico is implemented using the recommended parameter of 0.1 Hartree.⁵¹ Following a surface hop, nuclear velocities are rescaled along the full nuclear velocity vector to preserve the total energy. The total energy and changes in energy between consecutive timesteps is closely monitored for all trajectories; energy conservation violations falling outside a 0.2 eV threshold are flagged, and those trajectories are terminated at the timestep where the problem arises.

2.3 X-ray scattering simulations

For both sets of trajectory ensembles, ultrafast X-ray scattering (UXS) patterns are calculated using the independent atom model (IAM) approach.^30,52 The IAM assumes that the scattering patterns of a molecule can be approximated as a coherent sum of isotropic scattering amplitudes of all constituent atoms.⁵³ Any effects in the scattering signal which are caused by changes in electronic distribution resulting from bonding interactions within the molecule are thus not considered.^54–56 In short, the rotationally averaged IAM elastic scattering, I_el(q̃,R̄), for a molecule containing N atoms described by spatial coordinates R̄ = {R₁,…,R_N} is given by³⁰where q̃ = |q̃| is the magnitude of the momentum transfer vector (within the Waller–Hartree approximation),⁵⁷f^X_i(q̃) is the atomic form factor for atom i,⁵⁸ and R_AB = |R_A − R_B| is the distance between atoms A and B. The total X-ray scattering probability, I_tot(q̃, R̄), is then a sum of elastic and inelastic X-ray scattering components:where the second term approximates the inelastic component of the total scattering and involves an incoherent sum of tabulated coherent scattering functions, S_A(q̃).⁵⁸ Thus, by using eqn (2) to calculate the total scattering signal at every molecular geometry along a surface-hopping trajectory for all trajectories in an ensemble, the mean total scattering signals at each time – that is, the mean simulated UXS signal – can be calculated aswith equal weight given to each trajectory R̄_i(t). To evaluate how the UXS pattern varies across time it is convenient to plot the percent difference of the UXS signal, Δ%I(q̃,t), relative the signal at t₀ = 0:

We note that in the IAM, the inelastic component is invariant with respect to molecular geometry, meaning that only the elastic component affects the numerator in eqn (4).^56,59

3 Results and discussion

3.1 Electronic structure of ground and excited states

The LIIC potential energy curves, the electronic state characters and the corresponding CI coefficients at the S₀ (min) geometry, as well as the molecular structures at the S₀ (min), S₁ (min), and MECI (S₀,S₁) geometries, are all shown in Fig. 3. The ground state minimum, S₀ (min), adopts a chair-like conformation similar to cyclohexane. Previous literature suggests that there also exists a local minimum twist-boat structure at approximately 10–20 kJ mol⁻¹ above the chair minimum, which may be accessed via a half-chair transition state lying approximately 50 kJ mol⁻¹ above the ground state minimum.^60,61 In the present calculations, using SA(4|4)-CASSCF(6,4)/ANO-R1(t), the twist-boat energy was found to lie 23 kJ mol⁻¹ above the ground state. Assuming a two-state Boltzmann distribution, and accounting for the vibrational zero point energy of each conformer, only 0.5% of the molecules are expected to adopt the twist-boat geometry at 298 K. Thus, only initial conditions from the chair-like minimum are considered in the following, with further details on the twist-boat conformation provided in the SI.

In both the S₁ (min) and the MECI (S₀,S₁) geometries, the S–S bond is broken, and the two sulfur atoms have radical character. The S–S distance for the S₀ (min) geometry is 2.13 Å, while for the S₁ (min) and MECI (S₀,S₁) geometries, the S–S distances are 3.64 Å and 3.91 Å, respectively. The MECI geometry is thus quite similar to the S₁ minimum, with the root mean square deviation (RMSD) between the two Kabsch-rotated structures equal to 0.0914 Å. In contrast, the RMSD between S₁ (min) and S₀ (min) is 0.3529 Å.

The LIIC pathway shown in Fig. 3 is constructed by interpolation from S₀ (min) to S₁ (min), and then onwards to the MECI (S₀,S₁), i.e. S₀ (min) → S₁ (min) → MECI (S₀,S₁). The main change in molecular geometry along the LIIC pathway is a progressive stretching of the S–S bond; there are only minor changes in other internal coordinates, such as a small rotation in the carbon backbone dihedral angle. Looking at Fig. 3(a), it is evident that the excited S₁–S₃ and T₁–T₄ states are dissociative in the Franck–Condon (FC) region. As the S–S bond stretches, the potential energies of the excited states decrease smoothly, thus providing a barrierless path to the S₁ minimum. Accessing the MECI (S₀,S₁) geometry from S₁ (min) along the LIIC pathway requires surmounting a very small potential barrier—0.03 eV—on the S₁ state at this level of theory. From left to right along the LIIC pathway, the S₀ energy increases. As the S–S distance increases, all potential energy curves considered approach similar asymptotic energies corresponding to a broken S–S bond, with no or very little overlap between the orbitals on the two sulfur atoms.

The dominant electronic state characters for the S₀–S₃ and T₁–T₄ states are listed alongside their configuration interaction (CI) coefficients in Fig. 3(b).† The S₀ state in the FC region is defined predominantly by two electrons each in the σ_S–S, n^g_S, and n^u_S valence orbitals. Given that all excited states considered here populate the S–S σ* orbital, this qualitatively explains the repulsive nature of these states with respect to the S–S stretching coordinate, as seen in Fig. 3. By comparing the state characters in the FC region between the (6,4) and (10,8) active spaces (see SI), we see that the C–S orbitals are unlikely to influence the state characters at energetically accessible geometries considering the excitation window in the simulations, further supporting the use of the smaller (6,4) active space in the simulations.

Finally, all SOCMEs were calculated between the four lowest-energy singlet and four lowest-energy triplet states at geometries along the S₀–S₁ LIIC pathway. Across the LIIC pathway, strong spin–orbit coupling with SOCMEs greater than 100 cm⁻¹ (0.01 eV) is present between both singlet–triplet and triplet–triplet state pairs, as may be anticipated from the Z⁴ scaling of spin–orbit coupling in an atomistic picture. SOCMEs of this magnitude are also consistent with those reported for other disulfide systems.⁶² Such strong SOC effects should be expected to facilitate efficient population redistribution among the triplet states from the singlet manifold. This supports and justifies the need to investigate the effect of including the triplet states in the simulations.

3.2 UV absorption spectrum

The UV absorption spectrum for 1,2-dithiane, calculated at the SA(4|4)-CASSCF(6,4)/ANO-R1(t) level of theory, is shown in the top panel of Fig. 4. The absorption spectrum calculated with the (10,8) active space is presented in the bottom panel to aid with the assignment of the observed features in the experimental absorption spectrum.^46,47 Both the experimental and shifted computational spectra show a broad peak with a maximum centred at approximately 4.37 eV, which corresponds to direct excitation into the S₁ state. For both active spaces, much of the intensity of this peak can be ascribed to a transition. This peak energy is in reasonable agreement with the vertical excitation energy from S₀ to S₁ (ΔE_S1,S0 = 4.50 eV) as shown in Fig. 3.

Fig. 4 Experimental^46,47 (dotted line) and computational UV absorption spectra (shifted by −0.21 eV) for 1,2-dithiane calculated at the SA(4|4)-CASSCF(6,4)/ANO-R1(t) (top) and the SA(4|4)-CASSCF(10,8)/ANO-R1(t) (bottom) levels of theory, respectively. The proposed excitation window, 4.25–4.31 eV is indicated by the shaded light-grey region on both plots. For a discussion of the strong absorption at energies >4.5 eV in the experimental spectrum, see the text.

Noticeable deviations from the experimental gas-phase spectrum occur for both calculated absorption spectra at energies greater than 4.5 eV. A very strong absorption peak is seen in the experimental spectrum at energies exceeding 5.5 eV; a shoulder feature at approximately 5.2 eV is also seen in this spectrum. Excitation into the S₂ state, corresponding for the most part to a transition, only marginally accounts for absorption intensity for energies exceeding 4.5 eV when the (10,8) active space is used. This strongly suggests that excitation into states which are not included in the (6,4) and (10,8) active spaces, for example Rydberg states, must be responsible for the majority of the absorption intensity in this region.⁶³ Nonetheless, the excitation window used for the dynamics is well separated from this region, supporting the use of the smaller active space for the simulations.

3.3 Simulations

In total, 384 trajectories are propagated for both the triplet-inclusive and singlet-only trajectory ensembles, with all dynamics initiated on the S₁ state. For the triplet-inclusive ensemble, 371 trajectories (96.6%) were propagated successfully, while for the singlet-only ensemble, 346 trajectories (90.1%) were propagated successfully. A diagnostic analysis of the prematurely terminated trajectories is provided in the SI.

The classical populations, calculated within the molecular Coulomb Hamiltonian (MCH) representation of states,⁴² as well the changes in S–S bond distance, r_{S − S}, as a function of time are shown in Fig. 5. Both simulations exhibit a rapid decrease in the excited singlet-state population, S_n>0, dominated by loss of population from the S₁ state in the period 20 < t < 90 fs. This is mirrored by a commensurate increase in the S₀ population, in particular for the singlet-only simulation. In the simulation that also includes the triplet states, a rapid increase in the net triplet population dominates over the increase in S₀ population.

Fig. 5 Populations and S–S distances as a function of time (fs) for the (a) triplet-inclusive and (b) singlet-only simulations. (Top row) The classical molecular Coulomb Hamiltonian (MCH) state representation populations. The left panel shows the ground state population S₀, the net singlet excited-state population S_n (summing the populations for all singlet states with n > 0), and the net triplet population T (summing over all triplet states). The right panel is the same, but excluding the net triplet population. (Bottom row) The S–S bond length, r_S–S, for each trajectory with the mean bond length fitted to an exponentially damped sine function and shown as a thick black curve, analogously to the fitting in Rankine et al.¹⁴ For reference, the grey dot-dashed line indicates the S–S distance at the S₀ equilibrium geometry, 2.13 Å.

The mean times at which the first surface hopping events occur are 27.5 ± 6.4 fs and 32.9 ± 10.4 fs for the triplet-inclusive and singlet-only simulations, respectively. These values agree well with the timescales for depopulation of the initially excited S₁ state reported in previous studies.^13–15 The ultrafast nonadiabatic transfer out of S₁ is associated with stretching and subsequently cleaving the S–S bond, as seen in the plots in Fig. 5, which is also commensurate with the dissociative LIIC pathway shown in Fig. 3. The ground and excited states rapidly come together as the S–S bond stretches towards the MECI (S₀,S₁) facilitating effective population transfer. The S–S bond length is already greater than 3 Å at t = 25 fs, with most electronic states within a narrow 0.1 eV energy band at this point. The more rapid decay of the S₁ population in the triplet-inclusive case correlates well with the overall higher density of states in the triplet-inclusive case and thus the availability of more decay channels. Strikingly, the net triplet population, T_n, rises earlier than the S₀ population in the triplet-inclusive simulations, indicating that the proximity of singlet and triplet states already at short S–S distances leads to rapid and efficient ISC even before the onset of IC to the ground state. We also observe a moderate rise in the S₂ population—before any significant rise in the populations of any other electronic state—at approximately t = 20 fs in both trajectory ensembles; see SI for the variation in the populations of individual electronic states across the total trajectory time. This suggests a dominant S₂ ← S₁ pathway before population redistribution to other electronic states—especially triplet states—occurring immediately after. Reassuringly, Cao and Chen also make this observation in their investigation into the early-time population dynamics of 1,2-dithiane.¹³

Both simulations show a similar trend with respect to the S–S distance, r_S–S. The initial concerted ring-opening progresses from t = 0 to form linear dithiyl structures, reaching maximum values of r_S–S at approximately t = 150 fs before the thiyl termini return to close proximity at t ≈ 320 fs. The ring-closing dynamics during 150 < t < 320 fs is associated with a steady S₀ population, owing to the large energy gap between S₀ and all excited states as r_S–S approaches the equilibrium geometry S₀ bond length. The trajectories on S₀ during the ring-closing phase thus tend to remain on S₀ at least until the ring begins to open again. From approximately 320 fs, nuclear motion starts to significantly decohere as the vast majority of trajectories in both ensembles undergo a second ring-opening. Subsequently, most trajectories in both simulations remain in ring-open geometries for the remainder of the simulations. Fitting the time dependence of r_S–S to damped sine functions yields oscillatory time periods of 350.2 ± 1.1 fs and 351.1 ± 1.2 fs for the triplet-inclusive and singlet-only simulations, respectively. Remarkably, both time periods are within statistical error of each other and agree within error with the value of 349.6 fs reported by Rankine et al.¹⁴ Such quantitative agreement suggests that the nuclear dynamics is quite robust with respect to the number of electronic states considered and the electronic structure method used, with the main point for the latter being to reproduce the dissociative nature of the excited states.

The loss of concerted nuclear motions from approximately 400 fs correlates well with the equilibration of state populations seen in both simulations. Despite the fact that individual trajectories may still undergo oscillatory motions after the initial ring-opening and closing, the sulfur termini generally remain several Å ngströms apart for the majority of trajectories at these later times in both simulations. Thus, trajectories are almost always found in a region where the electronic state density is high—and thus also in a region where IC and ISC is efficient—resulting in a high probability of surface hopping. Indeed, the ubiquity of surface hops in this time range facilitates population equilibration of singlet and triplet states towards an almost statistical limit. At t = 1 ps, the total triplet population is 70.9%. Given the 12 triplet states and 4 singlet states included in the simulations, a statistical 75% triplet yield could be anticipated at longer times beyond 1 ps.

In addition, population transfer between electronic states is not localised to specific geometries, i.e. IC and ISC are not dominated by any specific conical intersections or minimum energy crossing points in this system. This is illustrated in Fig. 6, which shows the distribution of surface-hopping events in the triplet-inclusive ensemble for all IC and ISC transitions across the total trajectory time. The distribution of surface hops across a large region of molecular geometries in both simulations, for both IC and ISC, is a direct consequence of the comparatively flat potential energy landscape outside the ground state equilibrium well.

Fig. 6 Density heatmaps of crossing geometries—with respect to the disulfide distance, r_S–S, and the carbon backbone dihedral angle, ∠_CCCC—for all surface hopping events occurring in the triplet-inclusive trajectory ensemble. Surface hops are partitioned into subsets involving (a) singlet–singlet IC, (b) triplet–triplet IC, and (c) singlet–triplet ISC. The S₀ (min), S₁ (min), and MECI (S₀,S₁) geometries are indicated on each heatmap for reference. Details of additional electronic state minima and MECIs which were located—not shown here for clarity—can be found in the SI.

Returning to Fig. 5, both ensembles show very similar trends with respect to the variation of r_S–S across the trajectory time. However, several trajectories exhibit oscillations in r_S–S about the S₀ equilibrium disulfide bond length from approximately 300 fs in the singlet-only ensemble; these features are not observed in the triplet-inclusive ensemble. At 300 fs, 51/346 (14.7%) successful trajectories in the singlet-only ensemble propagate in the ground state potential well, whereas only 15/371 (4.0%) in the triplet-inclusive ensemble do so. Those in the former set are wholly responsible for the trajectories which show r_S–S becoming temporarily compressed below the S₀ equilibrium value (2.13 Å) at approximately 300 fs in the lower subplot in Fig. 5(b).

The presence of the triplet states effectively impedes trajectories from entering the S₀ potential well. All trajectories that enter the S₀ potential well at 300 fs in the triplet-inclusive ensemble eventually exit and return to ring-opened conformations. This reflection out of the ground state potential energy well results from the repulsive force between sulfur atoms as r_S–S is compressed beyond the S₀ equilibrium value. For the singlet-only simulations, out of the 51 trajectories in the S₀ potential well at t = 300 fs, 8 remain on the S₀ potential energy well for the rest of the simulation, giving rise to the oscillations in r_S–S around the equilibrium S₀ distance in Fig. 5(b); the other 43 trajectories return to the ring-open coupling region. It is possible that permanent S–S recombination may have been observed in the triplet-inclusive ensemble had an even greater number of trajectories been propagated. Nonetheless, it is clear that the recombination pathway is significantly disfavoured when ISC is accounted for when the total number of trajectories is fixed between both ensembles.

Trajectories that are not on S₀ are reflected by the dissociative nature of the excited states. This repulsive character of the excited state potentials with respect to r_S–S is the main factor behind the oscillatory Newton's cradle-type nature of most trajectories. Oscillations about the S₀ equilibrium bond length in Fig. 5(b)—which have a smaller amplitude and shorter period than oscillations resulting from the Newton's cradle-type behaviour—correspond to the eight trajectories which remain indefinitely on the ground state in this ensemble. The absence of such oscillations in the triplet-inclusive ensemble—with the exception of one trajectory which remains in the potential well from 870 fs onwards—is a consequence of the reduced probability of being in the S₀ state at any given time. This can be seen by comparing the trends in S₀ populations in Fig. 5. Notably, the proportion of trajectories remaining in the S₀ potential well after the first ring-closing event is very low (2.3%), even for the singlet-only ensemble.

Considering the singlet-only trajectories which enter the ground state potential well, the variation in r_S–S over the first 400 fs of simulations—separated into a subset which reflects out of the well and a subset which remains trapped for the duration of the simulation—is shown in Fig. 7. Whilst both trajectory subsets display coherent nuclear dynamics with respect to r_S–S over the first 200 fs, they deviate in how individual trajectories navigate through the potential well conformational landscape. Trajectories which ultimately leave the well generally exhibit turning points in r_S–S approximately 50 fs earlier than those which remain on S₀, thereby indicating a steeper descent into the well. This causes more rapid r_S–S bond compression below the ground-state equilibrium value, resulting in a more forceful collision of sulfur atoms. Conversely, trajectories which remain trapped on S₀ display a more gentle descent into the potential valley. This gives a greater period of time during which energy can be redistributed to other internal degrees of freedom, resulting in less available kinetic energy in the appropriate internal degrees of freedom for escaping the ground state.

Fig. 7 Variation in r_S–S over the first 400 fs of dynamics simulations for all trajectories which descend into the ground state potential well upon the first ring-closing event for the singlet-only ensemble. Trajectories which reflect out of the potential well are shown in grey, whilst those which remain trapped on S₀ for the rest of the simulation time are shown in purple.

Differential rates of vibrational relaxation and internal energy redistribution between individual trajectories also rationalises both the damped nature of subsequent oscillations in r_S–S for both individual trajectories and the fitted mean, and also the nuclear decoherence across both ensembles following the first ring-opening and ring-closing events. The latter observation results in ergodic ensemble dynamics over increasingly longer time periods as the excited state conformational landscape is explored more thoroughly. Finally we note that although the r_S–S coordinate has been found more than sufficient for the current analysis, it would nevertheless be interesting in future work to examine the trajectories using newly developed clustering tools.^64,65

3.4 X-ray scattering signals

Whilst slight differences in nuclear dynamics—namely the minor S–S recombination pathway in the singlet-only ensemble—are evident in the simulations, it is interesting to ascertain whether this difference may give rise to distinct and detectable signals in an experiment. Provided that there is sufficient spatial and temporal resolution, where the latter is ideally on the order of tens of femtoseconds or less, UXS experiments can be particularly sensitive to changes in molecular structure on ultrafast timescales.⁶⁶ Using the IAM to calculate UXS patterns is also straightforward and computationally inexpensive in that it only requires a set of molecular geometries and tabulated atomic form factors;⁵⁸ no prior knowledge of the electronic structure of the system is required. Despite the approximations inherent in the IAM, it is nonetheless a useful first-order approximation for calculating UXS patterns, particularly in light of the moderate molecular size of 1,2-dithiane and the relatively large number of intermediate timesteps involved in a single trajectory which would render more complex UXS simulations using ab initio electronic structure theory computationally expensive.³⁰

The mean percent difference UXS signals for both simulations, the difference in the UXS signals between both simulations, and the signal arising from the subset of trajectories which display permanent S–S recombination in the singlet-only simulations, are all shown in Fig. 8. Strikingly, the triplet-inclusive and singlet-only trajectory ensembles yield very similar UXS signals owing to almost identical mean nuclear dynamics. However, some differences persist, in particular for times t > 350 fs, as can be seen in Fig. 8(d). Whilst these differences are roughly an order of magnitude smaller than the overall scattering signal, they could in principle be detected, particularly considering the excellent signal-to-noise anticipated for new high-repetition rate experiments at X-ray free-electron laser (XFEL) facilities such as the upgraded linac coherent light source (LCLS-II).⁶⁷

Fig. 8 Rotationally averaged total X-ray scattering patterns for the simulations as a function of time. The scattering signal is expressed as a percent difference signal according to eqn (4), and is shown for the (a) triplet-inclusive and (b) singlet-only simulations, as well as (c) the subset of trajectories from the singlet-only simulation where S–S radical recombination occurs following the first ring-closing event at approximately t = 350 fs. The final panel (d) shows the difference between the signal shown in panels (a) and (b).

Further to this, the UXS signal arising from trajectories involving S–S recombination, shown in Fig. 8(c), is qualitatively similar to the total singlet-only signal in panel (b) for times t < 350 fs. This is due to the concerted Newton's cradle ring-opening/closing nuclear dynamics which is dominant during these initial stages of the reaction. However, Fig. 8(c) does deviate from the signals in panels (a) and (b) for times t > 350 fs. This is most clearly seen by the depletion in Fig. 8(c) of the strongly negative and strongly positive features centred at q = 0.8 Å⁻¹ and q = 1.7 Å⁻¹, respectively, which can be observed in panels (a) and (b) from 350 fs onwards. The recombination trajectories display weakly positive UXS bands for 2.5 Å⁻¹ < q < 4.0 Å⁻¹, whereas the ensemble averaged UXS signals in Fig. 8(a) and (b) are weakly negative in this region. It remains to be seen whether the scattering patterns in Fig. 8(a), which derive from the more physically realistic simulation that includes the triplet states, is in better agreement with UXS signals measured in future experiments, and if those experiments will require that the scattering is calculated using ab initio methods.⁵⁶

4 Conclusions

Simulations of the photoexcited dynamics in 1,2-dithiane on ultrafast timescales have demonstrated the characteristic oscillatory “Newton's cradle” molecular behaviour, which is shown to be present irrespective if triplet states are included or not. In both simulations, there is rapid population transfer and redistribution which commences approximately 30 fs after excitation. The rapid population transfer arises from the concerted ring-opening, which brings ground and excited electronic states into a narrow energy band as the S–S bond is elongated beyond dissociation. While both singlet and singlet–triplet simulations exhibit almost identical mean variation of the r_S–S distance, the inclusion of the triplet states noticeably deactivates the dithiyl radical recombination pathway observed from ∼300 fs onwards in the singlet-only simulations. This is a consequence of fewer trajectories descending into the ground-state S₀ potential well in the triplet-inclusive simulation, which can be seen as a statistical effect, with population distributed to the inactive triplet states. Nonetheless, we note that this study may somewhat underestimate the role of recombination in the singlet-only simulation as compared to the triplet-inclusive. The reason for this is that the electronic structure stability issues which prematurely terminate some trajectories occur when r_S–S is compressed and the molecule propagates on S₀ state. This issue disproportionally affects the singlet-only simulations, as explained in further detail in the SI. However, owing to the scarcity of S–S recombination events even in the singlet-only simulation, the effect on the predicted UXS signal is comparatively small, at least at the level of the IAM approximation. The differences in UXS signal between the two simulations nonetheless may still be experimentally resolvable. It is also worth noting that other experimental techniques, such as time-resolved photoelectron imaging,⁶⁸ may be better equipped to differentiate the populations on the singlet and triplet states. If anything, this strengthens the broad conclusion of this paper that the triplet states should be accounted for when considering the dynamics of this system and related molecular disulfides.

In terms of structural dynamics, the nuclear decoherence and equilibration of excited states after 400 fs are shared by both simulations and reflect the tendency for redistribution of kinetic energy and linear momentum along the S–S bond axis to other internal degrees of freedom at longer times. This behaviour is clearly quite robust and predicted irrespective of what electronic states are included in the simulations. Future ultrafast imaging experiments on 1,2-dithiane proposed at the European XFEL facility will hopefully eventually verify—or disprove—our predictions for the nuclear dynamics.

The hypothesis that the triplet states suppress the dithiyl recombination pathway in 1,2-dithiane may at first glance seem at odds with the experimental observation that disulfide linkages in proteins are particularly resilient to UV photolysis. However, there are differences between the currently examined small molecule gas-phase reaction and the S–S moiety in proteins. In proteins, polar or sterically bulky amino acid side chains in the vicinity of disulfide bridges, including those which may participate in hydrogen bonding, limit the conformational space available to the disulfide moiety upon S–S bond homolysis. This is likely to promote a more efficient and accessible radical recombination so as to retain the three-dimensional structure of the protein, and consequently its biological function, even on sub-picosecond timescales. Similar tendencies towards geminate recombination have been observed in solution.⁶⁹ By taking such effects into consideration in future studies, and by investigating how variations in disulfide molecular geometry and the carbon linkage affect the underlying electronic structure and dynamics,⁷⁰ a more holistic understanding of the impressive stability of disulfide bonds in proteins can hopefully be achieved.

Author contributions

James Merrick: conceptualisation (supporting); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); software (lead); validation (lead); visualisation (lead); writing – original draft (lead); writing – review and editing (lead). Lewis Hutton: software (supporting); writing – review and editing (supporting). Joseph C. Cooper: conceptualisation (supporting); writing – review and editing (supporting). Claire Vallance: conceptualisation (supporting); funding acquisition (equal); project administration (equal); supervision (equal); writing – review and editing (supporting). Adam Kirrander: conceptualisation (lead); funding acquisition (equal); project administration (equal); resources (lead); supervision (equal); writing – review and editing (supporting).

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting this work, including geometry optimisation and electronic structure benchmarking calculations in addition to auxiliary dynamics analyses, can be found in the SI. Supplementary information is available. See DOI:https://doi.org/10.1039/d5cp01776a

Acknowledgements

JM would like to thank Mats Simmermacher for helpful discussions and support in running the X-ray scattering calculations, and acknowledges a doctoral studentship from the Engineering and Physical Sciences Research Council (EPSRC UK) and the Carolyn and Franco Gianturco Theoretical Chemistry Scholarship from Linacre College at the University of Oxford. LH and AK acknowledge the support of the Leverhulme Trust grant RPG-2020-208, and JCC and AK acknowledge EPSRC grant EP/X026698/1. CV acknowledges funding from EPSRC Programme Grants EP/V026690/1 and EP/T021675/1. Finally, AK acknowledges funding from EPSRC grants EP/V006819/2, EP/V049240/2, and EP/X026973/1, as well as the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award number DE-SC0020276.

Notes and references

1. C. S. Sevier and C. A. Kaiser, Nat. Rev. Mol. Cell Biol., 2002, 3, 836–847.
2. R. W. Cowgill, Biochim. Biophys. Acta, Proteins Proteomics, 1967, 140, 37–44.
3. C. S. Sevier and C. A. Kaiser, Nat. Rev. Mol. Cell Biol., 2002, 3, 836–847.
4. W. Qiu, T. Li, L. Zhang, Y. Yang, Y.-T. Kao, L. Wang and D. Zhong, Chem. Phys., 2008, 350, 154–164.
5. W. E. Lyons, Nature, 1948, 162, 1004.
6. A. B. Callear and D. R. Dickson, Trans. Faraday Soc., 1970, 66, 1987–1995.
7. A. Rinker, C. D. Halleman and M. R. Wedlock, Chem. Phys. Lett., 2005, 414, 505–508.
8. C. Luo, W.-N. Du, X.-M. Duan, J.-Y. Liu and Z.-S. Li, Chem. Phys. Lett., 2009, 469, 242–246.
9. M. Ochmann, A. Hussain, I. von Ahnen, A. A. Cordones, K. Hong, J. H. Lee, R. Ma, K. Adamczyk, T. K. Kim, R. W. Schoenlein, O. Vendrell and N. Huse, J. Am. Chem. Soc., 2018, 140, 6554–6561.
10. A. B. Stephansen, M. A. B. Larsen, L. B. Klein and T. I. Sølling, Chem. Phys., 2014, 442, 77–80.
11. M. A. B. Larsen, A. B. Skov, C. M. Clausen, J. Ruddock, B. Stankus, P. M. Weber and T. I. Sølling, ChemPhysChem, 2018, 19, 2829–2834.
12. A. B. Stephansen, R. Y. Brogaard, T. S. Kuhlman, L. B. Klein, J. B. Christensen and T. I. Sølling, J. Am. Chem. Soc., 2012, 134, 20279–20281.
13. J. Cao and D.-C. Chen, Phys. Chem. Chem. Phys., 2019, 21, 4176–4183.
14. C. D. Rankine, J. P. F. Nunes, M. S. Robinson, P. D. Lane and D. A. Wann, Phys. Chem. Chem. Phys., 2016, 18, 27170–27174.
15. L. M. Ibele, Y. Lassmann, T. J. Martínez and B. F. E. Curchod, J. Chem. Phys., 2021, 154, 104110.
16. Y. Lassmann, D. Hollas and B. Curchod, J. Phys. Chem. Lett., 2022, 13, 12011–12018.
17. C. Middleton, C. Rankine and T. Penfold, Phys. Chem. Chem. Phys., 2023, 25, 13325–13334.
18. Z. Wang, R. Jing, Y. Li, D. Song, Y. Wan, N. Fukui, H. Shinokubo, Z. Kuang and A. Xia, J. Phys. Chem. Lett., 2023, 14, 8485–8492.
19. L. Ma, W. Du, H. Yong, B. Stankus, J. M. Ruddock, A. M. Carrascosa, N. Goff, Y. Chang, N. Zotev, D. Bellshaw, T. J. Lane, M. Liang, S. Boutet, S. Carbajo, J. S. Robinson, J. E. Koglin, M. P. Minitti, A. Kirrander, T. I. Sølling and P. M. Weber, Sci. Adv., 2025, 11, eadp9175.
20. P. A. Robertson, J. Merrick, D. Heathcote, M. S. Robinson, A. Butler, Y. Biddick, J. F. P. Nunes, C. Rankine, Z. Liu, S. F. Arrowsmith, J. O. F. Thompson, N. Murty Madugula, R. Chapman, E. Springate, E. A. Anderson, A. Kirrander and C. Vallance, Chem. Phys. Lett., 2025, 871, 142095.
21. W. Pauli, Z. Phys., 1927, 43, 601–623.
22. C. M. Marian, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 187–203.
23. A. D. Smith, E. M. Warne, D. Bellshaw, D. A. Horke, M. Tudorovskya, E. Springate, A. J. H. Jones, C. Cacho, R. T. Chapman, A. Kirrander and R. S. Minns, Phys. Rev. Lett., 2018, 120, 183003.
24. D. Bellshaw, R. S. Minns and A. Kirrander, Phys. Chem. Chem. Phys., 2019, 21, 14226–14237.
25. W. O. Razmus, K. Acheson, P. Bucksbaum, M. Centurion, E. Champenois, I. Gabalski, M. C. Hoffman, A. Howard, M.-F. Lin, Y. Liu, P. Nunes, S. Saha, X. Shen, M. Ware, E. M. Warne, T. Weinacht, K. Wilkin, J. Yang, T. J. A. Wolf, A. Kirrander, R. S. Minns and R. Forbes, Phys. Chem. Chem. Phys., 2022, 24, 15416–15427.
26. I. Gabalski, M. Sere, K. Acheson, F. Allum, S. Boutet, G. Dixit, R. Forbes, J. M. Glownia, N. Goff, K. Hegazy, A. J. Howard, M. Liang, M. P. Minitti, R. S. Minns, A. Natan, N. Peard, W. O. Rasmus, R. J. Sension, M. R. Ware, P. M. Weber, N. Werby, T. J. A. Wolf, A. Kirrander and P. H. Bucksbaum, J. Chem. Phys., 2022, 157, 164305.
27. F. Coppola, M. Hussain, J. Zhao, A. M. El-Zohry and M. Pastore, J. Phys. Chem. C, 2024, 128, 11998–12009.
28. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
29. A. Odate, A. Kirrander, P. M. Weber and M. P. Minitti, Adv. Phys.: X, 2023, 8, 2126796.
30. M. Simmermacher, P. M. Weber and A. Kirrander, Structural Dynamics with X-ray and Electron Scattering, Royal Society of Chemistry, United Kingdom, 1st edn, 2023, ch. 3, vol. 25, p. 85.
31. H. Yong, A. Kirrander and P. M. Weber, Structural Dynamics with X-ray and Electron Scattering, Royal Society of Chemistry, United Kingdom, 1st edn, 2024, ch. 3, vol. 25, p. 344.
32. P. M. Weber, B. Stankus and A. Kirrander, in Ultrafast X-Ray Scattering: New Views of Chemical Reaction Dynamics, ed. K. Ueda, Springer Nature, Singapore, 2024, pp. 195–227.
33. B. O. Roos, P. R. Taylor and P. E. Sigbahn, Chem. Phys., 1980, 48, 157–173.
34. G. Li Manni, I. Fdez Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. de Graaf, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov(Gulak), L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krosnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Mörchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. Olivucci, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023, 19, 6933–6991.
35. B. O. Roos and P.-Å. Malmqvist, Phys. Chem. Chem. Phys., 2004, 6, 2919–2927.
36. B. A. Heß, C. M. Marian, U. Wahlgren and O. Gropen, Chem. Phys. Lett., 1996, 251, 365–371.
37. D. Peng and M. Reiher, Theor. Chem. Acc., 2012, 131, 1081.
38. J. P. Zobel, P.-O. Widmark and V. Veryazov, J. Chem. Theory Comput., 2020, 16, 278–294.
39. N. Forsberg and P.-Å. Malmqvist, Chem. Phys. Lett., 1997, 274, 196–204.
40. T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-Structure Theory, 2013, pp. 183–186.
41. M. Oppel, sharc-md/sharc: SHARC Release 2.1.2, 2022 DOI:10.5281/zenodo.6077788.
42. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
43. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258.
44. E. Wigner, Phys. Rev., 1932, 40, 749–759.
45. M. Barbatti and K. Sen, Int. J. Quantum Chem., 2016, 116, 762–771.
46. J. P. F. Nunes, PhD thesis, Chemistry, University of York, 2017.
47. G. Bergson, G. Claeson, L. Schotte, A. Block-bolten, J. M. Toguri and H. Flood, Acta Chem. Scand., 1962, 16, 1159–1174.
48. R. Crespo-Otero and M. Barbatti, Theor. Chem. Acc., 2012, 131, 1237.
49. G. Granucci, M. Persico and G. Spighi, J. Chem. Phys., 2012, 137, 22A501.
50. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
51. G. Granucci, M. Persico and A. Zoccante, J. Chem. Phys., 2010, 133, 134111.
52. A. Kirrander, K. Saita and D. V. Shalashilin, J. Chem. Theory Comput., 2016, 12, 957–967.
53. M. Stefanou, K. Saita, D. V. Shalashilin and A. Kirrander, Chem. Phys. Lett., 2017, 683, 300–305.
54. T. Northey, N. Zotev and A. Kirrander, J. Chem. Theory Comput., 2014, 10, 4911.
55. T. Northey, A. M. Carrascosa, S. Schäfer and A. Kirrander, J. Chem. Phys., 2016, 145, 154304.
56. A. M. Carrascosa, H. Yong, D. L. Crittenden, P. M. Weber and A. Kirrander, J. Chem. Theory Comput., 2019, 15, 2836–2846.
57. I. Waller, D. R. Hartree and R. H. Fowler, Proc. R. Soc. London, Ser. A, 1929, 124, 119–142.
58. E. Prince and I. U. of Crystallography, International tables for crystallography: Volume C: Mathematical, Physical and Chemical Tables, Wiley, Dordrecht, 3rd edn, 2004.
59. N. Zotev, A. M. Carrascosa, M. Simmermacher and A. Kirrander, J. Chem. Theory Comput., 2020, 16, 2594–2605.
60. M. H. Issa Yavari and R. Amiri, Phosphorus, Sulfur Silicon Relat. Elem., 2004, 179, 2015–2023.
61. C. Rankine, PhD thesis, Chemistry, University of York, 2020.
62. V. Rishi, N. C. Cole-Filipiak, K. Ramasesha and L. M. McCaslin, Phys. Chem. Chem. Phys., 2024, 26, 23986–23997.
63. P. Kimber and F. Plasser, Comprehensive Computational Chemistry, Elsevier, Oxford, 1st edn, 2024, pp. 55–83.
64. K. Acheson and A. Kirrander, J. Chem. Theory Comput., 2023, 19, 6126–6138.
65. J. Kára, K. Acheson and A. Kirrander, J. Chem. Theory Comput., 2025, 5789–5802.
66. J. M. Budarz, M. P. Minitti, D. V. Cofer-Shabica, B. Stankus, A. Kirrander, J. B. Hastings and P. M. Weber, J. Phys. B, 2016, 49, 034001.
67. A. Moreno Carrascosa, J. P. Coe, M. Simmermacher, M. J. Paterson and A. Kirrander, Phys. Chem. Chem. Phys., 2022, 24, 24542–24552.
68. N. Kotsina and D. Townsend, Phys. Chem. Chem. Phys., 2021, 23, 10736–10755.
69. M. Ochmann, J. Harich, R. Ma, A. Freibert, Y. Kim, M. Gopannagari, D. H. Hong, D. Nam, S. Kim, M. Kim, I. Eom, J. H. Lee, B. A. Yorke, T. K. Kim and N. Huse, Nat. Commun., 2024, 15, 8838.
70. H. G. McGhee, R. Totani, O. Plekan, M. Coreno, M. de Simone, A. Mumtaz, S. Singh, B. C. Schroeder, B. F. E. Curchod and R. A. Ingle, J. Chem. Phys., 2024, 161, 134303.

---

Footnote

† The double-excitation character of S₃ and T₄ comprises two separate double excitations from the two sulfur lone-pair orbitals – n^u_s and n^g_s – displayed in Fig. 2.

---