Structure and Thermodynamics of Li⁺Ar_n Clusters beyond the Second Solvation Shell

Abstract

Small Li⁺Ar_n clusters are employed in this work as model systems to study microsolvation. Although first and second solvation shells are expected to be the most relevant ones for this type of atomic solvents, it is also interesting to explore larger clusters in order to identify the influence of external atoms on structural and thermodynamic properties. In this work, we perform a global geometry optimization for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters (with n = 41–100) and parallel tempering Monte Carlo (PTMC) simulations for some selected sizes. The results show that global minimum structures of large clusters always have 6 argon atoms in the first solvation shell while maintaining the number of 14 or 16 argon atoms in the second one. By contrast, third and fourth solvation shells vary significantly the number of argon atoms with the cluster size, and other shells can hardly be assigned due to the reduced influence of ${\mathrm{Li}}^{+}$ on the external argon atoms for large clusters. In turn, PTMC calculations show that the melting of the most external solvation shells of large microsolvation clusters occurs at $T\sim 50\mathrm{K}$, which is independent of cluster size. Structural transitions can be observed between quasi-degenerated structures at low temperatures. Moreover, the present results highlight the fluxional character of the external solvation shells of these large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters, which may be seen as typical “snowball” structures.

1. Introduction

The study of microsolvation is intimately related with the evaluation of the relative stability of clusters resulting from the stepwise addition of solvent to the solute. Energetic and thermodynamics properties of the clusters as well as their structural characterization are thus expected to offer detailed insight about the organization of the solvent around the solute. Since all those properties depend on the nature of the chemical entities of the microsolvation system, the development of a sufficiently accurate potential energy surface (PES) that encodes the relevant interactions among the particles of each cluster it is of primarily importance. Using the PES, the energy and low-energy structures of the clusters may be obtained by employing global optimization methods, while the effect of temperature on the stability of microsolvation aggregates can be assessed through parallel tempering Monte Carlo calculations.

The development of efficient global optimization algorithms [1,2,3,4,5,6,7,8] has been the endeavor of several research groups over the past decades, and some of those methods have been applied to the study microsolvation clusters [9,10,11,12,13,14,15,16,17]. In particular, genetic algorithms [18,19,20] have been applied to search for the global minimum structures of clusters resulting from the microsolvation of alkali metal ions with rare-gas atoms. This kind of systems present distinct solute–solvent and solvent–solvent interactions, because the former are established between the ion and the rare-gas atoms while the others are dominated by dispersion, and hence it is expected that they are weaker. Thus, global optimization results have shown the formation of a first solvation shell with a total number of argon atoms that depend on cation [21]. For instance, a six-atom octahedral structure constitutes the first solvation shell in the case of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters. Using the interaction potential based only on two-body terms [21], it has also been shown that the number of argon atoms in the first solvation shell increases to 8 and 12 when considering the microsolvation of ${\mathrm{Na}}^{+}$ and ${\mathrm{K}}^{+}$, respectively. These results are in agreement with those of a previous study on ${\mathrm{Na}}^{+}$${\mathrm{Ar}}_n$ using an analytical potential that includes three-body terms [14], and they are corroborated by a recent work on ${\mathrm{K}}^{+}$${\mathrm{Ar}}_n$ applying density functional theory [22]. Additionally, the closure of the first solvation shell in these systems is characterized by strong magic numbers [21,23], which may be attenuated by including three-body terms in the potential energy surface [23]. Another fingerprint for identifying structural features of the microsolvation clusters is heat capacity, which indicates [24] distinct temperatures for melting the first and the second solvation shells of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$. In fact, it has been shown that ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters show a rigid first solvation shell in comparison with the fluxional character of the second solvation shell [25].

In this work, we perform a global optimization search of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n(n$ = 41–100) using our genetic algorithm. The study is complemented for some specific cluster sizes with the analysis of the heat capacity as a function of temperature, which is calculated employing the parallel tempering Monte Carlo method. The potential energy surface that models the interactions arising in the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters is the so-called PES II from our previous work [23]. We aim to investigate whether other solvation shells beyond the second one can be observed for large clusters. In addition, PTMC calculations allow assessment of structural transitions and the temperature of dissociation of the clusters. Our main goal is to shed light on the organization of argon atoms around ${\mathrm{Li}}^{+}$ in terms of solvation shells and characterize melting and premelting phenomena that may occur in these large microsolvation clusters.

2. Methods

2.1. Potential Energy Surface

We assume that the energy of each cluster composed by ${\mathrm{Li}}^{+}$ and n argon atoms may be given by a sum over all two-body contributions, i.e.,

$$V_{cluster}=\sum _{i=1}^n{V}_{{\mathrm{Li}}^{+}-\mathrm{Ar}}\left(r_i\right)+\sum _{i=1}^{n(n-1)/2}{V}_{\mathrm{Ar}-\mathrm{Ar}}\left(r_i\right),$$

(1)

where $V_{{\mathrm{Li}}^{+}-\mathrm{Ar}}$ and $V_{\mathrm{Ar}-\mathrm{Ar}}$ are, respectively, the ${\mathrm{Li}}^{+}$-Ar and Ar-Ar pair potentials that depend on the corresponding internuclear distances ($r_i$). The analytical expressions for both pair potentials are the same as in our previous work [23], since they are accurate and also allow a fair comparative analysis of the present results; for completeness, the analytical expressions for the pair potentials are described in Supplementary Information. The ${\mathrm{Li}}^{+}$-Ar potential is fitted to ab initio points calculated with the CCSD(T) method and quadruple-zeta basis sets (i.e., the cc-pVQZ basis set for the ${\mathrm{Li}}^{+}$ ion [26] and the aug-cc-pVQZ for argon [27]); the ${\mathrm{Li}}^{+}$-Ar equilibrium distance is 4.484 ${\mathrm{a}}_0$ and the corresponding well depth is −10.6156 ${\mathrm{mE}}_{\mathrm{h}}$. For the $V_{\mathrm{Ar}-\mathrm{Ar}}$ term, the Rydberg–London potential proposed by Cahill and Parsegian [28] is employed, the parameters of which are obtained by fitting to the accurate HFDID1 model of Aziz [29] for ${\mathrm{Ar}}_2$. The minimum of the Ar-Ar potential is located at ∼$7.100{\mathrm{a}}_0$ and it has the value of about −0.4554 ${\mathrm{mE}}_{\mathrm{h}}$, which is more than 20 times weaker than that of the abovementioned ${\mathrm{Li}}^{+}$-Ar interaction.

Although we have shown [23,25] that three-body interactions are especially relevant for establishing the correct structural and energetic features of small-size ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters, the essential trends are expected to be reproduced by a PES based on pair potentials as the microsolvation system grows up. The reason for that relies on the obvious predominance of weak interactions among argon atoms for larger clusters when compared with the high-energy contribution of the terms involving ${\mathrm{Li}}^{+}$ which dominate for small systems. Thus, we believe that the PES employed in the present work is sufficiently accurate to capture the main structural and thermodynamic features of the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters beyond the second solvation shell. In fact, we observed, in a previous work [23], that global minimum structures are not overly affected by including three-body interactions in the PES for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters with $n\ge 30$. Further evidence supporting this assumption is revealed in Figure S1 of Supplementary Information, where we can observe a decreasing contribution of three-body interactions for the binding energy as the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters become large.

2.2. Evolutionary Algorithm for Structure Optimization

The evolutionary algorithm (EA) used for the geometry optimization of the clusters was developed in our group during more than one decade [18,19,20,30,31]. Across the years, it has been successfully applied to discover the global minima of several challenging systems involving different types of interactions. Among these, we especially refer to the works on Morse clusters [18,20], rare gases [20], transition metal alloys [21,32], charged colloids [31,33], and aggregates resulting from the microsolvation of alkali metal ions with rare-gas atoms [21,23,25,34,35].

The EA described below is applied for each cluster size. The initial step of the algorithm consists of generating a random population of structures; a set of 100 cluster structures is employed in this work. All the structures of the population are relaxed to the bottom of the corresponding attraction basin applying the local optimization method of Broyden–Fletcher–Goldfarb–Shanno (usually designated as the L-BFGS method [36,37]). Then, the global optimization cycle begins with the choice of the structures among the whole population that are the “parents” of the new generation. For that, the EA applies tournament selection with a tourney size of five. This means that each “parent” is chosen to be the best structure among a set of five randomly picked up from the population. The selected parent structures are modified by the application of genetic operators, i.e., crossover and mutation. Crossover combines pairs of structures to form offspring structures, and these can then undergo mutation, which consists of modifying the position of an atom of the cluster. For crossover, which is applied with a probability of 0.7, we rely on a generalized version [30] of the Cut and Splice proposed many years ago by Deaven and Ho [38]. In turn, the present calculations apply the Sigma mutation with a probability of 0.05, whereas 0.1 is the standard deviation ($\sigma$) of the Gaussian distribution from which the new location of the atom can be randomly sampled. Offspring structures are subject to local optimization with the L-BFGS method. In the replacement strategy of the EA, “offspring” and “parents” compete to constitute part of the population for the next generation: low-energy structures have the advantage to enter the pool, but a certain degree of structural diversity is enforced [19] to avoid a premature convergence to a local minimum. Finally, the whole cycle is repeated until the number of potential energy evaluations reaches $1\times {10}^7$. Since the EA is a stochastic algorithm, it is necessary to run it several times in order to obtain an adequate statistical treatment of the global optimization results. We accomplish 30 runs of the EA for each cluster size, and the lowest energy structure so obtained is the corresponding putative global minimum.

2.3. Parallel Tempering Monte Carlo Method

We employ our own implementation of the parallel tempering Monte Carlo (PTMC) method to calculate the heat capacity as a function of temperature for some specific cluster sizes. Since the method has been previously described [24,39], here, we just offer an overview of the main features and present some modifications introduced for this study.

The reduced heat capacity (or adimensional heat capacity per particle) can be calculated as

$$c_V=\frac{3}{2}+\frac{\langle V_{cluster}^2\rangle -{\langle V_{cluster}\rangle }^2}{N{k}^2{T}^2},$$

(2)

where $\langle V_{cluster}\rangle$ and $\langle V_{cluster}^2\rangle$ are the average values of the potential and the squared potential, respectively, computed by the Metropolis algorithm in a Monte Carlo simulation. In Equation (2), $N(=n+1)$, k, and T are, respectively, the number of atoms of the cluster, the Boltzmann constant, and the temperature.

The PTMC method considers M independent simulations for exploring the configuration space of the cluster system through random walks generated according to Boltzmann’s probability distribution. The simulations are run in parallel, each one at a different temperature, and exchanges between configurations of pairs of neighboring temperatures are attempted after a specified number of Monte Carlo steps ($N_{MC,exch}$); in this work, we consider $N_{MC,exch}=50$ and a total number of Monte Carlo moves of $4\times {10}^7$. The acceptance probability for exchanging configurations of different simulations is given by expression

$$A({\mathbf{R}}_i^{\left(m\right)}\to {\mathbf{R}}_i^{\left(p\right)})=\min \left\{1,e^{\left[\left(V_{cluster}\left({\mathbf{R}}_i^{\left(m\right)}\right)-V_{cluster}\left({\mathbf{R}}_i^{\left(p\right)}\right)\right)(1/k{T}_m-1/k{T}_p)\right]}\right\}.$$

(3)

In this equation, ${\mathbf{R}}_i^{\left(m\right)}$ and ${\mathbf{R}}_i^{\left(p\right)}$ are the configurations at the ith step of the Metropolis random walk for temperatures $T_m$ and $T_p$, respectively. The deterministic even-odd (DEO) procedure [40,41] was employed for exchanging structures of simulations at different temperatures. In turn, the step-size adjustment of Swendsen [42] was carried out on the fly to guarantee an acceptance rate of ∼$50\%$ for the MC moves in each simulation. In addition, we should also mention that the PTMC calculation begins with the random selection of a cluster geometry among the five lowest-energy structures obtained from global optimization, and a subsequent thermalization process is performed by employing $2\times {10}^7$ Monte Carlo moves.

Moreover, the set of temperatures ($\left\{T_m\right\}$) was chosen from constant steps in two predefined intervals that cover low- and high-temperature ranges. Typically, we use steps of 0.5 K in the low-temperature range and $\Delta T=10\mathrm{K}$ for temperatures above 200 K (or 150 K, for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{63}$ and ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{85}$).

In order to avoid prompt evaporation of the argon atoms that are less bounded, we add the following rigid-wall potential to the interaction energy:

$$V_{\mathrm{rw}}\left(\mathbf{R}\right)=\left\{\begin{array}{cc}\infty \hfill & |\mathbf{R}-{\mathbf{R}}_{CM}|\ge R_{cut}\hfill \\ 0\hfill & |\mathbf{R}-{\mathbf{R}}_{CM}|<R_{cut}\hfill \end{array}.\right.$$

(4)

In Equation (4), distance vector $\mathbf{R}$ is referred to the vector of the center of mass coordinates of the cluster (${\mathbf{R}}_{CM}$) and $R_{cut}$ is a parameter that establishes the volume of the available configurational space. As in our previous work [24], we define $R_{cut}=\alpha D_{max}$, where $\alpha =1.4$ and $D_{max}$ is the distance of the farthest atom from the center of mass of the initial structure. Since five low-energy minima are considered as possible departing structures for each simulation, one has five values for $D_{max}$. From these, we select the largest $D_{max}$ value for the calculation of $R_{cut}$, which guarantees that the relevant volume of the configurational space for all departing structures can be sampled in the simulation. The values of $R_{cut}$ are shown in Table S1, where we also represent other relevant parameters of PTMC calculations.

The main structures arising in the simulations are identified by storing hundreds of non-correlated geometries for all the temperatures. These geometries are then submitted to local optimization so that minimum-energy structures can be obtained.

3. Results and Discussion

3.1. Global Optimization

We employed the EA described in Section 2.2 to discover all the putative global minimum structures of the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n(n=41–100)$ clusters; the corresponding Cartesian coordinates are given in Supplementary Information. In a previous work [23], we already found out the putative global minima for the lower-size ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n(n=2–40)$ clusters. By using such previous results [23], we illustrate in Figure S2 how the different solvation shells are building up around the lithium ion as the size of the cluster increases. Although the definition of a solvation shell used here (which is solely based on the ${\mathrm{Li}}^{+}$-Ar distance) may be somehow arbitrary, Figure S2 offers at least an idea about the organization of the argon atoms around ${\mathrm{Li}}^{+}$; further details about the distribution of argon atoms by different solvation shells are given in Table S2. Whilst the first solvation shell closes at $n=6$ and the corresponding octahedral motif is kept for larger sizes, second and third solvation shells show more variable shapes among the clusters. For instance, a single argon atom appears in the third solvation shell at $n=19$ and 21, but none are shown for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{20}$, which has the greatest number of atoms in the second solvation shell. Also, structural changes in the second solvation shell can be observed in many cases, including the appearance of a perfect cubic shape for $n=38–40$. It is also interesting to notice the formation of an incipient fourth solvation shell at $n=36$ and 37, with one and three argon atoms, respectively. This disappears at $n=38$ to form the more compact octahedral structure and arises again for larger clusters.

Moreover, we show in Figure 1 how the different solvation shells are organized for eight selected sizes of the large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters studied in this work (i.e., $n\ge 41$); the distribution of argon atoms by different solvation shells are displayed in Table S3. The ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{41}$ global minimum maintains the shell structure observed for $n=38–40$ (by just adding a third atom to the fourth solvation shell), 8 and 24 atoms in the second and third solvation shells, respectively. A significant change is observed for larger clusters whose second and third solvation shells have, in general, 16 and 12 atoms, respectively. These values may slightly vary for some cluster sizes. For instance, ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{63}$ keeps 16 atoms in the second solvation shell, but it has only 10 atoms in the third one; as for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{85}$, the corresponding numbers of atoms are 14 and 15. Indeed, it is apparent from Figure 1 and Table S3 that the maximum numbers of argon atoms in the second and third solvation shells are 16 and 15, respectively, and such values do not change too much for large clusters (i.e., $n\ge 42$). Conversely, the number of argon atoms in the fourth solvation shell can rise up to 24 (and keeps at this value for $n\ge 89$), but it shows a great variation for specific cluster sizes, which is an indication of a very small influence of the ${\mathrm{Li}}^{+}$ ion. This means that one can hardly distinguish among solvation shells beyond the fourth one. Because of this, in the last column of Table S3, we collect the number of atoms that are more external (i.e., $d_{{\mathrm{Li}}^{+}-\mathrm{Ar}}>15.0{\mathrm{a}}_0$). It is worth noting that, for $n\ge 89$, the first four solvation shells do not change their structure, while the number of external argon atoms increases. Nonetheless, we can find sets of consecutive clusters with similar numbers of external argon atoms for $n<89$.

This behavior is perhaps better highlighted by the scatter plot of the ${\mathrm{Li}}^{+}$-Ar distances displayed in Figure 2. Clearly, the first solvation shell is defined by the set of ${\mathrm{Li}}^{+}$-Ar distances slightly above $4{\mathrm{a}}_0$ for all cluster sizes. Conversely, the ${\mathrm{Li}}^{+}$-Ar distances of atoms in the second solvation shell are shown to vary from ∼$8{\mathrm{a}}_0$ to ∼$10{\mathrm{a}}_0$. The variation of the ${\mathrm{Li}}^{+}$-Ar distances in the third solvation shell is smaller, but the separation between this shell and the second one is less than $2{\mathrm{a}}_0$ for most of the global minimum structures. The first argon atom in the fourth solvation shell arises for $n=36$ at a ${\mathrm{Li}}^{+}$-Ar distance visibly longer than those for the atoms in the third solvation shell, and such separation is maintained for larger cluster sizes. The construction of the fourth solvation shell is suddenly interrupted for $n=38$, the global minimum of which is a very compact fcc structure. This compact structure proceeds up to $n=41$, while the fourth solvation shell begins again to be built up. Furthermore, we may observe, in Figure 1, that an even more external argon atom appears for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{42}$. As the microsolvation cluster grows up, these external argon atoms tend to accumulate on one part of the surface (instead of spreading around the structure), mainly nearby the fourth solvation shell so that one can be hardly distinguished from the other. In a certain way, it is legitimate to consider such external argon atoms as forming a subcluster on the surface of the microsolvation cluster, which is essentially not affected by the ion. As shown in Figure 2b, such accumulation of external argon atoms contributes to the lithium ion becoming more off-center in the cluster for $n\ge 42$.

Another interesting feature of these large clusters that can be observed in Figure 2a is the quasi-continuous distribution of the ${\mathrm{Li}}^{+}$-Ar distances for the external argon atoms. Indeed, it is apparent from this figure for clusters above $n=88$ that the external atoms spawn ${\mathrm{Li}}^{+}$-Ar distances very close to each other for $d_{{\mathrm{Li}}^{+}-\mathrm{Ar}}>13.5{\mathrm{a}}_0$. Thus, it becomes impossible to assign a fourth solvation shell for such large clusters.

Now, we turn the discussion to the relative stability of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters. The maxima of the second energy difference, that is,

$${\Delta }_{2}E=-2E_n+E_{n-1}+E_{n+1},$$

(5)

are usually designated as “magic numbers”, since they correspond to clusters that are particularly stable in comparison with the corresponding neighbor sizes; in Equation (5), $E_n$, $E_{n-1}$, $E_{n+1}$ are the energies of clusters with n, $n-1$, and $n+1$ argon atoms, respectively. We observe in Figure 3 that the highest value of ${\Delta }_{2}E$ arises for closure of the first solvation shell at $n=6$. As previously reported [23], another significant magic number appears for the octahedral structure of the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{38}$ cluster. Whereas the neighbor $n=37$ and $n=39$ clusters have atoms in the fourth solvation shell (three and one, respectively), ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{38}$ shows a closed third solvation shell (cf. Figure S2 in Supplementary Information), which may justify its relative stability. It is also apparent from Figure 3 that other important magic numbers arise for larger clusters: ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{57}$, ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{81}$, ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{89}$ and ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{92}$. It is interesting to notice that the magic number for $n=81$ appears to be also present in the mass spectrum of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters reported by Froudakis et al. [43]. In contrast to ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_6$ and ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{38}$ where ${\mathrm{Li}}^{+}$ is at the geometric center of the cluster, the magic numbers of larger systems are assigned to structures with an off-center ion (in comparison with a corresponding neighbor size cluster); see also Figure 2b. Nonetheless, it is interesting to notice that, within the present approximation for PES, symmetry appears to play a role in relative stability, even for large clusters. In fact, all the magic number structures in Figure 3 are more symmetric than the neighbor counterparts (see also the corresponding symmetry point groups given along with the energies and coordinates in Supplementary Information).

3.2. PTMC Calculations

In order to gain insight about structural transitions and dissociation of large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters, we calculated the heat capacity at constant volume as a function of temperature by employing parallel tempering Monte Carlo simulations. The $c_V$ values were obtained by averaging over two runs of the PTMC program (as described in Section 2). The $c_V$ curves so obtained are represented in Figure 4 for six cluster sizes, i.e., $n=41$, 42, 44, 50, 52, 53, 63, and 85. We noticed from this figure that a strong peak arises at $T\sim 50\mathrm{K}$ for all cluster sizes, which may be attributed to the melting of the most external solvation shells. In fact, external argon atoms were loosely attracted by the ${\mathrm{Li}}^{+}$ ion and, as their number increased with cluster size, the corresponding $c_V$ peaks showed a larger amplitude for $n=63$ and $n=85$. A similar trend was observed in a PTMC calculation for rare-gas clusters (including ${\mathrm{Ar}}_n$), where the melting peak became sharper with the increasing size of the cluster [44].

In turn, this behavior in the $c_V$ curve was observed in a previous work [24], where the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters with $n=33$, 34 and 38 also presented a peak around $50\mathrm{K}$. In that paper [24], the peak was associated to the melting of the external solvation shell, which is likely to be similar for the clusters studied here. Furthermore, the amplitude of such peak for smaller cluster sizes in Figure 4 is similar to that previously reported [24] for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{38}$, whose global minimum structure has a completed third solvation shell with 24 argon atoms (see Table S2 and Figure S2). Indeed, the total number of argon atoms in the external solvation shells for $n=41$, 42, 44, 50, 52, and 53 did not differ too much from that value (cf. Table S2).

This prominent peak as well as other significant features that arise in the heat capacity curves can be rationalized by plotting the frequency of appearance of different structural motifs in the simulations as a function of temperature for each cluster size. Thus, we display in Figure 5 the frequency of appearance of various types of structures, namely the global minimum, the second-lowest minimum and, for $n=50$ and $n=53$, the third-lowest minimum detected in the PTMC simulations; such main configurations are shown in Figure 6 for each cluster size. Also included in each panel of Figure 5 are two additional lines, which report the cumulative values of frequency of the structures that appear in the simulation less than 2% of the time (blue curve) and those arising between 2% and 20% of the time at least at one temperature (magenta curve). The former corresponds to structures that present a certain degree of fragmentation (i.e., with atoms separated from the other ones). This curve reaches essentially 100% for all cluster sizes when the temperature approaches $50\mathrm{K}$, which is compatible with the position of the heat capacity prominent peak (Figure 4) already attributed to the melting of external solvation shells.

We observe, in Figure 5, that, in general, the global minimum is the dominant structure at low temperatures. An exception arises for $n=41$, where the two lowest-energy structures (A and B in Figure 6) have similar probabilities at $T=1\mathrm{K}$ and, after a small increase in temperature, structure B becomes the dominant one; a similar feature was previously observed for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{33}$ [24]. Despite an unambiguous structural difference between A and B, i.e., the root mean square deviation (RMSD) is 2.2 Å (see also Figure 6), it should be mentioned that the two structures are energetically quasi-degenerated ($\Delta E=0.01\mathrm{kcal}{\mathrm{mol}}^{-1}$). This is compatible with a solid–solid transition involving mainly the atoms on the surface, similar to what has been also reported for Lennard–Jones clusters [45]; by contrast with ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$, structural transitions in pure argon clusters, such as ${\mathrm{Ar}}_{131}$, ${\mathrm{Ar}}_{132}$, ${\mathrm{Ar}}_{196}$ and ${\mathrm{Ar}}_{309}$, have been associated to strong peaks in the corresponding heat capacity curves [46]. Then, structures A and B maintain approximately the same frequencies (∼$80\%$ for A and ∼$20\%$ for B) up to $T=20\mathrm{K}$, when a set of other structures with relatively small probability (represented by the magenta curve) begin to appear. At $T\sim 26\mathrm{K}$, such structures become dominant over A and B (whose frequencies fall down), which justifies the small shoulder that arises in the corresponding heat capacity curve. A similar structural transition is also perceived in the heat capacity at $T\sim 14\mathrm{K}$ for $n=42$, which corresponds to a dropping value of the frequency of structure A whereas structure B becomes the dominant one (cf. Figure 5). A different situation occurs for $n=44$, where no shoulder (or small maximum) is visible in the heat capacity curve prior the melting peak at around $50\mathrm{K}$. In contrast to $n=42$, the decrease in global minimum frequency at low temperatures is now more gradual and there is not a clear dominance of one type of structure over the others. In addition, the frequency of fragmented structures increases at a slow rate and is already non-negligible at ∼$22\mathrm{K}$, the temperature at which structure B and other low-frequency structures reach their maximum frequency of appearance.

Regarding cluster sizes $n=50,52$ and 53, a structural transition between A and B is only observed for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{52}$ (i.e., a very small shoulder in the heat capacity curve at $T\sim 14\mathrm{K}$). Indeed, the frequency of Structure B increases as the frequency of Structure A decreases, but its maximum value does not reach 50% due to the concomitant appearance of other structures (magenta and blue curves in Figure 5), which explains the very smooth shoulder for $n=52$ in comparison to those for $n=41$ and $n=42$ (cf. Figure 4). As for ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{50}$ and ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_{53}$, a third distinct structure that reaches frequencies above 20% (designated as C in Figure 6) is observed in PTMC simulations at low temperatures. However, none of Structures B and C become dominant as the frequency of Structure A drops down. In fact, low-frequency and dissociated structures (magenta and blue curves in Figure 5) begin to arise at very low temperatures, which leads the corresponding heat capacity curves to increase at a higher rate for $T>12\mathrm{K}$.

It is apparent from Figure 5 that the frequency of Structure A decays at a slower rate for $n=63$ than for $n=85$. Such rapid decrease in $n=85$ is accompanied by a prompt appearance of Structure B, which becomes dominant (i.e., with a frequency greater than 50%) at $T=8\mathrm{K}$, thus leading to a small peak in the corresponding heat capacity curve (see Figure 4). By inspection of Figure 6, we may conclude that the structural transitions reported here correspond to slight rearrangements of the most external argon atoms, which is compatible with small RMSD values and energy differences (less than $1\mathrm{kJ}{\mathrm{mol}}^{-1}$) among such structures. It is worth noting that similar premelting transitions, which are also associated to surface rearrangements, have been reported for Lennard–Jones clusters by several authors [47,48,49,50].

Finally, we notice in Figure 4 that the $c_V$ curve becomes broader after the prominent dissociation peak, which is independent from the cluster size. A similar feature has been observed for smaller ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters and it has been attributed to the melting of the first solvation shell, which is expected to occur at higher temperatures [24].

4. Conclusions

We performed a global optimization study of the clusters resulting from the microsolvation of ${\mathrm{Li}}^{+}$ by argon. Putative global minimum structures were obtained for cluster sizes between $n=41$ and $n=100$. In addition, we also employed the PTMC method to calculate the heat capacity as a function of temperature for some specific cluster sizes (namely $n=41,42,44,50,52,53,63$, and 85).

We concluded from the present study that global minimum structures of large clusters essentially maintain the number of argon atoms in the second solvation shell (i.e., 16 argon atoms). However, a different organization of the argon atoms around ${\mathrm{Li}}^{+}$ occurs for $n=41$ (with 8 argon atoms in the second solvation shell) and, mainly, for $n=85,88–100$ (with 14). Conversely, the third and especially the fourth solvation shells show significant variations in the number of argon atoms with cluster size, which is compatible with very fluxional structures. It is also clear that the third and fourth solvation shells are not totally fulfilled for clusters up to $n=88$. For large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters, the argon atoms are mainly located on specific regions of the cluster instead of spreading all around the structure, and the ion becomes quite off-centered. This effect results from the fact that, at large ${\mathrm{Li}}^{+}$-Ar distances, the influence of the ion is not relevant in comparison with the Ar-Ar interaction, and hence those argon atoms tend to accumulate as close as possible to each other. As expected, this becomes more acute for the most external argon atoms. Thus, we can hardly define a true fourth solvation shell as well as other more external shells. Nonetheless, magic number structures tend to be more symmetric than their neighbor-size clusters.

Regarding PTMC calculations, we noticed that the melting of the most external solvation shells (from the second shell forward) occurs at around 50 K for all the large cluster sizes analyzed in this work. Due to the larger number of external argon atoms, the intensity of the corresponding heat capacity peak was greater for $n=63$ and $n=85$ than for the smaller sizes. Moreover, such a prominent peak of the heat capacity became broader for high temperatures, which is a clear thermodynamic signature of a gradual fusion of the most internal solvation shells. In turn, it was possible to identify two or three (for $n=50$ and $n=53$) most prevalent structures in the simulations at very low temperatures, and structural transitions could be observed for some cluster sizes. Such structures are, however, quasi-degenerated, since the energy difference among them is less than $1\mathrm{kJ}{\mathrm{mol}}^{-1}$. In turn, tiny changes in the slope of the $c_V$ curve that occurred at pre-melting (low) temperatures for some ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters can be essentially associated to rearrangements on the surface of the clusters, which lead to the appearance of structures that mainly differ in the position of argon atoms from the most external solvation shells. In general, the global minimum structure is the prevalent one at very low temperatures, but it is essentially absent from the simulation for $T>30\mathrm{K}$, whereas dissociated clusters begin to appear and dominate the simulation. For low temperatures (typically, below 30 K), it was possible to observe a great amount of structures, each one with low prevalence in the simulation, which is another clear signature of the fluxional character of these large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters.

The results from this work corroborate the idea previously advanced [24] that the first solvation shell of the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters has a rigid-like character, while external solvation shells (including the second one) are essentially fluid-like. In other words, the ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_6$ cluster may be seen as a “snowball”, i.e., similarly to the solid-like sphere of helium atoms surrounding an ion that was proposed many years ago by Johnson and Glaberson [51]. Indeed, the concept of a “snowball” has been applied in the last years not only for the solvation of ions with helium atoms [52]; it has also been extended for the interpretation of ion solvation by molecular solvents [53].

As a final remark, we must say that the study of large ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters suggests the possibility of reaching a value of n from which the presence of the ion becomes irrelevant. The identification of such cluster size may be achieved comparing the melting temperatures of ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ with those obtained for pure ${\mathrm{Ar}}_n$ clusters. We observe, from this work and our previous calculations [24], that the melting temperature decreased from $T\sim 77\mathrm{K}$ for $n=7$ to $T\sim 50\mathrm{K}$ for $n=38$ and remains at this latter value for larger ${\mathrm{Li}}^{+}$${\mathrm{Ar}}_n$ clusters; $n=85$ was the largest cluster size for which the melting temperature was obtained. By contrast, the PTMC results of Pahl et al. [44] for ${\mathrm{Ar}}_n$ clusters showed increasing values of the melting temperature with cluster size (from $T\sim 41\mathrm{K}$ for $n=55$ up to $T\sim 68\mathrm{K}$ for $n=923$, with the experimental melting temperature of the bulk being $83.85,\mathrm{K}$ [54]). In addition, our $c_V$ melting peaks were not so sharp as those for ${\mathrm{Ar}}_n$ for the same cluster size. Thus, we expect that the presence of ${\mathrm{Li}}^{+}$ is irrelevant at sizes beyond the range studied in this work, i.e., for cluster sizes with sharper $c_V$ peaks and melting temperatures approaching those for ${\mathrm{Ar}}_n$ clusters with similar sizes.