Thermal stability of L-cysteine-protected Au₂₅ clusters: interplay between melting and ligand desorption

Abstract

We investigate the thermal stability and melting of a monolayer-protected Au₂₅Cys₁₈ cluster using classical reactive molecular dynamics simulations. While the enhanced thermal stability of thiol ligand-protected gold clusters compared to corresponding unprotected gold clusters is well known, the mechanism of melting of the protected clusters has not yet been studied in detail. Our results demonstrate that the covalent bonding of the thiol ligands in a Au₂₅Cys₁₈ cluster stabilises the gold core against thermally induced isomerisation and melting. The Au₂₅Cys₁₈ cluster undergoes a melting phase transition at temperatures of ∼580–760 K, which exceeds by approximately 400 K the melting temperature of a bare Au₂₅ cluster. The loss of thermal stability of the ligand-protected cluster occurs through an interplay of the metal core melting and the cascade evaporation of cysteine (Cys) ligands on the nanosecond timescale. The simulation results are validated by comparison with the results of ab initio calculations and relevant experimental data.

---

1. Introduction

Atomically precise ligand-protected gold clusters have attracted strong interest for their distinctive physicochemical properties^1–4 and their potential technological applications in various fields, including catalysis^5,6 and biological applications such as bioimaging, sensing, and targeted drug delivery.^7–9 The choice of specific ligands is known to determine the functional properties, stability and structure of gold clusters.^10,11 The gold–sulphur covalent bonding offers significant advantages in the stabilisation of gold surfaces, clusters and nanoparticles.¹² Some applications of these materials require their exposure to elevated temperatures, so the thermal stability is of both applied and fundamental interest.

Thiol-protected gold clusters have been experimentally synthesised with well-defined numbers of atoms,¹³ some of the most stable and widely studied being Au₂₅(SR)₁₈, Au₃₈(SR)₂₄, Au₁₀₂(SR)₄₄, and Au₁₄₄(SR)₆₀ (where SR stands for a thiolate ligand).^14–18 These systems consist of a metal core protected by a heavily reconstructed surface layer based on Au–S structural bonding. The most common motif is typically referred to as a “Au–S staple”.¹² Au₂₅(SR)₁₈ clusters have been studied experimentally with a variety of thiolate ligands^9,10,19 and have shown high resistance to thermal stress, core etching and spontaneous decomposition.^20–22 As a result, Au₂₅(SR)₁₈ clusters are versatile for biomedical applications, including biolabeling, photothermal cancer therapy and high-resolution in vivo imaging.^9,19 Successful biomedical investigations have focused mainly on glutathione-protected Au₂₅ clusters [Au₂₅(SG)₁₈].

Another attractive ligand with similar bio-functional capabilities is the amino acid L-cysteine (Cys; chemical formula C₃H₇NO₂S). Cysteine is the only amino acid with a terminal thiol group, which allows formation of covalent bonds with the gold clusters, while its amine and carboxyl functional groups enable conjugation to bioactive elements, such as tumour-targeting drugs, DNA, and proteins including antibodies. Au₂₅Cys₁₈ nanoclusters have been synthesized within 10 minutes via a simple protection–deprotection wet chemical method.²⁰ More recently, Au₂₅Cys₁₈ nanoclusters were synthesised even faster in a tube-in-tube membrane reactor.²³ In the later study, these clusters exhibited antibacterial activity and were shown to enhance the photobactericidal effects of crystal violet treated silicone.

Computational modelling greatly complements experimental analysis of the structure, dynamics, and properties of thiol-protected gold clusters. Most computational studies of such clusters have been performed using density functional theory (DFT), which has provided insight into the geometries and electronic structure of specific clusters of different sizes.^24–29 Much less attention has been paid to the study of the thermal stability and thermally-induced structural and phase transitions in such clusters, which is obviously relevant to some applications. While the thermodynamic properties and melting of bare metal clusters have been extensively studied both experimentally^30–34 and theoretically,^35–39 a much smaller number of studies have investigated the corresponding behaviour of functionalised metal clusters.

From an experimental point of view, thermal stability of thiolated gold clusters has been studied using several methods, such as differential scanning calorimetry, thermogravimetric analysis, EXAFS spectroscopy and transmission electron microscopy, which enabled to determine and characterise the cluster decomposition temperatures, the loss of ligands and the changes in cluster size and structure.^22,40,41 In the rare computational studies,^42,43 the thermodynamic stability and melting of Au₂₅(SH)₁₈ and Au₃₈(SR)₂₄ clusters have been investigated by ab initio molecular dynamics (MD). The inherent high computational cost of the ab initio MD method has limited the achievable temporal scales to several tens of picoseconds, resulting in the use of very high heating rates (∼40–50 K ps⁻¹) in the simulations.⁴³ However, the study of the thermal stability and thermally-induced phase and structural transitions in nanosystems requires the consideration of much longer time scales, which can be achieved using the classical MD method.³⁵

In this paper, we investigate the thermal stability and melting of a monolayer-protected Au₂₅Cys₁₈ cluster using classical reactive MD simulations. The thermal stability of the cluster is studied over a wide temperature range from 0 to 1000 K. It is found that the binding of the cysteine ligands to the gold cluster stabilizes the latter against isomer transitions and melting. As a result, the ligand-protected Au₂₅Cys₁₈ cluster undergoes the melting phase transition in the temperature range of ∼580–760 K, which is ∼400 K higher than the melting temperature of a bare/unprotected Au₂₅ cluster. The major thermal response of the Au₂₅Cys₁₈ cluster structure is the sequential desorption of the cysteine ligands, leading to rapid melting of the gold cluster. The desorption of ligands begins stochastically in the similar temperature range of ∼580–760 K and is not site-specific. The simulation results are validated by comparison with the results of DFT calculations and experimental data.

2. Computational methodology

The initial geometry of the Au₂₅Cys₁₈ cluster was determined through the DFT calculations performed using Quantum Espresso.^44,45 The calculations were carried out using the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional⁴⁶ and scalar relativistic projector augmented-wave (PAW) pseudopotentials with a non-linear core correction (NLCC).^47,48 The convergence threshold for self-consistency was set to 10⁻⁵ Ry. The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm was used for structural relaxations. The kinetic energy cutoffs for the wavefunctions and for the charge density and potential were 47 Ry and 364 Ry, respectively. We assumed an isolated system and employed the Martyna–Tuckerman correction to the total energy and self-consistent function potential.⁴⁹ van der Waals interactions were incorporated using Grimme's semi-empirical DFT-D3 correction.⁵⁰ The coordinates for the Au₂₅ cluster were taken from ref. 15. The atomic coordinates for the L-cysteine were taken from the PubChem database (PubChem CID: 5862). After the initial relaxation of the bare Au₂₅ cluster, 18 deprotonated cysteine molecules were attached to the outer-shell Au atoms in a symmetrical pattern, and the system was relaxed again.

Classical MD simulations have been performed using the MBN Explorer software package⁵¹ for multiscale modelling of the structure and dynamics of complex molecular and nanoscale systems.^35,52 MBN Studio,⁵³ a multitasking toolkit and a dedicated graphical user interface for MBN Explorer, has been used to prepare input files and analyse the simulation results. The initial geometry of the cluster was taken from the optimised structure obtained from the DFT calculations.

Interactions between all gold atoms have been described using the many-body Gupta potential⁵⁴ with the parameters taken from ref. 55. To enable the detachment of ligands from the metal core, Au–S bonds have been described using the reactive CHARMM (rCHARMM) force field⁵⁶ implemented in MBN Explorer. All other interatomic interactions within the ligands have been described using the standard (non-reactive) CHARMM27 force field for proteins.⁵⁷ The interaction parameters are listed in the SI, see Tables S1–S3.

To properly account for the Au–S bonding in the form of “staple” structures, two types of gold atoms corresponding to the inner Au₁₃ core and 12 surface atoms of the Au₂₅ cluster have been defined, as shown in Fig. 1. Atoms of the inner icosahedron (shown in orange colour) form a covalent bond with one sulphur atom each. 12 surface gold atoms (shown in green colour) form covalent bonds with two sulphur atoms each.

Fig. 1 Visualisation of different groups of interatomic distances and Au–S covalent bonds used to characterise the structure of the metal core and the Au–S “staples” in the Au₂₅Cys₁₈ cluster. (a) The (Au–Au)₁ group of distances involving the central Au atom and 12 atoms located at the vertices of the inner Au₁₃ icosahedron; (b) the (Au–Au)₂ group of distances between the nearest atoms of the Au₁₃ icosahedron; (c) the (Au–Au)₃ group of distances between two surface gold atoms in each Au–S “staple”; (d) the (Au–S)₁ bonds between the sulphur atoms (shown in yellow) and gold atoms of the Au₁₃ core; (e) the (Au–S)₂ bonds between the sulphur atoms and the surface gold atoms of Au₂₅; (f) optimised geometry of the entire Au₂₅Cys₁₈ cluster considered in this study.

The parameters of the bonded and non-bonded interactions between gold atoms and sulphur atoms were taken from ref. 58. These parameters were successfully used in the previous studies^59,60 for the structural characterisation of gold nanoparticles functionalized with poly(ethylene glycol) ligands. The molecular topology and CHARMM force field parameters for cysteine were obtained using the SwissParam web-service.⁶¹ Partial atomic charges in cysteine molecules were determined by the natural bond orbital (NBO) analysis using the Gaussian 16 software package.⁶² The hydrogen atom was then removed from the thiol group of each cysteine, and the cumulative charge sitting on all the removed hydrogen atoms (+0.179|e| per atom according to the performed DFT calculations) was redistributed equally among all the 25 gold atoms to maintain electrical neutrality of the entire cluster. Since 18 hydrogen atoms were removed in total, each gold atom has acquired a partial charge of (+0.179|e| × 18/25) ≈ +0.129|e|.

The optimised structure of the Au₂₅Cys₁₈ cluster was then used as the starting geometry for MD simulations of the cluster dynamics in a wide range of temperatures from 0 K to 1000 K. A series of successive simulations were performed to heat the cluster in a stepwise manner in steps of 20 K. For each temperature, the cluster dynamics was simulated for 4 ns using a Langevin thermostat, a damping time of 1 ps, and an integration time step of 1 fs. The total simulation time is thus equal to 200 ns, which is two orders of magnitude longer than the simulation times in previous ab initio MD simulations of the heating and melting of thiolated gold clusters.^42,43

Similar complementary simulations were performed to investigate the melting of a bare Au₂₅ cluster. The initial cluster structure was taken from the Cambridge Cluster Database⁶³ and corresponds to the global energy minimum for Au₂₅, as calculated using the Sutton–Chen interatomic potential.⁶⁴ This cluster structure was annealed to 300 K and optimised using the Gupta potential to generate the initial geometry for the present simulations. The bare Au₂₅ cluster was heated from 0 to 1000 K in 20 K steps, with the same simulation parameters as for the Au₂₅Cys₁₈ cluster. Additional benchmark simulations were also performed (see SI, Fig. S2), where the heating of the Au₂₅ cluster in 20 K steps was simulated via a series of consecutive 10 ns and 20 ns-long runs. Thus, the total simulation time for heating from 0 to 1000 K was equal to 500 ns and 1000 ns, corresponding to heat rates of 2 K ns⁻¹ and 1 K ns⁻¹, respectively.

3. Results and discussion

3.1. Structural characterisation of the Au₂₅Cys₁₈ cluster

Prior to analysing the thermal stability of Au₂₅Cys₁₈ cluster, we have validated the utilised parameters of the force fields and interatomic potentials by characterising the optimised structure of the cluster.

The structure of the ligand-protected Au₂₅ cluster is characterised by six [–S–Au–S–Au–S–] ‘staples’ that surround an icosahedral Au₁₃ core.^14,19 To characterise the structure of Au₂₅Cys₁₈, several distinct groups of distances between gold atoms and Au–S bond lengths were identified and quantitatively analysed, as shown in Fig. 1. The (Au–Au)₁ group refers to the distance between the central gold atom and 12 atoms of the icosahedral Au₁₃ core, see Fig. 1(a). The (Au–Au)₂ group defines the distances between two neighbouring atoms on the shell of the icosahedral Au₁₃, see Fig. 1(b). The (Au–Au)₃ group defines the distances between the gold atoms in each Au–S ‘staple’, see Fig. 1(c). The Au–S ‘staple’ structure is characterised by two groups of Au–S bond lengths: the (Au–S)₁ group includes the bonds between sulphur atoms of the ligands and the gold atoms of the Au₁₃ icosahedron (Fig. 1(d)), while the (Au–S)₂ group includes the sulphur atoms bonded to the 12 outer gold atoms (Fig. 1(e)).

The Au–Au interatomic distances and Au–S bond lengths for the optimised geometry of Au₂₅Cys₁₈, obtained using the classical interatomic potentials, are listed in Table 1 (see the column labelled “classical”). The calculated values are compared with the corresponding values obtained from the DFT-based geometry optimisations of Au₂₅Cys₁₈ (this work) and of a structurally similar Au₂₅(SH)₁₈⁻ cluster,⁴² as well as with the results of experimental characterisation of another structurally similar Au₂₅(SCH₂CH₂Ph)₁₈ cluster obtained by X-ray diffraction (XRD) methods.¹⁴

Table 1 Three groups of interatomic distances between gold atoms, (Au–Au)_1,2,3, and two groups of Au–S bond lengths, (Au–S)_1,2, in the Au₂₅Cys₁₈ cluster. These groups are visualised in Fig. 1(a)–(e). The values (in Å) are averaged over the entire cluster and are given with the corresponding standard errors

	Classical	DFT		Exp.^14
		This work	Ref. 42
(Au–Au)1	2.81 ± 0.01	2.83 ± 0.01	—	2.79 ± 0.01
(Au–Au)2	2.96 ± 0.02	2.98 ± 0.02	(2.82–3.03)	2.93 ± 0.06
(Au–Au)3	3.63 ± 0.01	3.55 ± 0.01	—	—
(Au–S)1	2.40 ± 0.01	2.39 ± 0.01	2.46	2.38 ± 0.01
(Au–S)2	2.38 ± 0.01	2.33 ± 0.01	2.36	2.32 ± 0.01

---

The interatomic Au–Au distances determined from classical and DFT optimisation calculations are in close agreement (within 0.05 Å) with the experimental values. The slightly larger distances obtained by DFT are to be expected, as the PBE functional is known to overestimate bonding.⁶⁵ The Au–S bond lengths in the structures obtained from DFT and classical optimisation calculations are also in agreement with each other, with the experimental data¹⁴ and with the results of previous DFT calculations⁴² to within ∼0.08 Å. The slight difference in the (Au–Au)₃ distances and the Au–S bond lengths may be due to the different ligands considered in this and the cited studies and to the different charge states of the cluster (neutral cluster in this study vs. singly charged anion in ref. 14 and 42). To complete the characterisation of the cluster structure, we also evaluated the angle between three sulphur atoms in each [–S–Au–S–Au–S–] ‘staple’. The resulting average value of (116.3 ± 1.2)° is close to the value of 109.3° obtained by DFT in ref. 42.

3.2. Dynamical response at elevated temperatures

In this study, the dynamics of the Au₂₅Cys₁₈ cluster has been studied over a wide temperature range from 0 to 1000 K. The thermal stability of the cluster is characterised by calculating the root-mean-square displacement (RMSD) averaged over all gold atoms. The RMSD over a period t is given by:where N is the number of considered atoms, r_i(t) is the position vector of atom i at time t and r^ref_i is the position vector of the same atom in the reference structure.

Fig. 2(a) shows the RMSD averaged over all gold atoms in the Au₂₅Cys₁₈ cluster (solid blue symbols) and in the bare Au₂₅ cluster (open green symbols) as a function of temperature. Each data point represents the mean value at a given temperature, computed from five independent MD simulations, with error bars indicating the standard error of the mean. For the ligand-protected Au₂₅Cys₁₈ cluster, the RMSD slowly increases from ∼0.1 Å to 0.5 Å at temperatures below ∼560 K, indicating the stability of the entire structure. As the temperature increases from 0 K to 560 K, the increased amplitude of thermal vibrations causes the cysteine ligands to reorient to different conformational states but the inner Au₁₃ core and the Au–S ‘staples’ do not undergo any major structural changes.

Fig. 2 The temperature dependence of the root-mean square displacement (a) and the Lindemann ratio (b) for gold atoms in the ligand-protected Au₂₅Cys₁₈ cluster and the bare Au₂₅ cluster. The results shown in panel (a) were averaged over five independent MD simulations. Error bars indicate the standard error of the mean.

At temperatures of ∼580–760 K, the RMSD for Au₂₅Cys₁₈ increases rapidly by an order of magnitude, indicating that the ligand-protected cluster has undergone a phase transition. Above 760 K, the averaged RMSD fluctuates around a value of ∼5 Å, suggesting that the gold cluster has melted and remains in a liquid-like state. In contrast, the RMSD of a bare Au₂₅ cluster (green symbols) increases rapidly at significantly lower temperatures. It begins to rise between ∼100 and 200 K, and then rises more quickly between ∼200 and 300 K. Above 300 K, the RMSD for the bare Au₂₅ cluster stabilises at a value of ∼5 Å, thus suggesting that the cluster has transformed into a liquid-like state. The results of five independent simulations also reveal the temperature range over which melting occurs for each cluster type. The RMSD values for the bare Au₂₅ cluster show minimal deviation from the mean RMSD. In contrast, the cysteine-protected Au₂₅Cys₁₈ cluster exhibits a broader melting temperature range of ∼580–760 K, as indicated by the respective standard error bars.

To verify the melting temperature values obtained from the analysis of RMSD, we also analysed the melting of ligand-protected and bare Au₂₅ clusters according to the Lindemann criterion.^66,67 This criterion states that a system undergoes a melting phase transition when the ratio of the root-mean-square displacement of atoms with respect to the average equilibrium interatomic distance exceeds a certain threshold, indicating a loss of long-range order. This threshold depends on several factors such as crystal structure, bonding, and system size, with typical values being in the range of ∼0.10–0.15.^67–69

Following the approach used in previous MD-based studies of the melting of nanoscale systems,^70–72 the melting of Au₂₅Cys₁₈ and Au₂₅ clusters was quantified via the Lindemann ratio, Δ:

where FWHM is the full width at half maximum of the first peak in the radial distribution function (RDF) for gold atoms and d₀ is the mean interatomic distance corresponding to the maximum of the first RDF peak. Thus, the Lindemann ratio Δ provides a measure of the fluctuation of atomic positions relative to the average interatomic distance for all gold atoms at different temperatures. At each temperature considered, the RDF was calculated by averaging the atomic positions over a 4-ns MD trajectory. Fig. 2(b) shows the temperature dependence of the Lindemann ratio for the Au₂₅Cys₁₈ and Au₂₅ clusters, obtained from one exemplar simulated MD trajectory. The horizontal dashed line indicates the typical threshold value of Δ = 0.1, above which the system melts according to the Lindemann criterion.

It is worth noting that the question of evaluating the Lindemann parameter for different bulk materials and finite-size systems has been widely discussed in the literature for several decades. A broad range of values for the Lindemann parameter, from ∼0.03 to 0.24, has been reported.^73,74 At the same time, a narrower range of values, from ∼0.08 to 0.15, has often been used to characterize the melting of metal clusters and nanoparticles.^69,75–77 Regarding finite-size gold systems, several studies have reported values of the Lindemann index of ∼0.08–0.10 corresponding to the onset of melting of the Au₅₆₁ and Au₅₅ clusters.^78,79 A similar range of values (approximately 0.07–0.10) has also been reported in an MD simulation of the melting of a larger-size gold nanoparticle containing 3538 atoms.⁸⁰ Finally, ref. 81 explored the idea that the Lindemann parameter could be linked to the periodic groups of the periodic table, with an exact value for each element belonging to a given periodic group. The study⁸¹ reported 12 distinctive values of the Lindemann coefficient corresponding to 12 groups of the periodic table containing solid elements with identifiable melting temperature. For group 11 (Cu, Ag, and Au), the reported value of the Lindemann parameter is 0.108, close to the above-mentioned values obtained from MD simulations for finite-size gold nanosystems. Based on the above considerations, we have concluded that a “universal” value of Δ = 0.10 would be suitable for the analysis presented in this study. The aim of this analysis is to complement and verify the evaluation of melting temperatures obtained from the temperature dependence of RMSD (see Fig. 2(a)). A more detailed quantitative analysis of the Lindemann criterion for the studied systems goes beyond the scope of this study.

The temperatures at which the calculated Lindemann ratio for each system exceeds the threshold value differ significantly, ∼240 K for the bare Au₂₅ cluster and ∼660 K for the Au₂₅Cys₁₈ cluster. Comparison of these values with Fig. 2(a) shows good agreement between the temperatures at which a sharp rise in the RMSD of gold atoms occurs and the temperatures at which the Lindemann ratio exceeds the chosen threshold value. This supports the interpretation that the melting of the Au₂₅ and Au₂₅Cys₁₈ clusters occurs at the temperature ranges of ∼200–300 K and ∼580–760 K, respectively.

The predicted melting temperature of Au₂₅ is close to the melting temperatures of slightly larger gold clusters (Au₃₀–Au₅₀), as determined by MD simulations using the Gupta potential.⁷⁹ Reported values of melting temperature vary from (260 ± 50) K for Au₃₅ to (330 ± 20) K for Au₅₀. Further discussion of the melting point of the Au₂₅ cluster is provided in the SI.

As shown in Fig. 2(b), the Lindemann ratio for the ligand-protected Au₂₅Cys₁₈ cluster exhibits a more irregular dependence on temperature as compared to the smooth, monotonic dependence observed for the bare Au₂₅ cluster. This irregularity reflects the different dynamics of the ligand-protected cluster at different temperatures. Between 0 K and ∼220 K, the Lindemann ratio for Au₂₅Cys₁₈ increases monotonically, suggesting that the gold atoms in the cluster vibrate thermally around their equilibrium positions. Above 220 K, the Lindemann ratio decreases, indicating a structural rearrangement of the gold cluster. This is further evidenced by changes in the RDFs, as shown in Fig. 3. Two narrow peaks at interatomic distances of ∼2.83 Å and 2.99 Å, which are present in the RDF at T = 200 K, merge into one broader peak centred at 2.87 Å at T = 300 K, indicating a structural rearrangement from the initial equilibrium geometry. The Lindemann ratio decreases until the temperature of ∼340 K, above which it begins to increase again. Between ∼400 K and 660 K, the Lindemann ratio for Au₂₅Cys₁₈ exhibits a similar trend to that observed in the bare Au₂₅ cluster. The RDFs at temperatures between 300 K and 600 K do not vary significantly, and the cluster structure is preserved.

Fig. 3 Radial distribution function (RDF) for gold atoms of the Au₂₅Cys₁₈ cluster at different temperatures.

In the selected MD simulation analysed in Fig. 2(b) and 3, the thermal vibrations of the gold atoms in Au₂₅Cys₁₈ become more pronounced above 600 K. According to the simulations performed, at 640 K, a sulphur atom in one of the Au–S ‘staples’ gains enough thermal energy to move to a different binding site on the Au₂₅ surface. This causes reorganisation of the entire Au–S “staple”. As the temperature increases to ∼680 K, the displaced ligand detaches from the gold cluster, accompanied by a significant reorganisation of the Au atoms within the cluster. The presence of an unprotected region on the cluster surface, created when one cysteine ligand detaches, causes the entire gold cluster to melt rapidly. A similar physical picture was observed in four other simulations of Au₂₅Cys₁₈ heating. As shown in Fig. 2(a) and discussed above, the temperature at which the gold cluster in Au₂₅Cys₁₈ starts to melt varies across the simulations performed, ranging from 580 to 760 K. The results in Fig. 2(b) and 3 therefore provide a characteristic quantitative analysis of the melting process.

3.3. Dynamics of thermally induced desorption of ligands

To gain a better understanding of the interplay between gold cluster melting and the detachment of cysteine ligands, we have analysed more closely the dynamics of the cluster in the temperature range from 640 to 680 K. The chosen temperature range corresponds to three subsequent 4 ns-long MD trajectories, which were combined into one 12 ns trajectory. Fig. 4(a) shows the evolution of the interatomic distances between the sulphur atom of the first detached ligand and the two gold atoms to which it is initially bonded: one atom in the Au₁₃ icosahedral core, d[(Au–S)₁], and one outer atom, d[(Au–S)₂]. The instant t = 0 ns corresponds to the beginning of the simulation at 640 K. Two vertical dashed lines indicate instants when the system temperature is equal to 660 and 680 K. Fig. 4(b) shows the RDFs for all gold atoms, averaged over 2 ns segments of the considered 12 ns combined trajectory. Four labelled vertical lines highlight specific peaks and features in the RDFs that evolve due to the desorption of cysteine ligands from the gold cluster.

Fig. 4 (a) Evolution of interatomic distances between the sulphur atom of the first detached cysteine ligand and two gold atoms to which this atom is initially bonded: one in the Au₁₃ icosahedral core, (Au–S)₁, and one outer atom, (Au–S)₂. The considered 12 ns segment of the trajectory corresponds to the temperature increase from 640 to 680 K. (b) RDFs for all gold atoms in the Au₂₅Cys₁₈ cluster averaged over 2 ns segments.

As can be seen in Fig. 4(a), the selected (Au–S)₁ and (Au–S)₂ bond distances practically do not change during the first ∼2 ns of the considered trajectory segment, with values of (2.40 ± 0.06) Å. This is followed by a sharp increase of the d[(Au–S)₁] distance to ∼5.1 Å, indicating the displacement of the cysteine ligand from its original binding site. The d[(Au–S)₂] distance remains unchanged during this displacement until the time instant of ∼10.7 ns, when the ligand gains sufficient energy to desorb from the cluster surface. This is indicated by a sharp increase in both d[(Au–S)₁] and d[(Au–S)₂] distances.

The ligand displacement and subsequent detachment events can be correlated with changes of the Lindemann ratio and the RDFs, shown in Fig. 2(b) and 4(b). The calculated Lindemann ratio for the gold atoms of the Au₂₅Cys₁₈ cluster exceeds the threshold value of Δ = 0.1 at ∼640 K, after which one cysteine ligand is displaced to a different binding site, causing a restructuring of the gold atoms. This restructuring is indicated by the broadening and merging of the RDF peaks labelled by the dashed lines 1 and 2 (see the RDFs for the segments 0–2 ns and 2–4 ns in Fig. 4(b)). Until the ligand detachment (∼10.7 ns), the RDF peaks labelled 2, 3 and 4 remain unchanged. Once the ligand has detached (see the RDF for the t = 10–12 ns segment), peak 3 broadens into peak 2, while peak 4 also broadens and becomes less distinct compared with the RDFs for other segments.

At temperatures above ∼700 K, the remaining ligands begin to sequentially detach from the Au₂₅ cluster surface. This process is not site-specific and the ligands desorb from the gold cluster randomly. The desorption of a larger number of ligands causes further structural rearrangement of the gold cluster, as seen by the variation of the Lindemann ratio in Fig. 2(b). Interestingly, in the simulations performed, the Au₂₅ cluster did not lose integrity due to ligand desorption. This observation is in contrast to the results of a previous experimental study of the desorption of dodecathiol ligands from a Au(111) surface.⁸² According to that study, dodecathiol ligands can weaken the bonding between the surface Au atoms and the Au atoms in the bulk region, causing surface Au atoms to detach together with the ligands. In the present study, however, we saw no indication of such behaviour. This can be explained by the twofold difference in the binding energy between an outer gold atom and other gold atoms in the Au₂₅ cluster (∼4.3 eV) and the Au–S bond dissociation energy (∼2.2 eV).

The results of this study can be discussed in connection to the previous experimental studies of the thermal stability of ligand protected gold clusters.^22,40,41 It was found in these studies that thermal decomposition of the ligand-protected clusters, including Au₂₅(SR)₁₈, occurs at temperatures of ∼470–520 K, with some small variations due to cluster size and its charge state. In the present simulations, the lowest temperature at which the first ligand detachment was observed is 580 K, which is reasonably close to the reported experimental values. Another experimental study⁸³ concluded that the threshold temperatures for the thermal stability of ligand-protected clusters may depend on the ligand type and can vary up to 80 K for different ligands. On this basis, one can conclude that the results of the present simulations agree reasonably well with the experimental data. As demonstrated for the bare Au₂₅ cluster (see SI, Fig. S2), increasing the simulation time by a factor of five does not change the melting temperature of the cluster. Similar behaviour can also be expected for the Au₂₅Cys₁₈ cluster.

4. Conclusion

In conclusion, we have studied the thermal stability and the mechanism of melting of ligand-protected gold clusters using classical reactive MD simulations. The Au₂₅Cys₁₈ cluster containing 18 cysteine amino acid ligands was considered as a case study. It was shown that the binding of the cysteine ligands to the gold cluster stabilizes the latter against isomer transitions and melting. The melting temperature of the Au₂₅Cys₁₈ cluster was found to be in the range of ∼580–760 K, which is approximately 400 K higher than that of a bare Au₂₅ cluster. We demonstrated that the loss of thermal stability of the ligand-protected Au₂₅Cys₁₈ cluster occurs through an interplay of the metal core melting and the cascade evaporation of cysteine ligands on the nanosecond timescale. Ligand desorption begins in the above indicated range of temperatures and is not site-specific. Understanding the mechanisms of ligand desorption from size-selected metal clusters is also relevant for the optimisation of these systems as catalysts and for the development of next-generation catalytic materials.^6,84 The parameters of the classical force field for describing the cysteine ligands and their interaction with the gold core have been validated by the comparison with the results of DFT calculations and relevant experimental data.

Author contributions

M. A. S.: investigation, formal analysis, visualization, writing – original draft, writing – review & editing; A. V. V.: investigation (MD simulations), formal analysis, data curation, methodology, visualization, writing – original draft, writing – review & editing; T. P.: investigation (DFT calculations), formal analysis, writing – original draft, writing – review & editing; R. E. P.: conceptualization, funding acquisition, supervision, writing – original draft, writing – review & editing; A. V. S.: conceptualization, funding acquisition, methodology, project administration, resources, supervision, writing – original draft, writing – review & editing.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: force field parameters, additional analysis and relevant benchmarks. See DOI: https://doi.org/10.1039/d5cp02394g.

Acknowledgements

This work has been supported by the RADON project (GA 872494) within the H2020-MSCA-RISE-2019 call and the COST Action CA20129 MultIChem supported by COST (European Cooperation in Science and Technology). The possibility of performing computer simulations at the Goethe-HLR cluster of the Frankfurt Center for Scientific Computing is gratefully acknowledged.

References

1. Frontiers of Nanoscience, Protected Metal Clusters: From Fundamentals to Applications, ed. T. Tsukuda and H. Häkkinen, Elsevier, Amsterdam, 2015, vol. 9.
2. Z. Wu and R. Jin, Atomically Precise Metal Nanoclusters, Springer, Cham, Switzerland, 2021.
3. S. Hossain, D. Hirayama, A. Ikeda, M. Ishimi, S. Funaki, A. Samanta, T. Kawawaki and Y. Negishi, Atomically precise thiolate-protected gold nanoclusters: Current status of designability of the structure and physicochemical properties, Aggregate, 2023, 4, e255.
4. M. F. Matus and H. Häkkinen, Understanding ligand-protected noble metal nanoclusters at work, Nat. Rev. Mater., 2023, 8, 372–389.
5. G. Li and R. Jin, Atomically precise gold nanoclusters as new model catalysts, Acc. Chem. Res., 2013, 46, 1749–1758.
6. Y. Du, H. Sheng, D. Astruc and M. Zhu, Atomically precise noble metal nanoclusters as efficient catalysts: A bridge between structure and properties, Chem. Rev., 2020, 120, 526–622.
7. S. M. van de Looij, E. R. Hebels, M. Viola, M. Hembury, S. Oliveira and T. Vermonden, Gold nanoclusters: Imaging, therapy, and theranostic roles in biomedical applications, Bioconjugate Chem., 2022, 33, 4–23.
8. S. M. Farkhani, P. Dehghankelishadi, A. Refaat, D. V. Gopal, A. Cifuentes-Rius and N. H. Voelcker, Tailoring gold nanocluster properties for biomedical applications: From sensing to bioimaging and theranostics, Prog. Mater. Sci., 2024, 142, 101229.
9. Q. Yan, Z. Yuan, Y. Wu, C. Zhou, Y. Dai, X. Wan, D. Yang, X. Liu, N. Xue, Y. Zhu and Y. Yang, Atomically precise water-soluble gold nanoclusters: Synthesis and biomedical application, Precis. Chem., 2023, 1, 468–479.
10. R. R. Nasaruddin, T. Chen, N. Yan and J. Xie, Roles of thiolate ligands in the synthesis, properties and catalytic application of gold nanoclusters, Coord. Chem. Rev., 2018, 368, 60–79.
11. E. Pensa, L. M. Azofra, R. C. Salvarezza and P. Carro, Effect of ligands on the stability of gold nanoclusters, J. Phys. Chem. Lett., 2022, 13, 6475–6480.
12. H. Häkkinen, The gold-sulfur interface at the nanoscale, Nat. Chem., 2012, 4, 443–455.
13. R. Jin, C. Zeng, M. Zhou and Y. Chen, Atomically precise colloidal metal nanoclusters and nanoparticles: Fundamentals and opportunities, Chem. Rev., 2016, 116, 10346–10413.
14. M. W. Heaven, A. Dass, P. S. White, K. M. Holt and R. W. Murray, Crystal structure of the gold nanoparticle [N(C₈H₁₇)₄][Au₂₅(SCH₂CH₂Ph)₁₈], J. Am. Chem. Soc., 2008, 130, 3754–3755.
15. M. Zhu, C. M. Aikens, F. J. Hollander, G. C. Schatz and R. Jin, Correlating the crystal structure of a thiol-protected Au₂₅ cluster and optical properties, J. Am. Chem. Soc., 2008, 130, 5883–5885.
16. H. Qian, W. T. Eckenhoff, Y. Zhu, T. Pintauer and R. Jin, Total structure determination of thiolate-protected Au₃₈ nanoparticles, J. Am. Chem. Soc., 2010, 132, 8280–8281.
17. P. D. Jadzinsky, G. Calero, C. J. Ackerson, D. A. Bushnell and R. D. Kornberg, Structure of a thiol monolayer-protected gold nanoparticle at 1.1 Å resolution, Science, 2007, 318, 430–433.
18. N. Yan, N. Xia, L. Liao, M. Zhu, F. Jin, R. Jin and Z. Wu, Unraveling the long-pursued Au₁₄₄ structure by x-ray crystallography, Sci. Adv., 2018, 4, eaat7259.
19. X. Kang, H. Chong and M. Zhu, Au₂₅(SR)₁₈: The captain of the great nanocluster ship, Nanoscale, 2018, 10, 10758–10834.
20. X. Yuan, Y. Yu, Q. Yao, Q. Zhang and J. Xie, Fast synthesis of thiolated Au₂₅ nanoclusters via protection–deprotection method, J. Phys. Chem. Lett., 2012, 3, 2310–2314.
21. Y. Shichibu, Y. Negishi, H. Tsunoyama, M. Kanehara, T. Teranishi and T. Tsukuda, Extremely high stability of glutathionate-protected Au₂₅ clusters against core etching, Small, 2007, 3, 835–839.
22. M. A. Tofanelli and C. J. Ackerson, Superatom electron configuration predicts thermal stability of Au₂₅(SR)₁₈ nanoclusters, J. Am. Chem. Soc., 2012, 134, 16937–16940.
23. G. B. Hwang, G. Wu, J. Shin, L. Panariello, V. Sebastian, K. Karu, E. Allan, A. Gavriilidis and I. P. Parkin, Continuous single-phase synthesis of [Au₂₅(Cys)₁₈] nanoclusters and their photobactericidal enhancement, ACS Appl. Mater. Interfaces, 2020, 12, 49021–49029.
24. H. Häkkinen, Ligand-protected gold nanoclusters as superatoms – insights from theory and computations, Metal Nanoparticles and Nanoalloys, Elsevier, Amsterdam, 2012, pp. 129–157.
25. Y. Li, G. Galli and F. Gygi, Electronic structure of thiolate-covered gold nanoparticles: Au₁₀₂(MBA)₄₄, ACS Nano, 2008, 2, 1896–1902.
26. Y. Pei and Z. Cheng, Investigating the structural evolution of thiolate protected gold clusters from first-principles, Nanoscale, 2012, 4, 4054–4072.
27. J. Zhong, X. Tang, J. Tang, J. Su and Y. Pei, Density functional theory studies on structure, ligand exchange, and optical properties of ligand-protected gold nanoclusters: Thiolate versus selenolate, J. Phys. Chem. C, 2015, 119, 9205–9214.
28. C. H. Chan, F. Poignant, M. Beuve, E. Dumont and D. Loffreda, Effect of the ligand binding strength on the morphology of functionalized gold nanoparticles, J. Phys. Chem. Lett., 2020, 11, 2717–2723.
29. W. Han, G. Wang, P. Liu, W. Li and W. W. Xu, Structural predictions of three medium-sized thiolate-protected gold nanoclusters Au₄₄(SR)₃₀, Au₅₆(SR)₃₂, and Au₆₀(SR)₃₄, Nanoscale Adv., 2023, 5, 4464–4469.
30. M. Schmidt, R. Kusche, B. von Issendorff and H. Haberland, Irregular variations in the melting point of size-selected atomic clusters, Nature, 1998, 393, 238–240.
31. H. Haberland, T. Hippler, J. Donges, O. Kostko, M. Schmidt and B. von Issendorff, Melting of sodium clusters: Where do the magic numbers come from?, Phys. Rev. Lett., 2005, 94, 035701.
32. A. Yalamanchali, K. L. Pyfer and M. F. Jarrold, Melting of size-selected aluminum clusters with 150–342 atoms: The transition to thermodynamic scaling, J. Phys. Chem. C, 2017, 121, 10242–10248.
33. C. Hock, C. Bartels, S. Straßburg, M. Schmidt, H. Haberland, B. von Issendorff and A. Aguado, Premelting and postmelting in clusters, Phys. Rev. Lett., 2009, 102, 043401.
34. D. M. Foster, T. Pavloudis, J. Kioseoglou and R. E. Palmer, Atomic-resolution imaging of surface and core melting in individual size-selected Au nanoclusters on carbon, Nat. Commun., 2019, 10, 2583.
35. I. A. Solov’yov, A. V. Verkhovtsev, A. V. Korol and A. V. Solov’yov, Dynamics of Systems on the Nanoscale, Springer International Publishing, Cham, 2022.
36. N. Gaston, Cluster melting: New, limiting, and liminal phenomena, Adv. Phys.:X, 2018, 3, 1401487.
37. Latest Advances in Atomic Clusters Collision: Fission, Fusion, Electron, Ion and Photon Impact, ed. J.-P. Connerade and A. V. Solov’yov, Imperial College Press, London, 2004.
38. J. A. Alonso, Structure and Properties of Atomic Nanoclusters, Imperial College Press, London, 2011.
39. A. Aguado and M. F. Jarrold, Melting and freezing of metal clusters, Annu. Rev. Phys. Chem., 2011, 62, 151–172.
40. Y. Zhu, H. Qian and R. Jin, An atomic-level strategy for unraveling gold nanocatalysis from the perspective of Au_n(SR)_m nanoclusters, Chem. – Eur. J., 2010, 16, 11455–11462.
41. S. Chen, M. Li, S. Yu, S. Louisia, W. Chuang, M. Gao, C. Chen, J. Jin, M. B. Salmeron and P. Yang, Ligand removal of Au₂₅ nanoclusters by thermal and electrochemical treatments for selective CO₂ electroreduction to CO, J. Chem. Phys., 2021, 155, 051101.
42. V. Mäkinen and H. Häkkinen, Density functional theory molecular dynamics study of the Au₂₅(SR)₁₈⁻ cluster, Eur. Phys. J. D, 2012, 66, 310.
43. R. Juarez-Mosqueda, S. Malola and H. Häkkinen, Ab initio molecular dynamics studies of Au₃₈(SR)₂₄ isomers under heating, Eur. Phys. J. D, 2019, 73, 62.
44. P. Giannozzi, et al., QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter, 2009, 21, 395502.
45. P. Giannozzi, et al., Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter, 2017, 29, 465901.
46. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
47. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
48. A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci., 2014, 95, 337–350.
49. G. J. Martyna and M. E. Tuckerman, A reciprocal space based method for treating long range interactions in ab initio and force-field-based calculations in clusters, J. Chem. Phys., 1999, 110, 2810–2821.
50. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132, 154104.
51. I. A. Solov’yov, A. V. Yakubovich, P. V. Nikolaev, I. Volkovets and A. V. Solov’yov, MesoBioNano Explorer – A universal program for multiscale computer simulations of complex molecular structure and dynamics, J. Comput. Chem., 2012, 33, 2412–2439.
52. I. A. Solov’yov, A. V. Korol and A. V. Solov’yov, Multiscale Modeling of Complex Molecular Structure and Dynamics with MBN Explorer, Springer International Publishing, Cham, 2017.
53. G. B. Sushko, I. A. Solov’yov and A. V. Solov’yov, Modeling MesoBioNano systems with MBN Studio made easy, J. Mol. Graphics Modell., 2019, 88, 247–260.
54. R. P. Gupta, Lattice relaxation at a metal surface, Phys. Rev. B:Condens. Matter Mater. Phys., 1981, 23, 6265–6270.
55. F. Cleri and V. Rosato, Tight-binding potentials for transition metals and alloys, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 22–33.
56. G. B. Sushko, I. A. Solov’yov, A. V. Verkhovtsev, S. N. Volkov and A. V. Solov’yov, Studying chemical reactions in biological systems with MBN Explorer: Implementation of molecular mechanics with dynamical topology, Eur. Phys. J. D, 2016, 70, 12.
57. A. D. Mackerell, et al., All-atom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B, 1998, 102, 3586–3616.
58. E. Pohjolainen, X. Chen, S. Malola, G. Groenhof and H. Häkkinen, A unified AMBER-compatible molecular mechanics force field for thiolate-protected gold nanoclusters, J. Chem. Theory Comput., 2016, 12, 1342–1350.
59. A. V. Verkhovtsev, A. Nichols, N. J. Mason and A. V. Solov’yov, Molecular dynamics characterization of radiosensitizing coated gold nanoparticles in aqueous environment, J. Phys. Chem. A, 2022, 126, 2170–2184.
60. M. D. Dickers, A. V. Verkhovtsev, N. J. Mason and A. V. Solov’yov, Atomistic modelling and structural characterisation of coated gold nanoparticles for biomedical applications, Eur. Phys. J. D, 2023, 77, 141.
61. V. Zoete, M. A. Cuendet, A. Grosdidier and O. Michielin, SwissParam: A fast force field generation tool for small organic molecules, J. Comput. Chem., 2011, 32, 2359–2368.
62. M. J. Frisch, et al., Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
63. D. J. Wales, J. P. K. Doye, A. Dullweber, M. P. Hodges, F. Y. Naumkin, F. Calvo, J. Hernández-Rojas and T. F. Middleton, The Cambridge Cluster Database, https://www-wales.ch.cam.ac.uk/CCD.html.
64. J. P. K. Doye and D. J. Wales, Global minima for transition metal clusters described by Sutton–Chen potentials, New J. Chem., 1998, 22, 733–744.
65. G.-X. Zhang, A. M. Reilly, A. Tkatchenko and M. Scheffler, Performance of various density-functional approximations for cohesive properties of 64 bulk solids, New J. Phys., 2018, 20, 063020.
66. F. A. Lindemann, Über die Berechnung molekularer Eigenfrequenzen, Phys. Z., 1910, 11, 609–612.
67. C. Chakravarty, P. Debenedetti and F. Stillinger, Lindemann measures for the solid-liquid phase transition, J. Chem. Phys., 2007, 126, 204508.
68. Y. Zhou, M. Karplus, K. D. Ball and R. S. Berry, The distance fluctuation criterion for melting: Comparison of square-well and Morse potential models for clusters and homopolymers, J. Chem. Phys., 2002, 116, 2323–2329.
69. E. C. Neyts and A. Bogaerts, Numerical study of the size-dependent melting mechanisms of nickel nanoclusters, J. Phys. Chem. C, 2009, 113, 2771–2776.
70. R. Gaillac, P. Pullumbi, K. A. Beyer, K. W. Chapman, D. A. Keen, T. D. Bennett and F.-X. Coudert, Liquid metal–organic frameworks, Nat. Mater., 2017, 16, 1149–1154.
71. R. Gaillac, P. Pullumbi and F.-X. Coudert, Melting of zeolitic imidazolate frameworks with different topologies: Insight from first-principles molecular dynamics, J. Phys. Chem. C, 2018, 122, 6730–6736.
72. N. Castel and F.-X. Coudert, Challenges in molecular dynamics of amorphous ZIFs using reactive force fields, J. Phys. Chem. C, 2022, 126, 19532–19541.
73. K. Zhang, G. M. Stocks and J. Zhong, Melting and premelting of carbon nanotubes, Nanotechnology, 2007, 18, 285703.
74. Y. Qi, T. Çağin, W. L. Johnson and W. A. Goddard, III, Melting and crystallization in Ni nanoclusters: The mesoscale regime, J. Chem. Phys., 2001, 115, 385–394.
75. Y. Shibuta and T. Suzuki, Melting and nucleation of iron nanoparticles: A molecular dynamics study, Chem. Phys. Lett., 2007, 445, 265–270.
76. F. Ding, K. Bolton and A. Rosén, Molecular dynamics study of the surface melting of iron clusters, Eur. Phys. J. D, 2005, 34, 275–277.
77. R. Radev and A. Proykova, Parallel J-W Monte Carlo simulations of thermal phase changes in finite-size systems, Comput. Phys. Commun., 2002, 147, 242–245.
78. Z. Yang, X. Yang and Z. Xu, Molecular dynamics simulation of the melting behavior of Pt-Au nanoparticles with core-shell structure, J. Phys. Chem. C, 2008, 112, 4937–4947.
79. D. Schebarchov, F. Baletto and D. J. Wales, Structure, thermodynamics, and rearrangement mechanisms in gold clusters – insights from the energy landscapes framework, Nanoscale, 2018, 10, 2004–2016.
80. R. Essajai, A. Rachadi, E. Feddi and N. Hassanain, MD simulation-based study on the thermodynamic, structural and liquid properties of gold nanostructures, Mater. Chem. Phys., 2018, 218, 116–121.
81. M. M. Vopson, N. Rogers and I. Hepburn, The generalized Lindemann melting coefficient, Solid State Commun., 2020, 318, 113977.
82. J. M. Keel, J. Yin, Q. Guo and R. E. Palmer, Layer by layer removal of Au atoms from passivated Au(111) surfaces using the scanning tunneling microscope: Nanoscale “paint stripping”, J. Chem. Phys., 2002, 116, 7151–7157.
83. S. Kumar and R. Jin, Water-soluble Au₂₅(Capt)₁₈ nanoclusters: synthesis, thermal stability, and optical properties, Nanoscale, 2012, 4, 4222–4227.
84. T. Kawawaki and Y. Negishi, Precise synthesis of ligand-protected metal nanoclusters for electrochemical/photoelectrochemical applications, Chem. Lett., 2024, 53, upae155.

---