Boltzmann Populations of the Fluxional Be₆B₁₁⁻ and Chiral Be₄B₈ Clusters at Finite Temperatures Computed by DFT and Statistical Thermodynamics

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.

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₁₉. 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₆₀ and B₈₀ 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₆B₁₁⁻, 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₆B₁₀⁻² and coaxial triple-layered anionic Be₆B₁₁ 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₃₀ 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₃₉ borospherene were reported as chiral due to their hexagonal and pentagonal holes [26]. Similarly, the B₄₄ 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₁₂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₁₃ 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₂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₄B₈ at a finite temperature. and reported that the fluxionality of the anionic Be₆B₁₁ 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₄B₈ and anionic Be₆B₁₁ 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₆B₁₁⁻ and chiral Be₄B₈ 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₂ and C₁ 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₁ 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)₂] complex, infrared spectroscopy on anionic Be₆B₁₁ 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).

$$Q\left(T\right)=\sum _i{g}_i{e}^{-{\mathit{\Delta E}}_i/k_{B}T}$$E1

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 ${\mathit{\Delta E}}_i$ is the total energy of a molecule [100, 138]. At high temperatures, all thermal states are accessible due to the term $-{\mathit{\Delta 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 $\psi ={\psi }_e{\psi }_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_{\mathit{\text{trans}}}{q}_{\mathit{rot}}{q}_{\mathit{vib}}{q}_{\mathit{\text{elec}}.}$$E2

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

Contribution	Partition function
Translational	$q_{\mathit{\text{trans}}}={\left(\frac{2{\mathit{\pi mk}}_{B}T}{h^2}\right)}^{\frac{3}{2}}\frac{k_{B}T}{P}$
Rotational linear	$q_{\mathit{rot}}^l=\frac{T}{{\mathit{\sigma \Theta }}_{\mathit{rot}}},{\Theta }_{\mathit{rot}}=\frac{{\mathrm{\hslash }}^2}{2{\mathit{Ik}}_B}$
Rotational nonlinear	$q_{\mathit{rot}}^{\mathit{nl}}=\frac{{\pi }^{1/2}}{\sigma }\left[\frac{T^{3/2}}{{\left({\Theta }_{\mathit{\text{rotA}}}{\Theta }_{\mathit{\text{rotB}}}{\Theta }_{\mathit{\text{rotC}}}\right)}^{1/2}}\right],{\Theta }_{\mathit{\text{rotj}}}=\frac{{\mathrm{\hslash }}^2}{2I_j{k}_B},j=A,B,C$
Vibrational	$q_{\mathit{vib}}^{\mathit{pol}}=\prod _{i=1}^{n_{{\upsilon }^{\left[a\right]}}}\frac{e^{-{\Theta }_{{\mathit{vib}}_i}/2T}}{1-e^{-{\Theta }_{{\mathit{vib}}_i}/T}},{\Theta }_{{\mathit{vib}}_i}=\frac{{\mathit{h\nu }}_i}{k_B}$
Electronic	$q_{\mathit{\text{elec}}}={\omega }_{\mathit{0}}$

Table 1.

Contributions 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_{\mathit{\text{trans}}}$, given in Table 1. In this study, the $q=q_{\mathit{\text{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, $\nu$, 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_{\mathit{\text{elect}}}$, is given by $q_{\mathit{\text{elect}}}={\omega }_{\mathit{0}}$, $q_{\mathit{rot}}$, $q_{\mathit{rot}}^{\mathit{nl}}$, $q=q_{\mathit{\text{trans}}}$ are used to compute the entropy contributions given in Table 2.

	Internal energy	Entropy
Translational	$U_{\mathit{\text{trans}}}=\frac{3}{2}\mathit{RT}$	$S_{\mathit{\text{trans}}}=R\left(ln{q}_{\mathit{\text{trans}}}+\frac{5}{2}\right)$
Rotational linear	$U_{\mathit{rot}}^l=\mathit{RT}$	$S_{\mathit{rot}}^l=R\left(ln{q}_{\mathit{rot}}^l+1\right)$
Rotational nonlinear	$U_{\mathit{rot}}^{\mathit{nl}}=\frac{3}{2}\mathit{RT}$	$S_{\mathit{rot}}^{\mathit{nl}}=R\left(ln{q}_{\mathit{rot}}^{\mathit{nl}}+\frac{3}{2}\right)$
Vibrational	$U_{\mathit{vib}}^{\mathit{pol}}=R\sum _i^{n_{{\upsilon }^{\left[a\right]}}}{\Theta }_{{\mathit{vib}}_i}\left(\frac{1}{2}+\frac{1}{e^{{\Theta }_{{\mathit{vib}}_i}/T}-1}\right)$	$S_{\mathit{vib}}^{\mathit{pol}}=R\sum _i^{\mathit{n\upsilon }}\left[\frac{{\Theta }_{{\mathit{vib}}_i}/T}{e^{{\Theta }_{{\mathit{vib}}_i}/T}-1}-ln\left(1-e^{-{\Theta }_{{\mathit{vib}}_i}/T}\right)\right]$
	${\Theta }_{{\mathit{vib}}_i}=\frac{{\mathit{h\nu }}_i}{k_B}$
Electronic	$U_{\mathit{\text{elec}}}=0$	$S_{\mathit{\text{elec}}}=Rln{q}_{\mathit{\text{elec}}}$

Table 2.

Contributions to internal energy and entropy.

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

$$H=U+\mathit{nRT},$$E3

$$G=H-\mathit{TS}.$$E4

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₆B₁₁⁻ and Be₄B₈ are computed through the probabilities defined in Eq. (5)

$$P\left(T\right)=\frac{e^{-{\mathit{\beta \Delta G}}^K}}{\sum e^{-{\mathit{\beta \Delta G}}^K}}$$E5

where $\beta =1/k_{B}T$, and $k_B$ is the Boltzmann constant, T is the temperature in Kelvin, $\mathit{\Delta 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)

$$\sum _{i=1}^{n}P\left(T\right)=1$$E6

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 $\mathit{\Delta 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).

2.4 Computational details

The global exploration of the potential and free energy surfaces of the Be₆B₁₁⁻ and Be₄ B₈ 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₉⁻ 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₆B₁₁⁻ isomers, and employing Orca code at the DLPNO-CCSD(T) theoretical level for the low-energy isomers of Be₄B₈ 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₆B₁₁⁻ 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₁ symmetry with a singlet electronic state ¹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₁ symmetry with a singlet electronic state ¹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, ¹A, respectively. This clearly, shows that the low-symmetry structure C₁ 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₃₂ [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₁ symmetry, and a singlet electronic state ¹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.

Figure 1.

The optimized geometries of Be₆B₁₁⁻ cluster. The most important energy isomers show in two orientations, front, and rotated 90 degrees up to plane paper. Relative Gibbs free energies in kcal/mol (in round parenthesis) and the relative population [in square parenthesis], at PBE0-D3/Def2-TZVP level of theory. The criterium to plot them is until the probability occupation is zero. The pink- and yellow-colored spheres represent the boron and beryllium atoms, respectively.

Be6B11	Level	i1	i2	i3	i4	i5	i6	i7	i8
	$\mathit{\text{CCSD}}\left(T\right)$	0.0	1.75	1.84	1.84	4.10	4.13	2.64	2.42
	$\mathit{\text{CCSD}}\left(T\right)+{\epsilon }_{\mathit{ZPE}}$	0.0	0.58	0.85	0.86	1.19	1.23	1.68	1.81
	$\mathit{\Delta G}$	0.0	−1.48	0.89	0.88	−0.63	−0.25	4.14	−0.87
Be6B11	${\epsilon }_0+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$	0.0	−0.29	1.51	1.52	2.41	2.42	5.0	−0.08
	${\epsilon }_0$	0.0	0.87	2.50	2.50	5.32	5.32	5.96	0.52
	Point Group Symmetry	C1	C1	C2	C2	CS	C2v	C1	C1
	Electronic ground state	^1A	^1A	^1A	^1A	^1A’	^1A1	^1A	^1A
	Frequencies (cm^−1)	230	119	102	100	46	43	161	151

Table 3.

The relative energies in kcal·mol⁻¹, coupled-cluster single-double and perturbative triple, $\mathit{\text{CCSD}}\left(T\right)$, $\mathit{\text{CCSD}}\left(T\right)$ with zero-point energy (${\epsilon }_{\mathit{ZPE}}$), ($\mathit{\text{CCSD}}\left(T\right)+{\epsilon }_{\mathit{ZPE}}$), Gibbs free energy ($\mathit{\Delta G}$) at 298.15 K, electronic energy with ${\epsilon }_{\mathit{ZPE}}$ ($\%{\mathit{\text{mcalEε}}}_{\mathit{0}}$), electronic energy (${\epsilon }_0$), point group symmetry, electronic ground state, and the lowest frequency in cm⁻¹ for eight low-energy isomers.

3.2 Boltzmann population of Be₆B₁₁⁻ 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₁ symmetry and this relative population is similar to -T⁻³ 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₁ 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₁ 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 2.

Panel (a) shows the Boltzmann population of the Be₆B₁₁⁻ (ensemble at thermal equilibrium) for the temperatures ranging from 10 to 1500 K computed at the PBE0-D3/def2-TZVP level of theory. Panel (b) shows the Boltzmann population for the temperatures ranging from 10 to 1500 K computed at the CCSDT/def2-TZVP//PBE0-GD3/def2-TZVP level of theory. At the temperature of 350 K, four structures co-exist with 20% of probability.

3.3 The lowest energy structures of Be₄B₈ clusters

Figure 3 shows the low-energy configurations of Be₄B₈ 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₈Be₄ 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₁ point group symmetry. Isomers 1 and 2 are also listed as i₁ and i₂ 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₁, C₂, and D₂, respectively, in good agreement with the bond length of the Be-Be in Be₂B₈ 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₁ to i₆ 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₇ to i₈ 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₇ 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₁ 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.

Figure 3.

Optimized geometries of a neutral Be₄B₈ cluster at the PBE0-D3/def2TZVP level of theory with zero-point correction energy. These are shown in front and side views. The first letter is the isomer label, the relative Gibbs free energies in kcal·mol⁻¹ (in round parenthesis) at 298.15 K, the relative population (in square brackets), and the group symmetry point (in red round parenthesis). The structures with labels (a and b), (c and d), (e and f), (i and j), (k and l) and (h) are chiral. The purple- and yellow-colored spheres represent the boron and beryllium atoms, respectively.

Level	i1	i2	i3	i4	i5	i6	i6	i8	i9	i10
$\mathrm{\Delta G}$	0.0	0.0	0.41	0.41	0.81	0.81	1.79	2.40	4.45	4.45
$\text{CCSD}\left(T\right)$	0.0	0.0	0.0	0.0	0.0	0.0	3.61	3.38	5.38	5.38
$\text{CCSD}\left(T\right)+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$	0.0	0.0	0.0	0.0	0.0	0.0	2.71	2.51	4.51	4.51
$\text{DLPNO}‐\text{CCSD}\left(T\right)$	0.0	0.0	0.0	0.0	0.0	0.0	0.75	1.37	5.0	5.0
$\text{DLPNO}‐\text{CCSD}\left(T\right)+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$	0.0	0.0	0.0	0.0	0.0	0.0	−0.20	0.50	4.10	4.10
${\epsilon }_0+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$	0.0	0.0	0.0	0.0	0.0	0.0	2.38	2.80	5.03	5.03
${\epsilon }_0$	0.0	0.0	0.0	0.0	0.0	0.0	3.28	3.68	5.90	3.28
Point Group Symmetry	C1	C1	C2	C2	C1	C1	C1	C1	C2	C2
$T_1$	0.019	0.019	0.019	0.019	0.019	0.019	0.019	0.019	0.019	0.019

Table 4.

Single-point relative energy calculations of the low-energy structures from i₁ to i₁₀ at different levels of theory: coupled cluster single-double and perturbative triple (CCSD(T)), CCSD(T) with zero-point energy ($\text{CCSD}\left(T\right)+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$, CCSD(T)) employing the domain-based local pair natural orbital coupled-cluster theory (DLPNO-CCSD(T)), with TightPNO setting, and with ${\epsilon }_{\mathrm{ZPE}}$ ($\text{DLPNO}‐\text{CCSD}\left(T\right)+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$), Gibbs free energy ($\mathrm{\Delta G}$) at 298.15 K, electronic energy with ${\epsilon }_{\mathrm{ZPE}}$ (${\epsilon }_0+{\mathrm{\epsilon }}_{\mathrm{ZPE}}$), electronic energy (${\epsilon }_0$), point group symmetry, and T₁ diagnostic. All relative energies are given in kcal·mol⁻¹.

3.4 Boltzmann population of Be₄B₈ 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₄B₈ 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₁, C_2, and D₂ 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₁ 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.

Figure 4.

Panel (a) shows the probability of occurrence for temperatures ranging from 20 to 1900 K at the PBE0-D3/def2-TZVP level of theory. Panel (b) shows the probability of occurrence for temperatures ranging from 20 to 1900 K at the CCSDT/def2-TZVP level of theory. In panel (a), the transition solid–solid point (T_ss1-g) is located at 739 K with 16.6% of probability, while in panel (b) the T_ss1-g is located at 780 K with 15% probability.

4. Conclusions

We computed the Boltzmann population of anionic Be₆B₁₁ and neutral Be₄B₈ 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₆B₁₁⁻ 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₆B₁₁⁻ 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₆B₁₁⁻ 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₁ 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₄B₈ 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.

Conflict of interest

The authors declare no conflict of interest.