Contributions to the partition function.

## Abstract

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 Be6B11− and chiral Be4B8 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 Be6B11− cluster is strongly dominant at hot temperatures, whereas the chiral Be4B8 cluster is dominant at room temperature. The methodology and results show the thermal effects in the relative population hence molecular properties.

### Keywords

- Global minimum
- infrared spectrum
- DFT
- boron cluster
- fluxional
- density functional theory
- temperature
- Boltzmann
- Gibbs free energy
- entropy
- CCSDT
- statistical thermodynamics

## 1. 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 temperature properties remain unexplored [61, 62]. Experimental studies are carried out in non-zero temperatures, then it is necessary to understand the effect of the temperature on the cluster properties and the lowest energy structure’s determination [61, 62, 63]. Herein, we investigate the effect of temperature-entropy term on the Boltzmann population, which needs the elucidation of the putative global minimum and its low-energy isomers [39, 64, 65, 66, 67, 68]. The properties observed in a molecule are statistical averages over the ensemble of geometrical conformations that are ruled by the Boltzmann distributions of isomers. So we need an efficiently sampling of the free energy surface to know the distribution of isomers at different temperatures [39, 68, 69, 70, 71]. A considerable change in the isomer distribution and the energetic separation among them is the first notable effect of temperature [39]. Useful materials work at finite temperatures; in that conditions, Gibbs free energy is minimized whereas, the entropy of the atomic cluster is maximized [39, 72]. and determines the putative global minimum at a finite temperature [39]. Although in the mid 1960’s, Mermin et al. [73] studied the thermal properties of the inhomogeneous electron gas, most DFT calculations are typically performed at zero temperature. Recently, over again, DFT was extended to finite temperature [74, 75, 76], but nowadays, as far as we know, it is not implemented in any public software and practical calculations are not possible. Taking temperature into account requires dealing with small systems’ thermodynamics; The Gibbs free energy of classical thermodynamics also applies for small systems, known as thermodynamics of small systems [77, 78, 79]. The thermodynamics of clusters have been studied by various theoretical and simulation tools [61, 68, 77, 80, 81, 82, 83, 84, 85, 86] like molecular-dynamics simulations. Previous reports investigated the behavior of Al_{12}C cluster at finite temperature employing Car-Pirinello molecular dynamics [61], and dynamical behavior of Borospherene in the framework of Born-Oppenheimer Molecular Dynamics [5, 10]. Under the harmonic superposition approximation, the temperature-entropy term can be computed with the vibrational frequencies on hand. The entropy and thermal effects have been considered for gold, copper, water, and sodium clusters [71, 87, 88, 89, 90, 91, 92, 93, 94, 95]. Franco-Perez et al. [96] reported the thermal corrections to the chemical reactivity at finite temperature, their piece of work validates the usage of reactivity indexes calculated at zero temperature to infer chemical behavior at room temperature. Gazquez et al. [97] presented a unified view of the temperature-dependent approach to the DFT of chemical reactivity. Recently, the effect of temperature was considered by Castillo-Quevedo et al. reported the reaction rate and the lowest energy structure of copper Cu_{13} clusters at finite temperature [98, 99]. Dzib et al. reported Eyringpy; A Python code able to compute the rate constants for reactions in the gas phase and in solution [100], Vargas-Caamal et al. computed the temperature-dependent dipole moments for the HCl(H_{2}O)_{n} clusters [101], Shkrebtii et al. computed the temperature-dependent linear optical properties of the Si(100) surface [102], several authors take into account the temperature in gold clusters [87, 88, 89], and thermochemical behavior study of the sorghum molecule [103], and more recently, Buelna-Garcia et al. [104] employing density functional theory and nanothermodynamics reported the lowest energy structure of neutral chiral Be_{4}B_{8} at a finite temperature. and reported that the fluxionality of the anionic Be_{6}B_{11} clusters depends strongly on temperature [39]. In this work, we employed density functional theory, statistical thermodynamics, and CCSD(T) to compute the Gibbs free energy and the Boltzmann population at absolute temperature T for each neutral chiral Be_{4}B_{8} and anionic Be_{6}B_{11} isomers. We think that this provides useful information about which isomers will be dominant at hot temperatures. No work has previously been attempted to investigate entropy-driven isomers in the fluxional Be_{6}B_{11}^{−} and chiral Be_{4}B_{8} cluster at CCSDT level of theory to the best of our knowledge. The remainder of the manuscript is organized as follows: Section 2 gives the computational details and a brief overview of the theory and algorithms 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.

## 2. Theoretical methods and computational details

### 2.1 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 nature-inspired 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].

### 2.2 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 temperature-dependent partition function Q shown in Eq. (1).

In Eq. (1), the

Table 1 shows the contributions of electronic, translational, vibrational, and rotational to the partition function.

Contribution | Partition function |
---|---|

Translational | |

Rotational linear | |

Rotational nonlinear | |

Vibrational | |

Electronic |

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,

Internal energy | Entropy | |
---|---|---|

Translational | ||

Rotational linear | ||

Rotational nonlinear | ||

Vibrational | ||

Electronic |

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).

### 2.3 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)

where ^{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

### 2.4 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 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.

## 3. Results and discussion

### 3.1 The lowest energy structures and energetics

Figure 1 shows the lowest energy structure of Be_{6}B_{11}^{−} clusters and seven low-energy 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 energy-ordering; therefore, the total molecular properties.

Be_{6}B_{11} | Level | i_{1} | i_{2} | i_{3} | i_{4} | i_{5} | i_{6} | i_{7} | i_{8} |
---|---|---|---|---|---|---|---|---|---|

0.0 | 1.75 | 1.84 | 1.84 | 4.10 | 4.13 | 2.64 | 2.42 | ||

0.0 | 0.58 | 0.85 | 0.86 | 1.19 | 1.23 | 1.68 | 1.81 | ||

0.0 | −1.48 | 0.89 | 0.88 | −0.63 | −0.25 | 4.14 | −0.87 | ||

Be_{6}B_{11} | 0.0 | −0.29 | 1.51 | 1.52 | 2.41 | 2.42 | 5.0 | −0.08 | |

0.0 | 0.87 | 2.50 | 2.50 | 5.32 | 5.32 | 5.96 | 0.52 | ||

Point Group Symmetry | C_{1} | C_{1} | C_{2} | C_{2} | C_{S} | C_{2v} | C_{1} | C_{1} | |

Electronic ground state | ^{1}A | ^{1}A | ^{1}A | ^{1}A | ^{1}A’ | ^{1}A_{1} | ^{1}A | ^{1}A | |

Frequencies (cm^{−1}) | 230 | 119 | 102 | 100 | 46 | 43 | 161 | 151 |

### 3.2 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.

### 3.3 The lowest energy structures of Be_{4}B_{8} clusters

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 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.

Level | i_{1} | i_{2} | i_{3} | i_{4} | i_{5} | i_{6} | i_{6} | i_{8} | i_{9} | i_{10} |
---|---|---|---|---|---|---|---|---|---|---|

0.0 | 0.0 | 0.41 | 0.41 | 0.81 | 0.81 | 1.79 | 2.40 | 4.45 | 4.45 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 3.61 | 3.38 | 5.38 | 5.38 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 2.71 | 2.51 | 4.51 | 4.51 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.75 | 1.37 | 5.0 | 5.0 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | −0.20 | 0.50 | 4.10 | 4.10 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 2.38 | 2.80 | 5.03 | 5.03 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 3.28 | 3.68 | 5.90 | 3.28 | |

Point Group Symmetry | C_{1} | C_{1} | C_{2} | C_{2} | C_{1} | C_{1} | C_{1} | C_{1} | C_{2} | C_{2} |

0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 |

### 3.4 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.

## 4. 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.

## Acknowledgments

C.E.B.-G. thanks Conacyt for a scholarship (860052). E.P.-S. thanks Conacyt for a scholarship (1008864). We are grateful to Dr. Carmen Heras, and L.C.C. Daniel Mendoza for granting us access to their clusters and computational support. Computational resources for this work were provided through the High-Performance Computing Area of the University of Sonora. We are also grateful to the computational chemistry laboratory for providing computational resources, ELBAKYAN, and PAKAL supercomputers.