Decoding ortho regiospecificity and high endo stereoselectivity in pyrazole synthesis via the activation/strain model

Abstract

A comprehensive density functional theory investigation was conducted to elucidate the regio- and stereoselectivity of the 1,3-dipolar cycloaddition reaction (13-DCR) between a diazo compound (1,3-dipole D-1) and an electron-deficient nitroethylene derivative (dipolarophile Dph-2), which serves as the key initiating step in a recently reported domino synthesis of pyrazole derivatives. Experimentally performed in dichloromethane at 80 °C, the reaction exhibits exclusive formation of the ortho regioisomer. Calculated activation free energies and rate constants quantitatively reproduce the experimentally observed complete ortho regioselectivity and high endo stereoselectivity. Activation strain model (ASM) analysis reveals that the regioselectivity is primarily governed by stabilizing interaction energies. Subsequent energy decomposition analysis using the recently developed sobEDA method identifies orbital interactions, particularly HOMO_D-1 to LUMO_Dph-2 charge transfer, as the main contributors. The Extended Transition State-Natural Orbitals for Chemical Valence (ETS-NOCV) analysis confirms this dominant orbital interaction, rationalizing the ortho selectivity. The stereoselectivity favoring the endo pathway is attributed to lower total strain energy, mainly originating from the flexibility of the D-1 fragment. This combined computational study offers a detailed mechanistic understanding of the 13-DCR regio- and stereoselectivity, providing results in excellent agreement with experimental outcomes.

---

1. Introduction

Heterocyclic compounds serve as essential building blocks and are widely found in numerous pharmaceuticals.¹ Among them, nitrogen-containing heterocycles hold exceptional industrial and commercial value, playing key roles in pharmaceuticals, agrochemicals, and natural products.² Within this group, pyrazoles—five-membered aromatic rings featuring two adjacent nitrogen atoms—stand out for their diverse biological activities, including anticonvulsant,³ antibacterial,⁴ antidepressant,⁵ and analgesic properties.⁶ Notably, several top-selling drugs such as celecoxib and sildenafil contain pyrazole moieties.² Given these valuable functions, the search for efficient methods to synthesize pyrazole derivatives has gained increasing attention.

Two main synthetic strategies are commonly employed to construct pyrazole scaffolds: cyclocondensation of hydrazines with 1,3-dicarbonyl compounds or their derivatives and the 1,3-dipolar cycloaddition reaction (13-DCR) of diazo compounds with alkenes or alkynes.^7,8 Among these, 13-DCRs are highly valued in organic synthesis due to their excellent atom economy, high yields, and remarkable regio- and stereoselectivity.⁹ Notably, the 13-DCR of diazo compounds as an established class of propargylic-type 1,3-dipoles toward alkenes, first reported by Buchner in 1889,¹⁰ offers an efficient route to pyrazole synthesis since no σ-bond cleavage occurs during the process.²

Nitroalkenes, characterized by their strong electrophilic nature due to the highly electron-withdrawing nitro group, are well-suited for various organic transformations, including 13-DCRs.¹¹ In this context, Kumar and co-workers¹² have recently reported the synthesis of a pyrazole derivative (E) through a domino sequence initiated by a 13-DCR between diazoacetate (A) and π-deficient nitroalkene (B), yielding the [3 + 2] cycloadduct (C). This intermediate then undergoes a 1,3-hydrogen shift to form intermediate (D), followed by a concomitant removal of HOCl and NaNO₂, producing the final pyrazole derivative (E) in a good yield of 82% ¹² (Scheme 1). Importantly, Kumar and co-workers¹² highlighted that E is obtained exclusively as the isomer bearing the Ph substituent ortho to the CO₂Et group, with no experimental evidence for the formation or isolation of the meta isomer.¹² Based on these experimental findings, it can be concluded that the 13-DCR step in this domino sequence proceeds with complete regioselectivity (regiospecificity), placing the Ph and CO₂Et substituents ortho to each other in intermediate C, which subsequently leads to intermediate D and the final product E (see Scheme 1).

Scheme 1 The domino reaction experimentally investigated by Kumar et al.¹² to synthesize pyrazole derivative E.

However, the sp²-hybridized nature of all carbons in the five-membered ring of compound E does not readily disclose the origin of the exo/endo stereoselectivity in the preceding 13-DCR step. This raises two fundamental questions regarding the 13-DCR between A and B in Scheme 1: (i) why does this reaction exhibit complete regioselectivity? (ii) Does it display any specific stereoselectivity, and if so, what factors govern this preference? To answer these questions, we carried out a detailed density functional theory (DFT) investigation on the 13-DCR between diazoacetate A (hereafter referred to as 1,3-dipole D-1) and nitroalkene B (hereafter referred to as dienophile Dph-2), which constitutes the key initiating step of the domino process shown in Scheme 1.

2. Computational details

All geometry optimizations followed by frequency calculations were performed in dichloromethane (DCM) at 80 °C (353.15 K) to fully replicate the experimental conditions,¹² using the M06-2X-D3/6-31G(d) level of theory as implemented in the Gaussian 16 software package.¹³ This computational level employs the M06-2X meta-generalized gradient approximation (meta-GGA) exchange–correlation functional,¹⁴ combined with Grimme's DFT-D3 dispersion correction¹⁵ (using zero-damping, denoted as D3) and the standard 6-31G(d) basis set.¹⁶

Solvent effects were modeled implicitly via the polarizable continuum model (PCM) based on the Tomasi solvation approach¹⁷ within the self-consistent reaction field (SCRF) framework.¹⁸ Frequency calculations provide not only thermochemical data at the desired temperature and pressure but also confirm the nature of each stationary point on the potential energy surface (PES). Specifically, reactants and products were identified by the absence of imaginary frequencies, while transition states (TSs), located as first-order saddle points, displayed exactly one imaginary frequency. Additionally, intrinsic reaction coordinate (IRC) calculations¹⁹ were performed for each TS using the second-order González–Schlegel integration method,²⁰ tracing both forward and backward paths to confirm that each TS correctly connects two relevant minima and also to check for the existence of any stable intermediates on the PES.

To refine electronic energies, single-point (SP) calculations were carried out using the double-hybrid rev-DSD-PBEP86 functional²¹ with Grimme's D4 dispersion correction²² and the ma-TZVPP basis set over the M06-2X-D3/6-31G(d)-optimized geometries, employing the ORCA 6.0.1 program package.^23,24 The ma-TZVPP basis set is a “minimally augmented” version of def2-TZVPP²⁵ with additional s and p diffuse functions on non-hydrogen atoms. It is also worth noting that ORCA employs the conductor-like polarizable continuum model (CPCM)¹⁷ to account for implicit solvent effects ensuring numerical stability and consistency across the dataset for the chosen DFT levels. Subsequently, Gibbs free energies were computed using Tian Lu's “Shermo” code,²⁶ incorporating the refined electronic energies. This code applies Grimme's quasi-rigid-rotor harmonic oscillator (RRHO) model to properly account for entropy contributions from low-frequency vibrations (below 150 cm⁻¹), yielding more accurate Gibbs free energy estimates. Furthermore, Shermo automatically adjusts Gibbs free energies with a correction factor of 2.36 kcal mol⁻¹ at 353.15 K to convert gas-phase standard-state values (1 atm) to solution-phase standard-state values (1 M) at the specified temperature.²⁷

The stability of the M06-2X-D3/6-31G(d) closed-shell wavefunctions for the 1,3-dipole D-1 and TSs involved in the studied 13-DCR pathways was verified using the “stable” keyword. Unless explicitly stated otherwise, no wavefunction instabilities were detected.

The temperature-dependent rate constant (k_T) for each reaction channel in the 13-DCR between D-1 and Dph-2 was calculated using the following expression²⁸

where K_B, T, h, C⁰, m, and R refer to Boltzmann's constant, absolute temperature, Planck's constant, standard molar concentration (1 mol L⁻¹), the molecularity of the reaction, and the universal gas constant, respectively. Here, ΔG°^‡ represents the standard Gibbs free energy of activation for the reaction under study.

To investigate the origins of regio- and stereoselectivity in the studied 13-DCR, the activation/strain model (ASM) analysis, originally developed by Bickelhaupt and co-workers,^29,30 was performed at the M06-2X-D3/6-31G(d) level along the IRC profiles of the TSs involved in the competing reaction pathways. This analysis was conducted using a custom Microsoft Excel 2016 spreadsheet developed in-house. In addition, A. M. Shariati (one of the co-authors) developed a highly efficient, flexible, and user-friendly code that significantly accelerates and automates the geometry preparation process along the IRC paths—a typically time-consuming step in ASM analyses.

Energy decomposition analysis (EDA) was performed along the entire IRC profiles of the TSs involved in the competing reaction pathways using the “sobEDA” method, a dispersion-corrected DFT-based approach recently introduced by Tian Lu.³¹ All these calculations were carried out in a fully automated manner using custom in-house scripts. It should be noted that the sobEDA method employs Gaussian software for performing the underlying quantum chemical computations.

The Extended Transition State-Natural Orbitals for Chemical Valence (ETS-NOCV)³² analysis was conducted using the Amsterdam Density Functional (ADF) program package, version 2020.102,^33–35 at the BLYP-D3(BJ)/TZ2P level of theory. This computational setup has been recommended for ETS-NOCV studies on 13-DCR systems.³⁶ Deformation density plots (Δρ) were visualized via the ADF graphical user interface (GUI).^33–35

3. Results and discussion

3.1. PES exploration of 13-DCR between D-1 and Dph-2 in the presence of DCM at 353.15 K

Due to the non-symmetric structures of both reactants, the 13-DCR between D-1 and Dph-2 can proceed via two regioisomeric pathways—meta and ortho—each of which allows for two possible stereochemical approaches, exo and endo. These four distinct reaction pathways are denoted as meta-exo (mx), meta-endo (mn), ortho-exo (ox), and ortho-endo (on). Along the meta (ortho) pathway, the C1 atom of D-1 interacts with the C5 (C4) atom of Dph-2, placing the CO₂Et and Ph substituents at the meta (ortho) position within the resulting [3 + 2] formal cycloadduct (CA). In contrast, an endo approach places these substituents on the same side of the corresponding TS, leading to a cis configuration in CA, whereas an exo approach positions them on opposite sides, resulting in a trans configuration. These four reaction channels along with the atom numbering are illustrated in Scheme 2. The calculated activation and reaction electronic energies (ΔE^‡ and ΔE_rxn, respectively) for each pathway in DCM are also provided in Scheme 2.

Scheme 2 Four competitive pathways involved in the 13-DCR between D-1 and Dph-2. Atom numbering is given as well (see the text for details). Numerical values indicate M06-2X-D3/6-31G(d) computed relative total electronic energy (kcal mol⁻¹) in the presence of DCM.

As shown in Scheme 2, the activation barrier (ΔE^‡) for TS-mx and TS-mn is 13.81 and 13.75 kcal mol⁻¹, respectively, relative to the separated reactants D-1 and Dph-2. The negligible difference of 0.06 kcal mol⁻¹ between these two values suggests no significant stereoselectivity along the meta regioisomeric pathway. In contrast, the activation barrier for TS-ox and TS-on is 10.35 and 8.68 kcal mol⁻¹, respectively. The energy difference of 2.00 kcal mol⁻¹ between these two values, although within the typical error range of DFT methods, clearly indicates a strong preference for endo stereoselectivity along the ortho regioisomeric pathway. Nevertheless, as far as total electronic energy is concerned, the formation of all [3 + 2] cycloadducts—CA-mx, CA-mn, CA-ox, and CA-on—proceeds through irreversible and exergonic pathways, with substantial energy releases ranging from 27.4 to 30.5 kcal mol⁻¹. These results suggest that the cycloadducts formed via the meta regioisomeric pathway are slightly more stable than those obtained through the ortho pathway.

It is well established that bimolecular reactions are typically associated with a substantial decrease in both activation and reaction entropy, resulting in strongly negative values for ΔS^‡ and ΔS_rxn. Consequently, the TΔS term acts as a highly unfavorable factor for both the formation of the TS and the corresponding CA, particularly at elevated temperatures. Therefore, to achieve a reliable and meaningful portrait of the PES of such reaction energetics, it is essential to compute relative Gibbs free energies, which account for both energetic and entropic contributions. To this end, the activation (ΔG^‡) and reaction (ΔG_rxn) Gibbs free energies, along with the corresponding rate constants for both the forward (k_f(T)) and reverse (k_b(T)) directions, were computed for the four reaction channels (Scheme 2) of the 13-DCR between D-1 and Dph-2. These calculations were performed at the rev-DSD-PBEP86-D4-CPCM(DCM)/ma-def2-TZVPP//M06-2X-D3-PCM(DCM)/6-31G(d) computational level and the results are summarized in Table 1.

Table 1 Activation Gibbs free energies for the forward (ΔG^‡_f) and backward (ΔG^‡_b) reactions, reaction Gibbs free energies (ΔG_rxn), in kcal mol⁻¹, forward second-order rate constants, k_f(T) in mol⁻¹ L s⁻¹, and backward first-order rate constants, k_b(T) in s⁻¹, calculated at the rev-DSD-PBEP86-D4-CPCM(DCM)/ma-def2-TZVPP//M06-2X-D3-PCM(DCM)/6-31G(d) level for the four reaction pathways of the 13-DCR between D-1 and Dph-2 in DCM at 353.15 K

Pathway^a	ΔG^‡f ^b	ΔG^‡b ^b	ΔGrxn ^b	k f(T)	k b(T)
a See Scheme 2 for details. b Computed relative to the separate fully optimized D-1 and Dph-2.
meta-exo (mx)	33.53	35.21	−1.68	1.31 × 10^−8	1.19 × 10^−9
meta-endo (mn)	32.08	34.76	−2.68	1.03 × 10^−7	2.27 × 10^−9
ortho-exo (ox)	28.99	28.76	0.23	8.44 × 10^−6	1.17 × 10^−5
ortho-endo (on)	28.03	29.19	−1.16	3.31 × 10^−5	6.35 × 10^−6

---

The ΔG^‡ values reported in Table 1 clearly show that among the four competing TSs located on the PES of the studied 13-DCR, TS-on exhibits the lowest forward activation Gibbs free energy ΔG^‡_f of 28.03 kcal mol⁻¹, corresponding to the highest second-order rate constant, k_f(T) = 3.31 × 10⁻⁵ mol⁻¹ L s⁻¹. Thus, the ortho-endo (on) pathway is identified as the kinetically most favorable channel. This pathway leads to the exergonic formation of the formal [3 + 2] cycloadduct CA-on, which slightly lowers the reaction PES by 1.16 kcal mol⁻¹. The next favorable pathway is the ortho-exo (ox) channel, with TS-ox located 28.99 kcal mol⁻¹ above the reactants—only 0.96 kcal mol⁻¹ higher than TS-on. The corresponding k_f(T) value for this pathway is 8.44 × 10⁻⁶ mol⁻¹ L s⁻¹. The remaining two channels, meta-endo (mn) and meta-exo (mx), are significantly less favorable, with k_f(T) values of 1.03 × 10⁻⁷ and 1.31 × 10⁻⁸ mol⁻¹ L s⁻¹, respectively. The last column of Table 1 lists the first-order backward rate constants, k_b(T), for the studied pathways. As is shown, there are no significant differences between the forward and backward rate constants across the channels. Consistent with this, the ΔG_rxn values reported in Table 1 indicate that (in contrast to the prediction obtained according to the values of relative electronic energies given in Scheme 2) the cycloadducts are not substantially more stable than the separate reactants, suggesting that the studied 13-DCR operates under kinetic control.

The Gibbs free energy profiles for the four reaction pathways, depicted in Fig. 1, provide a clear and intuitive comparison of the reaction energetics in DCM at 353.15 K.

Fig. 1 The rev-DSD-PBEP86-D4-CPCM(DCM)/ma-def2-TZVPP//M06-2X-D3-PCM(DCM)/6-31G(d) relative Gibbs free energy profile associated with the four competitive reaction pathways converting separate D-1 and Dph-2 into related CAs along a 13-DCR in the presence of DCM at 353.15 K (see the text for details). A considerable difference of 4.05 kcal mol⁻¹ in the relative Gibbs free energies between the two most favorable TSs of the competing regioisomeric pathways indicates an entirely regiospecific 13-DCR of D-1 toward Dph-2.

3.2. Kinetic evaluation of regio- and stereoselectivity in the 13-DCR of D-1 with Dph-2 in DCM at 353.15 K

In order to quantitatively assess the regio- and stereoselectivity of the studied 13-DCR, we analyzed the computed rate constants for the four competing pathways presented in Scheme 2. This kinetic evaluation allows for a direct and reliable estimation of the preferred reaction channels under the employed conditions.

According to fundamental kinetic principles, the ortho-to-meta regioselectivity can be quantitatively estimated using the following expression:³⁷

where the numerator and denominator represent the combined molar concentrations of the ortho and meta cycloadducts, respectively. Since the concentration of each cycloadduct is directly proportional to the rate constant of the reaction channel through which it is formed, this equation can be rewritten as:

By substituting the computed rate constants from Table 1 into this equation, we obtain:

This result indicates that over 99% of the products correspond to the ortho regioisomer, demonstrating a practically exclusive ortho regioselectivity (i.e., ortho regiospecificity). This finding is in excellent agreement with the experimental observation that only the ortho-substituted pyrazole derivative E was isolated, with no evidence for the formation of the meta isomer (Scheme 1).¹² Similarly, the endo-to-exo stereoselectivity within the ortho regiospecific pathway can be estimated as:³⁷

This corresponds to approximately 79.6% endo selectivity, indicating a strong preference for the endo approach along the kinetically dominant ortho pathway under the studied conditions (DCM, 353.15 K).

To complement this kinetic selectivity analysis based on computed rate constants, we further performed a time-dependent kinetic simulation using the Concvar software, aiming to visualize the dynamic evolution of the competing pathways over time. The Concvar program developed by Tian Lu³⁸ allows for the simultaneous solution of differential rate equations for a network of elementary reactions, providing time-dependent concentration profiles, C = f(t), for all species involved in a given reaction system. In this study, we employed Concvar by inputting the forward and reverse rate constants k_f(T) and k_b(T) listed in Table 1 to simulate the time evolution of D-1, Dph-2, and the four CAs generated in the studied 13-DCR (see Scheme 2). Fig. 2 presents the resulting concentration–time profiles for the reaction in DCM at 353.15 K, assuming initial concentrations of 0.2 mol L⁻¹ for D-1 and 0.1 mol L⁻¹ for Dph-2 over a total simulation time of 10 hours (36 000 s).

Fig. 2 Time-dependent concentration profiles, C = f(t), illustrating the evolution of all CAs formed during the 13-DCR of D-1 with Dph-2 in DCM at 353.15 K.

A closer inspection of Fig. 2 reveals that, even after a prolonged reaction time of 10 hours, the only detectable species in the reaction mixture is CA-on, present at a relatively low concentration of approximately 0.019 mol L⁻¹. This low concentration directly results from the small rate constant associated with its formation (see k_f(T) values in Table 1). Moreover, the concentrations of the meta-derived CAs (CA-mx and CA-mn) remain entirely negligible across the whole simulation.

Notably, this simulation clearly confirms that the ortho-endo (on) pathway overwhelmingly dominates the reaction course, fully consistent with both the kinetic rate analysis and the experimental observations.¹²

3.3. Exploration of the origin of ortho regiospecificity and high endo stereoselectivity in the 13-DCR investigated

Following the detailed energetic analysis, the origin of the complete ortho regioselectivity (ortho regiospecificity) and high endo stereoselectivity was investigated through activation/strain model (ASM) analysis.

To this end, comparing the ASM profiles of the most favorable TSs in each competing pathway provides key insights. Specifically, comparing the ASM profiles of TS-mn and TS-on clarifies the origin of the observed ortho regiospecificity, while comparing that of TS-on and TS-ox (within the ortho pathway) elucidates the underlying cause of the endo stereoselectivity preference in the studied 13-DCR.

ASM is a fragment-based theoretical model that systematically traces the PES from equilibrium structures toward TSs and even beyond, offering a clear and detailed understanding of chemical reactivity.²⁹ In this model, the reaction energy profile, ΔE(ξ), along the reaction coordinate ξ, is dissected into two key components: the total strain energy, ΔE^str_tot(ξ), and the interaction energy, ΔE^int(ξ), as shown below:

ΔE(ξ) = ΔE^str_tot(ξ) + ΔE^int(ξ). (4)

The total strain energy itself arises from the structural deformation of individual reactants as they progress toward the TS. For a bimolecular 13-DCR, it is expressed as:

ΔE^str_tot(ξ) = ΔE^str_D(ξ) + ΔE^str_Dph(ξ) (5)

where the subscripts D and Dph refer to the 1,3-dipole and the dipolarophile, respectively.

Solvent effects can also be incorporated into the ASM framework by including the solvation contribution as follows:

ΔE^solution(ξ) = ΔE^str_tot (ξ) + ΔE^int(ξ) + ΔΔE^solvation(ξ) (6)

where ΔE^solution(ξ) represents the total reaction energy in solution, ΔΔE^solvation(ξ) accounts for the change in solvation energy along the reaction coordinate, and ΔE^solute(ξ) = ΔE^str_tot(ξ) + ΔE^int(ξ) describes the intrinsic energy profile of the solute in the absence of solvent effects.

It is important to emphasize that, when including solvation effects viaeqn (6), single-point energy calculations must be performed in the gas phase using solution-phase optimized geometries—meaning the geometries are optimized with solvation included (e.g., PCM), but the subsequent single-point energies are computed without the “SCRF” keyword, if Gaussian is employed as the quantum calculating software. This approach properly captures the differential solvation effects (ΔΔE^solvation) along the reaction coordinate. In addition, Truhlar and co-workers³⁹ reported that solvation energies calculated using the M05-2X functional combined with the SMD solvation model show excellent agreement with experimental data. Therefore, in the present investigation, the M05-2X-SMD(DCM)/6-31G(d)//M06-2X-D3-PCM(DCM)/6-31G(d) computational protocol was employed to evaluate values of ΔE^solvation.

The ASM analysis was performed along the IRC profiles of the two competing TS-on and TS-mn, corresponding to the ortho-endo and meta-endo pathways, respectively, which connect the separate reactants D-1 and Dph-2 to the [3 + 2] cycloadducts CA-on and CA-mn (Scheme 2). This analysis aims to clarify why the ortho-endo pathway is entirely favored over the meta-endo one, as indicated by the energetic results (Table 1). Fig. 3 illustrates the activation/strain plots along the whole IRC profile of TS-on and TS-mn. In this 13-DCR, the reaction coordinate ξ was defined as the average of the two forming single bond distances, namely the C1–C4 and N3–C5 bonds, along the meta-endo (mn) and ortho-endo (on) pathways (see Scheme 2 for atom numbering):

Fig. 3 The magnified M06-2X-D3/6-31G(d) activation/strain plot along the whole IRC profile of TS-on (solid curves) and TS-mn (dashed curves) involved in the ortho-endo (on) and meta-endo (mn) channels of 13-DCR, respectively, between D-1 and Dph-2 in the presence of DCM. The blue circles characterize the positions of TS-on and TS-mn, while the black vertical line at ξ = 2.18 Å distinguishes the position of “consistent geometry point (CGP)” along the X-axis (see the text for details).

It is worth noting that the nature of the interacting atoms differs between these two forming bonds (C–C and C–N), and the bond formation rates are not identical. Thus, monitoring the change of only one bond distance would not provide a complete and accurate mechanistic picture for such reactions. This justifies the use of the average of both distances, as expressed in eqn (7). It is also worth noting that Fig. 3 is plotted using magnified X- and Y-axis scales to enhance the visual separation between the curves corresponding to the ortho-endo and meta-endo pathways considered. As clearly shown in Fig. 3, the ΔE^str_tot(ξ) curve of the energetically favored ortho-endo pathway lies above that of the meta-endo pathway throughout the entire reaction coordinate, indicating that the ortho-endo pathway suffers from a greater destabilizing strain energy. Meanwhile, the ΔE^solvation(ξ) curves are completely superimposed, indicating that both pathways experience virtually identical solvation effects throughout the reaction course. This results in ΔΔE^solvation(ξ) = 0, confirming that solvation has no contribution in determining the regioselectivity. In contrast, the ΔE^int(ξ) curve of the ortho-endo pathway consistently lies below that of the meta-endo pathway, with a larger separation than that observed between the strain energy curves. This indicates that the ortho-endo pathway benefits from substantially stronger stabilizing interaction energy, which not only fully compensates its slightly higher strain energy but also ultimately drives its energetic preference over the meta-endo pathway.

Given that the interaction energy governs the complete ortho regioselectivity over the meta pathway in the studied 13-DCR, a further energy decomposition analysis (EDA) was performed to break down the interaction energy into its fundamental components, aiming to identify which specific factor(s) primarily drives this regioselectivity. Very recently, Tian Lu introduced a new energy decomposition analysis (EDA) method, known as sobEDA, along with its variant sobEDAw, which is specifically designed to dissect weak interaction energies.³¹ The key advantages of sobEDA and sobEDAw include their ease of implementation, low computational cost, broad applicability, and high accuracy, with errors comparable to those of high-level CCSD(T)/CBS and DFT-SAPT reference data.³¹ According to the sobEDAw, the interaction energy ΔE^int(ξ) is decomposed into four distinct components as expressed in eqn (8):³¹

ΔE^int(ξ) = ΔE_els(ξ) + ΔE_orb(ξ) + ΔE_disp(ξ) + ΔE_xrep(ξ) (8)

where ΔE_els(ξ) represents classical electrostatic interaction between monomers. ΔE_orb(ξ) accounts for orbital interactions, including both polarization (the response of each monomer to the electric field of the other) and charge transfer between the monomers. ΔE_disp(ξ) denotes the dispersion contribution, arising from the Coulomb correlation between electrons in one monomer with those in another one. Finally, ΔE_xrep(ξ) captures the exchange-repulsion effects, which stem from the overlap of monomer wavefunctions and the antisymmetric requirements of the fermionic nature of electrons in the dimer.³¹ Tian Lu has also demonstrated that when the sobEDA or sobEDAw method is applied at the B3LYP-D3(BJ)/6-311+G(2d,p) level of theory, the results closely match those obtained from high-level CCSD(T)/CBS and DFT-SAPT methods.³¹ Accordingly, we employed sobEDA at the B3LYP-D3(BJ)/6-311+G(2d,p) computational level for this study. As previously mentioned, all sobEDA analyses were performed fully automatically over the entire IRC profiles of TS-on and TS-mn using our in-house developed scripts, and the resulting decomposed energy profiles are illustrated in Fig. 4. As clearly shown in Fig. 4, the destabilizing exchange-repulsion interaction, ΔE_xrep(ξ), along the ortho-endo (on) pathway is larger than that along the meta-endo (mn) pathway. This observation is consistent with the higher total strain energy, ΔE^str_tot(ξ), associated with the ortho-endo pathway, as previously discussed in Fig. 3.

Fig. 4 B3LYP-D3(BJ)/6-311+G(2d,p) sobEDA analysis along the M06-2X-D3/6-31G(d) whole IRC profile of TS-on (solid curves) and TS-mn (dashed curves) involved in the ortho-endo (on) and meta-endo (mn) channels of 13-DCR, respectively, between D-1 and Dph-2 in the presence of DCM.

Therefore, the origin of the energetic preference for the ortho-endo pathway over the meta-endo one must lie in one or more of the stabilizing energy components, namely the electrostatic, ΔE_els(ξ), orbital, ΔE_orb(ξ), or dispersion, ΔE_disp(ξ), terms. To further investigate this point, we focused exclusively on the stabilizing interaction components presented in Fig. 4 to obtain a clearer view of their relative differences. For this purpose, we re-plotted the stabilizing terms in Fig. 5 by adjusting the scales of both the X- and Y-axes to narrower ranges, thereby enhancing the visual separation between the curves. As shown in Fig. 5, all stabilizing components of the interaction energy along the ortho-endo (on) pathway lie below those of the meta-endo (mn) pathway throughout the entire reaction coordinate. Notably, the largest separation between the two pathways is observed for the orbital interaction component, ΔE_orb(ξ), clearly identifying this term as the dominant factor responsible for the energetic preference of the ortho-endo pathway over the meta-endo pathway.

Fig. 5 Magnified B3LYP-D3(BJ)/6-311+G(2d,p) sobEDA analysis along the M06-2X-D3/6-31G(d) whole IRC profile of TS-on (solid curves) and TS-mn (dashed curves) involved in the ortho-endo (on) and meta-endo (mn) channels of 13-DCR, respectively, between D-1 and Dph-2 in the presence of DCM.

It is also worth noting that electrostatic interactions, ΔE_els(ξ), contribute as the second most significant stabilizing factor, although to a lesser extent.

To deepen these qualitative conclusions, the sobEDA-derived numerical values of the interaction energy and its stabilizing components—used to construct the curves in Fig. 4—are compiled in Table 2 for the CGP along both the ortho-endo (on) and meta-endo (mn) pathways. It should be pointed out that the CGP along a given reaction path ensures us that the computed energy components reproduce the same relative trends observed at the TSs. This criterion is particularly critical when the TSs of competing pathways occur at different positions along the reaction coordinate, thereby avoiding bias in the comparison. For the 13-DCR studied, an appropriate CGP (indicated by the black vertical line in Fig. 3) may be considered at ξ = 2.180 Å.

Table 2 B3LYP-D3(BJ)/6-311+G(2d,p) values of interaction energy, ΔE^int, its stabilizing (constructive) components, ΔE_els, ΔE_orb, and ΔE_disp, in kcal mol⁻¹, together with the difference of these components, ΔΔE_x,^a calculated based on the sobEDA methodology at the CGP along the competitive ortho-endo (on) and meta-endo (mn) pathways involved in the 13-DCR studied

Pathway^b	ΔE^int	ΔEels	ΔEorb	ΔEdisp	ΔΔEels	ΔΔEorb	ΔΔEdisp
a ΔΔEx = ΔEx(on) − ΔEx(mn) where ΔEx could be any of ΔEels, ΔEorb, or ΔEdis contributors. b See Scheme 2 for details.
ortho-endo (on)	−17.33	−56.98	−61.10	−27.95	0.00	0.00	0.00
meta-endo (mn)	−14.92	−54.54	−57.30	−26.28	−2.44	−3.80	−1.67

---

As presented in the second column of Table 2, the ortho-endo (on) pathway exhibits a more stabilizing interaction energy (−17.33 kcal mol⁻¹) than that of the meta-endo (mn) pathway (−14.92 kcal mol⁻¹), with a difference of 2.41 kcal mol⁻¹. The relative contributions (ΔΔE_x) were evaluated by taking the ortho-endo pathway as the reference, such that ΔΔE_x = ΔE_x(on) − ΔE_x(mn), where ΔE_x represents any of the electrostatic (ΔE_els), orbital (ΔE_orb), or dispersion (ΔE_dis) terms. Accordingly, a negative ΔΔE_x signifies that the corresponding energy component is more stabilizing in the ortho-endo pathway than in the meta-endo pathway. Among the ΔΔE_x values listed in the last three columns of Table 2, the largest difference arises from orbital interactions (ΔΔE_orb = −3.80 kcal mol⁻¹), followed by electrostatic interactions (ΔΔE_els = −2.44 kcal mol⁻¹), while the smallest difference is associated with dispersion interactions (ΔΔE_dis = −1.67 kcal mol⁻¹). This quantitative analysis for both pathways, performed at the CGP, is in full agreement with, and reinforces, the qualitative interpretations discussed earlier.

Given that orbital interactions are the primary factor underlying the superiority of the ortho-endo (on) pathway over the meta-endo (mn) pathway, it is natural to ask which specific orbital interaction plays the dominant role. This question can be thoroughly addressed using the Extended Transition State-Natural Orbitals for Chemical Valence (ETS-NOCV) analysis. ETS-NOCV offers a powerful and intuitive approach for dissecting the nature of orbital interactions by visually decomposing the deformation density associated with inter-fragment interactions, thus providing deep insight into the bonding mechanism.³² Despite the growing popularity of ETS-NOCV in computational studies, our literature survey indicates that most published works employing this method only present its results, without offering a clear, detailed, and accessible explanation of its fundamental concepts and theoretical background. Therefore, while interested readers are referred to ref. 32 and 40–42 to reach a clear insight into the concepts and formulation of the ETS-NOCV method, we have taken the opportunity here to provide a concise yet thorough and easy-to-follow description of the basic principles and theoretical foundations of this powerful technique for the benefit of readers seeking a deeper understanding.

Based on the ETS-NOCV approach, orbital interactions between fragments bring about a change in the density matrix. Eigenvectors of the variation of the density matrix are called NOCVs. Each NOCV has a corresponding eigenvalue, and there is another NOCV with an eigenvalue of opposite sign. According to the absolute values of the eigenvalues, the NOCVs can be paired. Usually only a few NOCV pairs have noticeable eigenvalues and therefore are worth investigating. Each NOCV pair has its corresponding contribution to the density change caused by orbital interactions, which can be plotted to visually examine the nature of the inter-fragment interactions represented by the NOCV pairs. Moreover, each NOCV pair has a specific contribution to orbital interaction energy between fragments, and the contribution value directly reflects the strength of the interaction represented by the NOCV pair. Within the ETS-NOCV approach, the natural orbitals for chemical valence (NOCV) diagonalize the deformation density matrix. The change in the electron density, caused by the density re-organization as a result of inter-fragment interactions, is given as:

where ν_k and M stand for NOCV eigenvalues and the number of occupied orbitals or basis functions, respectively. Hence, with the formation of the AB complex, a fraction ν_k of electrons is transferred from the ψ_−k to the ψ_k orbital. It should however be noted that such an electron density transformation results in a depletion of electron density in some regions of the molecular complex AB for which Δρ < 0 (characterized with a red-colored isosurface) and an accumulation of electron density in some other regions for which Δρ > 0 (distinguished with a blue-colored isosurface). Only some NOCV pairs have ν_k significantly different from zero, and this subgroup is generally enough to describe the A⋯B interaction. For each value of k, an energy contribution associated with the k-th NOCV pair, ΔE_oi, is given. (ΔE_oi)_i is always negative and, hence, specifies how much stabilization is gained during the combination of NOCVs ψ_−k and ψ_k. (ΔE_oi)_i is calculated as the sum of the eigenvalue weighted energies of NOCVs ψ_−k and ψ_k; that is:

(ΔE_oi)_i = (Z × ν_k) + [W × (−ν_k)] (10)

in which Z and W denote the energy of NOCVs k and −k, respectively. Evidently, a sum of all (ΔE_oi)_i values should recover the value of the total orbital energy in eqn (8); i.e., . Thus, ETS-NOCV not only enables a quantitative partitioning of the orbital interaction energy but also provides direct visual insight into the nature of bonding interactions in molecular complexes.

Utilizing ADF software, the ETS-NOCV analysis at the BLYP-D3(BJ)/TZ2P computational level was carried out on the molecular complexes located at ξ = 2.180 Å along the IRC profiles of TS-on and TS-mn, and the key results from this analysis are illustrated in Fig. 6. The ETS-NOCV analysis reveals that the total energy associated with all orbital interactions—arising from both charge transfer and polarization—between D-1 and Dph-2 amounts to −50.83 kcal mol⁻¹ and −47.52 kcal mol⁻¹ at the CGP along the ortho-endo and meta-endo pathways, respectively. Among the identified NOCV pairs, only the first NOCV pair exhibits a significantly larger stabilizing contribution compared to the others and thus merits detailed consideration for both reaction pathways. As shown in Fig. 6, NOCV pair 1 alone accounts for 69.07% [−35.11/−50.83 × 100] of the total orbital stabilization along the ortho-endo pathway and 60.06% [−28.54/−47.52 × 100] along the meta-endo pathway. Therefore, the contributions from the remaining NOCV pairs can safely be neglected in the present analysis. A closer examination of the deformation density isosurfaces associated with NOCV pair 1 reveals that this dominant interaction arises from the overlap of the frontier molecular orbitals of the two fragments.

Fig. 6 BLYP-D3(BJ)/TZ2P//M06-2X-D3/6-31G(d) dominant NOCV-based deformation density contributions determined for the 13-DCR between D-1 and Dph-2 at the CGP along the ortho-endo (left) and meta-endo (right) channels. An electron density isovalue of 0.05 e au⁻³ was used through 3-D representation of the NOCV pairs (see the text for more details). ψ₂₂ and ψ₂₈₇ participating in the Δρ of NOCV pair 1 correspond to the wavefunction of NOCV orbitals 22 (HOMO of D-1, ψ_−k) and 287 (LUMO of Dph-2, ψ_k) mainly localized over D-1 and Dph-2 fragments, respectively.

Specifically, along both regioisomeric pathways, NOCV pair 1 corresponds to the interaction between the HOMO of D-1 and the LUMO of Dph-2. This is evident from the presence of a nodal plane between the carbon and nitrogen terminal atoms in the D-1 fragment (depicted as the lower fragment in Fig. 6), which clearly indicates its HOMO character. Conversely, the Dph-2 fragment exhibits a nodal plane between the two carbon atoms of the C=C double bond, characteristic of a LUMO-type orbital. The deformation density isosurfaces associated with NOCV pair 1 provide clear visual evidence of the direction of electron density flow. Specifically, the red-colored isosurfaces localized on the D-1 fragment (corresponding to the regions of electron density depletion, Δρ < 0) and the blue-colored isosurfaces on the Dph-2 fragment (corresponding to the regions of electron density accumulation, Δρ > 0) unambiguously indicate a net electron density transfer from the HOMO of D-1 toward the LUMO of Dph-2. The extent of this charge transfer is quantified by the Δq value, amounting to 0.693e and 0.638e at the CGP along the ortho-endo and meta-endo pathways, respectively (see Fig. 6). The significant magnitude of these values reflects the strong electron-withdrawing character of the NO₂ substituent on the Dph-2 fragment, which effectively facilitates this electron transfer process. Taken together, these findings from ETS-NOCV analysis assert that the more pronounced HOMO_D-1 → LUMOD_ph-2 electron density transfer along the ortho-endo pathway, relative to the meta-endo pathway, is the primary factor underlying the stronger orbital stabilization and, consequently, the energetic preference for the ortho-endo reaction channel.

As clarified in the light of the energetics exploration (Table 1, Fig. 1, and section 3.2), the 13-DCR of 1,3-dipole D-1 toward dipolarophile Dph-2 in the presence of DCM at 80 °C as the key step along the domino process experimentally followed by Kumar and co-workers¹² not only represents entirely ortho regioselectivity (ortho regiospecificity) but also displays a high degree of endo stereoselectivity. Therefore, the last section of the present investigation is assigned to explore the origin of such a high endo stereoselectivity. To this end, activation/strain analysis was performed over the whole IRC profile of TS-on and TS-ox involved in the endo and exo stereoselective approach modes, respectively, in the energetically preferred ortho regioselective channel. Fig. 7 depicts the corresponding plots. As shown in Fig. 7, although the interaction energy curve associated with the ortho-exo pathway lies slightly below that of the ortho-endo pathway—indicating a somewhat more stabilizing interaction energy in the former—the total strain energy curve of the ortho-exo pathway is significantly more destabilizing compared to that of the ortho-endo pathway. Notably, the separation between the strain energy curves is more pronounced than that between the interaction energy curves, clearly demonstrating that the higher strain energy along the ortho-exo pathway outweighs its marginally stronger interaction energy. Consequently, it becomes evident that the high endo stereoselectivity observed along the fully preferred ortho regioselective pathway in the 13-DCR under study primarily originates from the lower strain energy associated with the endo approach, which outweighs the slightly stronger interaction energy of the exo pathway. Given that the total strain energy comprises the individual strain energies of the D-1 and Dph-2 fragments, further decomposition of this term provides a clearer picture of which fragment predominantly contributes to the greater total strain energy observed along the exo pathway compared to the endo pathway.

Fig. 7 The magnified M06-2X-D3/6-31G(d) activation/strain plot along the whole IRC profile of TS-on (solid curves) and TS-ox (dashed curves) involved in the ortho-endo (on) and ortho-exo (ox) channels of 13-DCR, respectively, between D-1 and Dph-2 in the presence of DCM. The blue circles characterize the positions of TS-on and TS-ox.

Accordingly, as shown in Fig. 8, the total strain energy, ΔE^str_tot(ξ), was further decomposed into the strain energy of the D-1 fragment, ΔE^str_D-1(ξ), and that of the Dph-2 fragment, ΔE^str_Dph-2(ξ), along both the ortho-endo (on) and ortho-exo (ox) pathways. Based on Fig. 8, it can be concluded that the strain energy experienced by the D-1 fragment is evidently higher than that of the Dph-2 fragment. Additionally, based on the separation distance between the strain energy curves of D-1, ΔE^str_D-1, and the strain energy curves of Dph-2, ΔE^str_Dph-2, the D-1 fragment undergoes a greater destabilizing strain deformation. Therefore, the greater total strain energy observed along the ortho-exo pathway primarily arises from the higher strain energy of the D-1 fragment.

Fig. 8 Magnified M06-2X-D3/6-31G(d) plots showing the decomposition of the total strain energy, ΔE^str_tot(ξ), into the strain energy contributions from the D-1 fragment, ΔE^str_D-1(ξ), and the Dph-2 fragment, ΔE^str_Dph-2(ξ), along both the ortho-endo (solid curves) and ortho-exo (dashed curves) pathways.

It is also worth noting that this greater strain energy in the D-1 fragment, relative to Dph-2, likely stems from its higher structural flexibility and greater number of internal degrees of freedom.

4. Conclusion

Very recently, Kumar and co-workers synthesized a series of pyrazole derivatives via a domino process involving a 1,3-dipolar cycloaddition reaction (13-DCR), followed by a 1,3-hydrogen shift and subsequent HOCl/NaNO₂ extrusion, carried out in dichloromethane (DCM) at 80 °C. Experimentally, only the ortho regioisomeric pyrazole was isolated, despite the potential formation of two regioisomers. A closer examination of this domino process reveals that the initial 13-DCR between the diazo compound (serving as the 1,3-dipole, D-1) and an electron-deficient nitroethylene derivative (acting as the dipolarophile, Dph-2) dictates the relative positioning of the CO₂Et and Ph substituents on the resulting pyrazole skeleton.

To elucidate the underlying factors governing the energetics, regioselectivity, and stereoselectivity of this key 13-DCR step, we performed a comprehensive DFT study at the rev-DSD-PBEP86-D4-CPCM(DCM)/ma-def2-TZVPP//M06-2X-D3-PCM(DCM)/6-31G(d) level of theory, consistent with the experimental conditions. The computed kinetics parameters clearly predict an overwhelming ortho regioselectivity, with an ortho-to-meta ratio of 357.8, perfectly aligning with the experimental observation of exclusive ortho regioisomer formation. Furthermore, within the favored ortho regioisomeric pathway, the computed endo-to-exo stereoselectivity ratio of 3.92 indicates a pronounced preference for the endo approach.

Activation strain model (ASM) analysis across the full reaction pathways demonstrated that this complete ortho regioselectivity (ortho regiospecificity) is primarily driven by more favorable stabilizing interaction energies. Further decomposition of the interaction energy using the sobEDA method confirmed that orbital interactions, ΔE_orb(ξ), are the dominant contributors to this regioselectivity.

To pinpoint the specific orbital interactions responsible, we carried out ETS-NOCV analysis. This revealed that the HOMO of D-1 interacting with the LUMO of Dph-2 constitutes the key orbital interaction, with substantially stronger and more stabilizing HOMO → LUMO charge transfer along the ortho pathway compared to the meta pathway, thereby decisively favoring ortho regioselectivity.

Lastly, to rationalize the high endo stereoselectivity observed within the preferred ortho pathway, we applied ASM analysis to the competing endo and exo approaches. This analysis clearly showed that the endo pathway benefits from a somewhat lower total strain energy requirement for fragment deformation compared to the exo pathway, thus explaining the high endo selectivity of this 13-DCR.

The mechanistic analyses suggest that these principles could be extended to other electron-deficient alkene dipolarophiles and may inform future synthetic strategies for pyrazole library construction.

Conflicts of interest

The author declares that he has no conflict of interest.

Data availability

The supporting data for this article have been included in the supplementary information (SI). Supplementary information: the Cartesian coordinates of M06-2X-D3-PCM(DCM)/6-31G(d) optimized geometries of all stationary points located over the PES of the 13-DCR. See DOI: https://doi.org/10.1039/d5ob01142f.

All computational software used is commercially available and has been used with default settings. The programs are appropriately cited in the Theoretical methods section.

References

1. G. Bacolini, Syst. Synth. React. Prop., 1996, 1, 103.
2. K. Ruffell, F. R. Smith, M. T. Green, S. M. Nicolle, M. Inman, W. Lewis, C. J. Hayes and C. J. Moody, Chem. – Eur. J., 2021, 27, 13703–13708.
3. D. Kaushik, S. A. Khan, G. Chawla and S. Kumar, Eur. J. Med. Chem., 2010, 45, 3943–3949.
4. C. Brullo, D. Caviglia, A. Spallarossa, S. Alfei, S. G. Franzblau, B. Tasso and A. M. Schito, Pharmaceutics, 2022, 14, 1770.
5. S. S. Parmar, B. Pandey, C. Dwivedi and R. D. Harbison, J. Pharm. Sci., 1974, 63, 1152–1155.
6. A. T. Taher, M. T. M. Sarg, N. R. E.-S. Ali and N. H. Elnagdi, Bioorg. Chem., 2019, 89, 103023.
7. T. Baiju and I. N. Namboothiri, Chem. Rec., 2017, 17, 939–955.
8. Y. Wu, J.-Y. Chen, J. Ning, X. Jiang, J. Deng, Y. Deng, R. Xu and W.-M. He, Green Chem., 2021, 23, 3950–3954.
9. A. Padwa and W. H. Pearson, Synthetic Applications of 1, 3-Dipolar Cycloaddition Chemistry Toward Heterocycles and Natural Products, John Wiley & Sons, 2003, vol. 59.
10. E. Buchner, Ber. Dtsch. Chem. Ges., 1889, 22, 842–847.
11. S. E. Denmark and A. Thorarensen, Chem. Rev., 1996, 96, 137–166.
12. G. Kumar, S. Ray, S. Saha and M. Shankar Singh, Asian J. Org. Chem., 2024, 13, e202400077.
13. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian, Inc., Wallingford, CT, 2019.
14. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
15. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
16. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
17. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
18. V. Barone, M. Cossi and J. Tomasi, J. Comput. Chem., 1998, 19, 404–417.
19. K. Fukui, J. Phys. Chem., 1970, 74, 4161–4163.
20. C. Gonzalez and H. B. Schlegel, J. Chem. Phys., 1991, 95, 5853–5860.
21. G. Santra, N. Sylvetsky and J. M. Martin, J. Phys. Chem. A, 2019, 123, 5129–5143.
22. E. Caldeweyher, J.-M. Mewes, S. Ehlert and S. Grimme, Phys. Chem. Chem. Phys., 2020, 22, 8499–8512.
23. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 73–78.
24. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2025, 15, e70019.
25. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
26. T. Lu and Q. Chen, Comput. Theor. Chem., 2021, 1200, 113249.
27. J. Ho, A. Klamt and M. L. Coote, J. Phys. Chem. A, 2010, 114, 13442–13444.
28. K. J. Laidler, Theories of chemical reaction rates, McGraw-Hill, New York, 1969.
29. F. M. Bickelhaupt and K. N. Houk, Angew. Chem., Int. Ed., 2017, 56, 10070–10086.
30. F. M. Bickelhaupt, J. Comput. Chem., 1999, 20, 114–128.
31. T. Lu and Q. Chen, J. Phys. Chem. A, 2023, 127, 7023–7035.
32. M. P. Mitoraj, A. Michalak and T. Ziegler, J. Chem. Theory Comput., 2009, 5, 962–975.
33. G. T. Te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. van Gisbergen, J. G. Snijders and T. Ziegler, J. Comput. Chem., 2001, 22, 931–967.
34. C. Fonseca Guerra, J. Snijders, G. te Velde and E. J. Baerends, Theor. Chem. Acc., 1998, 99, 391–403.
35. Amsterdam Density Functional program, Scientific Computing & Modelling NV, Vrije Universiteit, Amsterdam, https://www.scm.com.
36. P. Vermeeren, Nat. Protoc., 2020, 15, 649.
37. S. Emamian, M. Soleymani and S. S. Moosavi, New J. Chem., 2019, 43, 4765–4776.
38. T. Lu, ChemRxiv, 2022, preprint,  DOI:10.26434/chemrxiv-2022-r6rh8-v2.
39. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
40. K. Dyduch, M. P. Mitoraj and A. Michalak, J. Mol. Model., 2013, 19, 2747–2758.
41. M. P. Mitoraj, M. Parafiniuk, M. Srebro, M. Handzlik, A. Buczek and A. Michalak, J. Mol. Model., 2011, 17, 2337–2352.
42. M. Mitoraj and A. Michalak, J. Mol. Model., 2007, 13, 347–355.

---