The 1,2-hydrogen shift reaction for monohalogenophosphanes PH2X and HPX (X = F, Cl)

ABSTRACT The aim of the present study was to perform a quantum chemical investigation in the 1,2-hydrogen shift reaction for the PH2X and HPX molecules (X = F,Cl). Several phosphorus–halogen-bearing molecules were studied, including PH2F, PH2Cl, HPF, HPCl, HPFH, HPClH, PFH and PClH. The energies of stationary and saddle points on the ground electronic potential energy surface were investigated with post-Hartree–Fock methods [CCSD(T), MP2, QCISD] and different DFT functionals. The PH2F 1,2-hydrogen shift energy barrier was 75 kcal mol−1 at the CCSD(T) level and only a small increase in this value was observed for the HPF isomerisation. In contrast, the HPCl 1,2-hydrogen shift barrier is higher than the PH2Cl one, which presented a barrier height of 69 kcal mol−1 among CCSD(T) and composite methods. The rate constants of these unimolecular rearrangements varied from 10−44 to 10−38 s−1, and these isomerisation channels exhibited large half-lives. In addition, the heat of formation of each monohalogenophosphane was also calculated. The Quantum Theory of Atoms in Molecules (QTAIM) and Natural Bond Orbital (NBO) analysis were also employed to characterise the differences between the phosphorous–halogen bonds.


Introduction
Phosphorus-bearing species have been detected in different interstellar media [1][2][3][4]. The phosphorus-bearing molecules are also important tracer gases in the atmosphere of giant planets [5,6] and the Titan moon [7] and may also have important implications for the Earth's troposphere [8,9] as well as in the phosphorus cycle [10]. Therefore, due to the importance of phosphorous species, several experimental [11][12][13][14] and quantum chemical investigations [15][16][17][18][19]  molecules with potential implications for atmospheric chemistry or the astrochemistry field. Among the phosphorous-bearing molecules, PH 2 F and PH 2 Cl have previously been the focus of different spectroscopic and quantum chemical studies [20][21][22][23][24][25]. The PH 2 F species was first detected by Andrew and Withnall in 1989 [20], where this molecule was produced from the reaction between phosphine and F 2 in an argon matrix and characterised by a high-resolution infrared analysis. In addition, Beckers [25] employed a different procedure for the gas-phase synthesis of PH 2 F and PH 2 Cl. Five years later, Beckers and co-authors [21,22] reassessed the vibrational modes of PH 2 F and PH 2 Cl as well as their rovibrational constants, and Smart and Ranking [23] performed a structure refinement from their rotational constants. In addition, Drean et al. [24] determined the millimeter-wave transitions of the PH 2 F and PH 2 Cl molecules and performed quantum chemical calculations to determine their electronic properties. Recently, several quantum chemical investigations have focused on analysing pnicogen-bonded complexes [26][27][28][29][30][31][32][33][34][35][36][37][38] using the PH 2 F and PH 2 Cl to understand the nature of these particular bonds. Conversely, other studies investigated analogous phosphorous halide molecules to understand the recoupled pair bonds [39][40][41].
Fourteen years later, Bramwell et al. [42] initiated a study of HPCl, which is considered to be a potential III-V semiconductor growth intermediate, and that study employed a low-resolution monochromator to assign the fluorescence transitions. After this investigation, Lee et al. [43] performed an ab initio investigation to describe the Franck-Condon profile of the HPCl emission spectrum, which was reevaluated by Baraille et al. [44] to include anharmonic corrections. Tackett et al. [45,46] reported the laser-induced fluorescence and single-vibronic-level emission spectra of HPCl under different experimental conditions. In addition, Chau and collaborators [47] performed a Franck-Condon simulation of the singlevibronic-level emission spectrum and chemiluminescence spectrum of HPCl including anharmonicity. In the case of the HPF species, there is any experimental investigation in the literature devoted to a spectroscopic characterisation. On the other hand, Jursic [48] described the ionisation potential and electron affinity of HPF by employing MP2 and density functional theory methods, and Schurmann et al. [49] reported the ab initio characterisation of its lowest electronic states.
The purpose of this investigation is to characterise the energy barriers associated with unimolecular arrangements of the PH 2 X and HPX molecules (X = F, Cl). To date, only analyses of the interconversion energy barriers linked to changes in the symmetry of PH 2 X molecules from a C s point group to a C 2V point group have been reported [50][51][52], and no descriptions of the 1,2-hydrogen shift reactions involved in these monohalogenophosphanes in the ground state have been described. Therefore, the objective of this study is to characterise the PH 2 X and HPX 1,2-hydrogen shift reactions: (1) PH 2 X → HPXH, and (2) HPX → PXH reactions (X = F, Cl). The energies of stationary and saddle points on the ground electronic potential energy surface were investigated with coupled-cluster and Møller-Plesset perturbation theories, as well as with quadratic configuration interaction calculations and 32 different density functional theory (DFT) methods. Furthermore, a Natural Bond Orbital (NBO) analysis and a topological analysis using the Quantum Theory of Atoms in Molecules (QTAIM) were performed to gain insight into the phosphorus-halogen bonds.

Methodology
All of the electronic structure calculations in this study were performed with the GAUSSIAN 09 quantum chemistry code [53]. The stationary points on the potential energy surface were characterised by the absence of imaginary frequencies. The transition state structures were confirmed by a harmonic frequency analysis, and all of the connections were evaluated by intrinsic reaction coordinate (IRC) calculations [54,55]. The stabilities of these Hartree-Fock wavefunctions were tested with respect to relaxing various constraints: allowing a restricted HF determinant to become unrestricted, allowing orbitals to become complex, and reducing the symmetry of the orbitals. If some instability was found, the wavefunction was reoptimised with the appropriate reduction in constraints and repeating the stability tests and optimising until a stable wavefunction was found. For a doublet system is expected a spin-squared operator value of 0.75, which was employed to evaluate the spin contamination in the unrestricted wave functions (the case of HPX isomers). The values ranged from 0.75 to 0.78 among the methodologies demonstrating that spin contamination was not significant for these structures.
The bond delocalisation index (bdi) [99], the Wiberg bond index (bo wi ) [100], natural resonance theory bond index (bo nrt ) [101][102][103], atom-atom overlap natural atomic orbital bond index (bo nao ) [104], atom-atom overlap natural localised molecular orbitals/natural population analysis bond index (bo nlmo ) [105] and the Mayer bond order (bo ma ) [106] were employed for the bond order analysis. For the atomic charge distribution, different population methods were employed as follows: Mulliken, Lowdin, QTAIM [107], electrostatic potential (Chelp) [108], electrostatic potential using a Grid base method (ChelpG) [109], the Merz-Singh-Kollman method (MK) [110,111] and natural population analysis (NPA) [104]. The NBO and QTAIM analyses were computed with the NBO6 [112] and AIMALL program [113], respectively. The bond order, atomic charge distribution, NBO and QTAIM analyses were predicted at B3LYP/cc-pVTZ level of theory. Moreover, in order to predict the rate constants were employed the canonical transition state theory (cTST, at 298.15 K) with W1BD due to the good performance of the Martin's W1 methods to thermochemical properties [114][115][116]. Figures 1 and 2 show a schematic representation of the H 2 PX and HPX 1,2-hydrogen shift reaction pathway (X = F, Cl) with the geometric parameters of the isomers as well as the transition state structure geometries determined with the QCISD/6-311G(2df,2pd), MP2/cc-pVQZ and M06/cc-pVQZ methods. An analysis of the results in Figure 1 for the PH 2 X structures indicates good agreement with those predicted by Drean et al. [24] who employed the MP2, CCSD and CCSD(T) methods using the cc-pVQZ basis sets. For example, Drean et al. [24] estimated the P-F and P-Cl bond lengths of PH 2 F and PH 2 Cl to be 1.6081 and 2.0735Å, respectively, with the CCSD(T)/cc-pVQZ method. These results are very similar to those demonstrated in Figure 1 for both structures. The estimated geometric parameters for HPCl (see Figure 2) are also in good agreement with the values computed by Chau et al. [47] with the CCSD(T) level employing aug-cc-pwCVQZ, cc-pV5Z, and aug-cc-pVQZ basis sets. Figure 3 demonstrates the mean absolute error (MAE) between the calculated harmonic frequencies and the experimental results for PH 2 F [21,22,25], PH 2 Cl [21,22,25] and HPCl [43,45]. Examining the MAE values among the monohalogenophosphanes, it is important to note that the methodologies are following almost the same behaviour for each molecule. In this aspect, B97D (15-18 cm −1 ), BP86 (7-8 cm −1 ) and BPW91 (10-11 cm −1 ) functionals provide the most accurate description of the vibrational frequencies versus experimental values. In the matter of PH 2 F, G4MP2 (15 cm −1 ) and B97D (16 cm −1 ) present almost the same calculated error in reproducing the observable vibrations. On the other hand, the largest MAE values were predicted by M06-HF (85-91 cm −1 ), LC-wPBE (61-64 cm −1 ) and MP2 (64-66 cm −1 ) methods. Furthermore, our data indicate that an increase in the percentage of Hartree-Fock exchange above 50% among the DFT functionals deteriorates the accuracy of the harmonic frequencies when compared with observed frequencies, leading to errors above 46 cm −1 . Table 1 shows the relative energy profile for the isomerisation pathway of PH 2 F → HPFH reaction. The barrier height for the PH 2 F 1,2-hydrogen shift route varied from 72.3 to 79.3 kcal mol −1 among the different methodologies. In comparison to analogous phosphorous 1,2-hydrogen shift reactions, H 2 PN has a similar energy barrier [17]. The H 2 PS [15] and H 2 P 2 [117] isomerisation routes exhibited barrier heights that were lower than those studied here. The smallest values for the PH 2 F → HPFH barrier energy were estimated by the functionals with the small percentage of exact Hartree-Fock exchange, while the post-Hartree-Fock and composite methods demonstrated values around 75 kcal mol −1 . Nevertheless, the reverse pathway gives a crucial picture to evaluate the performance among the different methods. The barrier energy for the HPFH → PH 2 F reaction is estimated as 1.4 and 1.8 kcal mol −1 when applied to the CCSD(T)/CBS and W1BD methods, respectively. In parallel with CCSD(T)/CBS and W1BD values, several DFT methods have also demonstrated a good performance as the popular M06 Truhlar's functionals (M06, M06-2X, M06-HF). In contrast, a significant barrier height was calculated with QCISD/6-311G(2df,2pd) (5.8 kcal mol −1 ) and BhandHLYP/cc-pVQZ (4.2 kcal mol −1 ) for the reverse pathway. Furthermore, among the 27 DFT methods employed in this reaction, 12 DFT functionals predicted that TS1 lies below HPFH. In the case of TPSSh, O3LYP, B97D, BP86 and BPW91 functionals, we could not obtain TS1 and HPFH optimised structures. Table 2 presents the energetics profile for the PH 2 Cl 1.2-hydrogen shift unimolecular rearrangement. The variations among the methodologies in the description of PH 2 Cl → HPClH barrier energy were very similar to what we predict for the fluorine analogous reaction. The barrier height for the PH 2 Cl isomerisation ranged from 66.4 to 72.6 kcal mol −1 . The lowest and highest values were estimated with M06-HF and M06L DFT functionals, respectively. LC-wPBE, wB97 and M05 functionals also overestimated the PH 2 Cl isomerisation barrier energy in comparison to the post-Hartree-Fock and composite methods. These show that the TS2 lies around 69 kcal mol −1 above PH 2 Cl, while MP2/cc-pVQZ estimates this barrier height to be 70.7 kcal mol −1 . In contrast, the HPClH → PH 2 Cl reverse route varied from 2.5 to 4.9 kcal mol −1 . The largest values predicted for the reverse mechanism were obtained by O3LYP, M06-HF, M06L, BP86 and BPW91 functionals.

Results and discussions
The rate constants and half-lifetimes (t 1/2 ) of each reaction route are presented in Table 3 (Table S4 in the supporting material contains the equilibrium constants of each reaction). The PH 2 F → HPFH and PH 2 Cl → HPClH rate constants are estimated in 9.57 × 10 −43 and 6.33 × 10 −38 s −1 , respectively. Based on these results, half-lifetimes (t 1/2 ) for the PH 2 F and PH 2 Cl isomerisation route are 1.04 × 10 42 and 1.58 × 10 37 s, respectively. Moreover, a glance at the adiabatic phosphorous-halogen dissociation of PH 2 F (varying from 101 to 115 kcal mol −1 ) shows that TS1 lies below the separated reactant energy level, which is a very similar picture for TS2 although with a small energy gap when compared to its respective reactants. Examining the phosphoroushalogen dissociation energies of the HPXH complexes [BDE(HP-FH) = −22.74 kcal mol −1 ; BDE(HP-ClH) = −15.43 kcal mol −1 ; with W1BD], these exothermic results demonstrate small stability, which would make the characterisation of these species difficult under experimental conditions because dissociation of these metastable isomers into PH(³ −1 ) and HX products is expected. Based on these results for the PH2+X reactions, a low energetic cost is required to overcome the rearrangement barriers and subsequent heterolytic dissociation of the HPXH exit complexes -especially for PH2F -because this species resides in a deep potential energy well with respect to the separated reactants.
The results obtained for the phosphorous-halogen bonds by different bond order indexes as well as the covalent/ionic character are listed in Table 4, while Table 5 shows the topological results for the phosphorous-halogen bond critical points (BCPs) obtained using a QTAIM analysis as the electronic charge density [ρ(r)] and its Laplacian [ 2 ρ(r)], the total energy density [H(r)], the ellipticity (ϵ), and the relationship between the modulus of the local potential energy and the local energy density [|V(r)|/G(r)]. The configuration changes for the conversion of PH 2 F to HPFH resulted in the phosphorous-fluorine bond becoming more delocalised on the fluorine, followed by an increase in its ionic character, which leads to a covalent interaction with a bond order of approximately 0.2 for HPFH. At this stage, it is important to note that although the different methodologies produce very small bond order values suggesting a covalent interaction, the bond delocalisation index and the NRT result confirmed the existence of a covalent bond. For HPFH, NRT also estimated a major resonance weight for the P-F σ bond (98%). From the perspective of the QTAIM analysis, in the PH 2 F → HPFH isomerisation process, a decrease in the electronic charge density on the P-F bond was predicted in addition to an increase on H(r) as well as in its ellipticity value, which indicates the instability of the P-F bond in HPFH. This result is consistent with  the calculated dissociation energy. In contrast, for the |V(r)|/G(r) relationship, a value between 2 and 1 was predicted for PH 2 F and HPFH, corresponding to partial covalent character in both cases, and a small difference between the values of each isomer was observed. In addition, a decreasing in the steric energy was estimated for the conversion of PH 2 F to HPFH based on the natural localised molecular orbital steric exchange energy. The main reason for this steric exchange energy is the pronounced steric repulsion difference between the donation from the fluorine lone pair electrons into P-H bonds [ lp (F) → σ P-H ] in PH 2 F and the lp (P) → σ F-H steric repulsion energy in HPFH. For the conversion of PH 2 Cl to HPClH, only a 30% decrease in the phosphorous-chlorine bond order was estimated by most of the methodologies with exception of the BDI and NRT results. Based on the NRT analysis, in PH 2 Cl, the P-Cl σ bond exhibited major covalent character compared to the P-F bond in PH 2 F. This result is also confirmed by the results provided by the QTAIM analysis, which provides a comparison between the |V(r)|/G(r) relationship results that indicated that the PH 2 Cl value is nearly two times larger than the PH 2 F value. The high covalent nature of the phosphoroushalogen bond in PH 2 Cl compared to that of the P-F bond is also consistent with the 2 ρ(r) shape of each bond in Figure 4. This nature is also enhanced by the more negative H(r) value of the P-Cl bond, indicating the large electronic charge concentration between the P and Cl basins. In addition, the second order perturbation theory analysis of the Fock matrix in the NBO basis provides useful insight to fully understand the stability of each isomer. For example, the main contribution to the stability of PH 2 Cl comes from the donation of lp (Cl) into the P-H σ anti-bond [ lp (Cl) → σ * P-H ] with an equal contribution from the charge transfer of lp (Cl) into the Rydberg phosphorous orbitals [ lp (Cl) → ry (P)]. A similar result was also observed for PH 2 F. However, the main interactions between HPFH and HPClH are different. In HPFH, the lp (P) → σ * F-H and lp (F) → ry (H) donations represent equal contributions to the stability of this species, and in HPClH, the lp (P) → σ * Cl-H donation is two times higher than the lp (Cl) → ry (P) one. Comparing between PH 2 Cl and HPClH, the steric energy exhibits a decrease of only 4.2 kcal mol −1 . Based on the pre-NLMO overlap interaction energies of HPClH, in addition to the steric repulsion from the lp (P) → σ Cl-H donation, it is also predicted an equal steric contribution associated with the lp (P) → lp (Cl) donation. Tables 6 and 7 show the potential energy profile of HPF and HPCl 1,2-hydrogen shift mechanisms. It was predicted that TS3 lies 77.1 kcal mol −1 below HPF with CCSD(T)/CBS. This is very similar to the results obtained with W1BD and CBS-QB3. G4MP2 estimated the largest barrier energy. With regard to the HPF phosphoroushalogen dissociation energy, it can be noted that TS3 lies below the reactants (PH+F) among all methodologies. Another important aspect is the reverse pathway. Examining the reverse pathway (PFH → HPF), most of the DFT methods give negative barrier energies as well as MP2. The exception among the functionals is TPSSh, BHandHLYP, M06-HF and B97D, although most of these functionals showed values larger than what was predicted by post-HF and composite methods. Other functionals predicted barrier heights lower than 0.2 kcal mol −1 as in the case of M06 and mPW2PLYP. This was also seen at the CCSD(T) level. In addition, the rate constant for the HPF isomerisation is predicted to be 2.89 × 10 −44 s −1 with a t 1/2 of 3.46 × 10 43 s. Of note, some DFT methods failed in the optimisation of the PFH structure, and thus their results are not presented in Table 6, including O3LYP, tHCTHhyb, BP86 and BPW91 functional results.
In contrast, most of the DFT methods yield similar barrier energies for the HPCl unimolecular arrangements. In the HPCl → PClH reaction, the barrier energy varied between 71.3 and 77.5 kcal mol −1 with the exception of BHandH (69.5 kcal mol −1 ) and M06-HF (68.5 kcal mol −1 ) that presented the lowest values. In contrast, the B98 functional showed the largest barrier energy (84.7 kcal mol −1 ) among the methods presented in Table 7. An inspection of the HPCl dissociation energy results shows that TS4 lies below the reactants with the exception when we applied QCISD/6-311G(2df,2pd). In addition, in the PClH → HPCl route, only the BHandH functional gives a negative barrier, while the post-HF methods have barrier energy ranges from 2.4 to 2.9 kcal mol −1 . The composite ones (W1BD, G4MP2, CBS-QB3) vary from 6.8 to 7.4 kcal mol −1 . In comparing the DFT methods with the results obtained among post-HF and composite methodologies, there is a good agreement in the description of the reverse mechanism (PClH → HPCl). In this case, the B98 functional also had the largest barrier energy with a value of 15.5 kcal mol −1 . Interestingly, there is little difference between unrestricted and restricted open-shell wave functions, indicating that the static correlation effects can be well described at CCSD(T) level, which is valid for HPF and HPCl isomerisations and also for the H-elimination routes. Furthermore, examining the kinetic rate constant of these reactions, the HPCl isomerisation is predicted to be 6.26 × 10 −41 s −1 (t 1/2 = 1.60 × 10 40 s), and the reverse pathway is 1.99 × 10 10 s −1 (t 1/2 = 5.03 × 10 −11 s).
Further analysis of the electronic properties of the 1,2-hydrogen shift pathways of the HPX structures indicated a decrease in the phosphorous polarisation and bond order values from HPX to the PXH species, which is comparable to that observed in the PH 2 X reactions. In addition, the overall result that can be extracted from these bond order analyses is that very small values were predicted by NLMO/NPA and NAO for PFH, which makes the existence of the covalent/ionic nature of this phosphorous-halogen bond more difficult as previously demonstrated in other studies [118,119]. Based on the topological analysis from the QTAIM results, the HPCl |V(r)|/G(r) value for the P-Cl σ bond indicates the strongest covalent nature among the HPX isomers, and further evidence is obtained by examination of the 2 ρ(r) shape of this bond, as shown in Figure 5. Figure 6 presents the transition state structures involving the hydrogen release pathways of PFH and PClH. Predictions of their relative energies are shown in Tables 6  and 7, respectively. Based on the CCSD(T) results, TS5 lies 27 kcal mol −1 above PFH, while the reverse route is accessible by overcoming an energy barrier of 32 kcal mol −1 . In the case of QCISD/6-311G(2df,2pd) and MP2/cc-pVQZ, the PFH → PF+H barrier heights are 32.2 and 41.5 kcal mol −1 , respectively. In the analogous chlorine process, MP2 and QCISD also give large values in comparison with G4MP2 and CCSD(T) results. Based on the computed energies from W1BD, the forward and reverse PClH H-release barrier energies are 8.4 and 6.7 kcal mol −1 , respectively. In addition, the forward and backward H-release channel of PFH leads to rate constants of 3.14 × 10 −10 and 3.44 × 10 −11 s −1 , respectively. For PClH H-loss rate constants, the forward and backward reactions are faster than the analogous fluorine reaction. The PXH exit complexes are also estimated to have an exothermic dissociation energy for the phosphorous-halogen bond [BDE(P-FH) = −31.56 kcal mol −1 ; BDE(P-ClH) = −24.81 kcal mol −1 ; with W1BD]. These values demonstrate their small stability and their metastable nature, which is similar to the literature for other exit complexes produced by a 1,2-hydrogen shift reaction [120][121][122]. Moreover, Lopez et al. [123] reported a hydrogen elimination mechanism for PFH + and PClH + . Lopez and co-workers [123] also demonstrated that the hydrogen elimination routes of PFH + and PClH + have barrier heights that are two-fold higher than those predicted for the neutral systems in this study. Nevertheless, Lopez et al. [123] predicted similar geometries for these Table . Relative energy profile (in kcal mol − ) for the HPCl ,-hydrogen shift mechanism and for the H-elimination reaction as well as the dissociation energy (in kcal mol − ) of each species by different methods. TS structures, and the major differences involved an elongation of the halogen-hydrogen distances and shorter phosphorous-halogen bonds. Further insights can be provided when we examine the results between LC-wPBE and CAM-B3LYP. Longrange exchange effects are reduced in CAM-B3LYP [82]. With regard to the reaction involving fluorine atoms, the LC-wPBE shows a better performance in the description of the barrier energies; however in the analogous chlorine reactions, the CAM-B3LYP functional demonstrated a better accuracy in comparison with post-HF and composite methods. These results highlight that long-range exchange effects with exchange functionals can better reproduce the PH 2 F and HPF reaction pathways. This was not applied to chlorine analogous reactions.
The influence of dispersion correlation effects in the use of double hybrid density functionals is important. The intramolecular dispersion correlation effects lead to a small influence on the description of these reaction pathways due to the addition of a dispersion correction in these double hybrid functionals, where only small decreases in the barrier energies and energy gaps are seen. Important information is achieved in comparing B3LYP and B2PLYP functionals because the Beck exchange and Lee-Yang-Parr correlation functionals are presented in both DFT methods [64,66,82]. In this case, the influence of MP2 correlation energy on Khon-Shan virtual orbitals were larger in fluorine 1,2-H shift reactions than in chlorine analogous ones. On the other hand, a look at the Hrelease reaction pathways shows a decreasing of almost 5 kcal mol −1 for the PClH → PCl+H barrier energy when  used in the B2PLYP functional instead of B3LYP. The perturbative correlation effect was smaller in the analogous fluorine process.
The B97 group of functionals used here also provides interesting information for B97D, B97-1, B97-2, B98, wB97, wB97X and wB97XD. In this aspect, the letter 'w' refers to the correction of the asymptotic behaviour of a pure exchange functional by the long-range Coulomb operator. The letter 'X' refers to the inclusion of the HF exact exchange in the short-range exchange, and 'D' indicates the addition of a dispersion correction [80,81]. The wB97, wB97X and wB97XD functionals only have a small effect on the barrier energies when the variations do not reach 3 kcal mol −1 . The dispersion correlation effects were more pronounced in the H-release routes. It is interesting to note that the B97D semiempirical dispersioncorrected functional only slightly underestimates the barrier height energies when compared to the composite methods, which has been seen in previous studies [124]. B97-1, B97-2 and B98 variants refer to hybrid GGA functionals [68][69][70]. The main difference between B97-1 and B97-2 is the fact that B97-2 presents an additional fitting of Zhao, Morrison, and Parr exchange-correlation potential provided from ab initio Brueckner doubles or MP2 electron densities [69]. In contrast, B98 presents the same formulas of B97-1 but was reoptimised with the extended G2 test set [70]. In this B97 hybrid GGA group, To elucidate the calculation of the heat of formation, a methodology that combines the ab initio and experimental values, which was proposed by Curtis et al. [125], was employed. In the total atomisation energy, the experimental heats of formation at 0 K for the elements were provided by NIST-JANAF [126]: P ( H f 0K = 75.42 kcal mol −1 ), F ( H f 0K = 18.47 kcal mol −1 ), Cl ( H f 0K = 28.59 kcal mol −1 ) and H ( H f 0K = 51.63 kcal mol −1 ). The fluorine (−0.39 kcal mol −1 ) and chlorine (−0.84 kcal mol −1 ) spin-orbit corrections were obtained from Moore [127]. Table 8 lists the heat of formation at 0 and 298 K calculated for each monohalogenophosphane molecule. It is important to note that Vasiliu et al. [128] employed a high-level electronic calculation for the heat of formation of H 2 PCl based on the CCSD(T)/CBS level of theory augmented by a number of additional corrections, and at 0 K, the authors estimated a value of −15.4 ± 0.3 kcal mol −1 . In this investigation, the H 2 PCl heat of formation was predicted to be −17.23 and −16.36 kcal mol −1 using W1BD and G4MP2, respectively. At 298 K, using G4MP2, the PH 2 Cl heat of formation was estimated to be −18.74 kcal mol −1 , which is only 1 kcal mol −1 higher than the result reported by Vasiliu et al. [128]. In general, for these phosphorous molecules, the difference among the three methods employed here does not reach 1 kcal mol −1 , but the difference among these methods can reach 2.5 kcal mol −1 in the case of PFH.
The molecular electrostatic potential (MEP) for each structure is shown in Figure 7. The electrostatic potential is an important tool for interpreting the potential sites for an electrophilic or nucleophilic attack where the negative regions are in red and the positive regions in blue. For the PH 2 X and HPX structures, positive charges are seen for the phosphorous atoms and negative ones for the halogen atoms. The HPXH structures have the opposite distribution. PFH has predicted negative charges for phosphorous and fluorine, and PClH is estimated to have a positive charge in chlorine and a negative one for phosphorous. Here, it is important to mention that MEP is a more reliable indicator by the effect of charge redistribution in forming molecules than atomic charges. This is because the atomic charges are not physical observables [129].
Nevertheless, due to the use of partial atomic charge in the characterisation of monohalogenophosphane structures to interpret the molecular aspects of the pnicogen bond [26][27][28][29][30][31][32][33][34][35][36][37][38], we now turn our attention to the atomic charge distribution among the several structures. In this context, it is interesting to note that there is a sensitivity to the population method employed to predict the correct nature of the partial atomic charge distribution [130][131][132][133][134]. In these analyses, eight different population methods were used including NPA, QTAIM, Chelp, ChelpG, MK, GAPT, Mulliken and Löwdin. Of note, QTAIM, Mulliken, NPA and GAPT population methods were good options for PH 2 X and HPX structures ( Table S5 in the supporting material contains phosphorous and halogen charge values of each structure). The electrostatic potential-derived-based methods estimate almost neutral charges for phosphorous atoms, while the Löwdin charge failed to describe the correct sign for the chlorine charge in HPCl. The QTAIM charges lead to a higher separation between the phosphorous and halogens charge, which is due to the QTAIM partition scheme. The good efficiency of the QTAIM methodology compared to other population methods was also confirmed for other systems [18,133,[135][136][137]. In HPFH, the high charge modulus is more concentrated in the fluorine than in the phosphorous atom [|q(F)| > |q(P)|] among all the population methods. For HPXH and PXH structures, QTAIM, GAPT and Löwdin failed to predict the correct nature of the charge distribution for these exit complexes.

Conclusion
This study reports the electronic properties that characterise the HPX and H2PX monohalogenophosphane 1,2-hydrogen shift reactions. For the PH 2 F 1,2-hydrogen shift reaction, a tight transition state structure with a barrier height ranging from 72 to 79 kcal mol −1 was calculated. The conversion between PH 2 Cl and HPClH occurs through a less tight transition state structure with a barrier height that ranged from 66 to 73 kcal mol −1 . Based on the cTST results, large half-lifetimes are expected for the PH 2 F (10 42 s) and PH 2 Cl (10 37 s) isomerisation routes. In the HPF isomerisation, the energy barrier was between 76.9 and 79.5 kcal mol −1 , and the HPCl 1,2-hydrogen shift reaction was accessible by overcoming an energy barrier of approximately 73 kcal mol −1 . In addition, the application of cTST to elucidate the H-release routes showed that the chlorine reaction is faster than the analogous fluorine process.