Ab initio molecular dynamics studies of Au38(SR)24 isomers under heating

Despite the great success in achieving monodispersity for a great number of monolayer-protected clusters, to date little is known about the dynamics of these ultra-small metal systems, their decomposition mechanisms, and the energy that separates their structural isomers. In this work, we use density functional theory (DFT) to calculate and compare the ground state energy and the Born-Oppenheimer molecular dynamics of two well-known Au38(SCH2CH2Ph)24 nanocluster isomers. The aim is to shed light on the energy difference between the two clusters isomers and analyze their decomposition mechanisms triggered by high temperatures. The results demonstrate that the energy that separates the two isomers is of the same order of magnitude as the energy difference between the fcc and hcp phases of bulk gold reported earlier. Moreover, the MD simulations show disordering and eventual fragmentation of the cluster structures at high temperature which seem to proceed via spontaneous formation of Aux(SR)y polymeric chains. Hence, these results greatly contribute to understanding the possible decomposition mechanism, stability and robustness of existing and new monolayer-protected clusters.


Introduction
Since the first reports on thiolate-passivated gold clusters [1][2][3][4], about hundred atomically precise monolayerprotected nanoclusters (MPC) have been synthesized with highly reproducible physical properties [5][6][7][8]. Despite this great success, however, over the last few years of research in the field the latent question regarding the stability and robustness of this type of ultra-small systems has become critical. Many of the existing MPC have been observed to decompose over a period of several days as consequence of atomic re-organizations or core-size conversion triggered during manipulations, synthesis, or by effect of the environmental conditions [9][10][11][12][13][14][15]. Moreover, highresolution electron microscopy studies have revealed that the exposure to normal electron dose inflicts significant morphological and structural changes in ultra-small metal nanoparticles and nanoclusters [16][17][18][19][20][21][22][23]. All these structural changes occurring in the structure of MPC, especially in solution phase, often lead to coexistence of different chemical species that might hinder the crystallization Contribution to the Topical Issue "Dynamics of Systems on the Nanoscale (2018)", edited by Ilko Bald, Ilia A. Solov'yov, Nigel J. Mason and Andrey V. Solov'yov.
Supplementary material in the form of one zip file available from the Journal web page at https://doi.org/10.1140/epjd/e2019-90441-5 a e-mail: rjuarezmos@gmail.com and the characterization of the main product, limiting thus their usage in practical applications [24,25].
In the particular case of Au 38 (SR) 24 nanocluster, it took over a decade to experimentally disentangle the atomic structure [26][27][28][29][30][31][32][33]. In 2010, Qian et al. found through X-ray crystallography experiments that the metal core of the Au 38 (SCH 2 CH 2 Ph) 24 nanocluster is a prolate biicosahedral with quasi-D 3h symmetry [33]. However, more recently Tian and co-workers have found, also by using single-crystal X ray crystallography, that the metal core of Au 38 (SCH 2 CH 2 Ph) 24 nanoclusters can also adopt an oblate structure with relatively high stability at low temperatures [34]. They observed that the prolate structure reported by Qian [33] and Lopez-Acevedo [31] can irreversibly be obtained from this less symmetric oblate Au 38 (SCH 2 CH 2 Ph) 24 cluster by heating the samples to 50 • C in toluene [34]. To date, this discovery is one in the large list of experimental data that underline the dynamical behavior of MPCs [9][10][11][12][13]. In spite of this, there is only a limited number of theoretical studies that focus on the dynamic behavior of MPC and that compare their energy with respect to structural isomers [35].
In this work, we used the Au 38 (SCH 2 CH 2 Ph) 24 structures of Qian and Tian (hereafter denoted as Au 38Q and Au 38T , respectively) as reference systems to compare the ground state energies, electronic structures, and dynamics of ligand-protected metal nanoclusters with same chemical composition. The aim is to shed light on the energy that separates these molecule-like isomers in their ground state and to study the dynamics of both systems with increasing the temperature. For this purpose, we calculated the energy difference between the two isomers by using different approximations to the exchange-correlation energy in DFT calculations and we performed DFT molecular dynamics simulations from 0-1100 K to analyze the structural changes induced by high temperatures. We found that, in accordance with the experimental observations, the Au 38T nanocluster is less stable than the Au 38Q , and that the lower stability of the structure of Tian can be largely attributed to a combined effect of high fluxionality of the Au atoms and the weakness of specific Au-S bonds when exposed to heat. Furthermore, the MD simulations provide a possible mechanism by which the clusters start to disorder and finally fragment.

Computational methods
Structure optimizations, energy calculations, and electronic structure analysis were carried out using the density functional theory (DFT) as implemented in GPAW [36,37] code package. The electron density was modeled using a real-space grid with spacing of 0.2Å. The structure optimizations were performed using a convergence criterion of 0.05 eV/Å for the residual forces on the atoms. The local-density approximation (LDA) [38], Perdew-Burke-Ernzerhof (PBE) [39], and van der Waals (vdW-DF-cx) [40,41] functionals were used to compare the total energy and the HOMO-LUMO energy gaps. In all-electron calculations, the Au(5d 10 6s 1 ), and S(3s 2 3p 4 ), C(2s 2 2p 2 ), and H(1s 1 ) electrons were treated as valence and the innermost ones were included as a frozen core. The structure optimizations were performed using the Au 38 (SCH 2 CH 2 Ph) 24 crystal structures of Tian [34] and Qian [31,33], as well as the corresponding simplified Au 38 (SCH 3 ) 24 models (hereafter denoted as Au 38t and Au 38q for Tian and Qian structures, respectively).
Born-Oppenheimer NVT molecular dynamics (MD) simulations were performed on the optimized Au 38t and Au 38q structures by using the PBE exchange-correlation functional. The atomic mass of the hydrogen atoms was replaced by the mass of deuterium in order to use simulation time step of 2.0 fs. For comparison, a short MD simulation on the Au 38t structure was also performed using a time step of 1.0 fs to validate the dynamics observed with larger time step (see Figs. S3-S4 in Supplementary material). To equilibrate the systems the Berendsen thermostat was used with a coupling constant to the bath of τ = 0.5 ps [42]. Initially, the MD simulations were performed at 500 K for 12 ps and 15 ps for Au 38q and Au 38t , respectively. However, to see more severe structural changes, the temperature was increased in steps of 100 K until it reached ∼1100 K. The total simulation time was 24.768 ps for both clusters.

Results and discussion
In Figure 1, we show the Au 38T and Au 38Q structures optimized with PBE functional [39]. Both clusters consist of a Au 23 core protected by a surface layer formed of several Au n (SCH 2 CH 2 Ph) m units (with n = 0, 1, 2, or 3, and m = 1, 2, 3, or 4). In the Au 38Q cluster, the Au 23 core is a biicosahedral structure with D 3 symmetry protected by six Au 2 (SCH 2 CH 2 Ph) 3 and three Au(SCH 2 CH 2 Ph) 2 units. In the Au 38T isomer, the Au 23 core is formed by one Au 13 icosahedron and one Au 12 cap (fused together via sharing two gold atoms), and the protecting layer consists of three Au 2 (SCH 2 CH 2 Ph) 3 , three Au(SCH 2 CH 2 Ph) 2 , two Au 3 (SCH 2 CH 2 Ph) 4 and one SCH 2 CH 2 Ph bridge. According to the superatom electron count formula, n * = N v -M -z [43], from which the number of delocalized electrons in a core of the Au 38 (SR) 24 clusters (n*) is calculated as the product of the number of core metal atoms (N = 38) and their atomic valence (v = 1), minus the number of electron-withdrawing ligands (M = 24) and the overall charge (z = 0), these clusters have 14 delocalized electrons [31].
To reduce the computational burden in simulations we replace the SCH 2 CH 2 Ph ligands in the Au 38Q and Au 38T structures with SCH 3 . To validate the choice of the Au 38 (SCH 3 ) 24 (Au 38q and Au 38t ) models, we calculate the density of states (DOS) of the Au 38Q and Au 38T clusters and compare them with the DOS of the simplified Au 38q and Au 38t structures. Figure 2 shows that the DOS of the Au 38 (SCH 2 CH 2 Ph) 24 clusters match the DOS of the Au 38 (SCH 3 ) 24 models within the energy range of −1.5 to 2.0 eV. This indicates that the states near the Fermi level are predominantly of the metal core and that the discrepancies outside the −1.5 to 2.0 eV energy window are due to the lack of the π-π * states of the phenyl rings in the ligands of the Au 38q and Au 38t structures.
Similarly, we investigate the effect of the ligands on the length of the chemical bonds by comparing the optimized Au-Au and Au-S bond distances between the Au 38Q and Au 38q structures. The results from the PBE calculations (shown in Tab. S1 in SI) demonstrate that the Au-Au bonds within the Au 23 core (Au core -Au core ), the Au-S bonds between Au atoms of the Au 23 core and S (Au core -S), and the Au-S bonds within the units Table 1. Relative energies, E r , HOMO-LUMO gaps (in parenthesis), and moments of inertia I x , I y and I z (in square brackets) of the optimized Au 38 (SR) 24 nanoclusters with R = CH 3 (Au 38t and Au 38q ) and R = CH 2 CH 2 Ph (Au 38T and Au 38Q ). Energies, given in eV per cluster, are calculated using the LDA, PBE and vdW-DF-cx functionals. The moments of inertia are in units of 10 3 amuÅ 2 .
In Table 1, we show the relative energies of the Au 38Q (Au 38q ) and Au 38T (Au 38t ) clusters calculated using different exchange-correlation functionals. The results confirm the highest stability of the Au 38Q structure over the Au 38T isomer by 1.8-3.1 eV (i.e., 3.7-6.7 meV per atom). The largest energy difference between these two crystal structures is obtained when using the vdW-DF-cx functional which accounts for long-range correlations that give rise to the vdW forces. These results importantly indicate that the effect of the van der Waals interactions is in the range of 1.0 to 1.4 eV. Moreover, when we calculate the ground-state energy of the simplified Au 38 (SCH 3 ) 24 structures the energy difference between the Qian and Tian clusters drops to 0.3-0.4 eV per cluster (i.e., 2.1-2.8 meV per atom), indicating also that ∼80% of the stabilizing energy of the Au 38Q nanocluster comes from the disposition of the SCH 2 CH 2 Ph ligands around the metal core.
Earlier, Wang et al. [44] reported that the energy difference between fcc and hcp phases of gold calculated using DFT/PW91 is 1.9 meV per atom. Comparing the results of Wang et al. with ours we find that the energy difference between the two forms of Au 38 (SR) 24 are of the same order of magnitude as the energy difference between the fcc phase and the metastable hcp phase of bulk gold. Moreover, the phase transition of square sheets of hcp gold, of length of 200-500 nm and thickness of ∼2.4 nm, to a fcc structure has been observed on exposure to an electron beam [45].
Also, in Table 1 we show the values of the HOMO-LUMO gap and the Cartesian moments of inertia, I x , I y , and I z , of the optimized Au 38T and Au 38Q , and Au 38t and Au 38q structures. The corresponding oblate and prolate shape of the Au 38T (and Au 38t ) and Au 38Q (and Au 38q ) is confirmed by the I z > I x ≈ I y and I z ≈ I y > I x values, respectively. In general, the calculated HOMO-LUMO energy gap of the Au 38Q (and Au 38q ) nanocluster is slightly larger than the one of the Au 38T (Au 38t ) isomer. A HOMO-LUMO gap of 0.9 eV has been reported for the Au 38 (SR) 24 species [31,[46][47][48].
To describe the structural changes induced by heat on the two clusters isomers, we performed MD simulations using the Au 38t and Au 38q models. As expressed above, the resulting energy values from the local optimizations demonstrate that the disposition of the aromatic ligands around the metal core and the vdW interactions are important in stabilizing the two isomers (particularly in the case of Au 38Q ). However, the use of the simplified models in the MD simulations also provides important insights regarding the stability and robustness of the thiolate-protected metal structures in these two geometric arrangements, while it significantly reduces the computational burden.
In Figure 3, we show the optimized Au 38t and Au 38q structures (Figs. 3a and 3e) and different snapshots taken during the MD simulations which were performed using a time step of 2 fs (Figs. 3b, 3c, 3d, 3f, 3g and  3h). The analysis of the MD trajectories shows that already at 500 K the oblate Au 38t structure completely loses its original shape; four Au-S bonds break (two between the Au 23 core and S of the Au 2 (SCH 3 ) 3 units, one between the Au 23 core and S of the SCH 3 bridge, and one within the Au 3 (SCH 3 ) 4 units (see Fig. S2 in Supplementary material), and a polymeric ring composed of Au 7−8 (SCH 3 ) 7 starts to form as an extension of the remaining Au 31−30 (SCH 3 ) 17 fragment (Fig. 3c). Consequently, the Au 7 (SCH 3 ) 7 ring detaches from the Au 31 (SCH 3 ) 17 fragment when the temperature reaches 1000 K (Fig. 3d). After detachment, only a weak bonding energy of 0.31 eV is seen between the Au 31 (SCH 3 ) 17 and Au 7 (SCH 3 ) 7 fragments, which was obtained from PBE calculations as E([ . On the other hand, the Au 38q cluster remains robust within a temperature range of 0-800 K (see Figs. 3f-3g) and Au-S bonds break only when the temperature exceeds 800 K (see Figs. 3g-3h and Fig. S2 in Supplementary material). Here it is important to point out that the temperature of fragmentation that could be measured experimentally might significantly differ from the values reported here since all the stabilizations arising from the ligands (such as ligand-ligand interactions and steric effects) have been omitted for practical purposes. Moreover, and despite the large number of observations of core-size conversion of MPC upon ligand exchange under heating [14,[49][50][51], to the best of our knowledge there are not reports on the fragmentation of this type of systems caused solely by the effect of temperature. In addition, modeling the time scales involve in the fragmentation of MPC at room temperature is unfeasible by using ab initio calculations. Figure 4 shows the evolution of the geometric structure of Au 38t and Au 38q clusters throughout the molecular dynamics simulation; the evolution of the three Cartesian moments of inertia (Figs. 4b and 4f), the evolution of the averaged Au-S bond distances (Figs. 4c and 4g), and the root mean square displacement (RMSD) of the Au atoms of the Au 23 core and units (Figs. 4d and 4h) as a function of time. Figures 4b and 4f show that while the deformation of the oblate Au 38t cluster into an almost-prolate one occurs already when the system reaches 500 K (∼3.70 ps), the Au 38q cluster remains close to the original prolate throughout the simulation. The averaged Au-S distances plotted along the simulations in Figures 4c and 4g are measured using as reference the original Au-S bonding as assigned in the optimized Au 38t and Au 38q structures displayed in Figures 3a and 3e. The evolution of the averaged Au-S distances shows that the bonds between Au atoms of the Au 23 core and S atoms (Au core -S) are more susceptible to break than the Au-S bonds within the protecting units (Au units -S). Particularly, the first bonds to break are the ones between the Au 23 core and S atoms of the Au 2 (SCH 3 ) 3 units. This was also observed in the simulations of Au 38t performed using a time step of 1 fs (see Fig. S4 in Supplementary material). Moreover, the RMSD plots of Au 38t show that the 15 Au atoms of the original the protecting Au 15 (SCH 3 ) 24 layer have larger RMSD values throughout the MD run than the 23 Au atoms in the original core. The same is true for the Au 38q structure but, in general, the RMSD values of all Au atoms in the Au 38q structure are significantly smaller than the ones calculated for Au 38t . The foregoing demonstrates the weakness of the Au core -S bonds and the higher fluxionality of the Au atoms in the Au 38t structure compared to Au 38q . Moreover, in both structures the almost-completed closure of the HOMO-LUMO energy gap was observed once the temperature reached 1000 K (see Fig. S5 in Supplementary material).

Concluding remarks
We present a theoretical analysis of stability and robustness of two well-known ligand protected metal nanoclusters synthesized independently by Tian [34] and Qian [33]. These nanoclusters, which chemical formula is Au 38 (SCH 2 CH 2 Ph) 24 , have been found to have either an oblate (Au 38T ) or a prolate (Au 38Q ) shape depending on the synthesis conditions. Using different approximations to the exchange-correlation functional in DFT calculations we found that these two isomers are separated by an energy of 1.8-3.1 eV per cluster, Au 38Q being the more stable. Moreover, the energy difference between the two isomers is the largest when the vdW interactions are considered. Hence, when the SCH 2 CH 2 Ph ligands are replaced by SCH 3 the energy difference between the two isomers significantly drops (by up to 80-85%), indicating that large part of the stability of the Au 38Q structure arises from the disposition of the ligands around the prolate core. Furthermore, we found that the energy difference between the two Au 38 (SR) 24 isomers (in eV per atom) is of the same order of magnitude as the energy difference between the fcc and hcp phases of bulk gold calculated earlier using the DFT/PW91 level of theory [44].
In addition, the results from the MD simulations performed on the Au 38 (SCH 3 ) 24 models demonstrate that while under the exposure of high temperatures (<800 K) the structure of Qian remains robust, the oblate structure of Tian undergoes drastic structural changes driven by a combined effect of high fluxionality of Au atoms and the weakness of the Au-S bonds between the Au 23 core and the S atoms of the units. Hence, these results offer a qualitative comparison of the robustness of the thiolate-protected Au 38 structure in the two geometric arrangements, at the same time that they propose a possible decomposition mechanism of passivated-gold systems under heating via formation and detachment of Au x (SR) y polymeric chains.
Open access funding provided by University of Jyväskylä (JYU). This work is supported by the Academy of Finland (Project 294217 and H.H. Academy Professorship). The computational resources were provided by CSC-the Finnish IT Center for Science in Espoo, Finland. We thank Satu Mustalahti for useful discussions.