Relative Populations and IR Spectra of Cu38 Cluster at Finite Temperature Based on DFT and Statistical Thermodynamics Calculations.

The relative populations of Cu38 isomers depend to a great extent on the temperature. Density functional theory and nanothermodynamics can be combined to compute the geometrical optimization of isomers and their spectroscopic properties in an approximate manner. In this article, we investigate entropy-driven isomer distributions of Cu38 clusters and the effect of temperature on their IR spectra. An extensive, systematic global search is performed on the potential and free energy surfaces of Cu38 using a two-stage strategy to identify the lowest-energy structure and its low-energy neighbors. The effects of temperature on the populations and IR spectra are considered via Boltzmann factors. The computed IR spectrum of each isomer is multiplied by its corresponding Boltzmann weight at finite temperature. Then, they are summed together to produce a final temperature-dependent, Boltzmann-weighted spectrum. Our results show that the disordered structure dominates at high temperatures and the overall Boltzmann-weighted spectrum is composed of a mixture of spectra from several individual isomers.

Introduction

Nanoclusters are of great interest since they allow us to study the transition from free atoms to bulk condensed systems (Wilcoxon and Abrams, 2006) by analyzing the size-dependent evolution of their properties (Ferrando et al., 2008). In particular, noble-metal nanoclusters (NMCs) have attracted attention in many fields of science due to their interesting plasmonic, catalytic (Mathew and Pradeep, 2014; Inwati et al., 2018), and photophysical properties at the nanoscale (Xavier et al., 2012). Specifically, Cu nanoclusters embedded in a dielectric matrix have attracted attention because of their tunable longitudinal surface plasmon resonance characteristics (Inwati et al., 2018). Copper is cheaper than gold and silver and has high photosensitivity, high thermal and electric conductivities, and optical properties (Zhang et al., 2019) that make it a good candidate for nanodevices (Liu et al., 2015) and nanoelectronics development (Jena and Castleman, 2006). Particularly, Cu38 has attracted attention because it has “magic number” structures (Tran and Johnston, 2009), which are defined in terms of geometric and energetic factors and related to the closing of electronic shells (Baletto et al., 2004), much like small sodium clusters (de Heer, 1993). The magicity of the Cu38 cluster is due only to energetic considerations (Doye and Wales, 1998; Baletto et al., 2004). In contrast, in small, packed barium clusters with magic numbers, stability is dominated by geometric rather than electronic effects (de Heer, 1993). It is believed that magic structures are the global minimum energy structures on the potential energy surface, and thus reflect the molecular properties of the system (Baletto et al., 2004).

From the experimental point of view, the Cu38 cluster has been widely studied via photoelectron spectroscopy (PES) in order to extract the electronic gap of the anionic Cu38 cluster resulting in a semiconductor with a 0.33 eV electronic gap (Pettiette et al., 1988; Zhang et al., 2019). PES has also been used to study the anionic Cu38 cluster, inferring that its putative global minimum should be an oblate structure instead of a highly symmetric structure (Kostko et al., 2005; Zhang et al., 2019), despite computational results for 38-noble metal atom clusters that frequently find a highly symmetric (cuboctahedral) structure (Fujima and Yamaguchi, 1989; Doye and Wales, 1998; Kostko et al., 2005).

From the theoretical point of view, several density functional theory (DFT)-based studies have been carried out, in which, for example, the thermodynamic properties of the Cu38 cluster have been reported (Taylor et al., 2008). The transition states and reaction energies of the water gas shift reaction on a Cu38 cluster and Cu slab have also been studied via DFT computations (Qi et al., 2020). In particular, the high-symmetry octahedral structure was reported as the lowest energy structure (Takagi et al., 2017) using the PW91 functional (Perdew et al., 1992), a plane-wave basis set, and the pseudopotential approximation (Itoh et al., 2009). The Cu38 cluster has also been investigated using a hybrid strategy (Hijazi and Park, 2010); in which the embedded atom potential method followed by DFT computations with the PBE functional and pseudopotential approximation was used. The authors reported a putative global minimum structure with octahedral (OH) symmetry. The second most stable structure was the incomplete-Mackay icosahedron (IMI) located 0.26 eV above the putative global minimum.

There has been some discussion about the lowest energy structure of the Cu38 cluster. Searches for the lowest-energy structure that employed many-body potentials identified a cuboctahedral structure (Doye and Wales, 1998; Grigoryan et al., 2005). Several authors used an empirical potential-energy function containing two-body atomic interactions and found that fivefold symmetry appears to be the putative global minimum in the Cu38 cluster (Erkoç, 1994; Erkoç and Shaltaf, 1999). In contrast, previous studies reported the cuboctahedron structure as the putative global minimum (Grigoryan et al., 2006) by employing empirical many-body Gupta and Sutton-Chen potentials; some other works also consider the Cu38 octahedron cluster to be the putative global minimum (Itoh et al., 2009; Hijazi and Park, 2010; Takagi et al., 2017), but several others found that the Cu38 cluster with truncated octahedron geometry is energetically more stable than other configurations (Darby et al., 2002; Itoh et al., 2009; Hijazi and Park, 2010; Núnez and Johnston, 2010; Park and Hijazi, 2012; Fernández et al., 2015; Zhao et al., 2017). Some of us pointed out that the energies computed via different methods such as DFT, second-order Møllet-Plesset approach (MP2), and Coupled cluster with single-double and perturbative triple (CCSDT) excitations yield different energetic ordering results (de la Puente et al., 1997; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Castillo-Quevedo et al., 2021). In the case of DFT, the functional and basis set employed, the zero-point energy correction (ZPE), the dispersion energy, and other parameters can change the energetic ordering of the low energy structures (Castillo-Quevedo et al., 2021). Moreover, practical molecular systems and materials must be studied at warm temperatures, so their molecular properties at finite temperature are dominated by Boltzmann distributions of isomers (Baletto and Ferrando, 2005; Li and Truhlar, 2014; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Castillo-Quevedo et al., 2021), and those of the associated materials are statistical averages over the ensemble of conformations (Mendoza-Wilson et al., 2020; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). Total energy computations using DFT methodology are typically carried out at absolute zero temperature, although the thermal properties of the inhomogeneous electron gas in the mid-1960s were studied (Mermin, 1965). Recently, DFT was extended to finite temperature (Pittalis et al., 2011; Gonis and Däne, 2018; Gazquezet al., 2019), but as far as we know, it has not been implemented in any software.

There are cases where the global minimum structure ceases to be the most likely at high temperatures, so other structures prevail. For instance, in small Ag clusters, the temperature leads to the transition from the initial Face-Centered Cubic (FCC) phase to other structures (Redel et al., 2015), thus temperature promotes face changes in materials. Interestingly, the molecular system minimizes the Gibbs free energy at temperatures other than zero and maximizes the entropy (Buelna-Garcia et al., 2021). Although the search for global and local minima is useful in understanding reactivities and catalytic efficiencies, such studies mostly neglect temperature-dependent entropic contributions to free energy when the temperature increases. Taking temperature into account requires dealing with nanothermodynamics (Hill, 1962; Baletto and Ferrando, 2005; Li et al., 2007; Li and Truhlar, 2014; Grigoryan and Springborg, 2019; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). The thermodynamics of clusters have been studied using various tools (Wales, 1996; Li et al., 2007; Li and Truhlar, 2014; Calvo, 2015; Buelna-Garcia et al., 2021) such as in the molecular-dynamics simulations of boron clusters (Martinez-Guajardo et al., 2015) and Cu38 clusters (Zhang et al., 2019).

Cluster properties depend heavily on the cluster structure, size, composition, and temperature. Therefore, the first step to understanding their molecular properties is to elucidate the lowest energy structure and its isomers close in energy (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Baletto and Ferrando, 2005; Darby et al., 2002; Ohno and Maeda, 2006), which is a complex task due to several factors (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). The second step relies on spectroscopy, which gives insight into the structure and has been proposed as a way of detecting structural transformations within clusters. The influence of temperature on Infrared spectroscopy (IR) has been computed before for a variety of clusters (Even et al., 1989; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). The present paper uses the statistical formulation of thermodynamics and nanothermodynamics (Li et al., 2007; Li and Truhlar, 2014; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021) to compute the thermodynamic properties of the neutral Cu38 cluster, define its putative global minimum at finite temperature, compute the relative populations among the isomers, and the IR spectra as Boltzmann-weighted spectral sums of individual spectra. Our findings show that an amorphous structure strongly dominates the putative global minimum at high temperatures, whereas the truncated octahedron dominates at low temperatures. The remainder of the manuscript is organized as follows: “Free Energy Surface Exploration Method and Computational Details” gives the computational details and a brief overview of the theory and algorithms. The results and discussion are presented in “Results and Discussion”. This includes the putative global minimum at room temperature, relative populations at temperatures from 20 to 1500 K, and IR spectra as functions of the temperature. Conclusions are given in “Conclusion”.

Free Energy Surface Exploration Method and Computational Details

The putative global minimum is determined by the enthalpy (at zero temperature) or the Gibbs free energy (at temperatures other than zero). A simple analysis of the Gibbs free energy given by ∆G = ∆H − ∆ST leads to the conclusion that the entropy must be maximized in order to minimize the Gibbs free energy (Sutton and Levchenko, 2020; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). From the theoretical point of view, and in order to understand molecular properties at finite temperature, the lowest Gibbs free energy structure (or the structure with the largest entropy), as well as all structures close in energy to the lowest energy structure (or all high-entropy structures close in entropy to the structure with the highest entropy) must be known (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021), considering that experiments are performed at finite temperature.

The search for the global minimum in atomic clusters is a complex task due to the number of possible combinations, which grows exponentially with the number of atoms, leading to a combinatorial explosion problem, among others (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). Despite the difficulty of this task, several algorithms have been successfully employed in a targeted way to explore the potential and free energy surfaces. These are coupled to a local optimizer of any electronic structure package. Examples include the Ab initio Random Structure Searching approach (Pickard and Needs, 2011), simulated annealing (Kirkpatrick et al., 1983; Metropolis et al., 1953; Xiang and Gong, 2000; Xiang et al., 2013; Vlachos et al., 1993; Granville et al., 1994), the kick methodology (Pan et al., 2014; Cui et al., 2015; Vargas-Caamal et al., 2016b,a; Cui et al., 2017; Vargas-Caamal et al., 2015; Flórez et al., 2016; Ravell et al., 2018; Hadad et al., 2014; Saunders, 2004, 1987; Grande-Aztatzi et al., 2014; Mendoza-Wilson et al., 2022), and genetic algorithms (Guo et al., 2017; Dong et al., 2018; Mondal et al., 2016; Ravell et al., 2018; Grande-Aztatzi et al., 2014; Rodríguez-Kessler et al., 2017; Alexandrova et al., 2004; Buelna-Garcia et al., 2021). Global minimum structure searches at the DFT level are too computationally expensive to be applied to intermediate and large cluster sizes since, as mentioned earlier, the number of candidates increases exponentially with the number of atoms. In this paper, we use a two-stage procedure to explore the potential energy surface efficiently. In the primary stage, we perform a global search using an empirical methodology. The Gupta interaction potential is used to describe the Cu-Cu interactions with default parameters taken from Refs. (Cleri and Rosato, 1993; Wilson and Johnston, 2002). It is coupled to the basin hopping global optimization algorithm implemented in Python code and part of the GALGOSON global search code (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). In the second stage, all of the lowest energy structures from the primary stage are symmetrized, which is followed by geometry optimization at the DFT level, using the Gaussian-suite code (Frisch et al., 2009). The calculations employ two exchange-correlation functionals, B3PW91 and PBE, and two basis sets, def2SVP and LANL2DZ, with and without considering the D3 version of Grimme’s dispersion corrections (Grimme et al., 2010), as implemented in the Gaussian 09 code (Frisch et al., 2009). Becke’s hybrid three-parameter (Becke 1993, 1988) exchange-correlation functional in combination with the Perdew and Wang GGA functional PW91 (Perdew et al., 1992; Perdew and Wang, 1992) is known as the B3PW91 exchange-correlation functional. The B3PW91 has been employed in other studies of reactivity in copper clusters, where it has provided good performance (Fernández et al., 2015). The PBE exchange-correlation functional (Perdew et al., 1996) has shown good performance with regarding thermochemical properties (del Campo et al., 2012). The LANL2DZ basis set (Dunning and Hay, 1977) has been used in previous computational studies of copper-based molecular properties with very good agreement with experimental values (Legge et al., 2001). In a previous DFT study, the def2-SVP basis set (Weigend and Ahlrichs, 2005) provided good results in the computation of the Cu-metal ligand bond lengths (Niu et al., 2014). The true minimum energy structures are validated via vibrational analysis.

Thermochemical Properties

All of the thermodynamic properties of an ensemble of molecules can be derived from the molecular partition function (Takadaet al., 2018; Dzib et al., 2019; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021). Similarly, the wave function contains all information about a molecular system (Buelna-Garcia et al., 2021). Previous theoretical studies used the partition function to compute thermodynamic properties of Cun clusters (n = 2, 150) as a function of the temperature and showed that the magic number structures are temperature dependent (Li et al., 2007; Grigoryan and Springborg 2019). The thermodynamics of unsupported neutral Aln (2 < n < 65) clusters have also been investigated by evaluating vibrational partition functions. They reported that the dominant cluster structure is temperature dependent (Li and Truhlar, 2014). Aditionally, the atomistic thermodynamics framework has been used to predict material behaviors at realistic temperatures (Sutton and Levchenko, 2020). More recently, the partition function was used by some of us to compute the temperature-dependent relative populations and IR spectra of neutral Be4B8 and anionic Be6B11 clusters (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021), and a similar procedure was employed to compute reaction rate constants in two representative hydrogen abstraction reactions (Dzib et al., 2019). Regarding the temperature-dependent entropic contributions, the [Fe(pmea)(NCS)2] complex was studied by Brehm et al. (Brehm et al., 2006). In this study, the thermodynamic properties are computed using the partition function Q given in Eq. 1 under the rigid rotor, harmonic oscillator, Born-Oppenheimer, ideal gas, and particle-in-a-box approximations.

In Eq. 1, gi is the degeneracy factor, KB is the Boltzmann constant, T is the temperature, and −∆Ei is the total energy of a cluster (Buelna-Garcia et al., 2021; Dzib et al., 2019; McQuarrie, 1975). Within Born Oppenheimer and rigid rotor harmonic oscillator approximations, the partition function Q(T) is factorized into electronic, translational, vibrational, and rotational contributions given by Eq. 2

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

TABLE 1

Contributions to the partition function.

Contribution	Partition function
Translational	q_trans=(2πmk_BTh^2)^32k_BTP
Rotational linear	qrot1=TσΘ_rot,Θ_rot=|h^22Ik_B
Rotational nonlinear	qrotnl=π^12σ[T^32(Θ_rotAΘ_rotBΘ_rotC)],Θ_rotj=h^22I_jk_B,j=A,B,C
Vibrational	qvibpol=∏i=1n01[a]e^−Θ_vib_i/2T1−e^−Θ_vib_i/T,Θ_vib_i=hv_ik_B
Electronic	q_elec=ω_0

We considered that the energy gap between the first and higher excited states is greater than ; consequently, the electronic partition function is given by . , , and were used to compute the internal energy (U), and entropy (S) contributions given in Table 2.

TABLE 2

Contributions to internal energy and entropy.

	Internal energy	Entropy
Translational	U_trans=32RT	S_trans=R(ln q_trans+52)
Rotational linear	Urot1=RT	Srotl=R(ln qrotl+1)
Rotational nonlinear	Urotnl=32RT	Srotnl=R(ln qrotnl+32)
Vibrational	Uvibpol=R∑in0[a]Θ_vib_i(12+1e^Θ_vib_i/T−1)	Svibpol=R∑in_o[Θ_vib_i/Te^Θ_vib_i/T−1−ln(1)−e^−Θ_vib_i/T]
	Θ_vib_i=hv_ik_B
Electronic	U_elec=0	S_elec=R ln q_elec

Equations shows in Tables 1 and 2 are the same as those employed in previous studies (Li et al., 2007; Grimme, 2012; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Dzib et al., 2019) and any standard thermodynamics textbook (McQuarrie, 1975; Hill, 1986).

The vibrational frequencies are computed employing Gaussian code. The Gibbs free energy (G) and the enthalpy (H) are computed employing Eqs. 3 and 4, respectively. In these Equations R is the ideal gas constant, n is the amount of substance, and T is the absoulte temperature.

To compute the probability of the occurrence of one particular Cu38 cluster in a Boltzmann ensemble at thermal equilibrium as a function of temperature, we employed the probability of occurrence (Shortle, 2003; Li et al., 2007; Grimme, 2012; Bhattacharya et al., 2017; Schebarchov et al., 2018; Dzib et al., 2019; Goldsmith et al., 2019; Grigoryan and Springborg, 2019; Mendoza-Wilson et al., 2020; Bhumla et al., 2021; Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021) given by Eq. 5: where , is the Boltzmann constant, T is the temperature, and k is the Gibbs free energy of the kth isomer. We point out that Gibbs free energies must be corrected based on the symmetry. Our previous work showed that the contribution of the rotational entropy to the Gibbs free energy depends on the symmetry, varies linearly with the temperature, and can be significant (Buelna-Garcia et al., 2021). Eq. 5 is restricted so that the sum of all occurrence probabilities at fixed temperature T, i.e., the sum of all Pi (T), is equal to 1. This is given by Eq. 6

In this study, the Boltzmann-weighted IR spectrum at a finite temperature is given by Eq. 7: where n is the total number of clusters in the ensemble, IRi is the spectrum of the ith isomer at temperature T = 0, and Pi (T) is the probability of the ith isomer, which is given by Eq. 5. We use the Boltzmann-optics-full-Ader code (BOFA) to compute the occurrence probability and the IR spectra (Buelna-Garcia et al., 2021).

Results and Discussion

Low-Energy Structures

The ball and stick models shown in Figure 1 depict the lowest-energy structures of neutral Cu38 clusters and some competing isomers. We use the B3PW91/def2SVP level of theory and consider the Grimme (DFT-D3) dispersion pairwise correction (Grimme et al., 2010), at room temperature and at 1 atm pressure. We found that the tetrakaidecahedron is the lowest energy structure. It has fourteen faces: six equivalent square FCC(100) faces and eight equivalent hexagons. This shape is obtained when cutting the corners off a 3D diamond shape. It is an FCC-like truncated octahedron (TO). The calculated structure belongs to the C1 symmetry point group and to the 1A electronic ground state. Its lowest IR active vibration frequency is 32.57 cm−1 and it is a semiconductor with an electronic gap of 0.623 eV. Previous works regarding the exploration of the potential energy surface of Cu38 using genetic algorithms via the Gupta potential have often found highly symmetric TO structures (Kostko et al., 2005; Darby et al., 2002), which have also been reported using the Sutton-Chen potential with Monte Carlo simulations (Zhao et al., 2017). The optimized Cu-Cu bond length is 2.4670 Å, which is in good agreement with the reported bond length in the Cu-Cu dimer via DFT calculations (2.248 Å) (Kabir et al., 2004; Guvelioglu et al., 2006) and is consistent with the experimental Cu-Cu bonding distance of 2.22 Å (Kabir et al., 2004). Our computed TO structure diameter is 7.8 Å, which is in good agreement with the 8 Å reported in previous DFT calculations (Guvelioglu et al., 2006).

FIGURE 1

(Color online) Optimized geometries in front and side views of neutral Cu38 clusters at the PBE-D3/def2-SVP level of theory. The calculations consider the D3 version of Grimme’s dispersion corrections and the ZPE correction energy. The first letter is the isomer label, relative Gibbs free energies are in kcal/mol (in round parenthesis) at 298.15 K, electronic groups and symmetry point groups are [in square parenthesis], and the probability of occurrence is shown (in round, red parenthesis) at 298.15 K.

The structure with the second-lowest energy lies at 0.16 kcal/mol at 298.15 K above the putative global minimum and is also a TO structure with C1 symmetry and 1A electronic ground state. Its lowest IR active vibration frequency is 32.13 cm−1 and it is a semiconductor with an electronic gap of 0.623 eV, thus, it is fairly similar to the putative global minimum. The next structure is slightly higher in energy. It is located 1.38 kcal/mol above the putative global minimum and is also a TO structure. It has D4h symmetry and 1A1g electronic ground state. Its lowest IR active vibration frequency is 33.44 cm−1. We also explore TO structures, starting with highly-symmetric OH and TH. After geometry optimization without constraints, the OH and TH symmetries become C1 and D4h. The perfect OH symmetry could be deformed due to the Jahn–Teller effect (Kabir et al., 2004; Guvelioglu et al., 2006). This effect must be considered when calculating the total energy (Zlatar et al., 2010; Zhang et al., 2018) because the computed optical properties could change due to relative population at finite temperatures (Opik and Pryce, 1957). In one of our recent works, we clarified the origin of Gibbs free energy differences between two similar structures with different symmetry point group due to rotational entropy, specifically the RTln(σ) factor (Buelna-García et al., 2021). In this work, the energy difference of 0.16 kcal/mol between the two isomers depicted in Figure 1A with C1 symmetries and a root-mean-square deviation (RMSD) of 0.08 is due to the Jahn-Teller effect. The structure located 1.38 kcal/mol above the putative global minimum with D4h symmetry exists due to rotational entropy. The next structure, shown in Figure 1D is located 5.65 kcal/mol above the putative global minimum and has point group symmetry C1 and electronic ground state 1A. Its lowest IR active vibration frequency is 24.16 cm−1. It is a distorted-structure semiconductor with an electronic gap of 1.0 eV, calculated Cu-Cu bond distance of 2.50 Å and molecular diameter of 9.1 Å. Its Cu-Cu bond distance and diameter are slightly larger than the global minimum. This structure possesses the smallest relative ZPE energy, as shown in Supplementary Figure S1, and the smallest vibrational frequency mode of all the isomers. The next two higher energy structures are shown in Figures 1E,F reach 5.8 kcal/mol. They are IMIs with C1 and C5V point group symmetries and electronic ground states 1 and 1A1, respectively. In both cases, the molecular diameter is 8.54 Å, the electronic gap is 0.97 eV, and the Cu-Cu bonding distance is 2.47 Å. Other, higher energy structures are shown in Figures 1G–J. These do not contribute to any molecular properties at zero and finite temperatures. Supplementary Figure S2 depicts lowest-energy-structure screening at the B3PW91/def2SVP level without considering the Grimme D3 atom-pairwise correction. The lowest energy structure is the IMI structure with point group symmetry C1 and electronic ground state 1A. The molecular diameter is 8.69 Å, which is slightly larger than that of the TO structure (7.8 Å). The average bond distance is 2.50 Å. We find the IMI structure to be the most stable at the PBE-D3/Def2SVP level. In contrast, at the PBE-D3/LANL2DZ level, we find the TO structure to be the putative global minimum. A complete description of the structures located at higher energies is presented in the Supplementary Material. For the Cu38 clusters, we point out that the energetic ordering of the isomers, the energy gaps among the isomers, and the putative global minimum interchange when we consider the dispersion interactions.

Energetics

Employing different methods to compute energies yields different results due to differences in the functional and basis set (Yanai et al., 2004), and the resulting energetic ordering changes (de la Puente et al., 1997; Buelna-Garcia et al., 2021). A comparison of total energy, among isomers, computed with two different exchange-correlation functionals and two basis sets, one which does and one which does not consider the Grimme D3 dispersion, is shown in Table 3. The optimizations performed at the B3PW91/PBE-def2TZVP level that consider the dispersion yield the same type of lowest-energy equilibrium geometries and similar energetic isomer ordering. From the energetic point of view, the inclusion of dispersion is more important than the type of functional or basis set. The first line of Table 3 shows the relative Gibbs free energies computed at the B3PW91-D3/def2TZVP level of theory. The isomer labeled ib in Table 3 and depicted in Figure 1B is located 0.16 kcal/mol above the putative global minimum. In contrast, the second line of Table 3 shows the relative Gibbs free energies computed at the B3PW91/def2TZVP level of theory. Here, the isomer ib in Table 3 that is depicted in Figure 1B is located 0.95 kcal/mol above the putative global minimum. For isomer with label ib, the inclusion of dispersion decreases the Gibbs free energy relative to the putative global minimum (from 0.95 to 0.16 kcal/mol). For isomer with label ic, considering dispersion decreases the Gibbs free energy relative to the putative global minimum from 2.0 to 1.38 kcal/mol. In contrast, the inclusion of dispersion for isomer with label id, the relative Gibbs free energy increases from 2.4 to 5.65 kcal/mol. In summary, considering dispersion reduces the Gibbs free energies of the lowest-energy structures where the Boltzmann factors are not zero. An overall comparison of free energies computed using functional B3PW91, the second line of Table 3, and PBE in four-line in Table 3, shows a reduction in the relative Gibbs free energies when the PBE functional is employed. The LANL2DZ basis set increases the relative Gibbs free energies of the low-energy isomers, as shown by comparing line six and one of Table 3.

TABLE 3

A comparison of the energetic isomer ordering as determined using the B3PW91/Def2SVP, PBE/Def2SVP, and PBE/LANL2DZ levels of theory. The Gibbs free energy is computed at room temperature. The electronic energy includes the ZPE energy correction.

Level of theory	Isomers (energy kcal/mol)
	Energy	ia	ib	ic	id	ie	if	ig	ih	ii
B B3PW91-D3/def2SVP	ΔG	0.0	0.16	1.38	5.65	5.79	5.81	6.76	8.85	9.81
	ε_0+ε_ZPE	0.	0.09	0.0	5.01	5.01	5.01	5.01	5.01	8.17
	ε_0	0.05	0.0	0.10	8.76	4.89	4.89	4.88	4.88	4.88
B B3PW91/def2SVP	ΔG	0.0	0.95	2.0	2.40	2.73	2.91	2.94	3.28	3.32
	ε_0+ε_ZPE	0.0	0.0	2.14	3.29	6.52	3.84	2.14	3.27	3.84
	ε_0	0.0	0.0	2.09	3.27	6.13	3.76	2.09	3.28	3.76
B PBE-D3/def2SVP	ΔG	0.0	0.86	0.92	5.02	7.23	7.59	7.81	8.88	12.27
	ε_0+ε_ZPE	0.0	0.89	0.0	8.70	7.94	7.92	8.61	7.94	14.24
	ε_0	0.0	0.90	0.0	9.14	7.92	7.93	8.84	7.94	12.47
P PBE/def2SVP	ΔG	0.0	0.34	0.92	1.37	1.38	1.77	2.82	5.79	8.70
	ε_0+ε_ZPE	0.0	0.0	0.37	0.41	1.77	1.77	1.77	9.08	15.75
	ε_0	0.0	0.0	0.38	0.39	1.73	1.73	1.74	9.47	9.86
P PBE-D3/LANL2DZ	ΔG	0.0	1.02	3.37	3.46	8.63	9.11	9.59	9.62	9.74
	ε_0+ε_ZPE	2.03	0.0	2.12	2.14	8.31	8.77	8.73	9.34	8.71
	ε_0	1.97	0.0	2.01	2.01	8.08	8.60	8.52	9.12	8.45
P PBE/LANL2DZP	ΔG	0.0	2.15	2.87	3.03	3.26	4.31	8.89	9.66	9.85
	ε_0+ε_ZPE	0.97	0.0	1.46	1.38	1.64	1.68	7.98	9.19	9.58
	ε_0	0.86	0.0	1.27	1.11	1.48	1.42	8.03	8.99	9

Occurrence Probabilities at Finite Temperature

Figure 2 shows the relative populations of neutral Cu38 clusters computed at different levels of theory for temperatures ranging from 20 to 1500 K. In order to gain insight into the effects of dispersion on the relative population, this is computed with and without the D3 Grimme dispersion. For ease of comparison, the results are displayed in side-by-side plots in Figure 2.

FIGURE 2

(Color online) The occurrence probabilities for temperatures ranging from 20 to 1500 K at six different levels of theory: (A) B3PW91-D3/def2-SVP, (B) B3PW91-D3/def2-SVP, (C) PBE/def2-SVP, (D) PBE-D3/def2-SVP, (E) PBE/LANL2DZ, and (F) PBE-D3/LANL2DZ. The cases (B,D), and (F) are computed considering the D3 Grimme dispersion. In all cases, the effect of the dispersion on the solid-solid transformation point in the temperature scale is large for the Cu38 cluster. At hot temperatures the dominant structure is an amorphous geometry depicted in Figure 1D, whereas, the TO structure depicted in Figure 1A is the strongly dominant structure at cold temperatures and at the B3PW91-D3/def2-SVP level of theory.

The occurrence probability of the TO structure with C1 symmetry at the B3PW91/def2-SVP level of theory is depicted in the solid, black line in Figure 2A. This occurrence probability of the TO structure strongly dominates from 0 to 300 K; thus, all of the molecular properties in this range of temperature are due only to this structure, and the occurrence probability starts to decay exponentially just before 300 K and nearly disappears at 900 K. The probability of finding the amorphous structure with point group symmetry C1 is depicted in the solid, violet line in Figure 2A. This probability starts to grow exponentially and becomes dominant at a temperature between 512.8 and 900 K. The TO and the amorphous structure co-exist at a solid-solid transition temperature of 512.8 K. The effect of dispersion can be seen Figure 2B. The relative population is computed at the B3PW91-D3/def2-SVP level of theory. The dispersion effect is dramatic; the solid-solid transformation point shifts from 512.8 to 824 K, an increase of 160%. In Figure 2B, one can see that the molecular properties below 600 K are due to only the TO structure. The probability of finding the amorphous structure, depicted via the solid, green line in Figure 2B, starts to increase exponentially just before 600 K. The TO and the amorphous structure co-exist at 824 K. The probability of finding the TO structure starts to decay exponentially at 600 K and is still around 20% at 900 K. The occurrence probabilities of various Cu38 isomers at the PBE/def2-SVP level of theory are displayed in Figure 2C. The dominant putative global minimum structure at T = 0 is the inverted incomplete-Mackay icosahedron (IIMI) structure depicted in Figure 3A with C1 symmetry.

FIGURE 3

(Color online) The inverted incomplete-Mackay icosahedron (IIMI) is labeled 1 and has C1 symmetry. The incomplete-Mackay icosahedron (IMI) is labeled 2, with symmetry Cs, and is located 0.34 kcal/mol energy above the putative minimum global at 298.15 K. The yellow, red, and blue colored spheres represent copper atoms. The IIMI structure is the result of interchanging the red Cu atom depicted in the IMI structure to the position of the blue atom in the IIMI structure. The IMI structure has been reported in reference Ref Zhang et al. (2019). as the low-energy structure. The Highest Occupied Molecular Orbital (HOMO)—Lowest Unoccupied Molecular Orbital (LUMO) gap of the IMI structure is 0.24 eV (0.356 eV reported in previous DFT studies (Zhang et al., 2019). In contrast, the HOMO-LUMO gap for the IIMI structure is 0.30 eV. This, provides a plausible explanation for the higher energetic stability of why the IIMI structure is energetically more stable.

The probability of finding the IIMI structure is shown via the black, solid line in Figure 2C. This probability decays almost linearly until a temperature of 1000 K, where the probability of occurrence almost disappears. At the solid-solid transformation point (759.7 K), the IIMI structure co-exists with an amorphous structure. The probability of finding the amorphous structure starts to increase at 600 K and starts to dominate heavily as the putative global minimum above the solid-solid transformation point. The probability of finding the IMI structure is depicted via the red, solid line in Figure 2C. The probability is maximized (30%) at room temperature. Interestingly, the probabilities of the IIMI and IMI structures do not cross at low temperatures. Previous work reported that the IMI structure can be highly competitive at finite temperature (Zhang et al., 2019). Still, our findings show that the amorphous structure with C1 symmetry is quite dominant at high temperatures, whereas the IIMI structure is strongly dominant at low temperatures.

For ease of comparison, Figure 3 displays the IIMI and the IMI structures side by side. The IIMI structure dominates at low temperatures. The dispersion effect shifts the solid-solid transformation point down from 759.7 to 654 K, as shown in Figure 2D. The probability of finding the IIMI structure is depicted via the black, solid line as a function of temperature. The probability decays approximately linearly from 50 to 500 K; after that, it decays exponentially until 900 K, where it disappears. At around 400 K, the probability of finding the amorphous structure (depicted via the green, solid line Figure 2D starts to grow exponentially. At 654 K, it co-exists with the IIMI structure. Above 654 K, the amorphous structure becomes energetically favorable.

IR Spectra at Finite Temperature

The properties observed in a molecule are statistical averages over the ensemble of geometrical conformations or isomers accessible to the cluster. Thus, the molecular properties are governed by the Boltzmann distributions of the isomers, which can change significantly with the temperature, primarily due to entropic effects (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Li et al., 2007). The many soft vibrational modes that the clusters possess are the major contributions to the entropy. The IR spectrum is related to vibrations or rotations that alter the dipole moment and it is observed in molecules with a dipole moment. The IR spectrum is also related to the curvature of the relationship between the potential and the interatomic distance. Complete information regarding molecular vibrations allows us to analyze catalytic chemical reactions (Tinnemans et al., 2006; Brandhorst et al., 2006; Hashimoto et al., 2019). IR spectra are used to identify functional groups and chemical bond information. However, assigning IR bands to vibrational molecular modes in measured spectra can be difficult and requires DFT calculations; as mentioned earlier, the temperature is not considered in these computations and discrepancies between experimental and computed IR spectra can result from finite temperatures, anharmonic effects, and the multi-photon nature of experiments. IR computations assume single-photon processes (Buelna-Garcia et al., 2021). The IR spectra of isolated metal clusters in the gas phase were measured for vanadium cluster cations and neutral and cationic niobium clusters (Fielicke et al., 2005). Even though Cu clusters are important in catalysis and were the first clusters produced experimentally (Powers et al., 1982), the available structural information is limited to photoelectron spectroscopy studies of anions, mass spectrometry, and visible-range photodissociation spectra (Lushchikova et al., 2019). Previous work determined the structures of small cationic copper clusters based on a combination of IR spectroscopy of Cun +-Arm complexes and DFT calculations (Lushchikova et al., 2019). In this work, the IR spectra of the isomers were computed using the Gaussian package under harmonic approximation at the PBPW91-D3 [106]/def2TZVP level and a full width at half maximum of 8 cm−1. The Grimme D3 dispersion was considered as implemented in the Gaussian code (Frisch et al., 2009). Imaginary frequencies were checked in all calculations to ensure that the resulting structures were not transition states. The computed frequencies were scaled by a factor of 0.98 to estimate the observed frequencies. Here, the total IR spectrum is computed as a weighted Boltzmann sum of the IR spectrum of each isomer in the distribution at a finite temperature (Buelna-Garcia et al., 2021; Buelna-Garcia et al., 2021; Lecoultre et al., 2011; Sieber et al., 2004). The spectrum is calculated using Eq. 7 and employing the occurrence probabilities displayed in Figure 2. We know of a few theoretical studies on the computation of IR spectra of metal clusters as weighted sums of the IR/UV spectra of the isomers (Sieber et al., 2004). The weighted Boltzmann IR spectra of Cu38 clusters at various temperatures are shown in Figure 4. The transition metal clusters are quite stable and their vibrational frequencies are found to be below 400 cm−1 (Lapoutre et al., 2013). This is in good agreement with our computed spectra displayed in Figure 4. In particular, the IR spectrum at a low temperature is displayed in Figure 4A. There are three dominant peaks at 125, 225, and 250 cm−1. The vibrational mode located at 125 cm−1 is a breathing mode that moves the atoms at the surface, whereas the mode at 250 cm−1 is a breathing mode where the core atoms move. The IR spectra in Figures 4A,B are similar in the 0–600 K temperature range because the relative populations in this temperature range are dominated strongly by TO structures with C1 and D4h structures, as shown in Figure 2B. The IR spectra starts to become small at 700 K. Large changes in the IR spectra happen at a temperature of approximately 800 K, where the solid-solid transition point is located, and where the amorphous structure depicted in Figure 2D and the TO structure depicted in Figure 2A coexist and contribute similarly to the overall IR Boltzmann weighted spectra. Figure 4C shows the IR Boltzmann weighted spectra at 800–1200 K. For temperatures up to 1300 K, the IR Boltzmann weighted spectra are displayed in Figure 4D. At 1200 K, the IR Boltzmann weighted spectra are similar to the spectrum of the amorphous structure depicted in Figure 2D. In contrast, the IR Boltzmann weighted spectra displayed in Figure 4A are similar to the individual IR spectrum of a TO structure. In general, the effect of temperature on the IR spectra is to extenuate the IR spectrum as the temperature increases.

FIGURE 4

(Color online) Temperature-dependent IR Boltzmann -spectra -weighted at room temperature of the neutral Cu38 cluster are shown in panels (A)–(D) for various temperatures. The computed IR spectrum of each isomer is multiplied by its corresponding Boltzmann weight at finite temperature; Then, they are summed together to produce a final Boltzmann-weighted IR spectrum. Each spectrum of each isomer is computed using density functional theory as implemented in Gaussian code at the B3PW91-D3/def2-TZVP level of theory. The large change in the IR spectra occurs at a temperature of 824 K, as we can see in panel (D), and this agrees with the relative occurrence displayed in Figure 2B. The bulk melting temperature of copper is 1358 K (Grigoryan and Springborg, 2019), considering this, our result below of this temperature are well-behaved.

Conclusion

The temperature and entropic effects produce several competing structures because energy separation between isomers on the free energy surface is small and changes the dominant structure. Thus, it is likely that various isomers interconvert at finite temperature. Our findings show that the amorphous structure with C1 symmetry is quite dominant at hot temperatures. These energetically competing structures provide various portions of the overall IR spectrum. In contrast, higher-energy structures with significant energy separation between isomers on the potential/free energy surface do not contribute to the overall IR spectrum. The main contribution to the molecular properties comes from low-energy structures that are close to the global minimum, where the temperature-dependent Boltzmann factor weights are not equal to zero. The Boltzmann-weights depend strongly on the energy separation; the IR spectrum is constant if the energy separation is significant. One motif is dominant in cold conditions and the other in hot conditions. In addition, the level of theory and dispersion influences the location of the TSS point in the temperature scale. Our computations at the six levels of theory clearly show (relative population) that low-symmetry isomers become more stable at high temperatures due to the entropic effect on a Boltzmann ensemble at thermal equilibrium. Our unbiased global search of the free energy surface shows that there is an amorphous structure that dominates at high temperatures. As far as we know, this is a novel high-temperature structure putative global minimum. Computations of the relative populations at a high level of theory is recommended as immediate future work.