First Principles Studies on Boron Sites in Zeolites

: A systematic computational investigation on protonated and nonprotonated boron-containing zeolites (boralites), performed by using different periodic density functional theory approximations, is presented. Both minimum energy structures and finite temperature behavior of model boron sodalites were analyzed. All of the adopted computational schemes agree in predicting an acid site composed of a silanol Si−OH group loosely linked to a planar BO3 structure in the protonated system and a BO4 tetrahedral site in the sodium-containing zeolite. Calculated structural and vibrational properties are in line with experimental data. Comparisons of the protonated boralite site with Al and Ga zeolitic acid sites are discussed as well. Results indicate that this class of mild acid catalysts is characterized by significant framework flexibility and pronounced thermal effects due to the loosely bound acid site.


INTRODUCTION
In the first part of the study, different basis sets and density functional approximations were compared by calculating minimum energy structures of B-SOD models containing either H + or Na + as counterions. The effect of framework relaxation on the B acid site structure was also investigated by optimizing the geometry of protonated B-SOD crystals in which atoms beyond the third coordination shell were kept fixed.
In the second part of the work, the room-temperature behavior of B-sites in a zeolite framework was studied by means of the Car-Parrinello approach. 35 Trajectories characterized by different simulation parameters were performed to systematically assess the reproducibility of results on B-zeolites. In the final part of the work, calculated physicochemical properties were compared with available experimental data on different Bzeolites and with results of previous computational studies on Al-and Ga-containing sodalites. [3][4][5]

Computational Methods
Sodalite, which has unit cell stoichiometry of [Si 12 O 24 ] in its all-silica form, 36 is characterized by six-membered rings as maximum pore openings. Its unit cell is formed by two cuboctahedral -cages. 37 The -cage is the building block of several larger zeolites and makes sodalite one of the choice model systems for carrying out accurate and extensive computational studies. Hydrothermal synthesis of B-SOD has been reported; 1 however, cell parameters and unit cell content were not available. Models of the B-SOD crystal were built from the all-silica sodalite structure by replacing with a B atom one of the 12 equivalent Si atoms in the unit cell. The model crystal is characterized by a B molar fraction of 0.0833, which is of the order of those commonly found in standard boralites. 1 The unit cell contraction due to the smaller size of the B cation with respect to Si was taken into account by applying Vegard's law 1 where V Si is the unit cell volume of the all-silica phase, d B  Two systems were considered: an acid B-sodalite model, H-B-SOD, and its Na-substituted counterpart, Na-B-SOD, where the extraframework species balancing the negative charge of the framework were H + and Na + , respectively. In both models, the starting configuration was built by assuming a tetrahedral geometry for the B-site. In H-B-SOD, the proton was placed at a bonding distance from one of the oxygen atoms bonded to B, while in Na-B-SOD the Na nucleus was placed close to the center of the six-membered ring containing B.
Geometry optimizations (GOs) and first principles molecular dynamics (FPMD) simulations were performed by using the plane waves (PWs) DFT code CPMD. 38 Wave functions were expanded in plane waves (Γ point only), and norm-conserving 39 semilocal 40 pseudopotentials were adopted for the ionic corevalence electron interactions. Nonlocality up to l ) 2 was adopted for all of the atomic species with the exception of H, for which a local pseudopotential was used. The pseudopotentials, in numerical form, were obtained via the Troullier and Martins approach. 41 GOs on Na-B-SOD and H-B-SOD were carried out by using quasi-Newton methods and adopting several gradient-corrected density functional approximations. In par-ticular, the following functionals or combinations of functionals were employed: Becke 42 -Perdew 43 (BP), Becke 42 -Lee-Yang-Parr 44 (BLYP), generalized gradient approximation (GGA), 45 Perdew-Burke-Ernzerhof 46 (PBE), and Hamprecht-Cohen-Tozer-Handy (HCTH/120). 47 Convergence of the results with respect to the kinetic energy cutoff in the PW expansion was also checked.
The performance of PW versus localized basis sets was investigated as well. Localized basis sets were used in the QUICKSTEP DFT-based approach. 48,49 QUICKSTEP (QS) is an electronic structure program, available in the CP2K code, 50 which computes energy and forces through the Gaussian-plane wave (GPW) method, a hybrid scheme of Gaussian and plane wave functions. In particular, the Kohn-Sham orbitals are expanded in a linear combination of atom-centered Gaussians while the electron density is expanded in a PW basis set. In analogy with PW-based methods, in the GPW approach only valence electrons are accounted for, and the ionic core-valence electron interactions are represented by pseudopotentials. A detailed review of the method can be found in ref 49. The PBE functional and pseudopotentials of the Goedecker-Teter-Hutter (GTH) 51 type optimized for PBE 52 were adopted in the GOs performed with the GPW method. The convergence of atomic basis set size was analyzed by comparing results obtained by using a DZVP and a TZV2P basis set. 50 A 300 Ry cutoff for the PW expansion of the density was adopted. A tight convergence criterion upon the wave function was enforced during the optimization runs (convergence criterion for the electronic gradient: e 10 -7 a.u.). The convergence criterion for nuclear gradients in GOs was δ e 1 × 10 -4 a.u. as in the PW calculations.
The effect of framework relaxation on the structure of the B acid site was studied by fixing the positions of the framework atoms beyond the third shell from the B-O(H)-Si oxygen in the H-B-SOD system to the all-silica sodalite X-ray positions. Relaxation effects were also investigated on three model sodalites H[T III Si 11 O 24 ], with T III ) Al, Ga, or Fe. Normconserving pseudopotentials (d-nonlocality) were used for these atoms as well, and in the case of Fe, a restricted open-shell formalism, with Fe III in the high-spin state, was adopted.
Car-Parrinello (CP) FPMD simulations on Na-B-SOD and H-B-SOD were performed at room-temperature conditions by using the PBE functional with a 70 Ry cutoff (PBE/70). Simulations were carried out also with the BP functional at a 60 Ry cutoff (BP/60) to (i) evaluate whether and up to what extent the choice of DFT approximation/PW cutoff might affect average structural properties and (ii) allow comparison with results previously obtained by FPMD on H-Al-SOD and H-Ga-SOD with the BP/60 computational scheme. 3,4 The validity of the BP/60 approach has been already proven in several examples dealing with many aspects of framework minerals and zeolite chemistry: structure, 53-55 reactivity, [56][57][58] and behavior at high pressure 59 and temperature 60 conditions. Here its accuracy is further tested by cross-checking the BP/60 results on H-Al-SOD 3 with those obtained on the same system with PBE/70. The equations of motion 35,38 were integrated by using a time step of 0.121 fs (5 a.u.) and adopting an inertia parameter of 500 a.u. for the electronic coefficients. Such values have been shown to ensure adiabatic separation of ionic and electronic degrees of freedom and have been adopted in all FPMD runs here presented. In the CP approach both ionic and electronic degrees of freedom are propagated by equations of motion derived from an extended Lagrangian 35 in which an inertia parameter (fictitious mass) is associated to the wave Both BP/60 and PBE/70 simulations were performed in the NVT ensemble with Nosé-Hoover chain thermostats and a target ionic temperature of 300 K. Different starting conditions, equilibration times, and thermostat parameters were chosen ( Table 1). Equilibrations were performed by alternating 0.1 ps bins in the NVE ensemble by 0.1 ps bins with velocity rescaling (target temperature rising from 100 to 300 K, with tolerance from 32 to 37 K). For some simulations, an additional equilibration period in the NVT ensemble at 300 K was carried out. For the H-B-SOD system, data were collected for periods longer than those typical of FPMD simulations to guarantee a converged description of the B-site properties at equilibrium conditions.

Optimized Structures of H-B-SOD and Na-B-SOD.
The boron site is found to be trigonal in H-B-SOD and tetrahedral in Na-B-SOD, independently of the adopted DFT functional and basis set type and size and in line with both available experimental data and results of previous calculations on model clusters. Relevant structural parameters are reported, for comparison, in Tables 2-5. Graphical representations of the optimized structures are shown in Figure 1.
In Na-B-SOD, the Na cation, which is coordinated to six oxygen atoms of the six-membered ring containing B, is found at a distance of 3.0 Å from B. The tetrahedral geometry of the B-site slightly deviates from ideality; in particular longer B-O bond distances are found for the two oxygen atoms coordinated to Na, namely, O* and O # (see Figure 1 for atom labels).
The H-B-SOD system is characterized by a silanol group Si*-O*H and a trigonal B-site in a nearly planar BO 3 geometry. These groups are loosely interacting, as indicated by the large separation between the silanol oxygen O* and the B atom (>2.5 Å). The B-O bonds in the BO 3 structure are shorter than those of the tetrahedral B-site in Na-B-SOD. Interatomic separations and B-O*-Si* angles are well-converged with respect to basis set size (Tables 2 and 3). In particular, the largest variation in the geometrical parameters of the B-site is found for the B-O* nonbonding distance in H-B-SOD and amounts to 0.7% in passing from 60 to 90 Ry with BP and to 0.5% from 60 to 110 Ry with PBE. The corresponding B-O*-Si* angle variations amount to 1.6% and 0.6%, respectively. Cutoff-dependent changes affecting the geometry of the B-site in Na-B-SOD are even smaller.
The dependency of the geometry on the DFT approximation and on the type of basis set is reported in Tables 4 and 5. It can be noticed how all of the adopted functionals provide structural descriptions of the B-site remarkably in agreement among each other. On this basis, the effect of localized versus plane wave basis sets on the calculated structural parameters was studied with the PBE functional only.
For both H-and Na-B-SOD, the GPW approach provides results well in line with PW calculations, confirming a planar trigonal and a tetrahedral B-site geometry, respectively. It is also worth pointing out that the nonbonding distances calculated with the GPW method are very close to the PW ones, thus demonstrating reproducibility of calculated structural parameters between the two approaches. In particular, the GPW B-O* and B-Si* distances in H-B-SOD (Table 4) indicate that the large separation between B and the silanol group found in the model acid B-zeolite is not dependent on the locality of the basis set. Basis set size slightly affects geometrical parameters; in particular the GPW/TZV2P geometries are closer to those obtained with PW at 70 Ry.

Framework Relaxation in H-B-SOD.
The minimum energy structure in H-B-SOD is characterized by a trigonal boron site weakly interacting with a silanol group, in keeping with previous calculations performed on cluster models of proton-exchanged B-zeolites. 29,[31][32][33][34] The nonbonding B-O* distance (2.57 Å with PBE/70) is however significantly larger than the calculated values reported in the literature, which range from 2.02 to 2.37 Å depending, beside the electronic structure treatment, on the size and extent of relaxation of the adopted cluster model. In this respect, it was shown that converged structural properties of zeolitic acid sites required cluster models built by positioning around O* at least 3-4 shells of atoms, the first two of which must be fully relaxed. 29 However, even though framework relaxation appeared to greatly affect the optimized geometry, extrapolation of these observations to the limiting case of a fully relaxed periodic crystal structure has not been performed yet. To shed light on this issue and establish whether differences among predicted B-O* distances should be attributed to geometrical constraints or to the calculation approach, GO results on the fully relaxed periodic H-B-SOD model were compared to those obtained on an identical model in which only the first three shells of atoms around O* were allowed to relax. The starting configuration for the constrained GOs was the same as that adopted in the full GOs on H-B-SOD, with the outer-shell atoms kept fixed to the SOD crystallographic positions. Results obtained with the PBE functional using both PW and GPW approaches are reported in Table 6. Comparison with data in Table 4 shows that the predicted B-site geometry is significantly altered by the imposed geometrical constraint, with a more pronounced effect in the case of localized basis sets. Both constrained optimizations led to structures characterized by significantly shorter B-O* and B-Si* separations and a larger Si*-O*-B angle with respect to the fully relaxed structure. Such a picture was also confirmed by constrained GOs performed with the other functionals. Interestingly, B-O* distances (2.06 and 1.96 Å for PW and GPW, respectively) compare well with the values of 2.020 and a tequil is the total equilibration time (in ps), tc is the data collection time (in ps), and ωi and ωe are the characteristic frequencies of the Nosé-Hoover thermostats on the ions and on the electrons, respectively (in cm -1 ). The Nosé-Hoover target ionic temperature and electronic kinetic energy were 300 K and 0.0035 a.u., respectively.
2.023 Å obtained recently from B3LYP/6-311+g(d,p) calculations on cluster models of B-MCM-22 34 and B-ZSM5 31 where a similar extent of framework atom relaxation was adopted.
The effect of the applied geometrical constraints on the resulting minimum energy structure is pictured in Figure 2, where the partially frozen PBE/70-optimized geometry is   overlapped with the corresponding fully relaxed structure. The fully optimized system is characterized by significant distortions of the four-and six-membered rings not containing B. For instance, the Si-Si separations in the four-membered ring, which were both equal to 4.39 Å in the constrained system, amounted to 4.40 and 4.46 Å in the relaxed structure. Larger relaxation effects are detected in a close six-membered ring. The Si-Si distances characterizing the circular ring opening (∼6.25 Å radius) are split to 6.42 and 5.46 Å in the relaxed system, bringing about an elliptical deformation of the ring. While the largest and most significant structural differences between the fully and the partially relaxed systems are detected in proximity to the B acid site, these results suggest that B brings about relevant long-range structural perturbations as well. In particular, larger deformations are detected in the more flexible six-membered rings. Relaxation effects also have been detected in other H-T III -SOD systems, where T III represents typical trivalent framework cations, namely, Al, Fe, or Ga. Relevant parameters of the different acid sites optimized at the PBE/70 level of theory, for the partially and fully relaxed H-T III -SOD framework, are reported in Table 7. Also in these cases, the framework relaxation has important effects on the structure of the acid site. By allowing full relaxation a significant lengthening of all of the T III -O distances and of the Si*-O* distance is observed, while the O*-H bond is slightly shortened. As in the case of the B-site, the Si*-O* and T III -O* bond lengths obtained with the partially frozen framework are in line with the corresponding results on cluster models. 31,34 The energetic stabilization associated with framework relaxation in the H-T III -SOD systems may be estimated from the energy differences between the fully and partially relaxed structures. The framework relaxation energy per unit cell amounts to 44.6, 52.1, 65.8, and 57.2 kcal/mol for H-B-SOD, H-Al-SOD, H-Ga-SOD, and H-Fe-SOD, respectively, indicating that full relaxation of the crystal structure has a large stabilizing effect on the acid site structure.

Equilibration Process in Acid B Sodalite.
Several simulations characterized by different parameters have been carried out on the H-B-SOD system ( Table 1). The first two simulations, namely, BP/60(1) and PBE/70(1), in which different electronic structure treatments, starting configurations, and equilibration procedures were adopted, provided a similar average structure of the zeolite matrix. However, a significant difference was found for the mean values of the nonbonding B-O* distance, which amounted to 2.59 and 2.71 Å in BP/ 60(1) and PBE/70(1), respectively. This difference, affecting the loosely interacting atoms, might depend on the different electronic structure treatment but could also be related to the simulation length. The equilibration period and/or the data collection time might have been insufficient to properly sample the system under study and specifically the weak interaction between the acid site atoms O* and B. That this could be the case was suggested by the two PBE/70 and BP/60 simulations a Results obtained by adopting the GPW method with the PBE functional as a function of the basis set size are also reported. Distances in are given in angstroms, and angles are given in degrees. a Results obtained by adopting the GPW method with the PBE functional as a function of the basis set size are also reported. Distances are given in angstroms, and angles are given in degrees. on Na-B-SOD that converged to the same structural description in simulation times on the order of those adopted in standard FPMD runs of zeolites (see section 3.4.2). Therefore, other simulations characterized by longer equilibration and sampling times were performed on H-B-SOD. Care was taken to check whether these simulations achieved adiabatic sampling. Indeed, the drifts of the fictitious electronic kinetic energy and of the conserved quantity in the longest trajectory (64 ps) were 5 × 10 -4 and 4.5 × 10 -5 a.u., respectively, thus confirming adiabatic sampling. As a further check, a BO simulation was carried out using the same initial configuration of the PBE/70(2) run.
Remarkably, the longer simulations converged to a closer description of the B acid site, which is characterized by a B-O* equilibrium separation of ∼2.7 Å and large thermal fluctuations. The differences in the results of the two shortest simulations are therefore not related to the choice of DFT approximation, basis set size, or thermostats but rather to the time required to equilibrate the weakly interacting units composing the B acid site, namely, the silanol group and the trigonal BO 3 moiety. By taking into account that the presence of the B-site perturbs the zeolite framework, a long thermalization period is needed by the zeolite lattice to reorganize and find the optimal structural arrangement (i.e., the equilibrium geometry) for simulations at ambient temperature conditions.

Finite Temperature Behavior of Boron Sites. 3.4.1. Protonated Site.
In view of the previous results, the behavior of the B acid site will be thoroughly discussed only in the case of the longest simulation (that is, PBE/70(2) in Table 1). However, comparison with results obtained with other DFT approximations and PW cutoffs, as well as with different sampling approaches, will be shortly presented. Indeed all of the trajectories longer than ∼20 ps, or with a long equilibration period (more than 10 ps), were found to provide well-converged results. Figure 3 shows the B-O running coordination numbers N(r) (top) and the radial distribution functions g(r) (bottom) calculated for the H-B-SOD system from FPMD simulations with different functionals and PW cutoffs (namely, PBE/70(2) and BP/60(2)). Results from BO simulations at the PBE/70 level are also shown for comparison.
The similarity of the curves indicates a very good agreement between the simulations. Moreover, the structural properties obtained from the two CP trajectories are consistent with the BO simulation as well, thus confirming that both the time step and the inertia parameter here selected for the CP equations (i.e., 5 and 500 a.u., respectively) are appropriate for the description of these systems.
The B-O g(r)'s are characterized by two maxima in the nearest neighbor region. The first peak is found at 1.40 Å; correspondingly, N(r) increases from 0 to 3 for distances comprised between 1.2 and 1.4 Å, remaining equal to 3 up to 2.4 Å. This is a clear indication of the fact that the B atom retains the trigonal geometry at room temperature. Relevant mean geometrical parameters along with their relative standard deviations are reported in Table 8. The distances between B and O1, O # , and O2 do not differ significantly (less than 0.2%) from those of the minimum energy structure, indicating that the temperature has a negligible effect on the internal geometry of the BO 3 moiety. Inspection of these values also reveals that both at 0 K and at room temperature the B-O2 distance is shorter than the B-O1 and B-O # ones. This can be due to the fact that, while the B-O1 and B-O # bonds belong to an unbroken four-membered ring, the B-O2 bond is part of the relaxed four-membered ring, which is broken due to the formation of the silanol group. Also the average O-B-O angles in the BO 3 site are practically identical to the minimum energy structure values and correspond to a planar trigonal geometry.
The second peak in the B-O g(r) is located at ∼2.7 Å, while N(r) reaches a value of 4 at 3.0 Å. This peak is much broader than the first one and corresponds to the B-O* contacts. Here, the thermal effect on the structural arrangement of the weakly interacting groups Si*-O*H and BO 3 is evident. In fact, by comparing the B-O* average distance (Table 8)    corresponding value obtained from GOs ( Table 2) an increase of about 0.15 Å (5%) is found. In addition, the trigonal B-O bond distances are characterized by standard deviations significantly smaller than that of the B-O* mean distance, indicating a large thermal fluctuation of this latter contact. Also the mean B-O*-Si* angle and the B-Si* distance are different from the minimum energy values; however their change in percent (2-3%) is lower than that found for the B-O* separation. These findings suggest that the silanol oxygen O* is the most mobile framework atom. Indeed, being less constrained by the framework, O* can explore at finite temperature a wide range of configurations characterized by different separations from the trigonal BO 3 unit.
These results suggested further analysis of the behavior of the silanol group during the longest run. Actually, after 11 ps, a structural rearrangement of the protonated B-site was observed. Such a reorientation, which was accompanied by concerted rigid rotations of the BO 3 unit and of SiO 4 tetrahedra connected to the site, led to the geometrical arrangement represented in Figure  4. Such a transition, detected in the 64 ps long simulation, brings the O*H silanol group, originally pointing to the center of a six-membered ring, in an adjacent six-membered ring. This reorientation of the silanol group, which occurs in less than 1 ps, does not affect the average structural properties of the system. Indeed, as clearly demonstrated by Figure 4, the final geometrical arrangement of the protonated B-SOD system may be considered as the mirror image of the starting structure, obtained by rigid unit rotations. 61 3.4.2. Unprotonated Site. Information on the average structure of the B-site in Na-B-SOD can be extracted from the B-O g(r) and N(r), shown in Figure 5, while average bond distances and angles are reported in Table 9. A first sharp peak in the g(r) is located at 1.49 Å, while the N(r) integrates to 4 at 1.65 Å for both BP/60 and PBE/70 simulations. This is a clear indication that in the Na-B-SOD system the B atom remains tetracoordinated also at room temperature. The thermally averaged structure of the B-site in Na-B-SOD is characterized by O-B-O angles typical of a tetrahedron and is very close to the minimum energy geometry. As found in the optimized structure, the longest B-O bond involves the O* atom, which is the one closest to the Na extraframework cation. The parameter affected by the largest oscillations around the  a Distances are given in angstroms, and angles are given in degrees. a Distances are given in angstroms, and angles are given in degrees. 〈O-B-O〉 indicates the average of the three angles in the trigonal unit BO3. equilibrium value is the O*-Na separation, as expected on the basis of the higher mobility of the extraframework species Na + . The calculated standard deviations of the four B-O bonds are 1 order of magnitude smaller than that found for the B-O* separation in H-B-SOD. Also, the Si*-B distance is shorter and its oscillation smaller with respect to the H-B-SOD system. This indicates a lower flexibility of a zeolite framework containing a tetrahedral B-site with respect to the case of the silanol-BO 3 -containing zeolite. Differences in the framework structure of the protonated and nonprotonated model boralites may be detected by inspection of the g(r)'s calculated for the Si-Si and B-Si atoms ( Figure  6). The Si-Si g(r)'s show a first maximum peaked at 3.18 Å, while the positions of the second peak are 4.4 and 4.3 Å for H-B-SOD and Na-B-SOD, respectively. As expected, differences are even more evident in the B-Si g(r)'s: The peak corresponding to the first neighbors is found at shorter distances in H-B-SOD because the B-O bonds in the BO 3 units are shorter than those in the BO 4 tetrahedron. Moreover, also the second and third peaks are displaced at shorter distances in H-B-SOD due to the breaking of the four-membered ring containing the B acid site. These results suggest therefore that the presence of different local geometries of the B-sites in Na-B-SOD and H-B-SOD affects the framework structure also at long range.

TABLE 7: Geometrical Parameters of the T III Site in H-T III -SOD (T III ) Al, Ga, or Fe) from Constrained and Fully Relaxed Geometry Optimizations with the PW Approach (PBE/70) a
3.5. Acid Site in B, Al, and Ga Zeolites. Significant flexibility of protonated boralites, due to the nonbonding interactions characterizing the B acid site, emerged from the previous discussion. In this respect, comparison with results obtained with the same approach on other zeolitic acid sites may be of relevance to gain further insight on this issue.
Al and Ga protonated model sodalites were studied by FPMD at the BP/60 level. [3][4][5] As confirmed by a PBE/70 simulation on the H-Al-SOD system (Figure 7), structural properties predicted by the BP/60 scheme are well in line with the PBE/ 70 description, thus allowing comparison between results obtained at finite temperature on T III -substituted zeolites (T III ) Al, Ga, or B).
The T III -O, Si-Si, and T III -Si g(r)'s for the B, Al, and Ga acid sodalites are shown in Figures 7 and 8 Figure 9.
In the system containing the strongest acid site, namely, H-Al-SOD, the Al-O g(r) is characterized by a first maximum with a peak at 1.77 Å and a shoulder at 1.95 Å. The peak    corresponds to the most probable bond distances of Al with the three unprotonated oxygens. The shoulder represents the Al-O* bond, which is elongated with respect to the other Al-O distances, but still is part of a tetrahedral AlO 4 unit. Such a tetrahedron is however distorted as indicated by the O-Al-O angle distribution, which is characterized by two maxima at 102°a nd 114°. 3 In the H-Ga-SOD case, the Ga-O g(r) first maximum is split into two peaks: The first one, which integrates to 3, is found at 1.85 Å, while the second is located at 2.20 Å and corresponds to the Ga-O* contacts. Such a value suggested that the Ga-O* bond is weaker than the Al-O* one and that the Ga acid site has silanol character. 4 Also, the O-Ga-O angle distribution indicates that the geometry of the Ga site is significantly distorted from an ideal tetrahedron and might be considered close to a trigonal arrangement. However, the limiting case of a planar trigonal geometry of the T III site and a nonbonded silanol group is achieved only in the case of the B acid site, as indicated by the O-B-O angle distribution. The peaks at 90°and 120°correspond respectively to the three angles involving O* and to the three angles of the BO 3 unit.
By comparing the distributions of the T III -O* distances, another interesting feature emerges. In the case of H-B-SOD, the full width at half-height (0.6 Å) is about twice that found for the Al-and Ga-containing systems. Such a result is a consequence of the weak interaction between O* and the BO 3 group and suggests that the high flexibility of the acid site, detected in the H-B-SOD simulations, is indeed a peculiarity of protonated boralites. In addition, this large local flexibility is also at the basis of significant structural deformations of the framework, which, as shown by the T III -Si g(r)'s, are much larger than in the Al and Ga zeolites. Moreover, the average T-O-T angle distribution is broader and tails toward larger values of the angles in the case of H-B-SOD.
The O-Si-O angle distributions calculated for all 11 Si sites show a single maximum at 109.5°. In the AlO 4 unit, the splitting of the O-Al-O distribution in two peaks indicates that the tetrahedral geometry is distorted to a trigonal pyramid structure. Indeed, while Si occupies the center of the SiO 4 tetrahedron, Al is found closer to the three unprotonated oxygens ( Figure  10). Such oxygens may be considered as forming a pyramid base, while the fourth and more distant protonated oxygen O* may be considered as the pyramid vertex. The pyramidal character becomes more pronounced in the case of GaO 4 , as indicated by a comparison of the T III -O* distance distributions and by the broadening and further splitting of the O-Ga-O distribution. However, the transition to a trigonal pyramid is fully accomplished only in the case of B, because the T III cation lies in the same plane of the three oxygen atoms forming the basis of the pyramid (Figure 10).   3.6. Simulated Vibrational Spectra. IR spectroscopy has been widely used for the characterization of boralites. Most studies were performed on B-MFI-type zeolites; therefore only a qualitative comparison with the present simulation results is possible.
In IR studies on evacuated B-MFI, a band at 1380 cm -1 was assigned to three-coordinate B, while bands in the 860-920 cm -1 range and at 670 cm -1 were assigned to B-O-Si symmetric stretching and bending modes, respectively. 14 A band at 700 cm -1 and a shoulder at 970 cm -1 were detected as well. 15 Further studies reported the presence of two bands at 1380 and 1405 cm -1 , attributed to trigonal B occupying different crystal sites. 16 Moreover, other authors reported bands at 1390 and 1150 cm -1 , ascribed to the simultaneous presence of trigonal and tetrahedral B-sites, respectively. 17 The latter band may be partially hidden by the crystalline silica matrix modes.
B-CHA has been recently synthesized and characterized by IR spectroscopy. 22 Upon calcination two main changes are reported: A typical signal of a trigonal BO 3 site appears at 1390 cm -1 , while bands at 800 and 890 cm -1 found in as-prepared samples are blue-shifted to 810 and 912 cm -1 , respectively. The 1390 cm -1 band disappears when B-CHA is contacted with water. 22 In general, while the presence of BO 3 is easily detectable from the ∼1390 cm -1 band, no definite attribution of a specific band to BO 4 modes has been possible by IR spectroscopy in evacuated boralites. This might also be due to the possible simultaneous presence of both kinds of sites in the same sample. 1 Moreover, preparation, heating, and dehydration procedures undergone by the samples may affect band position and intensity. 1,12,19 Vibrational spectra were obtained from the Fourier transform (FT) of T-O-T oscillation autocorrelation functions 62 calculated from the longest trajectory of each system. Partial spectra due to bending, symmetric, and asymmetric stretching modes of the T-O-T brigdes were also calculated by the FT of the relative correlation functions. Such an analysis allowed us to single out the contributions to the total spectrum of specific atoms or groups of atoms, thus enabling the assignments to the bands.
Due to the use of a fictitious inertia parameter in the CP equations, simulated vibrational frequencies are systematically red-shifted. Correction of this friction effect by application of a proper scaling factor is a well-established procedure, already applied to the study of other zeolitic systems. 58 In the present case, the scaling factor has been determined from the ratio of the H 2 O asymmetric stretching frequencies calculated from 30 ps long BO (3690 cm -1 ) and CP (3532 cm -1 ) trajectories of an isolated water molecule at the PBE/70 level.
Calculated vibrational spectra of Na-B-SOD and H-B-SOD are reported in Figure 11. The most significant differences between the two spectra are found in the 3600 and 1300 cm -1 regions. A peak at 3643 cm -1 is present only in the H-B-SOD spectrum and is due to the silanol hydroxyl group. The double band in the H-B-SOD spectrum, with peaks at 1338 and 1290 cm -1 , is missing in the Na-B-SOD one. Such a band has to be attributed to the trigonal B-site. The splitting into two peaks derives from the different B-O bond lengths in H-B-SOD (Table 8). Actually, by calculating the distinct contributions of the three B-O-Si subsystems to the total spectrum, it is found that the higher wavenumber component is due to the B-O-Si structure characterized by the shortest B-O bond.
In the T-O-T asymmetric stretching region, the Na-B-SOD spectrum presents a shoulder at 1184 cm -1 that is missing in the H-B-SOD system. This feature should be ascribed to B in tetrahedral sites and in particular to the B-O-Si asymmetric stretching component. At lower wavenumbers, a signal at 910 cm -1 is present in both systems and is associated with the symmetric B-O-Si stretching modes. In the so-called silica window region (800-1000 cm -1 ) of the Na-B-SOD spectrum, other bands are present with peaks at 850, 944, and 974 cm -1 , while a peak at 860 cm -1 is found in the case of H-B-SOD. Such peaks are associated with B-O-Si symmetric stretching modes. The presence of a multiple band is due to the different interactions of Na with the oxygen atoms of the BO 4 tetrahedron. Indeed only two of them are coordinated to Na with different distances (Table 9). It should be pointed out that, from the present analysis, such a Na-induced splitting seems to be more pronounced for the symmetric than for the asymmetric B-O-Si stretching modes.
In both systems, the B-O-Si bending modes give rise to a complex band in the 600-700 cm -1 region. A doublet can be identified in both H-B-SOD (678 and 702 cm -1 ) and Na-B-SOD (657 and 693 cm -1 ) spectra and associated with the 670 and 700 cm -1 bands experimentally detected.
Concerning the OH stretching mode, the one calculated for H-B-SOD is blue-shifted by 85 cm -1 with respect to the corresponding value found in H-Al-SOD (3558 cm -1 ). This difference can be considered a spectroscopic index of the different strengths of acid sites in Al-and B-zeolites. 1,2 It may be relevant to point out that, as opposed to both boroncontaining models, in the H-Al-SOD system no band is found either in the silica window or in the ∼1300 cm -1 region. Therefore, in the commonly occurring situation of the simultaneous presence of B and Al, framework B-sites can be experimentally detected by bands in the silica window, while the signature of trigonal B is found in the 1300 cm -1 region. However, the shoulder at 1184 cm -1 calculated for the tetrahedral B-site in Na-B-SOD may overlap with a similar signal calculated at 1179 cm -1 in the case of H-Al-SOD.

Summary and Conclusions
Extensive simulations on boron-containing sodalites have been performed, both at 0 K and at room temperature, by using periodic DFT level approaches. These studies allowed us to highlight a strict correlation between the framework flexibility and the local geometry of the B acid site. On the basis of the results presented here it emerges that, to achieve a proper structural description of these systems, the use of a computational scheme allowing the full relaxation of the zeolite framework is recommended. As opposed to Na + -exchanged B-SOD systems, characterized by tetrahedral B-sites, anhydrous protonated boralites are characterized by BO 4 polyhedra with a pyramidal geometry, where the base is a BO 3 planar structure and the vertex is an oxygen bearing a proton. Such a picture is confirmed by all of the different computational schemes here adopted. Comparison among T III -containing acid sites has shown a progressive transition of the T III O 4 polyhedron from a tetrahedral to a pyramidal structure. Such a geometrical transition implies an increasing silanol character in the Al f Ga f B sequence, which parallels the detected acid strength variation. 1,2 Calculated vibrational spectra, in line with experimental observations, allowed us to highlight features typical of both tetrahedral and pyramidal B-sites. Moreover, analysis of the finite temperature behavior of the loosely interacting BO 3 and silanol groups in the protonated system indicated a thermally induced B-O* distance elongation, from 2.57 Å at 0 K to 2.7 Å at 300 K, due to the silanol group mobility. Last but not least, the presence of a protonated B-site in a SiO 4 tetrahedra network induces a relevant framework flexibility, which may also allow long-range structural rearrangements through concerted rotations of rigid polyhedral units. Such a flexibility may be responsible for the peculiar behavior of acid boralites toward basic molecules and might also be related to the facile nature of the de-boronation process.