Boltzmann Populations of the Fluxional Be 6 B 11 (cid:1) and Chiral Be 4 B 8 Clusters at Finite Temperatures Computed by DFT and Statistical Thermodynamics

Total energy computations using density functional theory are typically carried out at a zero temperature; thus, entropic and thermic contributions to the total energy are neglected, even though functional materials work at finite temperatures. This book chapter investigates the Boltzmann populations of the fluxional Be 6 B 11 (cid:1) and chiral Be 4 B 8 isomers at finite temperature estimated within the framework of density functional theory, CCSD(T), and statistical thermodynamics. A couple of steps are taken into account to compute the Boltzmann populations. First, to identify a list of all possible low-energy chiral and achiral structures, an exhaustive and efficient exploration of the potential/free energy surfaces is carried out using a multi-level and multi-step global hybrid genetic algorithm search coupled with Gaussian code. Second, the thermal or so-called Boltzmann populations were computed in the framework of statistical thermodynamics for temperatures ranging from 20 to 1500 K at DFT and CCSD(T) theoretical levels. The results show the effects of temperature on the distribution of isomers define the putative global minimum at finite temperature due to the minimization of the Gibbs free energy and maximization of entropy. Additionally, we found that the fluxional Be 6 B 11 cluster is strongly dominant at hot temperatures, whereas the chiral Be 4 B 8 cluster is dominant at room temperature. The methodology and results show the thermal effects in the relative population hence molecular properties. procyanidin C 4 -C 8 and C 4 C 6 interflavan bonds reaction superoxide anion radical and H 2 O 2 forming NADH-oxidase flavoenzyme.


Introduction
Boron is the smallest and lightest semi-metal atom [1,2] and a neighbor of carbon in the periodic table. Moreover, it has high ionization energy of 344.2 kJ/mol [3], and an affinity for oxygen atoms, which is the basis of borates [3,4]. In recent years, the pure boron clusters, the metal, and non-metal doped boron clusters, have attracted considerable attention [1,[5][6][7][8][9][10][11][12][13] due to their unpredictable chemistry [14,15] and high potential to form novel structures [16]. The potential of boron atoms to form stable molecular networks [17] lies in the fact that they have three valence electrons and four available orbitals, which implies they are electron-deficient. Boron electron deficiency gives origin to a vast number of allotropic forms and uncommon geometries [6,16] such as nanotubes [13,18], borospherenes [19], borophene [16], cages [13,20], planar [21], quasi planar [22], rings [23,24], chiral [22, [25][26][27][28], boron-based helix clusters [25,29], and fluxional boron clusters [10,[29][30][31][32] that have recently attracted the interest of experimental and theoretical researchers. Aromaticity, antiaromaticity, and conflicting aromaticity dominate the chemical bonding in boron-based clusters [25,[33][34][35]. The two most-used indices for quantifying aromaticity are the harmonic oscillator model of aromaticity, based on the geometric structure, and the nucleus-independent chemical shift, based on the magnetic response. Aromaticity is not observable, cannot be directly measured [36], and correlates with electronic delocalization [37]. The fluxionality in boron and boron-doped-based molecular systems is highly relevant in terms of its catalytic activity [38] and is due to electronic delocalization [25]. Moreover, in boron-based nanoscale rotors, electronic localization o delocalization contributes significantly to stability, magnetic properties, and chemical reactivity [36], and it is a function of the atomic structure, size, bonding, charge, and temperature [39]. So far, doping a boron cluster with non-metals [40] dramatically affects its structure, stability, and reactivity, like shut-down the fluxionality of the boron-doped anion B 19 . In contrast, doping a boron cluster with metals [7,9,24,[41][42][43] like beryllium-doped boron clusters, exhibit remarkable properties such as fluxionality [16,29,32,[44][45][46], aromaticity [29,47], and characteristics similar to borophene [1]. Furthermore, previous theoretical studies showed that the boron fullerenes B 60 and B 80 can be stabilized by surrounding the boron clusters with beryllium atoms [48,49], which effectively compensates for boron electronic deficiency [49]. These effects make beryllium-doped boron clusters interesting, joined with the fact, nowadays, dynamic structural fluxionality in boron nanoclusters is a topic of interest in nanotechnology [19,50]. Particularly attractive are the chiral helices Be 6 B 11 À , reported by Guo [29], and Buelna-Garcia et al. [39] as one of the low-lying and fluxional isomers. Later, a chemical bonding and mechanism of formation study of the beryllium-doped boron chiral cluster Be 6 B 10 À2 and coaxial triple-layered anionic Be 6 B 11 sandwich structures were reported [25,46]. In these structures, the chirality arises due to the formation of a boron helix. Particularly, the chirality of nanoclusters has attracted attention due to their chiroptical properties, potential application in efficient chiral discrimination [51,52], nonlinear optics [53] and chiral materials with interesting properties [22,54], and of course, not to mention that chiral structures play a decisive role in biological activity [55]. Previous theoretical studies joint with experimental photoelectron spectroscopy reported the first pure boron chiral B 30 structure as the putative global minimum [22] at T = 0. In these pair of planar enantiomers, the chirality arises due to the hexagonal hole and its position. In the past years, the lowest energy structures of the B 39 borospherene were reported as chiral due to their hexagonal and pentagonal holes [26]. Similarly, the B 44 cluster was reported as a chiral structure due to its nonagonal holes [28]. That is, in these clusters, holes in the structure cause chirality. So far, the chirality depends on the geometry; In contrast, fluxionality strongly depends on temperature. A boron molecular Wankel motor [56][57][58] and sub nanoscale tank treads have been reported [59,60]; however, the temperature have not been considered. Nevertheless, most theoretical density functional studies assume that the temperature is zero and neglect temperature-dependent and entropic contributions; consequently, their finite used. The results and discussion are presented in Section 3. We discuss the effect of the symmetry in the energetic ordering and clarify the origin of the 0.41 kcal/mol difference in energy between two structures with symmetries C 2 and C 1 appear when we compute the Gibbs free energy. A comparison among the energies computed at a single point CCSD(T) against the DFT levels of theory and the T 1 diagnostic is presented. Conclusions are given in Section 4.

Global minimum search
Despite advances in computing power, the minimum global search in molecular and atomic clusters remains a complicated task due to several factors. The exploration should be systematic and unbiased [68,105]; a molecule's degrees of freedom increase with the number of atoms [68,[106][107][108][109]; a molecule composed of N number of atoms possesses 3 N degrees of freedom (i.e., a linear molecule has three degrees of translation, two of rotation, and [3 N-6] of vibrational modes); and, as a consequence, the potential/free energy surface depends on a large number of variables. The number of local minima increases exponentially as a function of the number of atoms in the molecule. Moreover, the total energy computation requires a quantum mechanical methodology to produce a realistic value for energy. In addition to that, there should be many initial structures. It is essential to sample a large region of the configuration space to ensure that we are not missing structures, making an incomplete sampling of the configurational space and introducing a significant problem to calculating the thermodynamic properties [64]. A complete sampling of the potential/free energy surface is impossible, but a systematic exploration of the potential energy surface is extremely useful. Although searching for a global minimum in molecular systems is challenging, the design and use of algorithms dedicated to the search for global minima, such as simulated annealing, [110][111][112][113][114][115] kick method [116,117], genetic algorithms [118][119][120], Gradient Embedded Genetic Algorithm [121][122][123] and basin hopping [124,125], has been accomplished over the years. In the past few years, one of us designed and employed genetic algorithms [12,13,29,39,98,99,104,126,127] and kick methodology [101,[127][128][129][130][131][132][133] coupled with density functional theory to explore atomic and molecular clusters' potential energy surfaces. They have led us to solve the minimum global search in a targeted way. In this chapter, our computational procedure to elucidate the low-energy structures employs a recently developed natureinspired hybrid strategy that combines a Cuckoo search [134] and genetic algorithms coupled to density functional theory that has been implemented in the GALGOSON code v1.0. Nature-inspired metaheuristic algorithms have been applied in almost all areas of science, engineering, and industry, work remarkably efficiently, and have many advantages over deterministic methods [135]. GALGOSON systematically and efficiently explores potential/free energy surfaces (PES/FES) of the atomic clusters to find the minimum energy structure. The methodology consists of a three-step search strategy where, in the first and second steps, we explore the PES, and in the third step, we explore the FES. First, the code builds a generation of random initial structures with an initial population of two hundred individuals per atom in the Be-B clusters using a kick methodology. The process to make 1D, 2D, and 3D structures is similar to that used in previous work [12] and are restricted by two conditions [12] that can be summarized as follows: (a) All the atoms are confined inside a sphere with a radius determined by adding all atoms' covalent radii and multiplied by a factor established by the user, typically 0.9. (b) The bond length between any two atoms is the sum of their covalent radii, modulated by a scale factor established by the user, typically close to 1.0; this allows us to compress/expand the bond length. These conditions avoid the high-energy local minima generated by poorly connected structures (too compact/ loose). Then, structures are optimized at the PBE0/3-21G level of theory employing Gaussian 09 code. As the second step, all energy structures lying in the energy range of 20 kcal/mol were re-optimized at the PBE0-GD3/LANL2DZ level of theory and joints with previously reported global minimum structures. Those structures comprised the initial population for the genetic algorithm. The optimization in this stage was at the PBE0-D3/LANL2DZ level of theory. The criterion to stop the generation is if the lowest energy structure persists for 10 generations. In the third step, structures lying in 10 kcal/mol found in the previous step comprised the initial population for the genetic algorithm that uses Gibbs free energy extracted from the local optimizations at the PBE0-D3/def2-TZVP, taking into account the zero-point energy (ZPE) corrections. The criterion to stop is similar to that used in the previous stage. In the final step, the lowest energy structures are evaluated at a single point energy at the CCSD(T)/def2-TZVP// PBE0-D3/def2-TZVP level of theory. All the calculations were done employing the Gaussian 09 code [136].

Thermochemistry properties
All the information about a quantum system is contained in the wave function; similarly, the partition function provides all the information need to compute the thermodynamic properties and it indicates the states accessible to the system at temperature T. Previous theoretical studies used the partition function to compute temperature-dependent entropic contributions [137] on [Fe(pmea)(NCS) 2 ] complex, infrared spectroscopy on anionic Be 6 B 11 cluster [39], and rate constant [100]. In this work, the thermodynamic functions are calculated using the temperaturedependent partition function Q shown in Eq. (1).
In Eq. (1), the g i is the degeneracy or multiplicity, using degeneracy numbers is equivalent to take into account all degenerate states and the sum runs overall energy levels, and k B is the Boltzmann constant, T is the temperature and ΔE i is the total energy of a molecule [100,138]. At high temperatures, all thermal states are accessible due to the term ÀΔE i =k B T tends to zero, and the partition tends to infinity. An exact calculation of Q could be complicated due to the coupling of the internal modes, a way to decouple the electronic and nuclei modes is through the use of Born-Oppenheimer approximation. (BOA) This approach says that the electron movement is faster than the nuclei and assumes that the molecular wave function is the electronic and nuclear wavefunction product ψ ¼ ψ e ψ n . The vibrations change the momentum of inertia as a consequence, affect the rotations; this fact tightly couple the vibrational and rotational degrees of freedom; The separation of rotational and vibrational modes is called the rigid rotor, harmonic oscillator (RRHO) approximation, under this approximation, the molecule is treated rigidly, this is generally good when vibrations are of small amplitude. Here the vibration will be modeled in terms of harmonic oscillator and rotations in terms of the rigid rotor. Within BOA and RRHO approximations, the partition function is factorized into electronic, translational, vibrational, and rotational energies. Consequently, the partition function, Q , can be given in Eq. (2) as a product of the corresponding contributions [100,139], and under the rigid rotor, harmonic oscillator, Born-Oppenheimer, ideal gas, and a particle-in-a-box approximations.
Q ¼ q trans q rot q vib q elec: (2) Table 1 shows the contributions of electronic, translational, vibrational, and rotational to the partition function.
We computed all partition functions at temperature T and a standard pressure of 1 atm. The equations are equivalent to those given in the Ref. [100], and any standard text of thermodynamics [138,139] and they apply for an ideal gas. The implemented translational partition function in the Gaussian code [136] is the partition function, q ¼ q trans , given in Table 1. In this study, the q ¼ q trans is computed as a function of T and is used to calculate the translational entropy. In addition to using vibrational modes to identify true lowest energy structures from transition states, we also used them to compute the vibrational partition function. In this study is considered vibrational modes, ν, under the harmonic oscillator approximation, and total vibrational energy consists of the sum of the energies of each vibrational mode. In computing the electronic partition, we considered that the energy gap between the first and higher excited states is more considerable than k B T, as a consequence electronic partition function, q ¼ q elect , is given by q elect ¼ ω 0 , q rot , q nl rot , q ¼ q trans are used to compute the entropy contributions given in Table 2. The vibrational frequencies are calculated employing the Gaussian code, and all the information needed to compute the total partition function is collected from the output. The Gibbs free energy and the enthalpy are computed employing the Eqs. (3) and (4).

Boltzmann population
The properties observed in a molecule are statistical averages over the ensemble of geometrical conformations or isomers accessible to the cluster [140]. So, the molecular properties are ruled by the Boltzmann distributions of isomers that can change due to temperature-entropic term [23, 71,101], and the soft vibrational modes that clusters possess make primary importance contributions to the entropy [93]. The Boltzmann populations of the low-energy isomers of the cluster Be 6 B 11 À and Be 4 B 8 are computed through the probabilities defined in Eq. (5)

Contribution Partition function
Translational where β ¼ 1=k B T, and k B is the Boltzmann constant, T is the temperature in Kelvin, ΔG is the Gibbs free energy of the k th isomer. Eq. (5) establishes that the distribution of molecules will be among energy levels as a function of the energy and temperature. Eq. (5) is restricted so that the sum of all probabilities of occurrence, at fixed temperature T, Pi (T) is equal to 1 and given by Eq. (6) It is worth mentioning that the energy difference among isomers is determinant in the computation of the solid-solid transition, T ss point. T ss occurs when two competing structures are energetically equaled, and there is simultaneous coexistence of isomers at T. in other words, the T ss point is a function of the energy difference between two isomers and the relative energy ΔG that the cluster possesses. Boltzmann distribution finds a lot of applications as like simulated method annealing applied to the search of structures of minimum energy, rate of chemical reaction [100], among others. For the calculation of the Boltzmann populations, we used a homemade Python/Fortran code called BOFA (Boltzmann-Optics-Full-Ader).

Computational details
The global exploration of the potential and free energy surfaces of the Be 6 B 11 À and Be 4 B 8 clusters were done with a hybrid Cuckoo-genetic algorithm written in Python. All local geometry optimization and vibrational frequencies were carried out employing the density functional theory (DFT) as implemented in the Gaussian 09 [136] suite of programs, and no restrictions in the optimizations were imposed. Final equilibrium geometries and relative energies are reported at PBE0 [141] /def2-TZVP [142] level of theory, taking into account the D3 version of Grimme's dispersion corrections [143]. and including the zero-point (ZPE) energy corrections. As Pan et al. [144] reported, the computed relative energies with PBE0 functional are very close to the CCSD(T) values in B 9 À boron cluster. The def2-TZVP basis from the Ahlrichs can improve computations accuracy and describe the Be-B clusters [29] To gain insight into its energetics, we evaluated the single point energy CCSD (T)/def2TZVP//PBEO-D3/def2-TZVP level of theory for the putative global

Internal energy Entropy
Translational minima and the low-energy Be 6 B 11 À isomers, and employing Orca code at the DLPNO-CCSD(T) theoretical level for the low-energy isomers of Be 4 B 8 cluster. Boltzmann-Optics-Full-Ader, (BOFA) is employed in the computation of the Boltzmann populations. The code is available with the corresponding author. Figure 1 shows the lowest energy structure of Be 6 B 11 À clusters and seven lowenergy competing isomers computed at the PBE0-D3/def2-TZVP basis set. For the putative global minimum at the PBE0-D3/def2-TZVP, the optimized average B-B bond length is 1.64 Å. In contrast, the optimized B-Be bond length is 2.01 Å. At the PBE0-D3/def2-TZVP and temperature of 298.15 K, the putative global minimum with 54% of the relative population has C 1 symmetry with a singlet electronic state 1 A. It is a distorted, oblate spheroid with three berylliums in one face and two in the other face. Nine-boron and one-beryllium atoms are forming a ring located around the spheroid's principal axes and the remaining two boron atoms are located close to the boron ring in one of its faces. The second higher energy structure, at 298.15 K, lies only 0.61 kcal/mol Gibbs free energy above the putative global minima, and it has C 1 symmetry with a singlet electronic state 1 A. It is a prolate spheroid with 19% of the relative population at a temperature of 298.15 K. The next two higher energy isomers, at 298.15 K, lies at 0.85 and 1.23 kcal/mol Gibbs energy above the putative global minimum. They are prolate, coaxial Triple-Layered structures with C s , and C 2v symmetries with singlet electronic states, 1 A, respectively. This clearly, shows that the low-symmetry structure C 1 become more energetically preferred than the C 2v symmetry by Gibbs free energy difference of 0.38 kcal/mol at 298.15 K, due to entropic effects and in agreement with a similar result found in Au 32 [105]. Indeed, and according to our computations, those structures are strongly dominating at temperatures higher than 377 K. The next structure is shown in Figure 1(5), is located at 1.48 kcal/mol above the global minimum; it is close to a spherical shape and correspond to a prolate structure with C 1 symmetry, and a singlet electronic state 1 A; this structure only has 4.4% of the relative population at 298.15 K. The next two structures, located at 2.37 kcal/mol Gibbs free energy above the global minimum, are the chiral helix-type structures, reported by Guo [29] as minimum global. They are prolate structures with C 2v symmetries, and their relative population is around only 1%. We must point out that those chiral-helix structures never become the lowest energy structures in all ranges of temperature. The relative population is zero for structures located at higher relative Gibbs free energy than 5.1 kcal/mol, and at 298.15 K, there is no contribution of these isomers to any total molecular property. A full understanding of the molecular properties requires the search of global minimum and all its closest low-energy structures [64]. The separation among isomers by energy-difference is an important and critical characteristic that influences the relative population and, consequently, the total molecular properties. We computed the global minima and the first seven low-energy to gain insight into how the energy-gap among isomers change and how the energy-ordering of the low-energy structures is affected at a single point CCSD(T)/def2-TZVP level of theory corrected with the zero-point energy computed at the PBE0-D3/def2-TZVP level of theory. At the CCSD(T) level of theory, the global minima, the seven lowest energy isomers, and the energy order agree with previous work [39], as seen in the first row of Table 3. The second row of Table 3 shows the corrected CCSDT+E ZPE energy. Interestingly, the energetic ordering does not change when we take into account the ZPE energy. Nevertheless, the energy difference among isomers was reduced drastically. we can deduce that the ZPE energy inclusion is essential in the isomers' energy ordering and molecular properties. The third row of Table 3 shows the energy-order considering the Gibbs free energy computed at 298.15 K; at this temperature, the isomers energy-ordering is changed, the second isomers take the putative global minima place, and the first isomers take the fifth place. Interestingly, this energy-ordering is at 298.15 K. This energy-ordering is a complete function of the temperature that we will discuss later in the relative population section. The fourth row in Table 3 shows the electronic energy taking into account the ZPE energy. It follows the same trend in energy-ordering when considering the Gibbs free energy, and it is the same putative global minima. The fifth row in Table 3 is just electronic energy. It almost follows the CCSD(T) energies trend, except the isomers number eight that take the second place located at 0.52 kcal/mol above the putative global minima. The sixth, seventh and eighth rows on Table 3 show the point group symmetry, electronic ground state, and the lowest vibrational frequency of each isomer. When we take the Gibbs free energy to energy-ordering structures, the second isomers interchange to the first place, becoming the lowest energy structure; The energy ordering change drastically, whereas the electronic energy almost follows the same trend CCSD(T) energy-ordering. This shows us that the level of theory and the inclusion of entropy and temperature change the energyordering; therefore, the total molecular properties.

Boltzmann population of Be 6 B 11
À cluster Figure 2a shows the most important and strongly dominating T ss1-g point that is located at 377 K temperature scale with a relative population of 33%. For temperatures ranging from 10 to 377 K, the relative population is strongly dominated by the putative global minima isomer distorted oblate spheroid with C 1 symmetry and this relative population is similar to -T À3 function with one point of inflection located at 180 K. After decreases monotonically up to 377 K. At the T ss1-g point, the distorted oblate spheroid with C 1 symmetry co-exist and compete with the coaxial Triple-Layered structures with C s symmetry; This implies that the distorted oblate spheroid will be replaced with the coaxial Triple-Layered structures. Above temperature 377 K, the relative population is strongly dominated by the coaxial Triple-Layered structures with C s symmetry, located at 0.85 kcal/mol above the global minima at temperature 298.15 K. This relative population depicted in blue-solid line in panel (a) has behavior as a sigmoid function, from temperatures ranging from 377 to 600 K, it grows rapidly and from temperatures ranging from 600 to 1500 K, it almost keeps constant with 60%. The second T ss2-g point is located at temperature 424 K with a relative population of 22.9%, and this point the global minima distorted oblate spheroid with C 1 symmetry co-exist, and compete with the coaxial Triple-Layered structures with C 2v symmetry, located at 1.23 kcal/mol above the global minima at 298.15 K. The relative population of the coaxial Triple-Layered C 2v symmetry depicted in green-solid line in panel (a) also has a behavior of a sigmoid function and up to 600 K it keeps constant with 32% of relative population. The T ss3-g , and T ss4-g points are located at 316.7 K, and 349 K axis temperature with relative populations 14% and 17%, respectively. These relative populations correspond to the second isomer located just 0.61 kcal/mol at 298.15 K above the global minima, and co-existing at the temperatures 316.7 K and 349 K with the coaxial Triple-Layered structures with C s , and C 2v symmetries, respectively. At low temperatures range, this isomer's relative population depicted in red-solid line of Figure 2a is around only 20%, and up to room temperature, it decreases exponentially to zero. At temperatures up to 600 K, the relative population is zero; hence, at high temperatures these isomers do not contribute to the molecular properties. The relative population lower than 10%, depicted in violet-solid line shows in Figure 2a, correspond to the isomers located at 1.48 kcal/mol above global minima at 298.15 K. Interesting, this structure is the putative minimum global when the CCSD(T) energy is employed in the ordering energetic, Despite that, this structure's relative population clearly shows that this structure does not contribute to molecular properties in all ranges of temperatures. Figure 3 shows the low-energy configurations of Be 4 B 8 clusters optimized at PBE0-D3/def2-TZVP level of theory taking into account ZPE energy correction. The optimized average B-B bond length of the putative chiral global minimum is 1.5867 Å, in good agreement with an experimental bond length of 1.57-1.59 Å [145,146], and also within agreement with others previous DFT calculations [39]. The most recurring motif within the lower energy isomers of B 8 Be 4 is a sandwich structure, (SSh) in which the boron atoms form a hollow distorted ellipsoid ring with each of the Be-Be dimers capping the top and bottom with C 1 point group symmetry. Isomers 1 and 2 are also listed as i 1 and i 2 in Table 4, are enantiomers differing in the orientation of the Be-Be dimers with respect to the boron skeleton. The Be-Be bond length for the six lowest energy enantiomers is 1.9874, 1.9876, and 1.9881 Å for symmetries C 1 , C 2 , and D 2 , respectively, in good agreement with the bond length of the Be-Be in Be 2 B 8 cluster 1.910 Å [44]. To gain insight into the energy hierarchy of isomers and validate our DFT calculations, relative energies were computed at different levels of theory, and differences between them are shown in Table 4. Energy computed at different methods yield different energies due mainly to the functional and basis-set employed, [39,147], so the energetic ordering change; consequently, the probability of occurrence and the molecular properties will change. The first line of Table 4 shows the relative Gibbs free energy computed at PBE0-D3/def2-TZVP and room temperature. The small relative Gibbs free energies (0.41, and 0.81 kcal/mol) differences among the six enantiomer structures i 1 to i 6 in Table 4 are caused by the rotational entropy being a function of the symmetry number that in turn depends on the point group symmetry. An increase/decrease in the value of rotational entropy changes the Gibbs free energy. The Gibbs free energy computed with and without symmetry will differ by a factor RTln(σ). Here, R is the universal gas constant, T, the temperature, and σ is the symmetry number. The computed factor at room temperature with σ = 2 is RTln (σ) = 0.41 kcal/mol, and it is RTln(σ) = 0.81 kcal/mol with σ = 4, in agreement with the values shown in the first line of Table 4. As the temperature increases, the energy differences between the factors RTln(σ) become larger. These small relative Gibbs free energies are responsible for different values of probability of occurrence at low temperatures for the similar isomers with different point group symmetry. This strongly suggests that there must be atomic clusters with low and high symmetries in the Boltzmann ensemble to compute the molecular properties correctly. The second line in Table 4 shows single point (SP) relative energies computed at the CCSD(T) [148], the energetic ordering of isomers listed in the first line of Table  4 follows almost the trend of energetic ordering at SP CCSD(T) level, notice that just the achiral isomers label i 7 to i 8 in Table 4 are interchanged in energetic ordering. The third line Table 4 shows single point relative energies computed at Point Group Symmetry the CCSD(T) [148]/def2-TZVP//PBE0-D3/def2-TZVP; the energetic ordering is similar to pure CCSD(T) energy. DLPNO-CCSD(T) relative energies, with and without ZPE correction, are shown in lines four and five of Table 4, the first follows the trend of pure CCSD(T) energy, and the second, the ZPE value, interchange the isomers, label i 7 in Table 4, to be the putative global minimum. Here we can say that the ZPE energy inclusion is essential in distributing isomers and molecular properties. The sixth and seventh lines of Table 4 show the electronic energy with and without ZPE correction, and both of them follow the trend of the Gibbs free energy given in line number one. Line number 8 in Table 4 shows the point group symmetry for each isomer. The T 1 diagnostic for each isomer is shown in line nine of Table 4, all of them are lower than the recommended value 0.02 [148] so the systems are appropriately characterized.

Boltzmann population of Be 4 B 8 clusters
As we mentioned earlier, the determination of the structure is the first step to study any property of a material. Moreover, we have to consider that an observed molecular property in a Boltzmann ensemble is a weighted sum of all individual contributions of each isomer that form the ensemble. At temperature 0 K, the electronic energy plus zero-point energy determine the putative global minimum and all nearby low-energy structures, whereas, at temperatures larger than 0 K, the Gibbs free energy defines the putative global minimum. Figure 4a shows the probability of occurrence computed at PBE0-D3/def2-TZVP level of theory for each particular chiral and achiral Be 4 B 8 isomers for temperatures ranging from 20 to 1900 K. Figure 4b shows the probability of occurrence computed at CCSD(T) level of theory. Notice, there is not a significant difference in the probabilities of occurrence between the two panels, thus the computation of probabilities at DFT level of theory is very similar to those computed at CCSDT level of theory. A closer examination of the panel (b) shown that in the temperature ranging from 20 to 300 K, all molecular properties are dominated by the chiral structure depicted in Figure 3a because its probability of occurrence is almost constant. We point out that in this range of temperature, the C 1 , C 2, and D 2 symmetries strongly dominate with different probabilities of occurrence of 28, 14 y 7% respectively. At this point, there is  a co-existence of chiral structures and achiral structures, shown in Figure 3, above this point the achiral structure (Figure 3g) becomes dominant. The second transformation solid-solid point located at 1017 K and 10% of probability also coexist the chiral putative global minimum with symmetry C 1 and achiral structure (Figure 3h) located at 2.51 kcal/mol CCSDT energy at above the putative global minimum. The Boltzmann population computed at PBE0-D3/def2-TZVP level of theory follows the trend of the Boltzamnn population computed at CCSD(T) level of theory.

Conclusions
We computed the Boltzmann population of anionic Be 6 B 11 and neutral Be 4 B 8 cluster at the SP CCSDT and DFT levels of theory. If one increases the system's temperature, entropic effects start to play an important role, and Gibbs free energy is minimized, and entropy is maximized. The fluxionality of the Be 6 B 11 À cluster is strongly dependent on temperature that is shown by its Boltzmann population. At the CCSDT level of theory, the Boltzmann population of the Be 6 B 11 À cluster indicate there are four competing structures, so a mixture of isomers co-exist at a specific temperature, so we expect that around a temperature of 350 K, four structures could be observed. The observed properties in a molecule are statistical averages over the ensemble of isomers. The molecular properties at cold temperatures are due to the lowest energy structure Be 6 B 11 À at CCSD(T) level of theory and zero temperature whereas at hot temperatures, the molecular properties are due to the coaxial Triple-Layered structure with C 1 symmetry. At room temperature the molecular properties are due to a mixture of spectra of the three systems that coexist at 350 K. Regarding Be 4 B 8 cluster, all molecular properties at cold and room temperatures are dominated by pair of enantiomers putative global minima. The computed Boltzmann populations at PBE0-D3/def2-TZVP level of theory is similar at the computed Boltzmann populations at CCSDT/def2-TZVP level of theory, so at the DFT level, the Boltzmann populations, hence the molecular properties are well calculated. As future work, the inclusion of anharmonic effects should be taken into account.  [23] An, W.; Bulusu, S.; Gao, Y.; Zeng, X.C. Relative stability of planar versus double-ring tubular isomers of neutral and anionic boron cluster B20 and B20-.