3.1. Stability and interior weak interaction of THE* and PHE*
At the PBE0-ma-def2TZVP level, two stable structures of hydrated electrons were obtained with an implicit water model (IEF-PCM). According to the shape, they were named "tetrahedral" and "prism" e−(aq) respectively (Fig. 1A and B). The structure of tetrahedral hydrated electron (THE*) and prism hydrated electron (PHE*) goes along with previous research (Kumar et al., 2015; Herbert, 2019). Especially, THE* model has also been observed both with the latest cutting-edge machine learning model (Lan et al., 2021) and the many-body wave-function theory (Wilhelm et al., 2019), therefore it is used as the default model in our later AIMD simulations. The extra electron in THE* is concentrated in the centroid and has a significant sigma characteristic (Fig. 1A). THE* prevents the formation of at least four hydrogen bonds, because hydrogen atoms pointed to the e−(aq) position, thus no longer to be the H bond donors. Therefore, THE* might dramatically change the local structure of water through an asymmetric electronic effect (Head-Gordon and Johnson, 2006). PHE* also has a sigma electronic characteristic, which prevented at least six hydrogen atoms from becoming the hydrogen bond donor. However, its water molecule orientation is significantly unlike that in THE*. A weak H-H interaction was suspected in PHE*, which goes along with Wilhelm’s results (Wilhelm et al., 2019). It is implied that the intermolecular interactions may have significantly different in THE* or PHE* structure.
The interaction among the water molecules was investigated with IGMH method (Fig. 2). In THE* system, there is almost no direct interaction between water molecules, which is manifested by the absence of the IGMH isosurface map between water (Fig. 2a). The prominent attractive interaction occurs between water and the central position of THE*, where a blue map stripped with green edges appears. This implies that THE* structure is dependent on the central electron attraction. THE* structure exhibits a higher energetically stability due to the absence of the prominent repulsive effects (Fig. 2a).
The intermolecular interaction in PHE* is divided into two parts. One part is the prominent interaction between the water molecules forming the bottom surface of the prism, and the other part is the relatively weak interaction in the central region. The water molecules on the bottom surface are subject to both the hydrogen bonding effect (blue region) and the steric effect (red region), forming a stable triangle-like structure (Fig. 2b). But the strength of the interaction between these two triangles is much less than that inside the triangle, showing the characteristics of van der Waals interaction (green region) (Fig. 2b). Therefore, PHE* structure is less stable than THE*, which is in agreement with the conclusion of vertical ionization potential (seen in Chap. 3.2).
The molecular orbitals of THE* and PHE* without an extra electron are shown in the Fig. 1a and b, the LUMO orbitals of THE0 (THE0 indicates a “THE*” w/o an extra electron but with an identical structure, same in PHE0) and PHE0 are not concentrated on a sole molecule or an atom. It is curled in a cavity composed of water molecules to form a spherical-like structure, implying that it is an ideal empty orbital to accept an extra electron. Spin density distribution has supported this hypothesis (Fig. 1A and B). In another hand, the extra electron has a low probability to "collapse" onto a sole water. In THE*, the extra electron mainly resides in the cavity and has a probability to collapse to one of the water molecules in the structure, as indicated by the spin density isosurface (Fig. 1A). Water molecules in the PHE* have a lower probability of gaining an extra electron (Fig. 3b, c), which was believed that is related to the number of water molecules, as the larger number of water molecules "dilutes" the probability to some extent (Fig. 1B). A partly reduced water molecule may be the starting point of the single electron transport (SET). To fully reveal the influence of the extra electron on a single water, a water anion is compared with THE* and PHE* by using the isosurface projection of the spin density on the water molecular plane (Fig. 3). Water molecule belongs to the C2v point group with a good symmetry along its angular bisector. Therefore, it is suspected that if an extra electron was injected, it would populate here in order to maintain symmetry. Contour line map of spin density in molecule plane suggested that the extra electron populated in the region between hydroxyl groups and induced a structure loosening effect (Fig. 3a). Compared with the neutral status, the structure of water anion changed insignificantly, an extra electron slightly increases the O-H bond length from 0.959 Å to 0.981 Å and causes ∠HOH to shrink from 104.4° to 101.1°. It is generally believed that the H-O covalent bond in the water molecule consists of the σ orbital of H and one of π orbitals of O, the diffusion of the extra electron leads to an increased electron on the H, enhancing the antibonding orbital, resulting in an increase in the bond length; H atom shifts in the direction of the extra electron due to electrostatic attraction, narrowing ∠HOH. The geometric of the extra electron spin density was also noticed to be especially close to the LUMO of neutral water doped with a little LUMO + 1 (Fig. S1 1, 2). The nature of e−(aq) reduction ability was investigated by vertical ionization potential (VIP, VIP = EN-EN+1). Calculations show that THE* has a higher VIP than PHE* (27.2 kcal mol− 1 vs 23.9 kcal mol− 1), indicating that PHE* is easier to release its extra electron.
All the above discussions suggested that e−(aq) was prefer to shrink to a cavity composed of water molecules, and the probability of the water molecules that make up the cavity getting additional electrons is relatively low and is inversely proportional to the number of water molecules; The extra electron inside the anionic water mainly occupies its "LUMO" orbital, exhibiting the mixed characteristics of σ-π electrons, of which the σ feature occupies a dominant position. THE* is more stable than PHE* from the energy aspect.
3.2. Organic halides reaction with hydrated electron (e-(aq))
Previous studies suggested two reaction pathways (Fig. 4) depending on pH for the reduction cleavage of the carbon-halogen bond with e−(aq) (Mohan et al., 1992; Andrieux et al., 1997; Yuan et al., 2015). The first one is the stepwise or concerted pathway in which the pollutants directly react with e−(aq) through a SET process before the anionic radical formation or carbon-halogen bond cleavage (Fig. 4a). Another pathway is that the SET process will firstly occur between the e−(aq) with the proton (Fig. 4c), and the formed hydrogen radicals (•H) could induce the cleavage through an induced polarization effect (Fig. 4b) (Calza and Pelizzetti, 2004; Zona et al., 2008; Yuan et al., 2015). Although the reaction rate of •H with halogenated organics (~ 107 M− 1 s− 1) (Calza and Pelizzetti, 2004) is slower than that of e−(aq) (~ 108 M− 1 s− 1) (Yuan et al., 2015). However, the reaction rate of e−(aq) with protons can reach ~ 2.2×1010 M− 1 s− 1 (Calza and Pelizzetti, 2004), so the •H contribution cannot be ignored in the acidic or proton-sufficient solution, e.g. water (Reitmeier et al., 2006). Therefore, two situations needed to be discussed. The first is the direct binding of e−(aq) to the pollutants, in this case, the vertical electron affinity (VEA) of pollutants and the carbon-halogen bond dissociation energy were discussed. The second is the reaction of hydrogen radicals with halogenated molecules, where the transition state and reaction energy barrier were discussed in the next chapter.
VEA is the energy difference between neutral molecules and monovalent anions with an identical geometry (VEA = EN-EN+1, N represents the number of electrons in the initial state). The results suggested that the VEA of halogenated molecules is significantly different in gas or in water (Table 1). For example, VEA of fluorobenzene is -0.573 eV in the gas phase, indicating that it is almost impossible to accept an extra electron, while it has a positive VEA at 0.839 eV in water, suggesting it could act as an electron acceptor. The VEA of the other aryl halides is all negative in the gas phase, indicating that their unstable binding to electron. Specifically, the VEA of chlorobenzene and bromobenzene is relatively close, is -1.183 eV and − 1.123 eV, respectively, while the VEA of fluorobenzene is -0.573 eV. These three molecules have very close VEA values in aqueous solutions, ~ 0.850 eV, in which the relatively largest VEA comes from chlorobenzene, the smallest is from fluorobenzene, shown in Table 1. It suggested that fluorobenzene is the hardest one to be reduced. For the alkane substitutes, VEA in the gaseous state has a smaller variance than in aqueous solutions. The VEA of fluoromethane, chloromethane and bromomethane in the gas phase are − 0.646 eV, -0.571 eV and − 0.518 eV, respectively. In water, they are 0.398 eV, 0.584 eV and 0.878 eV, respectively. Obviously, the results of VEA can only reveal the electron-accepting mechanism, it could not fully disclose the detail of the cleavage reaction. Bond dissociation energy (BDE) might provide a deeper insight as it is considered a thermodynamic determinant of bond cleavage.
Table 1
Vertical electron affinity (VEA) and bond dissociation energy (BDE) of halogenated molecules in gas or water
| | VEA (gas) eV | VEA (aq) eV | BDEN (aq) kcal mol-1 | BDEN+1(aq) kcal mol-1 |
| Fluorobenzene | -0.573 | 0.839 | 130.795 / 127.2 | -14.194 |
| Chlorobenzene | -1.183 | 0.880 | 99.654 / 97.1 | -26.219 |
| Bromobenzene | -1.123 | 0.872 | 91.027 / 84.0 | -19.674 |
| Fluoromethane | -0.646 | 0.398 | 122.430 / 115.0 | -32.737 |
| Chloromethane | -0.571 | 0.584 | 93.920 / 83.7 | -38.770 |
| Bromomethane | -0.518 | 0.878 | 94.943 / 72.1 | -25.630 |
| (CCSD(T) aug-cc-pVTZ tightSCF coupling with SMD water model). N in term BDEN indicates the number of electrons in a neutral state, hence N + 1 indicates a monoanionic state. Bold numbers are the experimental values (Blanksby and Ellison, 2003). |
For the aryl halides, the Carbon-halogen BDE can be arranged from large to small in the order of Fluorobenzene > Chlorobenzene > Bromobenzene, a trend consistent with experimental values (Daasbjerg, 1994; Blanksby and Ellison, 2003). For the alkyl halides, the calculated BDE also has a similar trend to the experimental values, namely the C-F bond is the strongest and the C-Cl and C-Br bonds are weaker. The calculation results of C-Cl and C-Br in this study are slightly different from the other report (Verevkin et al., 2003), in which the C-Cl is stronger than the C-Br. Perhaps related to the strong interaction between halogen anions and water molecules that cannot be correctly described by the implicit solvent model (Ayala et al., 2000).
The BDE calculated from the anionic radicals is negative (Table 1), indicating that the dissociation of the C-X bond is a thermodynamically preferred reaction in this respect (Isse et al., 2011), its reaction kinetics remains to be further investigated. The vertical electron affinity and BDE were believed jointly determine the reaction trend in the single electron transfer pathway.
3.3. Carbon-halogen bond cleavage with hydrogen radical (•H)
The possibility of e−(aq) reacting preferentially with protons is based on the following facts. First, the proton has a fully empty σ orbital that can accept an electron readily. Second, there is a natural electrostatic attraction between protons and e−(aq). Third, nuclear quantum effects may make it easier for protons to break free from the hydrogen bond networks (Rossi et al., 2016), making this reaction could occur at a higher barrier. Due to the complexity and divergent nature of hydrogen bonding frameworks’ calculations (de Almeida et al., 2021), we prefer to track the transition state (TS) between the formed hydrogen radicals and the halogenated molecules. The TS structures with one imaginary frequency are shown in the Fig. 5.
Unpaired electrons in the TS are mainly distributed in the region of carbon-halogen bond (Fig. 5), the electron distribution over the H nucleus has a sigma feature, while the electron distribution on the halogen atom appears to be related to its atomic orbital characteristics. For example, F atom has a half-filled orbital at its 2p orbital, Cl's half-filled orbital locates in its 3p orbital, and Br's half-filled orbital is at the 4p orbital and it is the only element here having a d orbital. From the perspective of the spin density distributions of F and Cl, the hybridization status of adjacent carbon seems to have no effect on spin density. However, the spin density of Br differs in the chemical environments of sp2 and sp3 carbons. Near the sp2 carbon, it is anisotropic, suggesting that it might be the p or d orbitals becoming the electron population destination. Hydrogen radicals are created by the reaction of protons with hydrated electrons. After accepting an extra electron, the hydrogen radical has the potential to give or take electrons from other molecules due to its half-filled S orbitals. The experimental results support these two pathways (Calza and Pelizzetti, 2004; Zona et al., 2008). In the former process, the electron in the hydrogen radical is transferred to the half-filled orbital of the halogen atom to form the halogen anion. In this case, the spin density of corresponding carbon, for example, is then contributed by the electron in its own half-occupied orbitals. Therefore, the spin density in the TS structure is not completely caused by the unpaired electrons of the hydrogen radical, it reflects the contribution of all the unpaired electrons during the intense change of covalent bond breaks.
Due to strong interaction with water, it is hard to calculate the enthalpy of production of individual halogen atoms accurately (Ayala et al., 2000). Therefore, electronic energy is an acceptable alternative and the change on the intrinsic reaction coordinate (IRC) step is sufficient to qualitatively and accurately reflect the energy barrier (Fig. S2). For benzene halogenates, fluorobenzene has the highest energy barrier and bromobenzene has the lowest. The same trendy is observed for halogenated methane, but the energy difference between the product of halogenated benzene and the reactants is smaller, while the energy difference between the product of halogenated methane and the reactants is larger (Fig. S2). We speculate that this is related to the deformation energy, since the product methyl radical has a planar structure while the reactant halogenated methane has a sp3 tetrahedral structure (Fig. S4) (Cole et al., 1958).
3.4. Simulation of hydrated electron migration and cleavage of carbon-halogen bond in AIMD
THE* is one of the most widely accepted models of e−(aq) in recent years (Shkrob, 2007; Kumar et al., 2015; Biswas et al., 2022), this model is considered to be the realities in bulk water and verifies many experiments (Karashima et al., 2016; Yamamoto et al., 2016). The C-F bond is generally considered to be the strongest carbon-halogen bond, and for the simplicity of the study, fluoromethane and fluorobenzene was selected as the representative pollutants for AIMD simulations. The evolution of e−(aq) with fluoromethane is presented in Fig. 6. At the beginning of the AIMD simulation, the single electron position is in the cavity between the four water molecules, proving the correctness of the initial guess (Fig. 6, 0 step). As the simulation proceeds, a vigorous expansion and contraction of the e−(aq) occurred. For example, from step 50 to step 200 (Fig. 6, 50–200 step), the electrons exhibit a partly dispersive nature. And at step 400, the extra electron undergoes a significant contraction (Fig. 6, 400 step). The induction of electron transfer by the LUMO + 1 orbital of fluoromethane can be observed from steps 600 to 2000 (Fig. 6, 600–2000 step), and by comparison with the molecular orbital of fluoromethane (Fig. S1), the electron enters the LUMO + 1 orbital of fluoromethane at the 2000 step. However, the expected C-F bond dissociation did not occur. In the subsequent simulations, electrons are continuously pulled back and forth by fluoromethane and the water cavity, demonstrating that the reduction of the hydrated electron in the real world is not absolutely controlled by thermodynamics. Under the influence of the environment, the formation of a fluoromethane anion does not necessarily lead to dissociation. The electron has a probability to return to the solution during our time scale and the same phenomenon was reported in other research (Saha, 2018; Yu et al., 2022). To ensure C-F bond dissociation, MTD was used in the AIMD simulation (Fig S3 B) (Laio and Parrinello, 2002). The free energy profile suggested that the cleavage of C-F bond in fluoromethane still has a high barrier (Fig S3 A) even with the extra electron populated (Fig S4).
In the simulation of fluorobenzene, three water molecules were added around the fluorine atom to further enhance the realism (Fig. 7). In addition to the accuracy of the initial guess was proved by the spin density, a trajectory in which C-F bond dissociation occurred without using the MTD technique was fortunately obtained among a total of 46 trajectories. Unlike fluoromethane, fluorobenzene easily traps hydrated electrons into its LUMO orbitals (Fig. 7 200 step and Fig. S1 10). Although this cleavage was partly related to the randomness, it did not violate the basic principles of quantum chemistry, we thought it could happen in the real world (Fedurco et al., 2001; Isse et al., 2011; Jiang et al., 2019). In this trajectory, the hydrogen atoms on the β carbon "attack" the α carbon after violent vibrations, resulting in the departure of halogen atoms, and the α carbon on the formed phenyl radicals is actually the previous β carbon, which may provide a new insight into the experimental results. The departing halogen formed a typical solvent layer structure after 200 steps. Due to the great uncertainty of ordinary AIMD for the success of rare event sampling, we used MTD enhanced sampling technology to study the C-F dissociation process, and unlike the AIMD trajectory, MTD observed C-F dissociation at about 50,000 steps, forming methyl radicals and halogen anion products, and methyl radicals vibrate near their own stable planar structure until the end of the simulation.
By comparing the shape of molecular orbitals and spin densities, hydrated electrons will migrate to the LUMO or LUMO + 1 orbitals of fluorobenzene or fluoromethane (Reitmeier et al., 2006). The additional electrons destroy the stable closed shell state, under the inducement of water molecules or hydrogen radical, the open shell molecules release halogen atoms, completing the conversion of organic halogens to inorganic halogen anions, which is crucial for the reduction of toxicity of environmental pollutants.
