## 3 Computational Methods  
All standard DFT calculations were performed using the Gaussian 09 program [17]. All structures were optimized and characterized using frequency analysis calculations to  
Figure 8: Isosurface maps of Fukui functions for Zr (\(f^{*}\)) in[Zr]-H and C\({}_{2}\), C\({}_{3}\) (\(f^{-}\)) in DVB-PP\({}_{1}\) and [BVPE]-PP\({}_{1}\). \begin{table}
\begin{tabular}{|c c c c c c c c c c c|} \hline  & \(\eta^{*}\) & \(\mu^{*}\) & \(\omega^{*}\) & \(N_{\text{Nu-local}}\) & \(N_{\text{Nu-local}}\) & \(f^{*}\) & \(f^{*}\) & \(f^{-}\) & \(h\) \\ \cline{5-10}  & & & & & & **Zr** & **C\({}_{2}\)** & **C\({}_{3}\)** & **Zr** & **C\({}_{2}\)** & **C\({}_{3}\)** \\ \hline
[Zr]-DVB-PP\({}_{1}\) & 4.994 & \(-\)6.838 & 5.088 & 1.064 & 0.842 & 0.040 & 0.012 & 0.166 & 0.088 & 0.011 \\
[Zr]-DVB-PP\({}_{1}\) & 4.295 & \(-\)6.637 & 5.128 & 1.619 & 0.854 & 0.031 & 0.028 & 0.167 & 0.075 & 0.039 \\
[Zr]-H & 4.945 & \(-\)7.212 & 5.260 & 0.588 & 1.298 & & 0.247 & 0.247 & 0.143 & 0.067 \\
[BVB-PP\({}_{1}\) & 8.748 & \(-\)3.271 & 0.611 & 3.311 & & 0.473 & 0.221 & 0.143 & 0.067 \\
[BVPE]-PP\({}_{1}\) & 8.145 & \(-\)3.247 & 0.647 & 3.264 & & 0.308 & 0.140 & 0.094 & 0.043 \\ \hline  & & & & & & & & & & \\ \end{tabular} \({}^{a}\) Chemical hardness (\(\eta\), in eV). \({}^{b}\) Electronic chemical potential (\(\mu\), in eV). \({}^{c}\) Global electrophilicity (\(\omega\), in eV). \({}^{d}\) Global nucleophilicity (\(N_{\text{Nu}}\), in eV). \({}^{e}\) Local electrophilicity (\(\omega_{local}\), in eV). \({}^{l}\) Local nucleophilicity (\(N_{\text{Nu-local}}\), in eV). \({}^{g}\) Nucleophilic attack Fukui function for Zr (\(f^{*}\), in e). \({}^{h}\) Electrophilic attack Fukui function for C\({}_{2}\), C\({}_{3}\) (\(f^{-}\), in e).  
\end{table} 
Table 3: The global reactivity index, local reactivity index, and Fukui function values for some molecules. be at a minimum (no virtual frequency) or a transition state (TS, with a unique virtual frequency) at the B3 LYP [18; 19]/BSI level. BSI denotes a basis set combining the SDD [20] for Zr and 6-31G (d, p) for other non-metal atoms. The pseudo-potential basis set was used for the Zr atom. To improve the energetic results, single-point energy calculations at the optimized structures were performed at the M06 [21; 22]/BSII, with solvent effects being simulated via the SMD [23] solvent model, using toluene as the solvent. BSII represents a basis set combining SDD for Zr, and 6-311++G (d, p) for other non-metal atoms. The basis set superposition errors (BSSEs) were accounted for at the same level by using the standard counterpoise method [24; 25], as implemented in Gaussian 09. The NBO charges [26] and Wiberg bond indices were obtained at the B3 LYP/BSI level. The gas phase B3 LYP/BSI harmonic frequency was used to modify the free energy by heat and entropy at 298.15 K and 1 atm pressure, respectively. BSSE-corrected free energies and enthalpies obtained at the M06 (SMD, solvent = toluene)/BSII level and free energies are discussed in the main text; the relative enthalpies are also given for reference. The reliability of the M06//B3 LYP combination is demonstrated by its successful application in explaining various transition metal catalytic reactions [27; 28; 29; 30; 31; 32; 33; 34; 35; 36]. The total energy and Cartesian coordinates of all the optimized structures are provided in the Supplementary Materials (Table S2). The global reactivity index, local reactivity index, and Fukui functions values were performed using the Multiwfn program [37; 38]. The calculations also used a mixed basis set (SDD for Zr and 6-31G(d) for other non-metal atoms). Vertical ionization potential: VIP \(=E_{N-1}-E_{N}\),  
Vertical electron affinity: \(\text{VEA}=E_{\text{N}}-E_{N-1}\), where N denotes the number of electrons for a stable system.