Phosphorus-centered ion–molecule reactions: benchmark ab initio characterization of the potential energy surfaces of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] systems

Abstract

In the present work we determine the benchmark relative energies and geometries of all the relevant stationary points of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] identity and non-identity reactions using state-of-the-art electronic-structure methods. These phosphorus-centered ion–molecule reactions follow two main reaction routes: bimolecular nucleophilic substitution (S_N2), leading to Y⁻ + PH₂X, and proton transfer, resulting in HX + PHY⁻ products. The S_N2 route can proceed through Walden-inversion, front-side-attack retention, and double-/multiple-inversion pathways. In addition, we also identify the following product channels: H⁻-formation, PH₂⁻- and PH₂-formation, ¹PH- and ³PH-formation, H₂-formation and HY + PHX⁻ formation. The benchmark classical relative energies are obtained by taking into account the core-correlation, scalar relativistic, and post-(T) corrections, which turn out to be necessary to reach subchemical (<1 kcal mol⁻¹) accuracy of the results. Classical relative energies are augmented with zero-point-energy contributions to gain the benchmark adiabatic energies.

---

Introduction

Bimolecular nucleophilic substitution (S_N2) is one of the most fundamental reactions in organic chemistry. Especially carbon-centered S_N2 reactions have been widely studied, both by experimental and theoretical means.^1–5 Initial kinetic measurements showed that nucleophilic substitution, where a reactant ion or other species that has a partial negative charge substitutes a similar part of the reactant molecule (leaving group), can occur in two different ways: in the unimolecular S_N1 reaction the leaving group first detaches from the molecule, and the new bond with the attacking nucleophile forms afterwards; while in S_N2 the breaking and the forming bonds are simultaneously present at the transition state (TS) of the reaction, which is why the reaction rate depends on the concentrations of both the molecule and the nucleophile in the latter case.⁶ Stereochemistry is also a key aspect of S_N2 reactions, whose stereoselectivity at C center was first measured in 1896 by Paul Walden;⁷ however, the underlying Walden-inversion mechanism was only described at the atomic level in the 1930s by Ingold and co-workers.^6,8 In the early work a so-called front-side-attack pathway was also mentioned, which leads to products with retained configuration, however, via a high-energy barrier.⁶ At the end of the 20th century the investigation of gas-phase S_N2 reactions reached a new level as a result of the rapid progress of quantum chemistry: the thorough theoretical mapping of the potential energy surfaces (PESs) of such reactions became achievable. Prototypic X⁻ + CH₃Y [X, Y = F, Cl, Br, I] S_N2 reactions turned out to mainly proceed through a central pentavalent Walden-inversion TS with a collinear X–C–Y arrangement, with pre- and post-reaction ion-dipole minima,^1–5 and front-side TS geometries with high energies relative to the reactants also being determined.^9–11 A comprehensive high-level ab initio study revealed so-called double-inversion¹² TSs in S_N2@C reactions,¹¹ while hydrogen- and halogen-bonded minima were also found to be present.^2,3,13–17

Central-atom effects in S_N2 reactions have also been investigated: experiments showed that reactivity doubled upon changing the central atom from carbon to nitrogen,¹⁸ and a classical TS structure was proposed for S_N2@N.¹⁹ After lower-level theoretical studies,^20–25 the benchmark energies and geometries of the stationary points of the X⁻ + NH₂Y [X, Y = F, Cl, Br, I] reactions have also been determined.²⁶ Silicon- and phosphorus-centered S_N2 reactions have also been the subjects of a number of theoretical studies.^27–41 Amongst Si@S_N2 reactions the identity X⁻ + SiH₃X [X = F, Cl] systems were first looked at,^{27–32,35,36} and the extensive studies of Bickelhaupt and co-workers revealed that replacing the central C atom with the valence-isoelectronic Si atom transforms the Walden-inversion transition state to a minimum on the PES, due to stronger central-atom–ligand interaction and smaller steric congestion.^10,35,39,41 The early investigations mainly used density functional theory methods; however, very recently a thorough high-level ab initio characterization of numerous stationary points of the X⁻ + SiH₃Y [X, Y = F, Cl, Br, I] reaction has been carried out.³⁷ Despite the biochemical relevance of phosphorus,^42–44 literature on P-centered substitution reactions is very scarce. In the case of S_N2@P a Walden-minimum was also observed,^35,39–41 and in both the P and Si-centered cases front-side TSs were usually found to be submerged below the energy of the reactants due to the larger size of the central atom with respect to carbon.^{30–32,37,38} S_N2 reactions have also been investigated in micro-solvated or aqueous environments, which were shown to increase their barrier heights stabilizing the nucleophile more than the TS geometries, where the negative charge is smeared.^45–47

In recent decades not only the stationary-points of a PES could be characterized, but studying the atomic-level dynamics of S_N2 reactions has also become a feasible task, mostly carried out for S_N2@C.^1–3 Cross-molecular beam experiments made it possible to detect the outcome of single ion–molecule collisions,⁴⁸ which could be first explained by the so-called direct dynamics approach, pioneered by Hase and co-workers describing various S_N2 mechanisms: direct rebound, stripping, complex formation, and roundabout.^1,2,48 Later, our group introduced a new approach to model S_N2 reactions,³ providing statistically accurate dynamics results on full-dimensional analytical PESs, fitted on high-quality ab initio energy points, and contributing to our knowledge of S_N2 reactions by discovering the double-inversion mechanism,¹² entrance-channel front-side complex formation,¹⁶ an unexpected leaving-group effect,¹⁵ and competition with bimolecular elimination,⁴⁹ in joint experimental-theoretical investigations in the latter two cases. Compromised stereospecificity has been recently revealed in S_N2 dynamics upon changing the central carbon atom by nitrogen due to the facilitated inversion around the N central atom.⁵⁰ A detailed S_N2@Si dynamics study,³⁸ as well as one for the F⁻ + PH₂Cl reaction⁵¹ have also appeared lately from our group.

To uncover the characteristics of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] P-centered ion−molecule reactions, here we aim to determine the benchmark geometries and relative energies of many relevant stationary points using the best-available electronic-structure methods, also taking into account core-correlation, scalar relativistic, finite basis-set, and zero-point-energy effects. Based on this comprehensive analysis we will be able to make predictions related to their dynamics as well.

Computational details

As a first round, we search for stationary-point structures at the MP2/aug-cc-pVDZ^52,53 level of theory. After identifying several stationary points, we further optimize their geometries at the CCSD(T)-F12b/cc-pVDZ-F12 and the CCSD(T)-F12b/cc-pVTZ-F12 levels of theory,^54,55 where the F12 basis sets, despite the lack of the aug prefix, do contain diffuse functions as well.⁵⁵ At the CCSD(T)-F12b/cc-pVTZ-F12 geometries we perform CCSD(T)-F12b/cc-pVQZ-F12 single-point-energy computations to minimize basis-set-convergence error. Using the CCSD(T)-F12b/cc-pVTZ-F12 geometries further auxiliary energy contributions are calculated to the ΔE relative energies: (1) core-correlation correction: the effect of the correlation of the sub-valence-shell electrons, calculated using all-electron (AE) and frozen-core (FC) CCSD(T) methods and the aug-cc-pwCVTZ⁵⁶ basis set:

Δ_core = ΔE(AE-CCSD(T)/aug-cc-pwCVTZ) − ΔE(FC-CCSD(T)/aug-cc-pwCVTZ), (1)

(2) scalar relativistic effects determined by using the second-order Douglas–Kroll (DK) Hamiltonian⁵⁷ and the appropriate basis sets:^56,58

Δ_rel = ΔE(DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK) − ΔE(AE-CCSD(T)/aug-cc-pwCVTZ), (2)

and (3) post-(T) energy contributions:

δ[CCSDT] = ΔE(CCSDT/aug-cc-pVDZ) − ΔE(CCSD(T)/aug-cc-pVDZ) (3)

δ[CCSDT(Q)] = ΔE[(CCSDT(Q)/aug-cc-pVDZ) − ΔE(CCSDT/aug-cc-pVDZ)]. (4)

Thus, the benchmark classical energies for the stationary points are determined as:

ΔE_classical = CCSD(T)-F12b/cc-pVQZ-F12 + δ[CCSDT]+ δ[CCSDT(Q)] + Δ_core + Δ_rel (5)

To obtain vibrationally adiabatic energies, we perform harmonic frequency analysis at the CCSD(T)-F12b/cc-pVTZ-F12 level of theory, and calculate the zero-point vibrational energy (ZPE) contribution to the relative energy of each stationary point, then determine the benchmark adiabatic energies as

ΔE_adiabatic = ΔE_classical + Δ_ZPE. (6)

In the case of all systems containing Br and I, we use effective core potentials (ECPs) and the corresponding pseudo-potential (PP) basis sets⁵⁹ to replace the inner-core 10 (1s² 2s² 2p⁶) and 28 (1s² 2s² 2p⁶ 3s² 3p⁶ 3d¹⁰) electrons of the halogen atoms, respectively, except for the DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK term of Δ_rel. For open-shell products ROHF-based DK-AE-UCCSD(T), AE-UCCSD(T), and UCCSD(T)-F12b as well as UHF-based UCCSD(T), UCCSDT, and UCCSDT(Q) methods are applied. All the electronic structure computations are carried out with the MOLPRO program package,⁶⁰ except for the CCSDT⁶¹ and CCSDT(Q)⁶² results, which are obtained using the MRCC program.^63,64

Results and discussion

1. The X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions

First we map the potential energy surfaces of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions. The benchmark relative energies, both classical and adiabatic, along with the geometries of the stationary points are shown in Fig. 1.

Fig. 1 Benchmark geometries of the stationary points of the potential energy surfaces of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions along with their classical and adiabatic (in parentheses) energies relative to the energy of the reactants. The underlined value is obtained at the MP2/aug-cc-pVDZ level of theory.

At first glance we can notice that the identity reaction with the F atom shows some marked differences with respect to the three other reactions. Amongst the proton-transfer (PT) channels only the F reaction is exothermic, although just 2.5 kcal mol⁻¹ below the reactant asymptote. The transition complex corresponding to Walden inversion is also the deepest in the case of F (−43.5 kcal mol⁻¹), and the remaining stationary points (except for the self-inversion (SI) barrier of PH₂F, which is the highest) are way lower in relative energy than those corresponding to the other halogens.

The traditionally very high-energy front-side TS (FS TS) of the S_N2 channels of nitrogen and carbon-centered ion–molecule reactions here is only of 9.1 kcal mol⁻¹ relative energy in the case of the F substituent, while it is still well above the reactants for Cl, Br, and I. An interesting FS transition-state structure (FS TS′) is also revealed much lower than the traditional FS TS, even submerged by 21.9 kcal mol⁻¹ below the reactants in the case of X = F, which might play an important role in the dynamics of these reactions. The non-reactively oriented FS premin is not present for the F reaction, but it turns out to be submerged below the reactants for all halogens. Taking a further look at the entrance channel of the F⁻ + PH₂F reaction, we find a deep H-bonded pre-reaction minimum (however, much less deep than the Walden minimum), and also a quite low-lying pre-DI TS. These two stationary points have recently been revealed to be responsible for the so-called multi-inversion mechanism in the dynamics of the F⁻ + NH₂Cl reaction, leading to almost a racemic mixture of the S_N2 products.⁵⁰ The TS and post-reaction minimum of the PT channel are also submerged when X = F, thus PT is also barrierless, like all the Walden-inversion reaction routes, in this case. We also find another transition state, H2 TS, at a similar height as FS TS with a somewhat distorted structure featuring two H atoms very close to each other. The imaginary vibrations in this geometry, with corresponding frequencies above 1000i cm⁻¹, involve a mixture of P–H and H₂-stretching motions and therefore suggest that this TS is related to the H₂ formation product channel (see below).

The stationary-point energies for the X = Cl, Br, and I reactions always follow a trend: their relative energies increase as the size of the halogen increases, except for the FS premin and the self-inversion barriers. SI barriers are generally high, especially with respect to those of the NH₂X [X = F, Cl, Br, I] molecules.²⁶ Moreover, entrance-channel inversion is also not expected to be very relevant in the X = Cl, Br, I reactions, because although H premin is quite deep, but the corresponding inversion TS (pre-DI TS) is rather high-lying (with a gap of 34.7 kcal mol⁻¹ for X = I), making this induced-inversion barrier difficult to surmount. In addition, the Walden minima for these halogens are about 15 kcal mol⁻¹ deeper than the H premins, suggesting a preference for the single-inversion mechanism. The PT product channel becomes more and more endothermic and its TS and postmin more and more high-lying with respect to the H-bonded premin and the Walden minimum, making it less and less competitive with the S_N2 channel as we go downwards in the periodic table with X. The front-side-attack S_N2 retention pathway is expected to be the most competitive with the Walden-inversion for X = F, owing to the submerged FS TS′ and the low-lying traditional FS TS.

The relative energies of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions obtained at different levels of theory, along with the auxiliary energy corrections and ZPE contributions are listed in Table 1. Basis-set-convergence data are also presented in Fig. 2; and core-correlation, scalar relativistic and post-(T) corrections are shown in Fig. 3. As can be seen from Table 1 and Fig. 2, the basis-set convergence of the CCSD(T)-F12b method is excellent, because while the difference between the cc-pVDZ-F12 and cc-pVQZ-F12 (the latter considered as a reference) results is usually above 1 kcal mol⁻¹ (and always positive), the cc-pVTZ-F12 and cc-pVQZ-F12 relative energies differ only by a few 0.1 kcal mol⁻¹, with mixed signs (all positive for X = F). The basis-set errors are the smallest for the X = Cl reaction.

Table 1 Energies (in kcal mol⁻¹) of the stationary points of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions relative to the energy of the reactants at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with the auxiliary energy corrections (also in kcal mol⁻¹) at CCSD(T)-F12b/cc-pVTZ-F12 geometries, as described in the computational details section

	QZ^a	δ[T]^b	δ[(Q)]^c	Δ core	Δ rel	Δ ZPE	Classical^g	Adiabatic^h
a CCSD(T)-F12b/cc-pVQZ-F12. b CCSDT/aug-cc-pVDZ − CCSD(T)/aug-cc-pVDZ. c CCSDT(Q)/aug-cc-pVDZ − CCSDT/aug-cc-pVDZ. d AE-CCSD(T)/aug-cc-pwCVTZ − FC-CCSD(T)/aug-cc-pwCVTZ. e Scalar relativistic effects calculated as DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK − AE-CCSD(T)/aug-cc-pwCVTZ. f Zero-point vibrational energy contribution. g Benchmark classical energy, the sum of QZ, δ[CCSDT], δ[CCSDT(Q)], Δcore, and Δrel. h Benchmark adiabatic energy, the sum of benchmark classical energy and ΔZPE.
F^− + PH2F
HF + PHF^−	−2.59	0.04	0.00	0.10	−0.04	−1.01	−2.50	−3.51
PH2F SI	50.52	−0.06	−0.12	0.16	0.56	−1.00	51.06	50.06
H premin	−25.83	0.02	−0.04	−0.01	−0.01	0.51	−25.88	−25.37
PT postmin	−20.76	0.07	−0.02	0.02	−0.06	1.05	−20.75	−19.70
PT TS	−20.76	0.07	−0.02	0.02	−0.06	1.05	−20.75	−19.70
Walden min	−43.65	0.06	0.01	0.05	0.02	1.26	−43.51	−42.24
FS TS	8.15	−0.01	−0.31	0.41	0.83	0.92	9.07	9.98
FS TS'	−22.42	0.08	−0.03	0.24	0.23	0.60	−21.90	−21.29
H2 TS	9.02	0.01	−0.22	−0.09	−0.22	−2.85	8.50	5.65
pre-DI TS	−18.40	0.04	−0.03	0.03	0.06	0.62	−18.31	−17.68
Cl^− + PH2Cl
HCl + PHCl^−	24.91	−0.03	−0.07	0.01	−0.04	−2.37	24.77	22.40
PH2Cl SI	43.55	−0.02	−0.09	0.14	0.57	−0.92	44.16	43.24
H premin	−11.06	−0.02	−0.06	−0.09	0.02	0.09	−11.20	−11.12
PT postmin	12.43	−0.03	−0.13	−0.11	−0.04	−1.29	12.13	10.84
PT TS	12.88	−0.02	−0.12	−0.11	−0.03	−1.60	12.59	11.00
Walden min	−26.18	−0.01	−0.14	−0.07	−0.16	0.65	−26.55	−25.90
FS premin	−2.29	−0.01	−0.06	−0.03	−0.08	0.16	−2.48	−2.32
FS TS	36.75	−0.14	−0.50	0.31	0.65	0.32	37.07	37.38
FS TS'	8.14	−0.07	−0.24	0.11	−0.06	0.31	7.88	8.19
H2 TS	37.18	−0.04	−0.41	−0.25	−0.49	−2.87	35.99	33.12
pre-DI TS	12.65	−0.02	−0.14	−0.10	0.06	−1.48	12.46	10.97
Br^− + PH2Br
HBr + PHBr^−	32.55	−0.04	−0.08	−0.10	−0.21	−2.75	32.12	29.37
PH2Br SI	42.34	−0.01	−0.09	0.23	0.59	−0.90	43.06	42.16
PT premin	−9.80	−0.02	−0.05	−0.19	−0.17	0.13	−10.24	−10.11
H postmin	21.12	−0.03	−0.15	−0.40	−0.41	−1.84	20.13	18.29
PT TS	21.59	−0.02	−0.14	−0.37	−0.37	−2.19	20.69	18.50
Walden min	−24.14	−0.01	−0.14	−0.30	−0.38	0.50	−24.97	−24.47
FS premin	−5.56	−0.01	−0.11	−0.12	0.08	0.05	−5.72	−5.68
FS TS	40.16	−0.18	−0.57	0.09	0.50	0.06	40.00	40.06
H2 TS	39.85	−0.04	−0.42	−0.49	−0.60	−2.94	38.30	35.36
pre-DI TS	20.73	−0.01	−0.18	−0.32	−0.24	−2.13	19.99	17.86
I^− + PH2I
HI + PHI^−	37.82	−0.05	−0.09	−0.18	−0.05	−3.06	37.45	34.40
PH2I SI	39.23	0.00	−0.07	0.36	0.53	−0.85	40.05	39.20
H premin	−8.52	−0.01	−0.05	−0.35	0.00	0.11	−8.93	−8.82
PT postmin	28.15	−0.04	−0.19	−0.75	−0.08	−2.45	27.09	24.64
PT TS	28.59	0.76	−0.95	−0.68	−0.08	−2.67	27.64	24.97
Walden min	−21.03	0.00	−0.17	−0.61	−0.13	0.38	−21.93	−21.55
FS premin	−9.41	−0.01	−0.17	−0.29	0.14	−0.07	−9.73	−9.80
FS TS	42.42	−0.22	−0.66	−0.19	0.77	−0.05	42.11	42.05
H2 TS	42.69	−0.02	−0.47	−0.82	−0.38	−3.03	41.01	37.98
pre-DI TS	26.20	−0.02	−0.28	−0.39	0.28	−2.47	25.80	23.33

---

Fig. 2 Basis-set convergence of the CCSD(T)-F12b method in the case of the relative energies of the stationary points of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions. ΔVTZ = ΔE(CCSD(T)-F12b/cc-pVQZ-F12) – ΔE(CCSD(T)-F12b/cc-pVTZ-F12) and ΔVDZ = ΔE(CCSD(T)-F12b/cc-pVQZ-F12) – ΔE(CCSD(T)-F12b/cc-pVDZ-F12).

Fig. 3 Post-(T), scalar relativistic and core-correlation corrections to the relative energies of the stationary points of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions. δ[CCSDT] = ΔE(CCSDT/aug-cc-pVDZ) – ΔE(CCSD(T)/aug-cc-pVDZ), δ[CCSDT(Q)] = ΔE(CCSDT(Q)/aug-cc-pVDZ) – ΔE(CCSDT/aug-cc-pVDZ), Δ_core = ΔE(AE-CCSD(T)/aug-cc-pwCVTZ) – ΔE(FC-CCSD(T)/aug-cc-pwCVTZ), and Δ_rel = ΔE(DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK) – ΔE(AE-CCSD(T)/aug-cc-pwCVTZ).

To take the effect of higher-order excitations into account, we calculate both the CCSDT and CCSDT(Q) contributions to the relative energies of the stationary points of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions, shown in Table 1 and Fig. 3. The difference in the CCSD(T) and CCSDT values is only 0.1–0.2 kcal mol⁻¹, sometimes even negligible, and using the CCSDT method mostly reduces the relative energies, except for the F⁻ + PH₂F reaction. In contrast, the (Q)–T corrections are larger, especially for the FS TS′ structures, and always negative, except for the Walden minimum when X = F. Both the δ[CCSDT] and δ[CCSDT(Q)] corrections have outstanding values, >0.5 kcal mol⁻¹, although of opposite signs, for the PT TS of the I⁻ + PH₂I reaction.

As to the effect of correlating the sub-valence-shell electrons, shown also in Fig. 3, one can observe mixed-sign corrections of a few 0.1 kcal mol⁻¹ for X = F and Cl, which become larger negative values for the X = Br and I reactions, suggesting a stronger contribution from the outer core electrons. Relativistic effects of the inner core of the Br and I atoms are partly considered by using ECPs, however, further relativistic energy corrections can be obtained by using the Douglas–Kroll approach. These latter contributions are also plotted in Fig. 3, and show quite large values (sometimes several 0.1 kcal mol⁻¹), both negative and positive, which are not usual for second-row central atoms in S_N2 reactions, therefore might emerge due to the presence of phosphorus.

The adiabatic ZPE contributions listed in Table 1 are more and more negative and significant values as we go down in the periodic table, especially for the PT channel, where these are the results of two effects: (1) the number of the vibrational degrees of freedom (DOFs) decreases from 6 to 4 from reactants to products, causing a decrease in the ZPE of the products, and (2) while the forming H–F bond is stronger, and therefore has higher frequency, than the breaking P–H bond, the H–X bond becomes weaker and weaker as the size of the halogen increases, resulting in smaller and smaller ZPE for the products. The ZPEs of the Walden and H-bonded entrance-channel minima are always larger than the ZPE of the corresponding reactant, and the FS premin also features larger ZPE contributions in the case of X = F, Cl, Br with respect to the reactants, however with smaller and smaller values when going from F to I(Br). The reason behind this, apart from the increasing number of DOFs, may be the lengthened P–X bonds in the minimum structures (see Fig. 4) due to the increase in the size of the halogen, and thus the lower frequencies in the complexes. A special case is the FS premin structure, where a halogen-halogen bond is present, which, however, follows the same tendency. The pre-DI TS geometries, where the H–X bonds lengthen with increasing halogen-atom size (Fig. 4), also show decreasing ZPE contribution relative to the reactants. The largest negative values are seen for H2 TS, where the ZPE of the TS does not contain the contribution of the imaginary frequency, which is >1000i cm⁻¹ in all four cases.

Fig. 4 Selected structural parameters of the stationary points of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions obtained at the CCSD(T)-F12b/cc-pVTZ-F12 geometries. Bond lengths are in Å, angles are in degrees. The underlined values are obtained at the MP2/aug-cc-pVDZ level of theory.

In Fig. 4 we plot some important structural parameters for the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions. Besides the obvious growth of the P–X and H–X bond lengths with increasing X size, and the practically constant values of the P–H bonds in the PH₂X and PHX⁻ geometries, we can notice some clearly reactant- (H premin and FS premin) and product-like (PT TS and PT postmin) geometries. In FS premin the P–X bonds are only slightly longer than in PH₂X, while the halogens are rather far from each other. The situation is similar for H premin, that is, for the P–X bonds and long H–X distances. PT postmin features short H–X bond lengths, increasingly longer as the X size grows, than those of the HX molecules, but shorter than in the PT TS geometry. The P–X distances in PT postmin are also increasingly longer than in PT TS with F → I (except for Cl). These observations for the PT channel are in accordance with its increasingly endothermic nature as going down in the periodic table, based on Hammond's postulate.⁶⁵ The X–P–X angles always increase with increasing halogen-atom size in both the FS and H2 TSs, in contrast to the X–P–H angles, which decrease as they involve larger halogen atoms, e.g., in the PHX⁻, PT TS, postmin, and FS TS′ structures. The P–H–(attacking)X angle is almost linear (only slightly less than 180°), and decreases from X = Cl to I, while this angle is slightly greater than 180° for X = F.

We show the relative energies of some alternative product channels of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions in Table 2. All of these product channels are endothermic, except H₂ + PF₂⁻ (−8.7 kcal mol⁻¹ reaction energy): the highest-energy products are the X₂ + PH₂⁻ for X = F and Cl, which correspond to excited electronic states as the two doublet species (X₂⁻ + PH₂) have deeper energies, and the H⁻ + PH₂X for X = Br and I, whereas the lowest-energy channel is always H₂ formation, followed by the one leading to triplet PH (except for X = I). The energy difference between the excited singlet-state PH and the triplet one is around 26.5 kcal mol⁻¹, but on a singlet PES the only way to produce a PH molecule is through the formation of ¹PH. The reaction energies of the H₂ + PX₂⁻, H⁻ + PH₂X and the XHX⁻ + PH (singlet or triplet) channels are larger and larger as we go from F to I. In contrast, the remaining two product channels, which differ only in the position of the negative charge, become less and less endothermic as the size of the halogen increases. The CCSDT corrections are the largest, around −3.8 kcal mol⁻¹ for the ¹PH formation channel, which is very unusual and may require further investigations, and the second largest for the X₂⁻ + PH₂ reactions, however, their absolute values are below 1 kcal mol⁻¹, except for X = F. The CCSDT(Q) contributions are usually negative values of a few 0.1 kcal mol⁻¹, except the largest correction, which corresponds to the F₂⁻ + PH₂ product channel. The core-correlation and scalar relativistic effects are also in the magnitude of the δ[(Q)] contributions, mostly with positive and negative signs, respectively. The ZPE-contributions are always negative values, between −2.5 and −5.2 kcal mol⁻¹, as in many of these alternative reaction routes multiple bonds are affected, along with the formation of some rather exotic species.

Table 2 Reaction energies (in kcal mol⁻¹) of alternative product channels of the X⁻ + PH₂X [X = F, Cl, Br, I] identity reactions relative to the energy of the reactants at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with the auxiliary energy corrections (also in kcal mol⁻¹) at CCSD(T)-F12b/cc-pVTZ-F12 geometries, as described in the computational details section

	QZ^a	δ[T]^b	δ[(Q)]^c	Δ core	Δ rel	Δ ZPE	Classical^g	Adiabatic^h
a CCSD(T)-F12b/cc-pVQZ-F12. b CCSDT/aug-cc-pVDZ − CCSD(T)/aug-cc-pVDZ. c CCSDT(Q)/aug-cc-pVDZ − CCSDT/aug-cc-pVDZ. d AE-CCSD(T)/aug-cc-pwCVTZ − FC-CCSD(T)/aug-cc-pwCVTZ. e Scalar relativistic effects calculated as DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK − AE-CCSD(T)/aug-cc-pwCVTZ. f Zero-point vibrational energy contribution. g Benchmark classical energy, the sum of QZ, δ[CCSDT], δ[CCSDT(Q)], Δcore, and Δrel. h Benchmark adiabatic energy, the sum of benchmark classical energy and ΔZPE.
F^− + PH2F
H^− + PHF2	25.26	0.22	0.05	−0.04	−0.05	−3.18	25.44	22.25
F2 + PH2^−	126.58	0.00	−0.73	0.40	−0.36	−2.92	125.90	122.97
F2^− + PH2	86.37	−1.53	−1.59	0.17	−0.36	−3.25	83.07	79.81
FHF^− + ^1PH	36.11	−3.72	−0.32	0.86	−0.28	−2.51	32.64	30.13
FHF^− + ^3PH	9.56	0.03	0.08	0.04	−0.28	−2.54	9.43	6.89
PF2^− + H2	−8.78	0.17	−0.08	0.11	−0.12	−2.49	−8.70	−11.20
Cl^− + PH2Cl
H^− + PHCl2	74.01	0.11	−0.10	−0.20	−0.27	−3.98	73.56	69.58
Cl2 + PH2^−	77.35	−0.09	−0.15	0.09	−0.23	−2.79	76.97	74.18
Cl2^− + PH2	50.78	−0.94	−0.62	0.12	−0.34	−2.86	49.00	46.14
ClHCl^− + ^1PH	57.54	−3.79	−0.34	0.82	−0.22	−4.92	54.01	49.09
ClHCl^− + ^3PH	30.99	−0.04	0.06	0.00	−0.21	−4.95	30.80	25.86
PCl2^− + H2	15.66	0.10	−0.22	−0.01	−0.25	−3.45	15.28	11.83
Br^− + PH2Br
H^− + PHBr2	83.63	0.11	−0.13	−0.52	−0.21	−4.27	82.89	78.62
Br2 + PH2^−	70.00	−0.07	−0.15	0.30	0.20	−2.83	70.28	67.45
Br2^− + PH2	39.68	−0.82	−0.45	0.25	0.06	−2.70	38.73	36.02
BrHBr^− + ^1PH	63.16	−3.79	−0.34	0.57	−0.54	−4.90	59.04	54.14
BrHBr^− + ^3PH	36.62	−0.04	0.05	−0.25	−0.54	−4.93	35.83	30.91
PBr2^− + H2	19.29	0.10	−0.24	−0.13	−0.28	−3.77	18.74	14.97
I^− + PH2I
H^− + PHI2	91.96	0.11	−0.17	−0.78	−0.05	−4.44	91.07	86.63
I2 + PH2^−	59.90	−0.06	−0.11	0.39	0.27	−2.73	60.39	57.66
I2^− + PH2	30.24	−0.68	−0.33	0.29	0.08	−2.53	29.59	27.06
IHI^− + ^1PH	70.18	−3.81	−0.35	0.24	−0.26	−5.16	66.00	60.84
IHI^− + ^3PH	43.64	−0.06	0.05	−0.58	−0.25	−5.19	42.79	37.60
PI2^− + H2	21.04	0.13	−0.29	−0.31	−0.21	−3.93	20.36	16.44

---

As is clear from the above observations, we must take into account all the above auxiliary energy corrections, if we want to reach subchemical, that is, <1 kcal mol⁻¹, accuracy for our benchmark relative energies. Moreover, based on such systematic analysis of the different energy contributions, we can assign an ‘error bar’ to our calculated results, which in the case of the title reactions is expected to be under ±0.5 kcal mol⁻¹.

2. The X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions

In Fig. 5 we plot the schematic energy diagram of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions, where we always consider the exothermic direction of the S_N2 reaction for an X-Y combination. As expected, similar stationary points are found here as for the identity reactions, only two additional structures are identified: FS min and FS postmin. The first thing we can notice is that there are always large energy gaps separating the X = F systems from the other half of the non-identity reactions. Both the PT and the S_N2 channels are exothermic when X = F, and substantial relative-energy differences can be observed for all the stationary points with respect to the X ≠ F cases. The most exothermic S_N2 reaction, −53.0 kcal mol⁻¹, corresponds to the X = F, Y = I pair, and the reaction energy increases, i.e., the exothermicity decreases, when going from Y = I to Y = Cl. The X ≠ F S_N2 reaction energies are also negative, however, the PT channel is endothermic for these reactions, and is 28.1 kcal mol⁻¹ in the highest-energy X = Br, Y = I case. The deep Walden-inversion complex, which is the global minimum for all three reactive systems, features relative energies between −27.0 and −66.9 kcal mol⁻¹, being around twice as deep for the X = F reactions as for the others, following the X, Y = F, I; F, Br; F, Cl; and about 25 kcal mol⁻¹ higher the X, Y = Cl, I; Cl, Br; Br, I increasing order. This energy order applies for all the other non-identity stationary points as well, except for H2 TS. Thus, the lowest energy stationary points can be found also for the F⁻ + PH₂I reaction. The largest differences in relative energy are seen for the H2 TS between the X = F and the X ≠ F reactions, which for X = F, Y = Br lies 47.3 kcal mol⁻¹ below, and for X = Br, Y = I 31.8 kcal mol⁻¹ above the reactant asymptote.

Fig. 5 Benchmark geometries of the stationary points of the potential energy surfaces of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions along with their classical and adiabatic (in parentheses) energies relative to the energy of the reactants. We consider the exothermic S_N2 reaction routes for a given X, Y pair. The underlined values are obtained at the MP2/aug-cc-pVDZ level of theory.

The H-bonded entrance-channel minima are also of significant depth, closer in energy to Walden min when X = F, with a low-lying (7 kcal mol⁻¹ higher) inversion TS (pre-DI TS), found only for the X = F, Y = Cl, Br reactions, which may promote pre-reaction inversion leading to undermined stereoselectivity, as seen for the F⁻ + NH₂Cl reaction.⁵⁰ Since the PT channel is exothermic when X = F and H premin has a quite PT-product-like geometry (see below), the pre-reaction inversion might take place during a proton-abstraction (induced inversion), which can then either be followed by S_N2 or PT reactions. The PT channel lacks both kinetic and thermodynamic barriers for the X = F reactions, and the former holds for the remaining three as well. The PT TS and postmin lie close to each other in the non-identity case as well, with around −30 kcal mol⁻¹ relative energies when X = F, while setting a substantial 17 kcal mol⁻¹ barrier for the X, Y = Br, I pair. These larger X ≠ F barrier heights can be attributed to the weaker Y···HX hydrogen bonds with respect to the Y···HF connection.

The front-side attack substitution route leading to a retained configuration around the central P atom has two kinds of TS structures here as well, as already seen for the identity reactions. However, here, while the X ≠ F geometries of FS TS′ are similarly pyramidal as those of the X = Y systems, this transition state features different shapes for the X = F case, that is, one of the H atoms is turned and placed between the two halogen atoms, made possible by less steric congestion. The pyramidal FS TS′ structures feature imaginary frequencies around 300–400i cm⁻¹ and their motion along the imaginary mode is the twist of the two hydrogens, whereas the non-pyramidal X ≠ Y ≠ F FS TS′ geometries have imaginary frequencies of only 100–200i cm⁻¹ with a corresponding motion of PH₂F bending. FS TS′ lies always much lower than the classic FS TS, with about 20 kcal mol⁻¹ for the X = F reactions, and about 30 kcal mol⁻¹ when X ≠ F, indicating a lower-barrier, or even barrierless in the X = F case, front-side-, or more appositely side-attack, retention mechanism for the P-centered S_N2 reactions, with respect to the C- and N-centered ones. FS postmin is a rather deep minimum found in the exit-channel, with a product-like (see below) geometry but similar to that of FS TS′, only a few kcal mol⁻¹ deeper than the S_N2 product asymptote, therefore it guides the post-reaction system effectively to form the products. The FS premin with non-reactive orientation is submerged below the reactants in all cases, giving rise to notable hindering effects on the reactivity. It has the deepest relative energy (−29.2 kcal mol⁻¹) in the P–I–F arrangement and the highest (−7.0 kcal mol⁻¹) for P–Br–Cl. In the former case FS premin is only 10 kcal mol⁻¹, in the latter it is only 5 kcal mol⁻¹ above H premin; while it has alternating depth with respect to FS min.

In Tables 3 and 4 we show the relative energies, determined at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with auxiliary energy corrections (see the computational details section), of the stationary points of the X⁻ + PH₂Y [X = F and Y = Cl, Br, I] and X⁻ + PH₂Y [X = Cl, Br and Y = Br, I] non-identity reactions. These corrections are also visualized in Fig. 6 (basis-set convergence) and 7 (core-correlation, scalar relativistic, and post-(T)). In Fig. 6 and in Tables 3 and 4 we can see that in the X = F reactions the deviations between the cc-pVDZ-F12 and cc-pVQZ-F12 basis sets are quite large and always positive, between 1.0–1.5 kcal mol⁻¹, larger than in the case of the identity reactions. The TZ–QZ differences are also somewhat larger than in the identity case, also with positive signs, however, they fall into an order-of-magnitude smaller range, only a few 0.1 kcal mol⁻¹, similarly to the identity reactions. These basis-set corrections are smaller and mostly of positive signs for the X⁻ + PH₂Y [X = Cl, Br and Y = Br, I] reactions. ΔVDZ is mostly below 1 kcal mol⁻¹, whereas ΔVTZ is only 0.1–0.2 kcal mol⁻¹. Some geometries could only be optimized at the MP2/aug-cc-pVDZ level of theory, thus no auxiliary energy corrections are calculated. Taken together, the excellent basis-set convergence properties of the CCSD(T)-F12b method are clearly reflected by these results, and we can also claim that using the cc-pVQZ-F12 basis set is necessary for achieving sub-chemical accuracy.

Table 3 Energies of the stationary points (in kcal mol⁻¹), determined at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with the auxiliary energy corrections (also in kcal mol⁻¹) at CCSD(T)-F12b/cc-pVTZ-F12 geometries, as described in the computational details section, of the X⁻ + PH₂Y [X = F and Y = Cl, Br, I] non-identity reactions, relative to the energy of the reactants

	QZ^a	δ[T]^b	δ[(Q)]^c	Δ core	Δ rel	Δ ZPE	Classical^g	Adiabatic^h
a CCSD(T)-F12b/cc-pVQZ-F12. b CCSDT/aug-cc-pVDZ − CCSD(T)/aug-cc-pVDZ. c CCSDT(Q)/aug-cc-pVDZ − CCSDT/aug-cc-pVDZ. d AE-CCSD(T)/aug-cc-pwCVTZ − FC-CCSD(T)/aug-cc-pwCVTZ. e Scalar relativistic effects calculated as DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK − AE-CCSD(T)/aug-cc-pwCVTZ. f Zero-point vibrational energy contribution. g Benchmark classical energy, the sum of QZ, δ[CCSDT], δ[CCSDT(Q)], Δcore, and Δrel. h Benchmark adiabatic energy, the sum of benchmark classical energy and ΔZPE.
F^− + PH2Cl
HF + PHCl^−	−14.55	0.07	0.00	0.13	0.00	−0.74	−14.35	−15.09
Cl^− + PH2F	−38.91	0.09	0.15	0.16	0.16	0.67	−38.35	−37.68
H premin	−35.76	0.06	−0.04	0.01	0.03	0.93	−35.71	−34.77
PT postmin	−29.27	0.07	−0.02	0.03	0.00	0.89	−29.20	−28.30
PT TS	−28.57	0.07	−0.03	0.02	0.00	0.56	−28.51	−27.95
Walden min	−58.51	0.07	0.05	0.07	0.06	1.35	−58.26	−56.92
FS premin	−7.24	−0.03	−0.12	0.01	−0.16	0.28	−7.54	−7.26
FS TS	2.19	−0.03	−0.32	0.45	0.84	0.93	3.13	4.05
FS TS'	−18.66	−0.02	−0.13	0.01	−0.19	0.31	−18.99	−18.68
H2 TS	−40.31	0.06	0.11	0.08	0.07	0.64	−39.98	−39.34
FS min	−26.72	0.07	−0.03	0.17	0.12	0.87	−26.39	−25.52
FS postmin	−46.51	0.07	0.10	0.11	0.12	1.19	−46.11	−44.92
pre-DI TS	−28.85	0.08	−0.03	0.05	0.10	0.83	−28.65	−27.83
F^− + PH2Br
HF + PHBr^−	−17.93	0.07	0.00	0.23	−0.05	−0.61	−17.69	−18.29
Br^− + PH2F	−47.13	0.10	0.17	0.41	0.18	0.96	−46.27	−45.31
H premin	−38.43	0.06	−0.03	0.09	−0.04	1.08	−38.36	−37.28
PT postmin	−31.37	0.07	−0.02	0.05	−0.12	0.91	−31.40	−30.49
PT TS	−30.89	0.07	−0.02	0.05	−0.08	0.64	−30.88	−30.24
Walden min	−63.45	0.08	0.08	0.15	−0.04	1.48	−63.18	−61.70
FS premin	−16.48	0.02	−0.24	0.18	0.17	0.21	−16.34	−16.12
FS TS	−0.46	−0.04	−0.33	0.48	0.72	0.93	0.38	1.31
FS TS'	−20.98	0.00	−0.14	0.05	−0.25	0.34	−21.33	−20.98
H2 TS	−47.87	0.08	0.14	0.28	0.05	0.85	−47.32	−46.47
FS postmin	−53.09	0.08	0.13	0.26	0.08	1.27	−52.54	−51.27
pre-DI TS	−31.61	0.08	−0.02	0.14	0.04	0.85	−31.37	−30.52
F^− + PH2I
HF + PHI^−	−21.96	0.08	0.00	0.29	−0.10	−0.46	−21.69	−22.15
I^− + PH2F	−53.92	0.12	0.19	0.61	0.01	1.21	−52.98	−51.78
H premin	−41.61	0.08	−0.04	0.11	−0.07	1.17	−41.52	−40.35
PT postmin	−33.91	0.08	−0.02	0.01	−0.10	0.95	−33.94	−32.99
PT TS	−33.38	0.08	−0.02	0.01	−0.09	0.70	−33.40	−32.69
Walden min	−67.23	0.10	0.11	0.18	−0.03	1.61	−66.88	−65.27
FS premin	−29.78	0.08	−0.22	0.55	0.14	0.21	−29.24	−29.03
FS TS	−2.82	−0.02	−0.35	0.46	0.75	0.96	−1.98	−1.02
FS TS'	−23.47	0.02	−0.15	0.05	−0.24	0.40	−23.78	−23.38
H2 TS	−7.33	−0.02	−0.11	−0.03	−0.36	−1.62	−7.83	−9.46
FS postmin	−56.95	0.10	0.16	0.36	−0.02	1.32	−56.35	−55.04

---

Table 4 Energies of the stationary points (in kcal mol⁻¹), determined at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with the auxiliary energy corrections (also in kcal mol⁻¹) at CCSD(T)-F12b/cc-pVTZ-F12 geometries, as described in the computational details section, of the Cl⁻ + PH₂Br, Cl⁻ + PH₂I, and Br⁻ + PH₂I reactions, relative to the energy of the reactants

	QZ^a	δ[T]^b	δ[(Q)]^c	Δ core	Δ rel	Δ ZPE	Classical^g	Adiabatic^h
a CCSD(T)-F12b/cc-pVQZ-F12. b CCSDT/aug-cc-pVDZ − CCSD(T)/aug-cc-pVDZ. c CCSDT(Q)/aug-cc-pVDZ − CCSDT/aug-cc-pVDZ. d AE-CCSD(T)/aug-cc-pwCVTZ − FC-CCSD(T)/aug-cc-pwCVTZ. e Scalar relativistic effects calculated as DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK − AE-CCSD(T)/aug-cc-pwCVTZ. f Zero-point vibrational energy contribution. g Benchmark classical energy, the sum of QZ, δ[CCSDT], δ[CCSDT(Q)], Δcore, and Δrel. h Benchmark adiabatic energy, the sum of benchmark classical energy and ΔZPE.
Cl^− + PH2Br
HCl + PHBr^−	21.53	−0.03	−0.07	0.12	−0.10	−2.25	21.44	19.19
Br^− + PH2Cl	−8.22	0.01	0.02	0.24	0.02	0.29	−7.93	−7.63
H premin	−11.66	−0.02	−0.06	−0.08	−0.03	0.06	−11.86	−11.80
PT postmin	10.11	−0.03	−0.12	−0.09	−0.16	−1.24	9.71	8.47
PT TS	10.37	−0.02	−0.12	−0.08	−0.13	−1.51	10.02	8.51
Walden min	−29.56	0.00	−0.12	−0.05	−0.26	0.74	−29.99	−29.25
FS premin	−6.96	−0.01	−0.13	−0.02	0.13	0.11	−6.98	−6.88
FS TS	34.33	−0.15	−0.52	0.32	0.56	0.33	34.54	34.87
FS TS'	5.40	−0.08	−0.25	0.09	−0.07	0.37	5.09	5.47
H2 TS	29.20	−0.04	−0.31	−0.21	−0.54	−2.33	28.09	25.76
FS min	−9.31	−0.01	−0.06	−0.08	−0.11	0.40	−9.56	−9.16
Cl^− + PH2I
HCl + PHI^−	17.50	−0.02	−0.07	0.18	−0.15	−2.10	17.44	15.34
I^− + PH2Cl	−15.01	0.03	0.04	0.45	−0.15	0.54	−14.64	−14.10
PT postmin	7.43	−0.02	−0.12	−0.15	−0.13	−1.17	7.01	5.84
PT TS	7.68	−0.01	−0.12	−0.14	−0.12	−1.39	7.29	5.91
Walden min	−32.07	0.01	−0.10	−0.06	−0.23	0.85	−32.45	−31.61
FS premin	−14.08	0.00	−0.19	0.09	0.08	0.03	−14.10	−14.07
FS TS	32.17	−0.15	−0.55	0.28	0.59	0.35	32.35	32.69
FS TS'	3.78	−0.09	−0.28	0.03	−0.08	0.43	3.36	3.79
H2 TS	26.75	−0.04	−0.29	−0.24	−0.52	−2.14	25.66	23.52
FS min	−10.26	0.00	−0.06	−0.11	−0.11	0.43	−10.53	−10.10
Br^− + PH2I
HBr + PHI^−	28.52	−0.02	−0.09	−0.04	−0.25	−2.60	28.12	25.52
I^− + PH2Br	−6.79	0.02	0.02	0.21	−0.17	0.24	−6.71	−6.47
PT postmin	18.39	−0.02	−0.15	−0.47	−0.37	−1.82	17.38	15.57
PT TS	18.81	−0.01	−0.14	−0.43	−0.35	−2.08	17.88	15.80
Walden min	−26.15	0.00	−0.14	−0.34	−0.34	0.57	−26.97	−26.40
FS premin	−11.76	0.00	−0.18	−0.08	0.02	−0.04	−12.01	−12.04
FS TS	37.83	−0.18	−0.60	0.06	0.54	0.11	37.65	37.76
FS TS'	8.93	−0.13	−0.31	−0.22	−0.07	0.24	8.20	8.44
H2 TS	33.24	−0.03	−0.32	−0.48	−0.56	−2.39	31.85	29.46

---

Fig. 6 Basis-set convergence of the CCSD(T)-F12b method in the case of the relative energies of the stationary points of the X⁻ + PH₂Y [X,Y = F, Cl, Br, I] non-identity reactions. ΔVTZ = ΔE(CCSD(T)-F12b/cc-pVQZ-F12) − ΔE(CCSD(T)-F12b/cc-pVTZ-F12) and ΔVDZ = ΔE(CCSD(T)-F12b/cc-pVQZ-F12) − ΔE(CCSD(T)-F12b/cc-pVDZ-F12).

In Fig. 7 the core-correlation, scalar relativistic and post-(T) corrections to the relative energies of the stationary points are plotted. The core-correlation contributions are of few to several 0.1 kcal mol⁻¹ in value with almost always positive signs when X = F, and become larger and more and more negative as the size of the halogen atoms increases. These corrections are usually larger positive values for the FS TS structures, while it is the largest for the proton-transfer product channel of the F⁻ + PH₂I reaction, indicating that the difference of the energy contribution of the core electrons here is the largest between the reactants and the products. The scalar relativistic effects are also outstanding in the case of FS TS resulting in large positive energy corrections. Interestingly the relativistic contribution is negligible for the I⁻ + PH₂F product channel. Otherwise, the relativistic corrections are of a few 0.1 kcal mol⁻¹ and of mixed signs, and are quite similar in the X = F stationary geometries, but become more and more negative and considerable in magnitude with increasing halogen-atom size. Among all three kinds of energy contributions highlighted in Fig. 7, the scalar relativistic effects are usually the most significant, in some cases exceeding 0.5 kcal mol⁻¹. The impact of higher-order excitations is reflected by the post-(T) corrections: the δ[CCSDT] values are small, only 0.1–0.2 kcal mol⁻¹, usually positive for the X = F reactions, and negative for the others; while the δ[CCSDT(Q)] contributions are a bit larger, a few 0.1 kcal mol⁻¹ with mixed signs when X = F, and mostly negative otherwise, reaching even −0.6 kcal mol⁻¹ in the case of the FS TS geometry of the Br⁻ + PH₂I reaction. The δ[(Q)] energy corrections are also generally large in the FS TS structures.

Fig. 7 Post-(T), scalar relativistic and core-correlation corrections to the relative energies of the stationary points of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions. δ[CCSDT] = ΔE(CCSDT/aug-cc-pVDZ) − ΔE(CCSD(T)/aug-cc-pVDZ), δ[CCSDT(Q)] = ΔE(CCSDT(Q)/aug-cc-pVDZ) − ΔE(CCSDT/aug-cc-pVDZ), Δ_core = ΔE(AE-CCSD(T)/aug-cc-pwCVTZ) − ΔE(FC-CCSD(T)/aug-cc-pwCVTZ), and Δ_rel = ΔE(DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK) − ΔE(AE-CCSD(T)/aug-cc-pwCVTZ).

The absolute value of the adiabatic ZPE corrections, listed in Tables 3 and 4, are usually between 0–1.5 kcal mol⁻¹, mostly with positive signs (more negative contributions are seen for the X ≠ F reactions). Taking into account ZPE always reduces the reaction energies of the PT channel, due to the same reasons as in the case of the identity reactions: the decreasing number of vibrational DOFs wins over the forming H–halogen bonds, featuring higher frequencies, than the breaking P–H bonds. The PT postmin shows larger positive ZPE corrections if X = F, which turn into larger negative values for the other three reactions. In the former case, the effect of the larger number of DOFs (from 6 to 9 from reactants to postmin) is supplemented by the stronger forming H–F bond, very close in length to the one in the HF molecule, with respect to the breaking P–H bond. In the X ≠ F case, the PT postmin is much less product-like, thus the forming H–halogen bonds are not that strong, and thus contribute less to the ZPE of the complex. Besides the PT channel, the H2 TS geometries also feature large negative ZPE contributions, due to the exclusion of the large (around 1000i cm⁻¹) imaginary frequencies. Exceptions are the F⁻ + PH₂Cl/Br reactions, where the imaginary frequency is only a few hundred cm⁻¹ and the corresponding motion is also different. The Walden min and FS postmin geometries have quite large positive ZPE corrections, especially when X = F, where, besides the increasing number of DOFs (from 6 to 9 from reactant to min), the forming P–F bond is much stronger than the other P–halogen bonds.

The structural parameters of the stationary points of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions, highlighted in Fig. 8, give us a chance to take a closer look at the steric and electronic-structure effects in these reactions. The Walden-inversion minimum geometry is very product-like, the P–Y bond is the closest to that of PH₂Y in the case of the Br⁻ + PH₂I reaction, therefore this complex can be imagined as the bottom of an (energetic and structural) slide, which turns the reactants immediately into products. Interestingly, in the case of the front-side-attack retention S_N2 pathway, the traditional FS TS geometry is more product-like, while the newly-identified FS TS′ has a rather reactant-like structure. Thus, if we think about Polanyi's rules of thumb,⁶⁶ the front-side-attack route could possibly be promoted both by translational and vibrational excitation. Moreover, in the X ≠ F reactions, which have a barrierless FS path, a longer reaction time can also be favorable. Just as in the case of the Walden minimum, since this deep minimum might be avoided at higher initial translational energies. The FS min is definitely a pre-reaction minimum with its long P–X bonds, especially in the X ≠ F reactions. FS premin, with its nonreactive orientation, shows decreasing X–Y and increasing P–Y bond lengths as the size of Y increases when X = F, while these values alternate in the other three cases. The geometry of FS postmin definitely resembles more to that of the product PH₂X molecule with long P–Y bonds. The H2 TS also features quite long P–Y bonds (Y is always the larger halogen atom of the two in these exothermic reactions considered), which might be a reason why the two H atoms can get closer to each other.

Fig. 8 Selected structural parameters of the stationary points of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions obtained at the CCSD(T)-F12b/cc-pVTZ-F12 geometries. Bond lengths are in Å, angles are in degrees. Underlined values are obtained at the MP2/aug-cc-pVDZ level of theory.

The proton-transfer TS and postmin have very similar structures in these non-identity reactions, as well, however, the TS has a more bent P–H–X angle than the postmin. These stationary points have clearly product-like geometries when X = F, with very similar H–X distances as in the HX molecules, but they become more and more reactant-like (much and only slightly longer H–X and P–Y bonds than in HX and PH₂Y, respectively). In PT TS the H–X distances are, of course, shorter than in PT postmin. H premin, which actually has a product-like geometry (however, it is still a pre-reaction minimum, which can be formed even before an S_N2 reaction, or can induce reactant inversion through pre-DI TS), features only a little longer H–X bonds than the PT postmin in the X = F reactions, and these H–X distances are much longer than in HX for the other three reactions. The P–Y bond lengths are usually larger at H premin than in PT postmin.

In Table 5 we list the reaction energies and the corresponding auxiliary energy corrections for several alternative product channels of the X⁻ + PH₂Y [X,Y = F, Cl, Br, I] non-identity reactions. Here, the lowest-energy channel, even exothermic for all the X = F reactions, is also the formation of a H₂ molecule. For the X = F case this channel also has submerged barriers, as we can see in Fig. 5 and in Table 3. It is followed by the formation of PH in its triplet ground state, and then by the HY + PHX⁻ channel, except for the Br⁻ + PH₂I reaction, where the second most stable products are IBr⁻ + PH₂. The thermodynamically least favorable reaction for X = F is the formation of PH₂⁻ (which is the excited-state version of the XY⁻ + PH₂ channel), while in the X ≠ F case it is the hydride ion formation with relative energies above 70 kcal mol⁻¹. H⁻- and PHF⁻-formations feature increasing reaction energies as we go down in Table 5, whereas the remaining product channels have decreasing reaction energies from Cl to I when X = F, and follow the increasing relative-energy order of X, Y = Cl, I; Cl, Br; and Br, I; except for the PH₂ channel, which is the most endothermic in the case of the X, Y = Cl, Br pair. The reaction energies of these alternative product channels are generally way higher for the X ≠ F reactions, however, e.g., dihalide (ion) production has quite high energy in the X, Y = F, Cl case. The auxiliary energy corrections are similar as in the case of the identity alternative product channels.

Table 5 Reaction energies (in kcal mol⁻¹) of alternative product channels of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] non-identity reactions relative to the energy of the reactants at the CCSD(T)-F12b/cc-pVQZ-F12 level of theory augmented with the auxiliary energy corrections (also in kcal mol⁻¹) at CCSD(T)-F12b/cc-pVTZ-F12 geometries, as described in the computational details section

	QZ^a	δ[T]^b	δ[(Q)]^c	Δ core	Δ rel	Δ ZPE	Classical^g	Adiabatic^h
a CCSD(T)-F12b/cc-pVQZ-F12. b CCSDT/aug-cc-pVDZ − CCSD(T)/aug-cc-pVDZ. c CCSDT(Q)/aug-cc-pVDZ − CCSDT/aug-cc-pVDZ. d AE-CCSD(T)/aug-cc-pwCVTZ − FC-CCSD(T)/aug-cc-pwCVTZ. e Scalar relativistic effects calculated as DK-AE-CCSD(T)/aug-cc-pwCVTZ-DK − AE-CCSD(T)/aug-cc-pwCVTZ. f Zero-point vibrational energy contribution. g Benchmark classical energy, the sum of QZ, δ[CCSDT], δ[CCSDT(Q)], Δcore, and Δrel. h Benchmark adiabatic energy, the sum of benchmark classical energy and ΔZPE.
F^− + PH2Cl
H^− + PHClF	31.42	0.19	0.06	−0.03	−0.08	−3.26	31.56	28.30
ClF + PH2^−	68.88	−0.04	−0.14	0.28	−0.15	−2.46	68.85	66.39
ClF^− + PH2	64.66	−0.47	−0.15	−	−	−2.71	64.04	61.33
ClHF^− + ^1PH	17.93	−3.70	−0.20	0.97	−0.14	−0.98	14.86	13.88
ClHF^− + ^3PH	−8.61	0.05	0.20	0.15	−0.14	−1.01	−8.35	−9.36
PClF^− + H2	−15.51	0.16	−0.07	0.14	−0.11	−3.53	−15.39	−18.92
HCl + PHF^−	−2.04	0.03	0.08	0.14	0.07	−1.98	−1.73	−3.70
F^− + PH2Br
H^− + PHBrF	33.05	0.18	0.06	−0.08	−0.05	−3.27	33.15	29.89
BrF + PH2^−	56.42	−0.02	−0.13	0.81	0.18	−2.32	57.27	54.95
BrF^− + PH2	33.43	−0.39	−0.05	0.54	0.05	−2.48	33.58	31.10
BrHF^− + ^1PH	13.35	−3.69	−0.18	1.08	−0.29	−0.75	10.27	9.52
BrHF^− + ^3PH	−13.19	0.06	0.22	0.26	−0.29	−0.78	−12.94	−13.72
PBrF^− + H2	−17.67	0.15	−0.07	0.20	−0.13	−3.39	−17.52	−20.91
HBr + PHF^−	0.76	0.03	0.08	0.16	−0.02	−2.18	1.03	−1.16
F^− + PH2I
H^− + PHIF	35.39	0.17	0.05	−0.13	−0.07	−3.25	35.42	32.16
IF + PH2^−	38.40	0.01	−0.06	1.35	0.07	−2.15	39.78	37.63
IF^− + PH2	19.63	−0.25	0.05	0.91	−0.08	−2.26	20.26	18.01
IHF^− + ^1PH	10.15	−3.68	−0.15	1.15	−0.29	−0.59	7.17	6.58
IHF^− + ^3PH	−16.39	0.07	0.24	0.33	−0.29	−0.62	−16.04	−16.66
PIF^− + H2	−19.80	0.15	−0.08	0.24	−0.20	−3.24	−19.69	−22.93
HI + PHF^−	3.27	0.03	0.10	0.24	0.02	−2.40	3.65	1.25
Cl^− + PH2Br
H^− + PHBrCl	74.83	0.11	−0.10	−0.24	−0.23	−3.98	74.37	70.39
BrCl + PH2^−	69.32	−0.08	−0.14	0.38	0.04	−2.66	69.51	66.85
BrCl^− + PH2	40.66	−0.30	0.00	0.33	−0.10	−2.63	40.59	37.96
BrHCl^− + ^1PH	54.22	−3.79	−0.29	0.91	−0.34	−3.35	50.70	47.35
BrHCl^− + ^3PH	27.68	−0.04	0.10	0.09	−0.34	−3.38	27.49	24.12
PClBr^− + H2	13.36	0.10	−0.22	0.05	−0.26	−3.85	13.04	9.19
HBr + PHCl^−	27.72	−0.03	−0.07	0.03	−0.13	−2.58	27.52	24.94
Cl^− + PH2I
H^− + PHICl	76.14	0.12	−0.11	−0.28	−0.24	−3.95	75.62	71.67
ICl + PH2^−	58.16	−0.07	−0.11	0.72	−0.04	−2.49	58.66	56.17
ICl^− + PH2	30.77	−0.25	0.03	0.59	−0.16	−2.42	30.97	28.55
IHCl^− + ^1PH	51.44	−3.78	−0.26	0.96	−0.32	−2.90	48.03	45.13
IHCl^− + ^3PH	24.89	−0.03	0.14	0.14	−0.32	−2.93	24.82	21.90
PClI^− + H2	11.00	0.12	−0.23	0.07	−0.32	−3.68	10.63	6.95
HI + PHCl^−	30.22	−0.03	−0.05	0.10	−0.09	−2.80	30.15	27.35
Br^− + PH2I
H^− + PHIBr	84.61	0.12	−0.14	−0.55	−0.21	−4.24	83.82	79.59
IBr + PH2^−	60.31	−0.06	−0.11	0.52	0.07	−2.66	60.73	58.08
IBr^− + PH2	30.85	−0.24	0.04	0.44	−0.02	−2.50	31.06	28.56
IHBr^− + ^1PH	61.70	−3.79	−0.30	0.60	−0.56	−3.82	57.65	53.83
IHBr^− + ^3PH	35.16	−0.04	0.10	−0.22	−0.56	−3.85	34.44	30.59
PBrI^− + H2	16.83	0.12	−0.25	−0.11	−0.33	−3.85	16.25	12.40
HI + PHBr^−	37.89	−0.04	−0.07	−0.03	−0.17	−2.97	37.58	34.62

---

Conclusions

Phosphorus-centered ion–molecule reactions have only been rarely studied so far, therefore in the present work we investigate the energetics of the X⁻ + PH₂Y [X, Y = F, Cl, Br, I] identity and non-identity reactions in detail by determining benchmark energies and geometries for their stationary points. These reactions follow the two main traditional reaction routes: bimolecular nucleophilic substitution (S_N2) and proton-transfer, leading to Y⁻ + PH₂X and HX + PHY⁻, respectively. The S_N2 reaction can proceed through three-different pathways: Walden-inversion, front-side-attack and double- or multiple-inversion, amongst which the first leads to inverted, the second to retained, and the third to inverted or retained, depending on the number of inversions during a reaction, configuration around the central P atom. We consider all relevant stationary points for these reaction paths. Besides the above-mentioned main products, we also identify several other product channels, as well: H⁻-formation, PH₂⁻- and PH₂-formation, ¹PH- and ³PH-formation, H₂-formation and HY + PHX⁻ formation, which reveal the rich chemistry of these ion–molecule reactive systems at the P-center. The benchmark geometries are optimized at the CCSD(T)-F12b/cc-pVTZ-F12 level of theory, where harmonic vibrational analysis is also carried out. Further CCSD(T)-F12b/cc-pVQZ-F12 single-point energy computations are performed at these geometries to minimize the already small basis-set error of the CCSD(T)-F12b method. The following auxiliary energy corrections are also calculated to gain benchmark classical relative energies of the stationary points: (1) core-correlation correction obtained by correlating the sub-valence shell electrons as well, (2) scalar relativistic effects are computed using the second-order Douglas–Kroll Hamiltonian, and (3) post-CCSD(T) energy corrections are also taken into account by considering the higher-order CCSDT and CCSDT(Q) terms. Benchmark adiabatic relative energies are calculated by adding zero-point-energy contributions to the final classical energy values. It turns out that the above energy corrections are absolutely necessary to achieve the sub-chemical (<1 kcal mol⁻¹) accuracy of the benchmark energies, while they also provide additional insight into the electronic-structure properties of these P-centered reactive systems.

Regarding the energetics of the two main reaction channels, in the identity (X⁻ + PH₂X) case of the title reactions, the one with the F ligand shows notable differences with respect to the other halogens: an exothermic PT channel, low-lying front-side-attack transition-state, a very high self-inversion barrier, as well as an extremely deep Walden-well and a H-bonded pre-reaction minimum with a submerged inversional TS are found for X = F. These findings suggest that proton transfer, the front-side retention S_N2 route, and the multi-inversion mechanism are the most competitive with the traditional Walden-inversion in the X = F reaction. For the other halogens we observe increasing relative energies of the stationary points and enhanced preference for Walden-inversion with increasing halogen-atom size. The (exothermic) non-identity (X⁻ + PH₂Y) reactions also differ significantly when F⁻ is the attacking ion compared to the remaining three reactions. We find similar discrepancies here as in the identity case, predicting more indirect dynamics for all the X = F reactions, with deep minima, exothermic PT, competing S_N2 routes, and favorable entrance-channel inversion. Stereoselectivity of the S_N2 reactions, just as in the N-centered case, is certainly questioned in the title systems, which may encourage future dynamical investigations of central-atom effects in S_N2 reactions.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was supported by the National Research, Development and Innovation Office-NKFIH, K-125317; project no. TKP2021-NVA-19, provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021-NVA funding scheme; and the Momentum (Lendület) Program of the Hungarian Academy of Sciences.

References

1. J. Xie, R. Otto, J. Mikosch, J. Zhang, R. Wester and W. L. Hase, Acc. Chem. Res., 2014, 47, 2960.
2. J. Xie and W. L. Hase, Science, 2016, 352, 32.
3. I. Szabó and G. Czakó, J. Phys. Chem. A, 2017, 121, 9005.
4. G. Czakó, T. Győri, B. Olasz, D. Papp, I. Szabó, V. Tajti and D. A. Tasi, Phys. Chem. Chem. Phys., 2020, 22, 4298.
5. R. Wester, Mass Spectrom. Rev., 2022, 41, 627.
6. C. K. Ingold, Structure and Mechanisms in Organic Chemistry, Cornell University Press, Ithaca, NY, 1953.
7. P. Walden, Ber. Dtsch. Chem. Ges., 1896, 29, 133.
8. W. A. Cowdrey, E. D. Hughes, C. K. Ingold, S. Masterman and A. D. Scott, J. Chem. Soc., 1937, 1937, 1252.
9. M. N. Glukhovtsev, A. Pross, H. B. Schlegel, R. D. Bach and L. Radom, J. Am. Chem. Soc., 1996, 118, 11258.
10. A. P. Bento and F. M. Bickelhaupt, J. Org. Chem., 2008, 73, 7290.
11. I. Szabó and G. Czakó, J. Phys. Chem. A, 2015, 119, 3134.
12. I. Szabó and G. Czakó, Nat. Commun., 2015, 6, 5972.
13. J. Zhang, J. Mikosch, S. Trippel, R. Otto, M. Weidemüller, R. Wester and W. L. Hase, J. Phys. Chem. Lett., 2010, 1, 2747.
14. J. Xie, J. Zhang and W. L. Hase, Int. J. Mass Spectrom., 2015, 378, 14.
15. M. Stei, E. Carrascosa, M. A. Kainz, A. H. Kelkar, J. Meyer, I. Szabó, G. Czakó and R. Wester, Nat. Chem., 2016, 8, 151.
16. I. Szabó, B. Olasz and G. Czakó, J. Phys. Chem. Lett., 2017, 8, 2917.
17. X. Ji, C. Zhao and J. Xie, Phys. Chem. Chem. Phys., 2021, 23, 6349.
18. R. Gareyev, S. Kato and V. M. Bierbaum, J. Am. Soc. Mass Spectrom., 2001, 12, 139.
19. P. Beak and J. Li, J. Am. Chem. Soc., 1991, 113, 2796.
20. M. Bühl and H. F. Schaefer, J. Am. Chem. Soc., 1993, 115, 364.
21. M. N. Glukhovtsev, A. Pross and L. Radom, J. Am. Chem. Soc., 1995, 117, 9012.
22. Y. Ren and H. Zhu, J. Am. Soc. Mass Spectrom., 2004, 15, 673.
23. J. Yang, Y. Ren, H. Zhu and S.-Y. Chu, Int. J. Mass Spectrom., 2003, 229, 199.
24. Y.-M. Xing, X.-F. Xu, Z.-S. Cai and X.-Z. Zhao, J. Mol. Struct.: THEOCHEM, 2004, 671, 27.
25. X. Liu, J. Zhang, L. Yang and R. Sun, J. Phys. Chem. A, 2016, 120, 3740.
26. B. Hajdu and G. Czakó, J. Phys. Chem. A, 2018, 122, 1886.
27. A. P. Bento, M. Solà and F. M. Bickelhaupt, J. Comput. Chem., 2005, 26, 1497.
28. S. C. A. H. Pierrefixe, C. Fonseca Guerra and F. M. Bickelhaupt, Chem. – Eur. J., 2008, 14, 819.
29. Y. Ren, X. Wang, S.-Y. Chu and N.-B. Wong, Theor. Chem. Acc., 2008, 119, 407.
30. Z.-Z. Yang, Y.-L. Ding and D.-X. Zhao, J. Phys. Chem. A, 2009, 113, 5432.
31. Y.-L. Ding, J.-R. Mu and L.-D. Gong, J. Chin. Chem. Soc., 2013, 60, 327.
32. T. Matsubara and T. Ito, J. Phys. Chem. A, 2016, 120, 2636.
33. E. Hupf, M. Olaru, C. I. Rat, M. Fugel, C. B. Hübschle, E. Lork, S. Grabowsky, S. Mebs and J. Beckmann, Chem. – Eur. J., 2017, 23, 10568.
34. M. V. J. Rocha, N. W. G. Smith, L. P. Wolters, A. de Cózar, C. Fonseca Guerra, T. C. Ramalho and F. M. Bickelhaupt, Int. J. Mass Spectrom., 2017, 413, 85.
35. T. A. Hamlin, M. Swart and F. M. Bickelhaupt, ChemPhysChem, 2018, 19, 1315.
36. M. Fugel, A. Dittmer, F. Kleemiss and S. Grabowsky, J. Phys. Chem. A, 2021, 125, 4070.
37. A. Á. Dékány, G. Z. Kovács and G. Czakó, J. Phys. Chem. A, 2021, 125, 9645.
38. A. Á. Dékány and G. Czakó, J. Chem. Phys., 2023, 158, 224303.
39. M. A. van Bochove, M. Swart and F. M. Bickelhaupt, J. Am. Chem. Soc., 2006, 128, 10738.
40. O. I. Kolodiazhnyi and A. Kolodiazhna, Tetrahedron: Asymmetry, 2017, 28, 1651.
41. M. A. van Bochove, M. Swart and F. M. Bickelhaupt, ChemPhysChem, 2007, 8, 2452.
42. J. M. Denu, D. L. Lohse, J. Vijayalakshmi, M. A. Saper and J. E. Dixon, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 2493.
43. M. Oivanen, S. Kuusela and H. Lönnberg, Chem. Rev., 1998, 98, 961.
44. S. D. Lahiri, G. Zhang, D. Dunaway-Mariano and K. N. Allen, Science, 2003, 299, 2067.
45. J. Zhang, Y. Xu, J. Chen and D. Wang, Phys. Chem. Chem. Phys., 2014, 16, 7611.
46. X. Niu, P. Liu and D. Wang, J. Phys. Chem. A, 2020, 124, 141.
47. X. Wu, S. Zhang and J. Xie, Phys. Chem. Chem. Phys., 2022, 24, 12993.
48. J. Mikosch, S. Trippel, C. Eichhorn, R. Otto, U. Lourderaj, J.-X. Zhang, W. L. Hase, M. Weidemüller and R. Wester, Science, 2008, 319, 183.
49. J. Meyer, V. Tajti, E. Carrascosa, T. Győri, M. Stei, T. Michaelsen, B. Bastian, G. Czakó and R. Wester, Nat. Chem., 2021, 13, 977.
50. D. Papp and G. Czakó, Chem. Sci., 2021, 12, 5410.
51. A. Giricz, G. Czakó and D. Papp, Chem. – Eur. J., 2023 DOI:10.1002/chem.202302113.
52. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618.
53. T. H. Dunning, Jr., J. Chem. Phys., 1989, 90, 1007.
54. T. B. Adler, G. Knizia and H. J. Werner, J. Chem. Phys., 2007, 127, 221106.
55. K. A. Peterson, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2008, 128, 084102.
56. K. A. Peterson and T. H. Dunning, Jr., J. Chem. Phys., 2002, 117, 10548.
57. M. Douglas and N. M. Kroll, Ann. Phys., 1974, 82, 89.
58. W. A. de Jong, R. J. Harrison and D. A. Dixon, J. Chem. Phys., 2001, 114, 48.
59. K. A. Peterson, D. Figgen, E. Goll, H. Stoll and M. Dolg, J. Chem. Phys., 2003, 119, 11113.
60. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al., Molpro, version 2015.1, a package of ab initio programs, see https://www.molpro.net.
61. J. Noga and R. J. Bartlett, J. Chem. Phys., 1987, 86, 7041.
62. M. Kállay and J. Gauss, J. Chem. Phys., 2005, 123, 214105.
63. M. Kállay, P. R. Nagy, D. Mester, Z. Rolik, G. Samu, J. Csontos, J. Csóka, B. P. Szabó, L. Gyevi-Nagy, B. Hégely, et al., Mrcc, a quantum chemical program suite, see https://www.mrcc.hu.
64. M. Kállay, P. R. Nagy, D. Mester, Z. Rolik, G. Samu, J. Csontos, J. Csóka, P. B. Szabó, L. Gyevi-Nagy and B. Hégely, et al. , J. Chem. Phys., 2020, 152, 074107.
65. G. S. Hammond, J. Am. Chem. Soc., 1955, 77, 334.
66. J. C. Polanyi, Science, 1987, 236, 680.

---