Mapping the Finite-Temperature Behavior of Conformations to Their Potential Energy Barriers: Case Studies on Si6B and Si5B Clusters.

Dynamical simulations of molecules and materials have been the route to understand the rearrangement of atoms within them at different temperatures. Born-Oppenheimer molecular dynamical simulations have further helped to comprehend the reaction dynamics at various finite temperatures. We take a case study of Si6B and Si5B clusters and demonstrate that their finite-temperature behavior is rather mapped to the potential energy surface. The study further brings forth the fact that an accurate description of the dynamics is rather coupled with the accuracy of the method in defining the potential energy surface. A more precise potential energy surface generated through the coupled cluster method is finally used to identify the most accurate description of the potential energy surface and the interconnected finite-temperature behavior.

Introduction

It has been more than a decade since the fact that the physical and chemical properties of aggregates with 3–2000 atoms evolve in a complex and irregular pattern as compared to their bulk counterparts has been established.[1,2] The optical, magnetic, and reactivity patterns in these aggregates have been studied in an far-reaching manner by both experimental and computational methods.[3−5] One of the most notable outcomes of these exhaustive studies is an unexpected thermal stability trend within small-sized molecular clusters of 10–100 atoms. Some of the size-selective clusters of Ga, Al, and Na have been, experimentally, found to have a thermal stability surpassing that of their bulk counterparts in spite of the actuality that the former systems have higher surface to volume ratios.[6−11]

Based on this surprising revelation, a flurry of computational simulations based on Born–Oppenheimer molecular dynamics followed in attempting to identify the origin of this unanticipated finite-temperature behavior.[12−22] Simulations that followed over the years attempted to discern the conformational dynamics followed by the various-sized clusters and the origin of their enhanced thermal stabilities. These computational studies also further revealed that the shape- and size-specific clusters of Au and Sn have enhanced thermal stability as compared to their bulk counterparts, thereby indicating it to be a more general property than thought of among the atomic clusters.[22,23] Along with these newer aspects, detailed simulations also highlighted interesting aspects such as the fluxionality of clusters and/or role of factors such as charge in enhancing the stability of a given conformation at finite temperatures.[24−28] Several studies continue to follow where the thermal stability of a given conformation or alloy is explored using first-principles-based molecular dynamics.[29−35]

The origin of thermal stability of the cluster has been so far attributed to a more generic aspect, viz., its innate molecular structure. How consistent the predictions obtained by various computational methods are, their sensitivity with respect to the method used etc. were so far overseen aspects. In an attempt to identify these aspects, we have essayed upon the fact that it is the potential energy surface that leads to the thermal stability of the finite-sized system. In the following sections, we bring forth this understanding using a simple case of the Si6B atomic cluster.

The choice of case study is based on the fact that the conformational space should not be very large so as to follow it with relative ease. Hence, a Si6B cluster is taken. Recently combined experimental and theoretical studies have delineated five distinct conformations of the Si6B cluster as shown in Figure .[36] The nomenclature used for each of these conformations is mentioned in the figure.

Figure 1

Various Si6B conformations considered in the present study. The nomenclature used for the various conformations is given below. The relative energies of these conformations as obtained from BLYP and PBE exchange correlation functionals are given in Supporting Table ST1. Si6B_0, Si6B_1, Si6B_2, Si6B_3, and Si6B_4 have Cs, C1, C2, C1, and C2 symmetry, respectively.

Results and Discussion

The finite-temperature behavior of the Si6B_0 cluster as obtained with two different basis sets, viz. DZVP and TZVP, for PBE and BLYP functionals is quantified in terms of the root-mean-square bond-length displacements and given in Figure . A notable observation is that the cluster maintains its original Cs conformation up to 900 K for the PBE functional, and this is quantified by the δrms values remaining less than 0.05 up to 900 K for both DZVP and TZVP basis sets. On the other hand, the Cs conformation is retained only up to 500 and 600 K, respectively, for DZVP and TZVP basis sets when the BLYP functional is considered. Coming down to the more intricate features for each method, it is noted that the δrms values for the BLYP functional are quite similar for DZVP and TZVP functionals, with the main variance being the precise temperature at which the values increase; the values for the DZVP basis-based study increase 100 K earlier as compared to the ones obtained from the TZVP basis set. On the other hand, the finite-temperature behavior for the PBE method is identical for both DZVP and TZVP basis sets, with a deviation only at 1100 K. Analysis of the ionic motions reveals that the onset of conformational transitions after 500, 600, and 900 K for BLYP-DZVP, BLYP-TZVP, and PBE-DZVP/PBE-TZVP, respectively, is marked primarily by the presence of Si6B_1 in each method. The variations are corroborated by the mean-square displacement plots at various temperatures given in Supporting Figures SF1–SF4, where it is noted that the mean-square displacement of atoms for the PBE functional is substantially low up to 900 K. It is also perceived from the mean-square displacement plots that the boron atom more or less remains in the center, with the silicon atoms diffusing around it at higher temperatures of 1100, 1300, and 1500 K independent of the functional and basis sets considered for the simulations. Snapshots of BOMD simulation trajectories for Si6B_0 at 1100 K using the TZVP basis and PBE functional are demonstrated in Supporting Figure SF5. In summary, the precise temperature up to which the conformation Si6B_0 remains stable varies between 500 and 900 K depending upon the functional and the basis set used. The atomic displacements on the onset of the conformational distortions are primarily on the surface involving only the silicon atoms. Analysis of ionic motion reveals that after the onset of conformational distortion, the cluster undergoes conformational changes between Si6B_0, Si6B_1, and Si6B_2 between 600/900 and 1100 K. The four conformations for Si6B are observed only beyond 1100 K. It is pertinent to note that a δrms value (commonly referred to as Lindemann measure) of 0.15 is generally considered to be the onset of a disordered state of the cluster.[37,38]

Figure 2

Root-mean-square bond-length fluctuations (δrms) for Si6B_0 as obtained with different functionals and basis sets. The color coding used for different methods (viz. basis set and functional) is given in the inset.

However, a notable point is the variation in the onset temperature of the conformational dynamics. In order to identify the underlying factors contributing to the variation in thermal stability as a function of the computational methodology applied, we have analyzed the electronic properties as obtained by the respective methods that have been tabulated in Table . Analyses of the averaged-out Si–Si bond distance values, presented in Table , show a small correlation with the temperature up to which the Si6B_0 conformation is retained. The average Si–Si bond lengths are smaller for the PBE method, which indicates the higher thermal stability of the Si6B_0 cluster as compared to the BLYP. Similarly, the averaged-out Si–B bond distance values also demonstrate a soft correlation with the thermal stability values of the conformation. The charge on the boron atom as obtained by Löwdin’s partitioning scheme consistently shows a negatively charged boron atom (which means the rest of the silicon atoms have a net positive charge). However, the trend in the values does not correlate to the finite-temperature behavior predicted by the different methods as seen from Table . Frontier orbital analysis of the electronic structure as done by the BLYP-DZVP and PBE-DZVP methods (given in Supporting Figure SF6) indicates that the compositions of the frontier orbitals are more or less similar. Next, we reassess the differences in relative energies between Si6B_0 and the next highest-energy conformation, which is noted to be Si6B_1 for all of the four sets of calculations (see Table ). The relative energy differences between the lowest-energy conformation and the next highest conformation, given in Table , correlate to the thermal stability trends. Based on this observation, we next compute the transition-state barriers for the reaction dynamics from Si6B_0 to Si6B_1, as presented in Figure . The transition-state barriers for the BLYP functional are found to be lower than those obtained from the PBE functional, as seen from Figure . The transition-state barrier between the TZVP- and DZVP-based calculations for BLYP are more notably different as compared to the corresponding values from the PBE method. Thus, the transition-state barrier trends correlate to the thermal stability temperatures as obtained from each method.

Table 1

Electronic Properties of Si6B_0 as Obtained by the Different Methods Studied

		electronic properties
computational method (basis set functional)	thermal stability temperature (K)	average Si–Si bond distance (Å)	average Si–B bond distance (Å)	charge on the boron atom as obtained by the Löwdin partitioning scheme (e)	relative energy of Si6B_1 with respect to Si6B_0 (kcal/mol)
BLYP-DZVP	500	2.34	2.10	–0.90	5.12
BLYP-TZVP	600	2.33	2.09	–1.10	5.46
PBE-DZVP	900	2.31	2.09	–0.90	8.53
PBE-TZVP	900	2.31	2.08	–1.1	8.85

Figure 3

Transition-state barriers (kcal/mol) from Si6B_0 to Si6B_1 as obtained from different basis sets and DFT functionals. The conformational stabilities noted for BLYP-DZVP, BLYP-TZVP, PBE-DZVP, and PBE-TZVP are 500, 600, 900, and 900 K, respectively. The transition-state barriers for Si6B_0 to Si6B_1 as obtained by BLYP-DZVP, BLYP-TZVP, PBE-DZVP, and PBE-TZVP methods are 8.16, 8.54, 12.05, and 11.92 kcal/mol, respectively.

Next, we examine the finite-temperature behavior of Si6B_1, the next high-energy conformation of Si6B, with BYLP and PBE functionals using both DZVP and TZVP basis sets for each functional. The finite-temperature simulations were carried out with Si6B_1 as the starting conformation at 200, 300, 400, 500, 700, 900, 1100, 1300, and 1500 K. The mean-square displacement of atoms as obtained for these simulations is given in Figure . Analysis of the ionic motion of Si6B_1 trajectories reveals that the starting conformation of Si6B_1 undergoes a structural transition to Si6B_0 in each method. Following the transition to Si6B_0, the cluster follows the same path as discussed for each method, i.e., the cluster is retained up to 900 K for the PBE-based studies and at lower temperatures for the BLYP-based methods. The exact temperature at which this transition occurs (quantified in terms of a small peak noted in the root-mean-square bond-length fluctuations (δrms) as seen from Figure ) varies between 300 and 400 K for the different methods. The transition-state barriers from Si6B_1 to Si6B_0 are 3.1 kcal/mol for BLYP-DZVP, BLYP-TZVP, and PBE-TZVP methods and 3.5 kcal/mol for PBE-DZVP. As the transition-state barriers are more or less similar, the transition temperatures for Si6B_1 to Si6B_0 are nearly identical for all of the methods. Following the transition to Si6B_0 from Si6B_1, the δrms values more or less begin to follow similar trajectories as discussed in earlier.

Figure 4

Root-mean-square bond-length fluctuations (δrms) for Si6B_1 as obtained with different functionals and basis sets. The color coding used for different methods (viz. basis set and functional) is given in the inset.

The correlation between the thermal stability of a given conformation and the potential energy barriers to be overcome to enter the next local minima is further established by expanding the study to other basis sets and functionals for the case of Si6B_0. The relative energies between the lowest-energy conformation and next two high-energy conformations as obtained from various functionals and different basis sets are demonstrated in Figure a. The black point in Figure a is indicative of the relative energy of Si6B_0 given by 0.0. The blue and red points are indicative of the next two high-energy conformations. From the values next to the dotted lines in Figure a, it is understood that the relative energy of the first high-energy conformation from Si6B_0 is smaller for both the hybrid functionals, viz., BLYP and B3LYP, where the next high-energy conformation lies between 5.3 and 6.25 kcal/mol. The GGA functional, viz., PBE, as discussed earlier predicts a higher energy difference between the Si6B_0 and Si6B_1 ranging from 8.53 to 8.85 kcal/mol. The meta-GGA functional, viz., M06, predicts Si6B_2 to be the next highest-energy conformer to Si6B_0 for both DZVP and TZVP functionals. As noted from Figure b, the transition-state barrier for Si6B_0 to Si6B_2 as predicted by the meta-GGA functionals is higher than those predicted by the BLYP and B3LYP methods for the DZVP basis set. Based on the examination of the values given by these studies, the Si6B_0 cluster should have lower thermal stability for the B3LYP functional than obtained for PBE, but slightly higher than those obtained for the BLYP functional. Similarly, the meta-GGA functional M06 with DZVP basis set should predict a higher thermal stability for Si6B_0 than all of the methods studied so far. In Figure , we plot the root-mean-square bond-length fluctuations (δrms) for Si6B_0 as obtained with the B3LYP and M06 functionals for the DZVP basis set. It can be observed that the M06 functional shows a higher thermal stability when compared to the B3LYP functional for the DZVP basis. Contrastingly, upon comparing with the δrms plot in Figure , B3LYP is seen to predict a higher thermal stability of Si6B_0 conformation as compared to that of the BLYP method and lower than that of PBE method. As anticipated from the transition-state barrier studies, the thermal stability of Si6B_0 as obtained from the M06 method with the DZVP basis set is the highest, i.e., around 1000 K. Interestingly, the transition from Si6B_1 to Si6B_2 (or Si6B_2 to Si6B_1) appears to be a barrierless transformation as understood from Figure b; this is ratified by the fact that the Si6B cluster undergoes interconversion between Si6B_0, Si6B_1, and Si6B_2 in all methods up to 1100 K. Thus, a clear and distinctive correlation is indicated between the finite-temperature behavior of a conformation and its potential energy surface.

Figure 5

(a) Relative energy of different conformations with respect to the ground-state Si6B_0 conformation as obtained by different methods and basis. The black dot is representative of the relative energy of Si6B_0, which is set as 0.0. The blue and red dots in the graph represent the next two high-energy conformations (Si6B_1 or Si6B_2) as the case may be. These values are also given in the Supporting Table ST2. (b) Transition-state barriers for transition from Si6B_0 to Si6B_1 and Si6B_1 to Si6B_2 as obtained from the different basis sets and DFT functionals.

Figure 6

Root-mean-square bond-length fluctuations (δrms) for Si6B_0 as obtained with the B3LYP and M06 functionals at the DZVP basis set. The color coding used for different methods is given in the inset.

On the other hand, the electronic properties, viz., average Si–Si bond distances, average Si–B bond distances, and charge on the boron atom as obtained from different methods (given in Supporting Figures SF7 and SF8) do not show conclusive correlation to the variations in the conformational stability.

Finally, a more relevant question that may arise out of this discussion is what may be the exact temperature to which the Si6B_0 conformation is stable given the fact that different methods give different transition-state barriers and consequently different conformational stability trends. Hence, we performed single-point coupled cluster theory-based calculations (CCSD) on the final energy conformations obtained from the DZVP-BLYP and DZVP-PBE methods. These values are given in Supporting Table ST2. The coupled cluster calculations reveal that the relative energy between Si6B_0 and Si6B_1 conformations is 8.68 and 9.32 kcal/mol, respectively, and is more closer to the values obtained by the meta-GGA method. Similarly, the transition-state barrier for Si6B_0 to Si6B_1 is 12.11 kcal/mol and from Si6B_1 to Si6B_2 is 5.12 kcal/mol, respectively, for the single-point calculations carried out on the final geometries obtained from BLYP methods. The corresponding values when single-point calculations are carried out on the final optimized geometries of the PBE method are 12.60 and 4.66 kcal/mol, respectively. These values are closer to the values obtained from the meta-GGA method, viz., M06. We would like to point that we use geometries optimized at DFT for the coupled cluster level. While the absolute values are likely to be slightly different on optimization at the CCSD level, we expect that the trend we observe will not change. Finally, we conclude that the thermal stability of the Si6B_0 cluster is around 1000 K based on the above discussions and trends.

The results obtained on Si6B_0 are ratified by performing BOMD simulations using BLYP, PBE, and M06 functionals with the DZVP basis set on the lowest-energy conformation of Si5B36 (referred to as Si5B_0). The δrms values obtained for the Si5B_0 are given in Supporting Figure SF9a. It is clearly noted from the figure that once again the thermal stability of the Si5B_0 as predicted by various functionals follows the order M06 > PBE > BLYP. The thermal stabilities of Si5B_0 as predicted by BLYP, PBE, and M06 functionals are 600, 700, and 800 K, respectively. These thermal stability values corroborate once again the transition-state barriers between Si5B_0 and Si5B_1 (which is the next highest-energy conformation of the Si5B cluster as obtained from all of the three functionals). The transition-state barrier between Si5B_0 and Si5B_1 is calculated to be 7.15, 10.29, and 12.36 kcal/mol, respectively. Thus, in principle significantly, BLYP functionals underestimate the thermal stability (or the transition-state barrier), at least in the case of two studied Si5B clusters. A salient understanding that can be noted by comparing the stability trends between Si5B and Si6B clusters is that while BLYP suggests that Si5B is more stable than Si6B, GGA and meta-GGA functionals suggest that Si6B is thermally stable as compared to Si5B, ratifying the fact that employment of M06 functionals for understanding the thermodynamic properties is more rational.

Conclusions

The present work, through a finite-temperature study on the Si6B cluster, demonstrates that hybrid functionals predict a lower thermal stability for the Si6B cluster ranging between 500 and 700 K. GGA functionals and meta-GGA functionals describe that the Si6B ground-state conformation is stable up to at least 1000 K. The study also identifies meta-GGA functionals to have a more closer description of the potential energy surface by the coupled cluster method.

Computational Details

All of the five conformations are optimized using the density functional theory-based methodology as implemented in deMon 2K code.[39] The optimizations were carried out using Becke exchange with Lee, Yang and Parr correlation functional (BLYP)[40,41] and Perdew–Burke–Ernzerhof (PBE)[42] exchange and correlation functionals. The basis set was also varied to be DZVP and TZVP for both the functionals. No additional polarization functions are added. The GEN-A2* auxiliary functions are set to fit the charge density.[43] The convergence of the geometries is based on gradient and displacement criteria with a threshold value of 10–3 au, and the criterion for convergence of SCF cycles was set to 10–6 au. Only the lowest spin state is considered for all of the Si6B clusters in the study.

Two of the clusters have odd number of electrons. For the considered multiplicity of 2, Restricted Open shell Kohn Sham (ROKS) formalism is used for the calculations. Analyses of relative stabilities bring out a consistent lowest-energy trend for the Si6B_0 conformation, which is quantified in Supporting Table ST1. In other words, all of the four calculations describe Si6B_0 to be the lowest-energy conformation. Following the geometry optimization calculations, finite-temperature simulations are carried out on the Si6B_0 cluster using Born–Oppenheimer Molecular Dynamics (BOMD) as implemented in deMon 2K for an NVT ensemble. We hold the total angular momentum of the cluster at zero by suppressing the cluster rotation. The simulations are carried out between 200 and 1500 K at temperatures of 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1300 and 1500 K. For each temperature, the cluster is equilibrated for a time period of 70 ps followed by a simulation time of 80 ps. The temperature of the cluster is maintained using the Nosé–Hoover thermostat[44,45] as implemented in deMon 2K. The nuclear positions are updated with the velocity Verlet algorithm with a time step of 1 fs. The atomic positions and bond-length punctuations of the atoms are analyzed using traditional parameters such as root-mean-square bond-length fluctuations (δrms) and the mean-square displacements of atoms (MSDs). δrms is defined aswhere N is the number of particles in the system, r is the distance between the i-th and j-th particles in the system, and ⟨...⟩ denotes a time average over the entire trajectory. The MSD of an individual atom is defined aswhere N is the number of atoms in the system and R(t0) is the instantaneous position of atom i at t0 and R(t0 + t) is the corresponding position of atom i after a time interval t. Thus, we average over M different time origins t0 spanning the entire trajectory. The average interval between consecutive t0 values was taken to be about 0.0125 ps. The MSD indicates the displacement of atoms in the cluster as a function of time. In the solid-like region, all atoms perform oscillatory motion about fixed positions, resulting in a negligible MSD of individual atoms. In a liquid-like state, on the other hand, atoms diffuse throughout the cluster and the MSD eventually reaches a saturated value of the order of the square of the cluster radius.

Following the analysis of the finite-temperature behavior of the Si6B_0 cluster, we performed BOMD simulations on Si6B_1 and Si6B_2 clusters also in order to validate our surmises. The simulations on these clusters were carried out at the aforementioned temperatures and conditions for Si6B_0.

Single-point calculations coupled cluster theory-based calculations (CCSD) were carried out on the final energy conformations obtained from the DZVP-BLYP and DZVP-PBE methods using the domain-based local pair natural orbital (DLPNO)[46] approach with the correlation consistent triple zeta (cc-pVTZ) basis set. These calculations are performed using orca 5.0 software package.[47]