Radicals in aqueous solution: assessment of density-corrected SCAN functional †

We study self-interaction eﬀects in solvated and strongly-correlated cationic molecular clusters, with a focus on the solvated hydroxyl radical. To address the self-interaction issue, we apply the DC-r 2 SCAN method, with the auxiliary density matrix approach. Validating our method through simulations of bulk liquid water, we demonstrate that DC-r 2 SCAN maintains the structural accuracy of r 2 SCAN while eﬀectively addressing spin density localization issues. Extending our analysis to solvated cationic molecular clusters, we find that the hemibonded motif in the [CH 3 S ‘ CH 3 SH] + cluster is disrupted in the DC-r 2 SCAN simulation, in contrast to r 2 SCAN that preserves the (three-electron-two-center)-bonded motif. Similarly, for the [SH ‘ SH 2 ] + cluster, r 2 SCAN restores the hemibonded motif through spin leakage, while DC-r 2 SCAN predicts a weaker hemibond formation influenced by solvent–solute interactions. Our findings demonstrate the potential of DC-r 2 SCAN combined with the auxiliary density matrix method to improve electronic structure calculations, providing insights into the properties of solvated cationic molecular clusters. This work contributes to the advancement of self-interaction corrected electronic structure theory and offers a computational framework for modeling condensed phase systems with intricate correlation effects.


Introduction
Density-corrected density functional theory (DC-DFT) has a rich history dating back to the early days of practical DFT simulations. Initially, researchers sought to enhance the Hartree-Fock (HF) density with density functionals containing correlation effects in a non-self-consistent manner. 1,2 However, the rapid advancements in Kohn-Sham DFT (KS-DFT) and advanced exchange-correlation functionals have made these approaches less relevant, as Kohn-Sham densities often provide excellent approximations to the exact density. 3 Moreover, self-consistent densities enable straightforward access to various molecular properties of interest, such as forces and stress tensors, crucial for molecular dynamics simulations.
The core concept behind DC-DFT 4,5 is the separation of the total energetic error in Kohn-Sham DFT calculations into functional-driven and density-driven components. Based on this rigorously defined framework, the field of DC-DFT has undergone a renaissance in recent years. [6][7][8][9][10][11][12] In contrast to the initial interest in density-corrected methods, current attention has shifted towards the limitations of exchange-correlation functionals based on the generalized gradient approximation (GGA), which can lead to density-driven errors and enhanced electron density delocalization. 13,14 In such cases, the selfinteraction-free Hartree-Fock density has emerged as a viable surrogate for the exact density. When combined with a Kohn-Sham energy functional employing the strongly-constrained and appropriately normed SCAN functional, 15 this approach, referred to as DC-SCAN, 8,9 yields exceptional results, demonstrating chemical accuracy in the study of water clusters.
The SCAN functional, 15 known for satisfying all 17 exact constraints for meta-GGA functionals, has garnered significant interest in the field. It has demonstrated superior performance compared to other functionals in aqueous systems. 16 However, like many other functionals, SCAN is susceptible to densitydriven errors, which can result in the artificial over-stabilization of hydrogen-bonded systems. 3,9 To address numerical stability and grid convergence issues of the original functional, we used the regularized-restored r 2 SCAN 17 in this study.
In certain cases, the utilization of the Hartree-Fock density instead of the self-consistently obtained Kohn-Sham density has shown promise in mitigating density-driven errors. 6,7,11 The application of DC-DFT, 5 when appropriately employed, has demonstrated improved simulation accuracy, particularly when density-driven errors overshadow functional-driven errors. Notably, DC-SCAN has shown favorable performance, even surpassing the sophisticated deep-learned local-hybrid functional DM21 developed by DeepMind. 10 The research by Dasgupta et al. 8,9 has repeatedly demonstrated that the addition of dispersion corrections, such as D3, 18 negatively impacts the energetics of SCAN and DC-SCAN in both charged and neutral water clusters. Based on these finding, we chose not to incorporate dispersion corrections in our study. Although, Song et al. 12 recently emphasized the importance of adding a reparametrized D4 dispersion correction 19 to HF-r 2 SCAN, for accurately describing long-range dispersion interactions.
The usability of density-corrected DFT methods is constrained by two primary issues. Firstly, the evaluation of the DC-DFT energy functional relies on a density, typically the Hartree-Fock density, which was variationally optimized for a different energy functional. This non-variational nature of the DC-DFT energy functional with respect to the reference density poses limits on its applicability. This becomes evident in the challenges associated with accessing forces and performing geometry optimizations, which are not trivial tasks. 20 Secondly, the computation of the reference density through Hartree-Fock calculations poses a significant bottleneck. Unlike KS-DFT, which exhibits cubic scaling O(N 3 ) with system size N, obtaining the Hartree-Fock exchange involves solving the computationally demanding, formally quartic scaling O(N 4 ) four-center twoelectron repulsion integrals.
Overcoming the challenge of accessing self-consistent analytical nuclear gradients in DC-SCAN has in the past been addressed through the development of data-driven potential energy functions within a many-body formalism. 8 While the theoretical framework for a variational formulation of DC-DFT has been proposed for some time, 21 its widespread implementation in quantum chemistry codes is still limited. 22 However, the availability of analytical forces in DC-DFT would offer significant advantages, as it not only addresses the challenge of self-consistency but also improves standard DFT forces and geometries in density-sensitive calculations. The wider adoption of such implementations would enable comprehensive investigations into the performance of DC-DFT in terms of both geometries and energies. 20 Here, we present a variational formulation of densitycorrected DFT, implemented in the CP2K software package. 23,24 The implementation follows the framework proposed by Verma et al., 21 and which extends the recently developed variational Harris functional energy correction. 25 Our approach utilizes the auxiliary density matrix method (ADMM) 26 and truncated interaction potentials, 27 available in CP2K, for efficient computation of the Hartree-Fock density. Effectively, the DC-DFT method employed here should be named ADMM-DC-r 2 SCAN, but for sake of readability the ADMM label will be omitted henceforth.
We aim to explore the potential of this advanced DFT method in addressing self-interaction errors that arise in stronglycorrelated materials. Specifically, we focus on the aqueous hydroxyl radical as a challenging system to simulate accurately using conventional exchange-correlation functionals. In DFT simulations, a two-center, three-electron (2c-3e) hemibond is observed in the first solvation shell of the hydroxyl radical due to favorable overlap between the singly-occupied molecular orbital of the radical and the Frontier lone pair of a neighboring water molecule. However, this hemibond disappears when the same structures are relaxed at the MP2 level, leaving behind a ''pseudo hemibond'' with a close contact between oxygen atoms. 28 Previous investigations into the hemibonded motif in theoretical studies have attributed its origin to factors such as selfinteraction errors in GGA-based DFT, 29,30 size effects, 31 or even an intricate interplay between these factors. 28,32 In line with these studies, our investigation aims to contribute to the understanding of the hemibonded motif by offering a new perspective to this extensively researched and complex puzzle.
In addition to revisiting lingering questions regarding the prevalence of the hydroxyl-water hemibond in theoretical studies, we also investigate the radical cation clusters [CH 3 S'CH 3 SH] +33 and [HS'SH 2 ] + . [34][35][36] The sulfur-sulfur hemibonded configuration in these clusters has been predicted to be more stable than the hydrogen-bonded motif in the gas phase. It has been suggested that the hemibonded ion core in these clusters remains stable upon microhydration with a single water molecule, while the (2c-3e) bond breaks upon proton transfer in the presence of a single methanol or ethanol molecule. 37,38 Here, we expand on these findings by investigating the behavior of these systems in a fully solvated environment, aiming to determine how the stability of the hemibonded motif is influenced by the solvent.
The article is structured as follows: Section 2 covers the development and implementation of the DC-r 2 SCAN method. Section 3 provides details regarding the computational setup for the MD simulations. Section 4 validates the DC-r 2 SCAN implementation by comparing its performance in reproducing the structure of bulk liquid water. Section 5 investigates hemibonding in solvated radicals and radical cation clusters. Finally, Section 6 presents the conclusion and outlines future research possibilities.

The density-corrected energy functional
The density corrected DFT energy functional is the Fock/Kohn-Sham-DFT functional that is evaluated using the electronic density obtained from a converged reference calculation. It comprises the one-electron core Hamiltonian h, the Coulomb J operator, as well as exact exchange K and/or contributions of the exchange-correlation (XC) energy E xc [n]. The scaling factor a x assumes a value of 1 in the context of Hartree-Fock theory. In the case of Kohn-Sham theory, this factor is set to zero unless a hybrid functional is employed, requiring an additional contribution from the exchange-correlation potential F xc .
The density matrix P is defined as a function of the occupied molecular orbital (MO) coefficients C and is expanded in the atomic orbital (AO) basis. This expansion is expressed by the equation: where {m, n, . . .} correspond to the atomic orbitals (AOs), and {i, j, . . .} refer to the occupied MOs. The index s corresponds to the spin component.

The variational DC-DFT Lagrangian
Nuclear gradients and the stress tensor are indispensable to compute molecular properties. To avoid the computational expense associated with determining the first-order response of the MO coefficients, which takes into account the geometry dependence of the AOs, the DC-DFT energy functional can be extended by means of two additional constraints. This extension leads to the formulation of the DC-DFT Lagrangian, which can be expressed as: The first additional constraint, (eqn (3b)), ensures the stationarity of the ground-state equations, whether they are based on Kohn-Sham or Hartree-Fock theory. By rearranging the equations and utilizing the definition of the projection operator Q the equivalence with MO-based formulations becomes apparent: where e represents the orbital energies. The overlap matrix S is defined in terms of the AOs j m (r) as: Only the virtual-occupied blocks of the Z vector, denoted as Z, contribute, while the occupied-occupied block is implicitly assumed to be zero: The Fock/Kohn-Sham operator F[P] appearing in the Brillouin constraint is defined as: In the case of DC-r 2 SCAN, F simplifies to the Fock operator, with a x = 1 and without the exchange-correlation potential F xc . Efficient treatment of the exact Hartree-Fock exchange (HFX) in the reference and linear response equations is achieved through the ADMM approximation in CP2K.
The second additional constraint ensures the orthogonality of the occupied MOs. Analogous to Z, this constraint is assumed to be symmetric. The determination of the Lagrange multipliers W and Z involves taking the derivative of L DC-DFT with respect to the MO coefficients C and projecting onto either the virtual or occupied space: The resulting expressions for the W multiplier are as follows: The matrix W exclusively includes occupied-occupied contributions. The intermediate H is derived from the KS matrix contributions of the kernel matrix and varies depending on the chosen reference functional: À a x d ss 0 ðkisjljs 0 Þ þ ðkjsjlisÞ ð Þ þ 2f xc kls 0 ;ijs Þ; with the exchange-correlation kernel f xc mnkl , defined as the second-order functional derivative of E xc [n], and where the two-electron repulsion integrals are defined as ðmnjklÞ ¼ ðð j m ðrÞj n ðrÞ 1 jr À r 0 j j k ðr 0 Þj l ðr 0 Þdr dr 0 : The linear Z-vector equation, derived from eqn (9), can be expressed as: The terms on the left-hand side of eqn (15) originate from the chosen ground-state energy functional and are independent of the density correction on the right-hand side. The right-hand side B represents the difference between the DC-DFT and the reference energy functional. Specifically, in the case of DCr 2 SCAN, it reduces to the difference between the XC-potential of the r 2 SCAN functional and the exact exchange contribution to the Fock operator of the reference calculation: The response density matrix P z is defined as: where symmetrization is employed to ensure numerical stability by making the response density matrix symmetric and positive semi-definite.
The response equations (eqn (15) are solved iteratively using the conjugate gradient method. 39 A trial solution X is constructed and updated until the matrix-vector product P ais X ais A ais converges to the right-hand-side property gradient B. The orbital Hessian consists of a Coulomb and exchange part, both of which must be consistently updated with the calculation of Fock matrix contributions during the SCF procedure.

Auxiliary density matrix method
The addition of non-local Hartree-Fock exchange (HFX) to density-corrected DFT methods introduces a significant computational burden. The computational cost of HFX scales with the fourth power of the system size due to the evaluation of 4-electron 2-center electron repulsion integrals (ERIs) required to obtain the Hartree-Fock exchange energy. The HFX energy is given by: Although techniques such as integral screening 40 and the utilization of short-range exchange operators 41,42 have partially addressed this issue, the scaling with respect to basis set quality remains poor in periodic systems. To overcome this challenge, the auxiliary density matrix method (ADMM) 26,27 was developed and implemented in the CP2K software package.
The ADMM approach accelerates the evaluation of HFX by introducing an auxiliary density matrix P that is either smaller in size or decays more rapidly than the original density matrix. The HFX energy is then approximated using the following expression: EE HFX By introducing a smaller and/or less diffuse auxiliary basis, ADMM provides an efficient approximation for computing HFX integrals. The method assumes that the basis set incompleteness error in the exchange energy between the primary and auxiliary density matrices can be well captured by a GGA, even in the presence of qualitative differences between GGA exchange and HFX. By introducing a GGA functional in ADMM, the method is no longer self-interaction-free, unlike the exact exchange that is approximated. Nevertheless, ADMM has been shown to maintain intermediate degrees of localization, depending on the employed exchange fraction. 26 Among the many variants of the ADMM method, 26,43,44 two commonly studied approaches are purified wave function fitting (ADMM1) and non-purified wave function fitting (ADMM2). In ADMM2, the auxiliary MO coefficients minimize the square difference between occupied wave functions in the primary and auxiliary bases, while the auxiliary density matrix is obtained through a simple projection. Purified wave function fitting, in addition to this, enforces orthonormality, ensuring that the auxiliary density matrix satisfies the idempotency conditions. 45 Consequently, an eigenvalue correction step is necessary with purified wave function fitting before using the resulting Kohn-Sham matrix for DC-DFT.
In the following discussion, we focus on ADMM2, which does not require any post-SCF treatment. The ADMM2-approximated results will be referred to using the general acronym ADMM. The speedup achieved in the evaluation of the Hartree-Fock exchange with the ADMM method depends on the choice of primary and auxiliary basis sets and can often reach orders of magnitude. The ADMM method has been employed before in the context of linear response and time-dependent density functional theory (TDDFT), where the linear response equations have been extended to include the ADMM terms. 46,47

Nuclear gradients
Finally, let us consider the gradient of the DC-DFT Lagrangian with respect to the nuclear coordinate R. The expression for the nuclear gradient is given by: where the intermediate quantity L mns is defined as: These equations provide a means to compute the gradient of the energy with respect to the nuclear coordinates. Analogous equations can be derived to obtain the gradient with respect to the strain tensor, which allows for the calculation of the stress tensor.

Computational details
All simulations were performed using the QUICKSTEP module in the CP2K simulation software package, 23,24 in the Gaussian and plane waves (GPW) framework 48 with 3D periodic boundary conditions. The atomic core electrons were described using Goedecker-Teter-Hutter (GTH) pseudopotentials, 49,50 while the valence electron molecular orbitals were expanded in triple-z double-polarized basis sets (TZV2P), 51 which were optimized with the numerical ATOM code of CP2K for the SCAN functional. The kinetic energy cutoff for the plane-wave expansion of the density was set to 1200 Ry. 52 To mitigate the computational cost associated with the quartic scaling (O(N 4 )) of the four-center two-electron integrals, the exchange calculations of the DC-SCAN simulations employed ADMM. 26 A contracted and polarized auxiliary basis set consisting of three Gaussian exponents for each valence orbital (cpFIT3) was used to construct the smaller and rapidly decaying auxiliary density matrix. The PBE exchange functional was employed to correct the basis set incompleteness error introduced by the auxiliary basis. The truncated Coulomb operator 27 had a cutoff radius of 6.0 Å, approximately half the length of the smallest edge of the simulation cell for all systems. The Schwarz integral screening threshold was set to 10 À6 atomic units.
The Born-Oppenheimer ab initio molecular dynamics simulations were performed in the canonical (NVT) ensemble with the CSVR thermostat 53 at a temperature of 350 K and a time step of 0.5 fs. The solvated radical and cluster systems were generated using Packmol. 54 The systems were initially equilibrated through a 1 ps MD simulation with a thermostat attached to each degree of freedom and a time constant of 10 fs, followed by 2.5 ps of equilibration with a single thermostat controlling all degrees of freedom and a time constant of 100 fs. The production runs were conducted for a duration of 30 ps. Evolution of the conserved energy and temperature fluctuations are listed in the ESI. † The radial distribution functions were calculated using VMD 55 over the entire 30 ps trajectory of the production runs, with a bin width of 0.01 Å. Hydrogen bond analysis was performed using MDanalysis. 56

Validation
In order to validate our implementation, we first examined a system consisting of 64 water molecules. The objective was to compare the simulated structure with those of bulk liquid water under ambient conditions, specifically at a temperature of 300 K and volumetric mass density of 1 g cm À3 . To achieve the desired volumetric mass density, we arranged the molecules within a cubic simulation cell, with a side-length of 12.4138 Å.

Structural properties of bulk liquid water
The radial distribution functions of bulk liquid water (H 2 O) 64 are presented in Fig. 1. We found that r 2 SCAN and DC-r 2 SCAN accurately reproduce the structure of bulk liquid water at ambient temperatures, 57 specifically at 350 K. This temperature was carefully selected to avoid freezing and ensure a high level of agreement between the theoretical RDFs and experimental data. Previous simulations using the SCAN functional at 300 K exhibited excessive structuring compared to experimental observations, 52 while at 330 K, SCAN resulted in an augmented first solvation shell. 58,59 Notably, a recent study employing a machine-learned KS-SCAN potential predicted the freezing temperature of water to be 312.5 K. 60 In this work we do not consider the explicit treatment of nuclear quantum effects (NQEs). It has been debated in the past, that elevated temperatures may also account for NQEs, 59,61,62 which aligns with our cautious approach to prevent freezing. The necessity of incorporating NQEs into the analysis remains a topic of ongoing discussion. 63,64 Our results confirm the suitability of 350 K as an appropriate temperature for DFT simulations of bulk liquid water using the r 2 SCAN and DC-r 2 SCAN functionals, offering an accurate description of the water structure that aligns with experimental measurements.
In Table 1, we report the positions and intensities of the first maximum, first minimum and coordination number between the oxygen atoms, and compare them with the experimental values for pure liquid water at ambient conditions. The peak positions and intensities, obtained for both r 2 SCAN and DCr 2 SCAN, align well with the experimental values. The coordination number n OO in Table 1 is defined 57 as with a cutoff radius r cut of 3.3 Å and the volumetric mass density of the simulation cell r (1 g cm À3 ).

Solvated hydroxyl radical
Density functional theory commonly relies on GGA functionals, which strike a favorable balance between computational cost and accuracy, making them suitable for large-scale molecular dynamics simulations. 29 However, certain systems, such as the solvated hydroxyl radical, display exceptional behavior that challenges conventional GGA simulations. In particular, a distinctive bonding arrangement known as a two-center threeelectron hemibond between the hydroxyl radical and the solvent has been observed and extensively discussed over the past two decades. The origin of this hemibonded structure has been the subject of ongoing research efforts. Recently, Rana and Herbert 32 provided a comprehensive summary of these investigations.
In this study, we reexamine the solvated hydroxyl radical using the DC-r 2 SCAN functional and compare our results to those obtained from the r 2 SCAN functional, aiming to contribute new insights to this complex phenomenon. 5.1.1 Radial distribution function. In Fig. 2 we plot the radial distribution functions (RDFs) between the atoms of the hydroxyl radical and the solvent. The hemibonded configuration can be identified in r 2 SCAN calculations as well defined shoulder with onset at 2.1 Å. Despite being among the most sophisticated generalised gradient approximation based functionals, the r 2 SCAN exchange-correlation functional predicts the formation of a hemibond between the hydroxyl radical and the aqueous solvent. In contrast, the minimum approach distance with DC-r 2 SCAN is 2.5 Å and immediately resolves into the hydrogen bonding peak located at 2.815 Å. No hemibonded configuration can be identified from the O-O* radial distribution function. However, in recent PBE0 + D3 simulations of OH(H 2 O) 63 at 370 K, Rana and Herbert 28 showed that the hemibonded distance can be enlarged and thus be hidden in the contiguous hydrogen bonding peak. We observe that the intensity of the first peak is lower for the DC-r 2 SCAN method (1.92) than for r 2 SCAN (2.23). The first solvation shell of DC-r 2 SCAN is much broader than that of r 2 SCAN, resulting in a much further positioned first minimum from 4.07 Å to 4.98 Å. This suggests a lower coordination number for the first solvation shell in the simulation with DC-r 2 SCAN, indicating a weaker solvation effect and reduced water coordination around the radical compared to the simulation with the r 2 SCAN functional.
The various O-H pair correlation functions in Fig. 2 We find for both methods that the hydrogen bond is shorter when the hydroxyl acts as donor. Considering the overall Table 1 Positions and intensities of the 1st maximum (r 1 , g 1 ), 1st minimum (r 2 , g 2 ), and oxygen-oxygen coordination number (n O*-O* )  difference in peak height, there is a higher probability of finding a donor rather than an acceptor bond. In particular, the DC-r 2 SCAN donor peak is higher than that for r 2 SCAN, indicating a higher probability of occurrence and increased strength of the donor hydrogen bond in DC-r 2 SCAN. This difference in donor peak height then also predicts a less stable donor hydrogen bond with the r 2 SCAN functional. Integrating the first donor hydrogen bond peak yields 0.96 for r 2 SCAN and 1.00 for DC-r 2 SCAN. The radical can donate only one hydrogen bond, so the donor bond exists always in DC-r 2 SCAN. On the other hand, the r 2 SCAN acceptor (O*-H) peak is higher than the peak observed for DC-r 2 SCAN, suggesting a higher probability of acceptor bond occurrence. This might be due to the enhanced delocalization of the electron density predicted by the r 2 SCAN functional, leading to more acceptor hydrogen bonds. This can be seen by integrating the first solvation shell, which yields a value of 1.4 for r 2 SCAN and 0.99 for DC-r 2 SCAN.
In Table 2, we report the positions and intensities of the first maximum, first minimum and coordination number between the radical and solvent oxygen atoms, and compare them with the experimental values for pure liquid water at ambient conditions. The lower g 1 peak in DC-r 2 SCAN in comparison to r 2 SCAN is likely linked to the much more diffuse first solvation shell, as can be seen from the far location of the first minimum at 4.985 Å. In comparable hybrid functional studies of systems with 64 molecules the first minimum is located in the range of 3.2 Å to 3.7 Å, 28,32 and at 4.4 Å with the revAPBEK functional. 65 5.1.2 Hydrogen bond analysis. The information on the hydrogen bond network derived from the radial distribution functions in Fig. 2 may be extended by performing a hydrogen bond analysis. The hydrogen bonds were determined by using a geometric criterion. 66 Two molecules are considered hydrogen bonded if the donor-acceptor oxygen distance d O*-O is less than 3.5 Å, and the donor-hydrogen-acceptor angle y O*-H-O is between 1351 and 1801.
In Fig. 3, we show the distribution of the number of hydrogen bonds between the solvent and the hydroxyl radical. In r 2 SCAN calculations, the number of hydrogen bonds between the solvent and radical is on average three. The likelihood of three hydrogen bonds on the hydroxyl radical is reduced by 25.8% in DC-r 2 SCAN, in comparison to r 2 SCAN. Conversely, configurations with two hydrogen bonds are 25% more likely with DC-r 2 SCAN, compared to r 2 SCAN. DC-r 2 SCAN reproduces the trend of self-interaction corrected BLYP simulations. 29 Furthermore, in cases where three hydrogen bonds are formed between the radical and the solvent, the radical acts twice as acceptor and once as donor, confirming the results reported by VandeVondele et al. 29 The hydrogen bonding between the solvent molecules in the system was analyzed, excluding the hydroxyl radical. The average number of hydrogen bonds per water molecule is 3.59 AE 0.11 and 3.65 AE 0.11 for r 2 SCAN and DC-r 2 SCAN, respectively. These values are in agreement with the inferred experimental value of 3.58 67 and previous hybrid functional results. 68 Notably, the higher average number of hydrogen bonds in DC-r 2 SCAN indicates that the solvent is more strongly bound to itself compared to r 2 SCAN. This observation is consistent with the lower average number of hydrogen bonds between the radical and the solvent in DC-r 2 SCAN compared to r 2 SCAN, as shown in Fig. 3.
It needs to be stressed that the hydrogen bond analysis strongly depends on the criteria that are chosen to define a hydrogen bonding event. Since it is not an observable in liquid water experiments that can be measured or calculated, definitions of hydrogen bonding are manifold. 69,70 In particular, within the simple geometric criterion we apply here, adjusting the cutoff criteria has drastic effects on the hydrogen bond analysis. Tightening the criterion on the hydrogen bond angle from 1351 to 1501 reduces the average number of hydrogen bonds per water molecule to 2.83 AE 0.17 and 2.981 AE 0.16 for r 2 SCAN and DCr 2 SCAN, respectively.
The hydrogen bonding network in liquid water DFT simulations is strongly affected, not only by the extensive liberty of choosing the hydrogen bonding parameters, but also by dispersion correction methods such as the D3 correction by Grimme, 18 which can significantly enhance the hydrogen bonding network. However, the use of the D3 method in its original parametrization, has been shown to introduce significant overbinding, leading to an overall deterioration of the energetics for SCAN and DC-SCAN. 8,9 Therefore, these corrections were omitted in our study.
In a recent study by Song et al., 12 it has been argued that despite delivering high accuracy for pure water simulations, DC-r 2 SCAN without a dispersion correction cannot accurately describe longrange dispersion interactions in other non-covalent systems. Table 2 Positions and intensities of the first maximum (r 1 , g 1 ), first minimum (r 2 , g 2 ), and coordination number (n O*O ) in the radical oxygen-solvent oxygen RDF for r 2   Therefore, the use of dispersion corrections may still be necessary to accurately describe non-covalent interactions of different nature, and a reparametrized D4 dispersion correction 19 was required to achieve this accuracy in standard non-covalent datasets.
In Fig. 4 we plot the distribution of hydrogen bond angles between the radical and the solvent when the distance criterion is satisfied, considering the radical as either a donor (upper panels) or an acceptor (lower panels). For both r 2 SCAN and DC-r 2 SCAN, the distributions of donor hydrogen bonds exhibit a prominent peak centered around 1651. In particular with the latter method, a well-defined peak in the donor hydrogen bond distribution can be identified. Together with the previously mentioned increase in hydrogen bonds involving the radical as a donor (see Fig. 2), this signifies a preferential structuring of the solute-solvent system in DC-r 2 SCAN functional simulations.
Furthermore, the observed slower decrease in the radial distribution function (RDF) and the extended first solvation shell in DC-r 2 SCAN simulations suggest the presence of hydrophobic effects. This implies that water molecules exhibit a more ordered structure around the hydroxyl radical, leading to a slower decay in the RDF and an enhanced solvation shell, as shown by the prominent peaks in the hydrogen bond angle distribution, compared to r 2 SCAN simulations.
Conversely, the distribution of hydrogen bond angles where the radical acts as acceptor is much more uniform. Although the r 2 SCAN functional simulations showed a preference for the radical to act as a hydrogen bond acceptor, no specific orientation can be identified. Consequently, the more flexible and weaker acceptor hydrogen bonds do not contribute to a more structured and enhanced hydrogen bonding network.

Mulliken population analysis. In 2005, VandeVondele and Sprick29
proposed an electronic definition of the hemibond based on spin delocalization (or ''spin leakage'') from the radical to the solvent. For this, a Mulliken population analysis was performed every 4 fs for the 30 ps trajectory of the production run. The ''spin charge'', which distinguishes between hemibonding and hydrogen bonding, is calculated from the difference in Mulliken spin densities obtained from the a and b density matrix.
It is worth mentioning that for DC-r 2 SCAN, the spin charge q A on atom A is determined using the gradient consistent Mulliken analysis In this modified formulation of the Mulliken charge, 71 the density matrix is replaced by the sum of the reference density and the response density, given by P = P + P z . This simple metric plays a crucial role in characterizing the obstructed hemibond structure and has recently been used in the hybrid functional analysis of the aqueous hydroxyl radical by Rana and Herbert. 28 Fig. 4 Distribution of the hydrogen bond angles between solvent and the OH radical, where the radical acts as donor (upper panels), and acceptor (lower panels).
In Fig. 5 we plot the distribution of the spin density for r 2 SCAN and DC-r 2 SCAN. The r 2 SCAN distribution exhibits a stretched-out tail reaching approximately 0.74 in electron charge units, indicating significant spin leakage from the radical to the solvent. Conversely, the spin density distribution in DC-r 2 SCAN displays a narrow peak at 1.0, confirming the absence of any hemibonded configuration in the O*-O radial distribution function in Fig. 2.
Furthermore, Fig. 6 shows the net average charge of the radical, which are also calculated from summing the Mulliken charges on the atoms of the radical. The net average charge is around 0.057 for r 2 SCAN and approximately 0.017 for DCr 2 SCAN. This outcome is consistent with the reduced likelihood of hydrogen bond formation in DC-r 2 SCAN, which maintains neutrality in the radical species.
The charge and spin charge outcomes align with the selfinteraction corrected BLYP results reported by Vandevondele and Sprick, 29 reinforcing the possibility that the source of the hemibonded configuration can be attributed to self-interaction errors. The utilization of the Hartree-Fock density instead of the self-consistent Kohn-Sham r 2 SCAN density helps alleviate density-induced errors responsible for the formation of a (2c-3e)bonded interaction between the radical and water.

Dipole moment.
To explore the behavior of the hydroxyl radical and the solvent molecules, we calculated their molecular dipole moments using Voronoi tessellation. 72 The trajectory was sampled every 4 fs during the 30 ps production run. Similar to the Mulliken charge analysis above, it is important to note that for the DC-r 2 SCAN method, the dipole moment is calculated from the gradient consistent density, which incorporates the reference density and the response density.
For r 2 SCAN and DC-r 2 SCAN, the average dipole moment of the radical is 2.09 AE 0.2 D and 2.05 AE 0.2 D, respectively. Notably, we observe an induced dipole moment by the solvent on the radical of 0.423 D and 0.382 D for r 2 SCAN and DC-r 2 SCAN, respectively, based on the reported experimental dipole moment of 1.668 D for the radical in the gas phase. 73 The average dipole moment of the solvent water molecules are 2.249 AE 0.2 D and 2.282 AE 0.2 D for r 2 SCAN and DC-r 2 SCAN, respectively. However, these outcomes underestimate the experimental dipole moment  of 2.95 D. 74 The dipole moment distribution of the solvated hydroxyl radical is illustrated in the ESI. † To address the observed underestimation of the experimental dipole moment, it is important to examine the potential impact of dispersion corrections. The absence of these corrections in our simulations, which account for long-range van der Waals interactions, may contribute to the deviation from the experimental value. By enhancing the attractive forces between water molecules, dispersion corrections can lead to a more ordered structure, increasing the alignment of dipole moments and potentially yielding larger average dipole values. Moreover, the overbinding introduced by dispersion corrections can influence the hydrogen bond network, altering the distribution of dipole moments and consequently impacting the observed average dipole moment.
However, it is important to note that incorporating dispersion corrections, particularly in SCAN and DC-SCAN, has been shown to result in overbinding and compromise energetics. 8,9 This trade-off suggests that accurately capturing the dipole moment may come at the expense of deteriorating the system's structural representation and hydrogen bonding network. Exploring alternative approaches, such as utilizing maximally localized Wannier centers, may offer a solution for determining molecular dipole moments, although incorporating the response density in this formalism presents its own challenges.

Solvated sulfanyl radical
The analysis of the solvated mercapto radical (sulfanyl) SH follows the foodsteps of the solvated hydroxyl radical analysis. Given the lower electronegativity of the Sulfur atom (2.58) compared to oxygen (3.50), we expect -by design -fewer interactions with the solvent. 5.2.1 Radial distribution function. In Fig. 7 we plot the radial distribution functions between the atoms of the radical and the solvent. In the r 2 SCAN case, a pronounced shoulder with onset at 2.4 Å can be observed. The identification of a similar shoulder in DC-r 2 SCAN is not straightforward. An early onset of the first solvation shell can be identified at 2.78 Å, without forming a shoulder.
In analogy to the solvated hydroxyl radical, the various O-H pair correlation functions in Fig. 7, reveal a first insight on the hydrogen bonding network, and show its, in comparison, much reduced complexity.
The RDF where the SH radical acts as donor (O-H*) does not feature a clearly defined first solvation shell, but instead reaches a peak at 2.82 Å and 3.31 Å for r 2 SCAN and DC-r 2 SCAN, respectively. In the case where SH acts as acceptor (S*-H), a weak hydrogen bond can be identified at 3.55 Å and 3.62 Å. We observe that the hydrogen bond length is shorter, and therefore stronger, when the radical is acting as the donor, as opposed to when it is the acceptor.

Hydrogen bond analysis.
We investigated the distribution of hydrogen bonds between SH and the solvent to gain insights into their behavior. Both r 2 SCAN and DC-r 2 SCAN analyses consistently revealed a clear dominance of the case where the SH radical forms a single hydrogen bond, accounting for 83% of occurrences, while the twofold case constituted a smaller fraction at 17%. Notably, in both scenarios, the sulfanyl radical solely acts as the acceptor in these hydrogen bond interactions. For a visual representation of this distribution, we refer readers to Fig. S1 in the ESI. † The absence of hydrogen bonds where the sulfanyl radical acts as donor can also be identified from the (O-H*) radial distribution function. In gas phase CCSD(T)/aug-cc-pV5Z simulations, the hydrogen bond length of the donor bond S*-H*Á Á ÁO is 2.19 Å. 75 While the onset of the RDF is at approximately 1.9 Å for both methods, the increase is linear until r = 3.29 Å and no peak around 2.19 Å can be identified.
In gas phase simulations at the CCSD(T)/aug-cc-pV5Z level, the hydrogen bond length, where sulfanyl acts as the acceptor, is 2.62 Å. While no distinct peak can be identified around this value, the RDF shows a shoulder in both methods, ranging from 2.2 Å to 2.8 Å. The rate of increase is slower than the curve leading to the peak at 3.78 Å and 3.51 Å.
The hydrogen bond peak between the solvent molecules (O-H) is located at 1.85 Å and 1.88 Å for r 2 SCAN and DC-r 2 SCAN, respectively. This result reproduces the solvent in the hydroxyl system. Similarly, the average number of hydrogen bonds between water molecules are 3.56 and 3.53 for r 2 SCAN and DC-r 2 SCAN, respectively.

Mulliken population analysis.
We employ again the electronic definition of the hemibond via the spin charge distribution to identify the hemibonded configuration S'O between the radical and the solvent. The distribution of the spin charge is presented in Fig. 8 for r 2 SCAN and DC-r 2 SCAN, respectively. For the r 2 SCAN method, the spin charge distribution exhibits a tail extending to 0.8. This amounts to a similar, but less pronounced, spin leakage effect that was observed for the solvated OH radical. Similarly, the net average charge of the radical, displayed in the ESI † is 0.057 for r 2 SCAN, whereas it is approximately 0.017 for DC-r 2 SCAN. The neutrality of the radicals within the DC-r 2 SCAN scheme is consistent with its reduced tendency to form hydrogen bonds.

Dipole moment.
The distribution of the dipole moment on the SH radical, shown in the ESI, † reveal that the average dipole moment of the radical for r 2 SCAN and DC-r 2 SCAN are 0.968 AE 0.2 D and 0.818 AE 0.2 D, respectively. Given the reported experimental dipole moment of the SH radical of 0.758 D in the gas phase, 73 we find an induced dipole moment by the solvent on the radical of 0.210 D and 0.060 D for r 2 SCAN and DC-r 2 SCAN, respectively. The origin of the smaller predicted molecular dipole moment on the sulfanyl radical with DC-r 2 SCAN is not entirely clear. The average molecular dipole moment on the water molecules is 2.24 D and 2.27 D with r 2 SCAN and DC-r 2 SCAN, respectively. This reproduces the average values encountered for the solvent.
The origin of the smaller predicted molecular dipole moment on the sulfanyl radical with DC-r 2 SCAN is not entirely clear. The average molecular dipole moment on the water molecules is 2.24 D and 2.27 D with r 2 SCAN and DC-r 2 SCAN, respectively. This reproduces the average values encountered for the solvent.

Hemibonding in radical cation clusters
Sulfur-containing radical cation clusters have recently received significant attention due to their potential relevance in various biological processes, such as redox reactions, enzyme catalysis, and oxidative stress, where sulfur-centered reactive species play essential roles. 76 Investigating these clusters provides insights into their reactivity, stability, and implications for biological systems. Notably, clusters such as [CH 3 S'CH 3 SH] + , (CH 3 SH) 2 + , and (SH 2 ) n + exhibit a strong tendency to form (3e-2c) hemibonds, as demonstrated in prior studies. 35,77,78 The high reactivity of these hemibonded species limits however direct experimental characterization. In 2018, Xie et al. 33 explored the possible structures of [CH 3 S'CH 3 SH] + in the gas phase. They identified two hemibonded (1a and 1b) and one hydrogen-bonded structure (1c). The geometries of these structures were optimized using B3LYP-D3(BJ)/aug-cc-pVTZ calculations. The two hemibonded structures were found to have similar energies, with (1b) being 0.024 kcal mol À1 higher in energy compared to (1a). On the other hand, the hydrogen-bonded structure exhibited a significantly higher relative energy of 1.936 kcal mol À1 compared to (1a).
To verify the ability of r 2 SCAN and DC-r 2 SCAN to distinguish between hydrogen and hemibonded configurations in the radical cation clusters (CH 3 SH) 2 + and (SH 2 ) n + , we performed geometry optimizations of these structures. The computational details can be found in the ESI. † For (CH 3 SH) 2 + , the hemibonded structure (1a) remained the most stable for r 2 SCAN and DC-r 2 SCAN. The structures (1a) and (1b) converged to the same geometry. The hydrogen-bonded equilibrium structure (1c) transitioned to the hemibonded motif during optimization with the r 2 SCAN functional. In DCr 2 SCAN, the relative energy of the hydrogen-bonded structure was 7.72 kcal mol À1 higher in comparison to the hemibonded structure. The same procedure was applied to the hydrogenbonded and hemi-bonded geometries of (SH 2 ) + 2. 34 The hydrogen-bonded r 2 SCAN and DC-r 2 SCAN optimized geometries were 15.87 kcal mol À1 and 15.99 kcal mol À1 higher in energy than the hemi-bonded dimers, respectively. The results for both radical cation clusters confirm the ability of DCr 2 SCAN to differentiate between hydrogen-bonded and hemibonded configurations and predict the relative stability of the hemibonded motif in gas-phase calculations. The identification of hemi-bonded and hydrogen-bonded motifs becomes significantly more intricate when solvation effects need to be considered, since these configurations are already in competition in the gas-phase. 34,79,80 Notably, neutral hemibonded radicals have been detected in solvated systems, 33 which suggests a potential stabilization mechanism for these (3e-2c)-bonded entities upon solvation. The introduction of solvent molecules introduces novel intermolecular patterns for the hemibonded species, as they can participate in interactions involving hydrogen bonding or even hemibonding with the surrounding solvent molecules.
So far, only the effects of microsolvation have been studied for the closely related CH 3 SH molecule in heterodimers, 38 and the radical cation dimers (CH 3 SH) 2 + . 81 For the latter, it was found that the hemibonded (CH 3 SH) 2 + ion core was maintained when a single molecule of H 2 O, (CH 3 ) 2 CO, or CH 3 SH was bound, whereas the hemibond was broken by a NH 3 molecule, because the proton transfer led to a more stable hydrogenbonded structure. While solvation effects have been studied for lone SH 2 molecules, 82-84 the effects of solvation effects on (SH 2 ) 2 + have not been investigated yet.
To gain a comprehensive understanding of the behavior of radical cation clusters under realistic solvation conditions, we placed the [CH 3 S'CH 3 SH] + , and [HS'SH 2 ] + clusters in a solvation environment consisting of 64 water molecules, using the same simulation protocol as for the hydrated OH and SH radicals.
5.3.1 Radial distribution function. In Fig. 9 we plot the radial distribution functions of the [CH 3 S'CH 3 SH] + cluster.
The most prominent observation is the large peak centered at 3.02 Å in the S*-S* RDF of the r 2 SCAN simulation. The peak position slightly underestimates previously reported hemibond lengths of 3.12 Å and 3.28 Å between the sulfur atoms of the [CH 3 S'CH 3 SH] + cluster in the gas phase. 33 The comparably large peak height indicates a high probability of finding this configuration throughout the r 2 SCAN simulation.
Interestingly, the reported peak is absent in the DC-r 2 SCAN simulations. Since the DC-r 2 SCAN MD simulation was started after 15 ps of the r 2 SCAN AIMD run, where the hemibonded motif existed, further investigation is required. To examine the time-dependent changes, we plot the incremental S-S radial distribution function (RDF) in Fig. 10, with each data point corresponding to a 5 ps interval. We are able to confirm that with r 2 SCAN, the hemibond continuously exists. In the simulation with DC-r 2 SCAN, the hemibonded motif briefly exists at the beginning of the AIMD run as a remnant of the starting r 2 SCAN geometry. Within a few picoseconds, the hemibonded structure is lost and is not recovered during the remainder of the 30 ps trajectory.
The other RDFs in Fig. 9 describe the solvent interaction with the SH and CH 3 S moiety. The S*-O RDF provides insights into the solvation behavior of the [CH 3 S'CH 3 SH] + cluster in water. While the onset of the first solvation shell is similar in both simulations, notable differences can be observed in the peak positions and intensities of the RDF. In the r 2 SCAN simulation, the first peak occurred at 3.28 Å, extending into a plateau. This indicates a poorly defined second solvation shell. In contrast, the DC-r 2 SCAN simulation shows a much broader and more intense first peak at 3.59 Å. These differences suggest that in DC-r 2 SCAN, the cluster molecules are slightly more surrounded by water than with r 2 SCAN, potentially leading to the destabilization of the hemibond.
Given its electronic configuration, only the CH 3 SH species is capable of acting as a hydrogen bond donor and thus appears in the O-H* RDF. The results from r 2 SCAN and DC-r 2 SCAN are comparable, as no distinct peaks can be identified. However, both the SH moiety of CH 3 SH and the CH 3 S moiety can act as hydrogen bond acceptors. In the corresponding RDF (S*-H), r 2 SCAN shows a peak at 2.24 Å, indicating a slight preference for forming acceptor bonds. This preference may contribute to maintaining the stability of the hemibond, resulting in a more well-defined hydrogen bond network between the solvent and solute. In contrast, DC-r 2 SCAN does not exhibit this peak, suggesting a lower likelihood of forming acceptor hydrogen bonds. This finding is intriguing, considering the seemingly stronger solvent-solute interactions observed in the S*-O RDF.
In Fig. 11 we show the radial distribution functions of [HS'SH 2 ] + . In the S*-O RDF, r 2 SCAN has a first peak at 3.705 Å and DC-r 2 SCAN at 3.525 Å. This length corresponds to the intermolecular bond lengths between the sulfur atoms of the solvated cluster molecules and H 2 O, when they are hydrogen bonded and H 2 S acts as donor 3.56 Å or as acceptor 3.49 Å. 85 The radial distribution function (O-H*) and (S*-H), where the radical cation cluster may act as hydrogen bond donor or acceptor, respectively, show the same trends for r 2 SCAN and DC-r 2 SCAN. With either method, no distinct hydrogen bond donor or acceptor peak can be identified. However, the (O-H*) distribution reaches a maximum at 3.65 Å, while the peaks for (S*-H) is located at 3.82 Å.
In the RDF for S*-S*, a large peak is located at 2.915 Å for the r 2 SCAN functional. Although the DC-r 2 SCAN method also displays a peak at 3.036 Å, its amplitude is considerably smaller. Notably, Ghanty et al. 79 reported a hemibond distance of the sulfur-sulfur ion core of 2.899 Å using the B3LYP method. Additionally, Do and Besley 34 found the hemibond distance in the closely related (H 2 S) 2 + system to be 2.83 Å, whereas the hydrogen-bonded distance between S*-S* was  To initiate the DC-r 2 SCAN production run, we used a geometry obtained after 15 ps of the r 2 SCAN trajectory, following the approach employed for [CH 3 S'CH 3 SH] + . By plotting the incremental RDF in 5 ps intervals in Fig. 12, we observe that the hemibonded motif does not exist after the thermalization in r 2 SCAN, and is absent during the first 15 ps. However, the peak gradually reappears after 15 ps and grows in strength, becoming the dominant configuration from 20 ps onward. Conversely, the hemibonded motif appears to be the prevailing configuration between the sulfur ionic core in DC-r 2 SCAN in the first 5 ps of the simulation but subsequently loses in strength as the simulation progresses before it disappears in the interval ranging from 10 to 15 ps. Notably, in the final   To gain further insights into the initial 10 ps following the 15 ps mark of the r 2 SCAN simulation, we examined the incremental RDF at 1 ps intervals, starting from 15 ps for r 2 SCAN and 0 ps for DC-r 2 SCAN. In r 2 SCAN, the peak at 20 Å gradually strengthens, taking 5 ps to develop the hemibonded motif. Similarly, in DC-r 2 SCAN, the hemibond peak emerges after 3 ps, reaches its maximum amplitude at 6 ps, and subsequently diminishes. In contrast, for [CH 3 S'CH 3 SH] + , the hemibonded motif in the DC-r 2 SCAN simulation, a remnant of the r 2 SCAN thermalization, disappeared within the initial 5 ps and made no reappearance.
These findings suggest several intriguing possibilities. First, the continuous presence of the hemibond with r 2 SCAN in [CH 3 S'CH 3 SH] + and its immediate disappearance in DC-r 2 SCAN may be attributed to the strong localization of DC-r 2 SCAN, promoted by the underlying Hartree-Fock density. Second, the formation of the hemibond with both r 2 SCAN and DC-r 2 SCAN in [HS'SH 2 ] + and its spontaneous appearance imply that the hemibond does not appear to be stabilized by solvation; instead, it competes strongly against other interactions with the surrounding solvent. Finally, the hemibond is likely to be formed and destroyed repeatedly, with a frequency that extends beyond the available simulation time frame of 30 ps, emphasizing the dynamic and transient nature of the hemibond formation in solvated radical cation clusters.
5.3.2 Mulliken population analysis. Following the spin charge analysis of the solvated radicals, we plot in Fig. 14 the distribution of the spin density on the sulfur atoms of CH 3 S and the SH moiety of CH 3 SH. We observe significant spin leakage in the r 2 SCAN functional simulations from the sulfur atom of CH 3 S (average of 0.78) to the SH moiety of CH 3 SH (average of 0.25). The corresponding spin density distribution confirms the absence of the hemibonded ion core, as anticipated from the radial distribution function. Notably, the spin density primarily resides on the sulfur atom of CH 3 S, exhibiting an average value of 1.03, while the SH moiety of CH 3 SH remains devoid of any spin charge.
In Fig. 15 we show the temporal evolution of the Mulliken spin charge, which remains constant throughout the entire simulation for both r 2 SCAN and DC-r 2 SCAN simulations. The absence of any spin leakage in DC-r 2 SCAN is in good agreement with the rapid disappearance of the hemibonded motif depicted in Fig. 10.
The spin charge distribution of the solvated [HS'SH 2 ] + in Fig. 16 reveals the appearance of the hemibonded sulfur-sulfur ion core in the r 2 SCAN simulations. The spin leakage is however much less pronounced in comparison to the one observed for the [CH 3 S'CH 3 SH] + cluster with the r 2 SCAN functional. This is because in contrast to [CH 3 S'CH 3 SH] + , the hemibonded motif only emerges slowly after 15 ps, instead of being present during the entire length of the production run.
On the other hand, the DC-r 2 SCAN spin charge appears to be located entirely on the sulfur atom of the SH moiety, in Fig. 16, despite the short-lived emergence of the hemibond after 18 ps, that could be seen in the incremental RDF Fig. 12 and 13. As for r 2 SCAN, we suspect that the averaged spin charge over the entire production run obscures the dynamical evolution. For this reason, we plot the spin charge as a function of time in Fig. 17. Between 19 ps and 20 ps, the simulation performed with the r 2 SCAN functional shows the characteristic occurrence of spin leakage, confirming the existence of a hemibond between the sulfur atoms. Even after its formation, the hemibond exhibits strong competition with the hydrogen-bonded motif. This competition is evident through the localization of spin charge on the SH moiety between 23 and 24 ps, as well as from 25.5 to 27.5 ps. The impact of this competition can be observed in Fig. 12 by the reduced peak height of the RDF between 25 and 30 ps, indicating a weakening of the hemibond during that timeframe.
Remarkably, the incremental RDF plot of DC-r 2 SCAN (Fig. 13) exhibits the appearance of a peak around 3 Å, which reaches its maximum intensity between 4 and 6 ps before vanishing. Correspondingly, in Fig. 17, there is a slight yet discernible deviation of the spin charge from 1.06 during the    The hydrogen bond analysis of the cation clusters follows the same geometric criterion introduced earlier for the solvated radicals. We limit our analysis here on the [CH 3 S'CH 3 SH] + cluster. In Fig. 18 we show the number of hydrogen bonds formed with the sulfur atom of the CH 3 S and SH moieties in the left and right panels, respectively. It is evident that in all cases at least one hydrogen bond is formed.
Notably, the sulfur atom of the CH 3 S moiety can only act as a hydrogen bond acceptor. The r 2 SCAN functional exhibits a 20% higher likelihood of encountering two acceptor hydrogen bonds compared to DC-r 2 SCAN, supporting the presence of the observed peak in the S*-H radial distribution function (see Fig. 9). The peak can now be linked to the hydrogen bonding network, characterized by multiple acceptor interactions involving the CH 3 S moiety. However, the sulfur atom in CH 3 SH also acts solely as a hydrogen bond acceptor, as shown by the absence of any peak in the O-H* radial distribution function, indicating a lack of donor hydrogen bonding interactions.
Moreover, it is crucial to emphasize that there is an absence of hydrogen bond formation between the CH 3 S moiety and CH 3 SH. This observation aligns with previous theoretical investigations 77,78 as well as experimental studies 86 conducted on analogous methanethiol clusters (CH 3 SH) n + with n = 2-5.
These studies have revealed that S-HÁ Á ÁS hydrogen bonding between clusters occurs exclusively for n = 3-5, while no such bonding has been observed in radical cation dimers.

Conclusions
This work provides a theoretical study of self-interaction effects in fully solvated and strongly-correlated cationic molecular clusters. The molecular dynamics application of the variational formulation of ADMM accelerated density-corrected DFT has been presented, specifically using the DC-r 2 SCAN method. The method was successfully validated for bulk liquid water, demonstrating its ability to accurately describe the structure. The DC-r 2 SCAN method not only reproduces the results  obtained with r 2 SCAN, but also improves upon them by addressing the shortcomings of r 2 SCAN, particularly related to spin density. This nearly self-interaction free method allows for the observation of hydroxyl radical solvation beyond the GGA approach and accurately reproduces results obtained from self-interaction corrected electronic structure theory methods.
Our study primarily focused on revisiting the solvated hydroxyl radical, and the solvated sulfanyl radical. Subsequently, we extended our investigation to solvated cationic molecular clusters, namely [CH 3 S'CH 3 SH] + and [HS'SH 2 ] + . In the case of [CH 3 S'CH 3 SH] + , the DC-r 2 SCAN simulation revealed the destruction of the initially existing hemibond, which was maintained in the r 2 SCAN simulation. The observed spin leakage and strong spin charge localization during the DC-r 2 SCAN simulation supported this finding.
For the solvated [HS'SH 2 ] + cluster, the hemibonded motif was lost during thermalization, but r 2 SCAN was able to restore it, evidenced by the observation of spin leakage. DC-r 2 SCAN also predicted the formation of a hemibond, although with weaker structural and electronic signatures, indicating a stronger competition from solvent-solute interactions on the survivability of the hemibonded motif.
Future directions for this research include extending the DC-DFT framework by incorporating long-range dispersion corrections and exploring the use of hybrid functionals to improve accuracy. The computational framework developed in this work holds promise for efficiently modeling condensed phase systems with complex correlation effects.
In summary, we demonstrate the effectiveness of selfinteraction corrected electronic structure theory, particularly the DC-r 2 SCAN method, in accurately describing solvated cationic molecular clusters.

Data availability
The CP2K software package can be found at https://github.com/ cp2k/cp2k. Data for this paper, including the AIMD trajectories of the solvated radicals and solvated cationic molecular clusters are available at the Materials Cloud Archive at https://archive. materialscloud.org/record/2023.85.

Conflicts of interest
There are no conflicts to declare.