Nuclear quantum effects on the stability of cationic neon clusters

Chemical Physics Letters 551 (2012) 38–41

Contents lists available at SciVerse ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Nuclear quantum effects on the stability of cationic neon clusters F. Calvo a,⇑, F.Y. Naumkin b, D.J. Wales c a

LASIM, CNRS UMR 5579, Université de Lyon, 43 Bd. du 11 Novembre 1918, F69622 Villeurbanne Cedex, France Faculty of Science, UOIT, Oshawa, Canada ON L1H 7K4 c University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK b

a r t i c l e

i n f o

Article history: Received 4 August 2012 In final form 6 September 2012 Available online 14 September 2012

a b s t r a c t The stable structures of cationic neon clusters containing up to 57 atoms have been located using a diatomic-in-molecules potential energy surface and basin-hopping hierarchical optimization. The effects of vibrational delocalization were included either in the harmonic approximation, or by performing Langevin molecular dynamics simulations coupled to a quantum thermal bath at T ¼ 0. For most clusters, zeropoint motion is sufficiently high to blur the picture of a single well-defined structure. However, structural diversity of the ground state wavefunction is found to be lower at sizes 14, 21, and 56, which correspond to special stabilities in experimental mass spectra. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Among the rare gases, neon stands in an important position between the highly quantum helium isotopes and the rather classical behavior of the heavier elements argon and krypton. Quantum nuclear effects in bulk neon are usually considered as a minor correction, accounting for some small deviations in the thermodynamic [1–6] or structural [6–8] properties. However, quantum effects in finite size Ne clusters could be more significant, as suggested by various experimental [9–11] and theoretical [7,12–16] works. Zero-point motion is predicted to be large enough to cause structural transitions in (neutral) neon clusters [7,15]. It has also been invoked to explain the unexpectedly high stability of multiply charged Neqþ n clusters [11], although this mechanism is not fully understood [16]. Nuclear quantum effects have also been discussed in the context of evaporation and nucleation [17,18], but could not be tested experimentally. Singly-charged Neþ n clusters are probably the best documented thanks to mass spectrometry investigations [9,10]. They can be experimentally produced by electron bombardment [9] or photoionization [10], and magic numbers were identified for sizes 14, 21, and 56, at variance with the special stabilities for the heavier rare gases resulting from icosahedral packing [19]. The magic numbers of Neþ n clusters can be explained from the formation of a cationic dimer or trimer surrounded by essentially neutral atoms [19– 22]. However, this interpretation relies on a fully classical picture for the nuclei, which could be inconsistent with the known importance of zero-point energy (ZPE) [7,15,16]. In the case of cationic

⇑ Corresponding author. Fax: +33 4 72 43 15 07. E-mail address: fl[email protected] (F. Calvo). 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.09.013

clusters, the covalently bound cationic core is expected to exhibit much higher ZPE than arises from the remaining van der Waals or polarization interactions. We thus anticipate that nuclear quantum effects could be even stronger in Neþ n clusters than in their neutral counterparts. The purpose of this Letter is to evaluate the magnitude of nuclear quantum effects and their role on the stability of cationic neon clusters. We have performed a survey of the energy landscapes for Neþ n clusters, 3 6 n 6 57, analysing large numbers of local minima based on a diatomic-in-molecules (DIM) model for the potential energy surfaces. Nuclear quantum effects were first incorporated in the simplest harmonic approximation. Delocalization of the vibrational ground state wavefunction was then assessed by performing Langevin molecular dynamics (MD) simulations using the recent quantum thermal bath (QTB) approach [23,24]. The QTB scheme offers a simple but efficient way of describing quantum delocalization in many-body systems, which scales much more favorably with system size and temperature than quantum Monte Carlo, path-integral MD, and wavepacket-based methods. As will be shown below, zero-point motion turns out to be quite significant in Neþ n clusters, to the extent that most systems change structure upon inclusion of the ZPE contribution. More importantly, delocalization usually occurs over a set of underlying minima, and the dispersion of these structures exhibits size dependence that correlates with the abundances in experimental mass spectra [9,10]. The Letter is organized as follows. In the next section, we briefly describe the methods used, including some elements of the DIM model and the basics of the Langevin MD scheme with QTB thermostatting. The results are presented and discussed in Section 3, before some concluding remarks are finally outlined in Section 4.

39

F. Calvo et al. / Chemical Physics Letters 551 (2012) 38–41

generated for each cluster size, starting from the quantum (ZPE included) global minimum.

Cationic neon clusters have previously been studied by various computational schemes, ranging from fully explicit treatments of electronic structure [20,21] to empirical models [16,22]. The diatomic-in-molecules [25] approach chosen in the present Letter is a good compromise in terms of chemical accuracy and computational efficiency [21]. In brief, the DIM approach is a configuration-interaction type method providing the electronic ground and excited states of the system from a Hamiltonian exactly divided into the Hamiltonians of the diatomic and atomic fragments [25]. For Neþ n the possibility to locate the positive charge on any of the n atoms and the associated half-occupied p orbital lead to 3n valence-bond type basis functions constructed as simple products of atomic wavefunctions. The DIM matrix also contains some corrections for the interactions between induced dipoles. As in Ref. [21], the diatomic input curves for Neþ 2 were obtained at the coupled-cluster level, but to better deal with large clusters we have attached a long-range polarization correction C 4 =r4 to all Neþ 2 energy curves, with C 4 ¼ aNe =2 ¼ 1:34 Å3, aNe being the atomic polarizability of neon. The Aziz–Slaman potential was used for Ne2 [26]. Evaluating the energy and gradient on the DIM potential energy surfaces involves a 3n  3n DIM matrix whose diagonalization scales as n3 with the number of atoms, making sampling significantly more expensive for large sizes. For n 6 30, favorable structures for the clusters were located by unbiased basin-hopping global optimization, employing five series of 10 000 local minimizations each. Above n ¼ 30, databases of putative global minima were generated using a simpler model inspired by Sebastianelli et al. [22] who assumed a fixed charge localization on Neþ 2 , neglecting charge resonance. Nuclear quantum effects were first included in the simplest harmonic approximation, with the ZPE correction expressed as

DE ¼

X hxi =2;

ð1Þ

i

where fxi g are the 3n  6 vibrational frequencies calculated at the local minimum geometry. The harmonic approximation is useful for detecting possible transitions due to favorable ZPE, but does not give any real insight into the extent of vibrational delocalization. For this reason we have carried out Langevin MD simulations coupled to a quantum thermal bath, a technique recently proposed by several groups [23,24] to describe quantum nuclear effects in many-particle systems at thermal equilibrium. We will not describe the method in detail, but simply emphasize the main concepts. The QTB approach is based on the numerical solution of the Langevin equation in which the random force component is taken from a correlated (colored) noise hðtÞ instead of the usual white noise employed in conventional classical Langevin MD. The correlated noise is designed as a time series in such a way that its power spectral density HðxÞ fulfils the quantum fluctuation–dissipation theorem [24]

HðxÞ ¼ hx

  1 1 þ : 2 expðhx=kB TÞ  1

ð2Þ

We follow the procedure of Barrat and Rodney [27] to generate the noise series hðtÞ on the fly, rather than the original method of Dammak and coworkers [24] which lacks scalability with simulation time. This procedure builds upon the introduction of a filter, HðxÞ, defined in a fixed frequency range Xmax 6 x 6 Xmax over N f bins, and in the present Letter we chose Xmax ¼ 1000 cm1 and N f ¼ 1000. The Langevin equations of motion were solved numeri1 cally using a time step of 1 fs and a damping rate of 104 fs . 10 5 4 trajectories of 10 steps following 5  10 equilibration steps were

3. Results and discussion The total energies per atom, EðnÞ, relative to the atomic energies are shown in Figure 1 for the putative global minima of Neþ n clusters in the size range 3 6 n 6 57. For sizes n 6 25 we did not find differences in the global minimum with respect to the previous work of Ref. [21], which indicates that the long-range potential correction added to the DIM potential does not play a major role in these relatively small clusters. One notable exception is Neþ 13 , for which we found a new global minimum with symmetry D3h lower in energy by 3.4 meV than the original structure. In fact, this structure is also the most stable with the original uncorrected DIM model. The variations of EðnÞ show two main regimes, with a steeper polarization dominated slope until the shell of neutral atoms is completed around the ionic core, and a more gentle slope above n ¼ 14 signaling the onset of the van der Waals dominated regime. Upon including harmonic ZPE, many clusters change structure due to a favorable set of lower vibrational frequencies in some of their higher lying minima. Such structural transitions are increasingly likely with increasing cluster size, which reflects the larger structural diversity in higher dimensional energy landscapes. When the ZPE is included the binding energy shifts by a significant amount, but the two bonding regimes remain clear. Anharmonic ZPEs were also estimated from the Langevin MD runs, using the QTB thermostat in the limit of T ¼ 0 where HðxÞ ¼ hx=2. In such simulations the ZPE is simply the averaged sum of kinetic and potential energies [24]. After further averaging over ten independent trajectories, the binding energies including the anharmonic ZPE contribution are essentially similar to the harmonically corrected values, only slightly higher. A pure static interpretation of relative stabilities based on first ½DE ¼ Eðn þ 1Þ  EðnÞ or second ½D2 E ¼ Eðn þ 1Þ þ Eðn  1Þ  2EðnÞ energy derivatives with the classical or harmonic ZPE indicates some special stabilities at n ¼ 14 and, to a lesser extent, n ¼ 21 and n ¼ 9. However, due to the significant dispersion (errorbars) in the anharmonic ZPE values these derivatives cannot be reliably determined.

0.12

Binding energy (eV)

2. Methods

0.1

Anharmonic ZPE Classical GM Harmonic ZPE Harmonic ZPE (same minimum)

0.08

0.06

0

10

20

30

40

50

60

Cluster size N Figure 1. Binding energy of Neþ n clusters (3 6 n 6 57) predicted by the DIM model, from the classical global minima (black circles) or the minima including the harmonic zero-point energy (red squares, empty squares indicating that the quantum and classical minima differ). The blue diamonds are the results of the Langevin MD simulations with QTB at T ¼ 0, with error bars obtained from ten independent trajectories. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

40

F. Calvo et al. / Chemical Physics Letters 551 (2012) 38–41

-1.85

Probability

Energy (eV)

-1.90

-1.95

Energy fluctuations (meV)

20

0.6 PIMD QTB

0.4 0.2 0.0

-2.03

-2.02

Energy (eV) -2.00

15

56

21

14

10

5

6 0

200

400

600

800

1000

0 0

Quench number

30

40

50

60

Figure 3. Fluctuation of the energies for local minima (inherent structures), dEIS obtained after quenching the Langevin MD trajectories with QTB at T ¼ 0. The red circles highlight the clusters with sizes n ¼ 6, 14, 21, and 56, which exhibit higher vibrational localization, with their (quantum) global minimum also depicted. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

nonzero fluctuations, indicating that delocalization is relatively strong, and cannot be accurately represented as a zero-point energy shift for a single structure. The small fluctuations in dEIS correlate with particularly stable sizes in experimental mass spectra [9,10]. This correlation suggests that nuclear quantum effects are partly responsible for the magic numbers observed in the experiments, in the sense that the most delocalized clusters are also more prone to dissociation. We have finally investigated the specific effects of the charge on the extent of delocalization, taking the 55-atom clusters as an example. In neutral form, this cluster adopts a highly symmetric two-shell Mackay icosahedron [29], which is stable against ZPE corrections. The cationic cluster has a short dimeric core exerting some significant strain on the first surrounding shell and contributing to destabilization of the outer shell. For this cluster, a structural transition due to ZPE is predicted to occur. The average energies of the two clusters were evaluated from Langevin MD simulations with quantum and classical ( h ! 0 leading to a white random noise) thermal baths in the temperature range 0 6 T 6 10

-0.028

Energy (eV)

The stability of individual clusters is better addressed by monitoring the structural diversity within the nuclear wavefunction, as approximated by the Langevin trajectories with QTB at T ¼ 0. More precisely, we have performed periodic quenches along the MD trajectory to unravel the inherent structures burried in our approximate representation of the wavefunction. Figure 2 illustrates the results obtained for a typical trajectory of Neþ 15 , by showing the time variations of the locally minimized energy against the cumulated average of the total energy giving the anharmonic ZPE. Three minima are found by local minimization, having different shell arrangements around a cationic, mostly dimeric core. The presence of several minima shows that the nuclear wavefunction is strongly anharmonic; the simple correction of Eq. (1) to EðnÞ is unable to detect such coexistence. From a more quantitative perspective, the isomer probabilities were obtained from the Langevin simulations, and for comparison we have performed additional path-integral MD simulations of the same Neþ 15 cluster, using the methodology as described in a previous paper [28], at T ¼ 2 K with a Trotter discretization number of P ¼ 32. Quenches were then initiated from the centroid positions. The quenched statistics for the three isomers, depicted as an inset in Figure 2, show good agreement with the results obtained for the Langevin trajectories with QTB. This result validates our computational strategy as a convenient way to incorporate vibrational delocalization. The same quenching analysis carried out for Neþ 14 indicates that only the global minimum structure is visited along the MD trajectory. The slightly higher energetic stability of this cluster is thus also reflected on the absence of structural diversity for the nuclear wavefunction. More generally, we have calculated the fluctuations dEIS of the energies for the local minima from the quenched MD trajectories, and the variations of dEIS with cluster size are represented in Figure 3. The most striking results from dEIS are the low values reached for a few special sizes, and the generally nonzero values for most clusters except n ¼ 3 and n ¼ 14. The trimer has two isomers, but it appears that only the most stable linear global minimum is actually visited. Already for the tetramer, the most weakly bound atom exhibits a significant degree of delocalization around the more tightly bound trimer core. Smaller fluctuations are also found for the 14-, 21-, and 56-atom clusters. However, except for Neþ 14 , the other clusters still show

20

Cluster size

-0.030 -0.032

Ne 55

-0.034 -0.036 -0.090

Energy (eV)

Figure 2. (Black line) Cumulative average of the energy of the Neþ 15 clusters along a Langevin MD trajectory with QTB at T ¼ 0. The empty red circles are the energies of the inherent structure obtained after local minimization, with the three isomers obtained being depicted. The inset shows the probabilities of finding each isomer as a function of its energy, as obtained with the QTB method and by alternative pathintegral MD at T ¼ 2 K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

10

-0.095 +

Ne 55

-0.100 -0.105 -0.110

0

2

4

6

8

10

Temperature (K) Figure 4. Average energy of Ne55 (upper panel) and Neþ 55 (lower panel), as obtained from Langevin MD simulations with quantum (circles) or classical (squares) thermal baths. The dashed lines show the corresponding harmonic results obtained from the classical or quantum global minima.

F. Calvo et al. / Chemical Physics Letters 551 (2012) 38–41

K. The results of these simulations are shown in Figure 4, together with the predictions of the corresponding harmonic result, as obtained from the most stable structures only. Those results illustrate the much larger magnitude of zero-point energy in the cationic cluster, which originates from the tighter covalent interactions in the charged core. The average energy of Neþ 55 obtained from the MD simulations show some shift above the harmonic value, especially in the quantum regime. This shift is a manifestation of higher lying minima in the vibrational wavefunction, whereas Ne55 only exhibits some intrinsic anharmonicity within its single icosahedral basin above about 10 K. 4. Conclusions In this Letter we have explored the interplay between electronic and vibrational delocalization in cationic neon clusters. Although geometric arguments based on icosahedral packing around a dimeric core contributes to explaining the enhanced stability of þ Neþ 14 and Ne56 , we have shown that vibrational delocalization is important and somewhat nullifies the picture of well-defined minima, except for the trimer and 14-mer. Using a refined diatomic-in-molecules quantum model, global optimization and an approximate but efficient simulation technique based on Langevin molecular dynamics with a quantum thermal bath, we obtained insight into the structural diversity of the vibrational wavefunction by quantifying the fluctuations among the local minima (inherent structures). It generally appears that the special stabilities are more connected to these fluctuations than to the binding energies of specific structures. The quantum thermal bath technique used here for neon clusters appears suitable for treating other cryogenic environments, and it would be useful to test it on systems exhibiting even stronger quantum effects, such as parahydrogen. As a relatively inexpensive alternative to path-integral methods, it could also be useful for spectroscopic studies involving quantized modes, such as hydrogen stretches, and perhaps for processes such as proton transfer.

41

Acknowledgments FC wishes to thank Jean-Louis Barrat for useful discussions about the implementation of the quantum thermal bath method. Support from the Pôle Scientifique de Modélisation Numérique (PSMN) in Lyon is gratefully acknowledged. References [1] L.M. Ses’, Mol. Phys. 78 (1993) 1167. [2] M. Asger, Q.N. Usmani, Phys. Rev. B 49 (1994) 12262. [3] J. Solca, A.J. Dyson, G. Steinebrunner, B. Kirchner, H. Huber, J. Chem. Phys. 108 (1998) 4107. [4] M. Venkatraj, C. Bratschi, H. Huber, R.J. Gdanitz, Fluid Phase Eq. 218 (2004) 285. [5] N. Tchouar, F. Ould-Kaddour, D. Levesque, J. Chem. Phys. 121 (2004) 7326. [6] C.P. Herrero, R. Ramirez, Phys. Rev. B 71 (2005) 17411. [7] F. Calvo, J.P.K. Doye, D.J. Wales, J. Chem. Phys. 114 (2001) 7312. [8] M. Neumann, M. Zoppi, Phys. Rev. E 65 (2002) 031203. [9] T.D. Märk, P. Scheier, Chem. Phys. Lett. 137 (1987) 245. [10] M. Fieber, G. Bröker, A. Ding, Z. Phys. D: At. Mol. Clusters 20 (1991) 21. [11] I. Mähr, F. Zappa, S. Denifl, D. Kubala, O. Echt, T.D. Märk, P. Scheier, Phys. Rev. Lett. 98 (2007) 023401. [12] E. Pahl, F. Calvo, L. Koçi, P. Schwerdtfeger, Angew. Chem. Int. Ed. 47 (2008) 8207. [13] C. Chakravarti, J. Chem. Phys. 102 (1995) 956. [14] J.P. Neirotti, D.L. Freeman, J.D. Doll, J. Chem. Phys. 112 (2000) 3990. [15] P.A. Frantsuzov, D. Meluzzi, V.A. Mandelshtam, Phys. Rev. Lett. 96 (2006) 113401. [16] F. Calvo, J. Phys. Chem. Lett. 1 (2010) 2637. [17] F. Calvo, P. Parneix, J. Phys. Chem. A 113 (2009) 14352. [18] W. Unn-Toc, N. Halberstadt, C. Meier, M. Mella, J. Chem. Phys. 137 (2012) 014304. [19] J.M. Soler, N. Saenz, N. Garcia, O. Echt, Chem. Phys. Lett. 109 (1984) 71. [20] M. Fieber, A.M.G. Ding, P.J. Kuntz, Z. Phys. D: At. Mol. Clusters 23 (1992) 171. [21] F.Y. Naumkin, D.J. Wales, Mol. Phys. 93 (1998) 633. [22] F. Sebastianelli, F.A. Gianturco, E. Yurtsever, Chem. Phys. 290 (2003) 279. [23] M. Ceriotti, G. Bussi, M. Parrinello, Phys. Rev. Lett. 102 (2009) 020601. [24] H. Dammak, Y. Chalopin, M. Laroche, M. Hayoun, J.-J. Greffet, Phys. Rev. Lett. 103 (2009) 190601. [25] F.O. Ellison, J. Am. Chem. Soc. 85 (1963) 3540. [26] R.A. Aziz, M.J. Slaman, Chem. Phys. 130 (1989) 137. [27] J.-L. Barrat, D. Rodney, J. Stat. Phys. 144 (2011) 679. [28] F. Calvo, F.Y. Naumkin, D.J. Wales, J. Chem. Phys. 135 (2011) 124308. [29] A.L. Mackay, Acta Cryst. 15 (1962) 916.