Ground state energies and chemical potentials of parahydrogen clusters are
calculated from 3 to 40 molecules using the diffusion Monte Carlo technique
with two different pH2-pH2 interactions. This calculation improves a previous
one by the inclusion of three-body correlations in the importance sampling, by
the time step adjustement and by a better estimation of the statistical errors.
Apart from the cluster with 13 molecules, no other magic clusters are
predicted, in contrast with path integral Monte Carlo results

F. Baletto

I.F. Silvera

J. Tang

J. Vrbik

K.E. Schmidt

P. Sindzingre

S. Baroni

S. Grebenev

S.A. Chin

U. Buck

Open Physics

English

arXiv

Ground state energies and chemical potentials of parahydrogen clusters are calculated from 3 to 40 molecules using the diffusion Monte Carlo technique with two different pH2-pH2 interactions. This calculation improves a previous one by the inclusion of three-body correlations in the importance sampling, by the time step adjustement and by a better estimation of the statistical errors. Apart from the cluster with 13 molecules, no other magic clusters are predicted, in contrast with path integral Monte Carlo results

Guardiola Barcena, Rafael

Navarro Faus, Jesús

Repositori d'Objectes Digitals per a l'Ensenyament la Recerca i la Cultura

arXiv:0709.4557v1  [physics.atm-clus]  28 Sep 2007A diffusion Monte Carlo study of small para-Hydrogen clustersRafael Guardiola1 and Jesu´s Navarro21 Departamento de F´ısica Ato´mica y Nuclear, Facultad de F´ısica, 46100 Burjassot, Spain2 IFIC (CSIC-Universidad de Valencia), Apdo. 22085, 46071 Valencia, Spain(Dated: January 21, 2014)Ground state energies and chemical potentials of parahydrogen clusters are calculated from 3 to 40molecules using the diffusion Monte Carlo technique with two different p-H2–p-H2 interactions. Thiscalculation improves a previous one by the inclusion of three-body correlations in the importancesampling, by the time step adjustement and by a better estimation of the statistical errors. Apartfrom the cluster with 13 molecules, no other magic clusters are predicted, in contrast with pathintegral Monte Carlo results.PACS numbers: 67.40.Db, 36.40.-c, 61.46.BcTheoretical studies of parahydrogen clusters have at-tracted a growing interest in the past years, partly mo-tivated by a recent experiment [1] in which Raman scat-tering was used in cryogenic free jets of the pure gas.Small changes in the frequency near the Q1(0) line of themonomer, were observed and interpreted as intermolecu-lar effects on the intramolecular potential. (p-H2)N clus-ters with N = 2 − 8 were clearly identified through fre-quency shifts ranging from ∆ν = −0.40cm−1 for N = 2to ∆ν = −2.35cm−1 for N = 8. The experiment alsoshowed a bump at N = 13, N = 33 and N = 55, whichwere interpreted as a signal of magical clusters. However,it is worth mentioning that these three values are actu-ally extrapolations from smaller clusters and presumablyare approximate.Magical numbers appear in classical Lennard-Jonesclusters, related to geometrical shapes [2]. Several papershave appeared in the last year with the main objective ofchecking the magical numbers found in Ref. [1], and/orstudying possible superfluidity effects in parahydrogenclusters. Indeed, Path Integral Monte Carlo (PIMC) cal-culations [3] have found a large superfluid fraction inclusters with N=13 and 18 at temperatures T ≤ 2 K.A superfluid response has been observed in small clus-ters consisting of a carbonyl sulfide cromophore sur-rounded by 15-17 p-H2 molecules, all within a large he-lium droplet [4]. This has been confirmed by several MCsimulations [5, 6, 7, 8, 9] of doped p-H2 clusters.Systematic studies of (p-H2)N clusters, covering therange from N = 3 to N = 50 molecules, have beendone based on powerful many-body techniques, as di-fusion Monte Carlo (DMC) [10], PIMC [12, 13, 14], andPIMC adapted to the ground state (PIGS) [11]. Whereasup to N ≃ 22 all these calculations are substantially inagreement, for heavier clusters there are noticeable dif-ferences between DMC and PIMC results, particularlyfor N ≥ 26. PIMC chemical potentials show very promi-nent peaks at N=26, 29, 34 and 39, in contrast with thesmooth behavior obtained with DMC.In this work we present new DMC calculations, im-proving our previous ones [10] so as to get very preciseresults within our computational capacity. Specifically,we consider three aspects: the importance sampling func-tion, the time step adjustement and the estimation of thestatistical errors.The DMC procedure is significantly improved whenusing a good importance sampling wave function, themain effect being the reduction of the variance of thestochastic procedure. We have used a Jastrow functionwith two- and three-body correlations:ΦT = exp(u¯2 + u3), (1)whereu¯2 =∑i<j[u2(rij) + λT ξ2(rij)r2ij], (2)withu2(r) = −∑i<j[p5r5ij+p1rij], (3)ξ(r) = exp(−(r − sT )2w2T), (4)andu3 = −λT2∑ℓGℓGℓ, (5)withGℓ =∑i6=ℓξ(rli)rli. (6)Indices i, j, l run over the number of molecules in thecluster. This function is described in terms of five vari-ational parameters, p5, p1, sT , ωT , and λT . In our pre-vious calculations [10] we used the standard two-bodyJastrow function ΦT = exp(u2). The present trial func-tion includes an enlarged the two-body variational spaceand also three-body correlations in the form suggested inRef. [19], which still has O(N2) computational complex-ity.DMC is based in a short-time approximation ofthe Green’s function related to the imaginary time2Schro¨dinger equation. In this way, an initial wave func-tion ΦT (t = 0) evolves to the exact ground state wavefunction Ψ at large t after many short-time steps τ . Wehave used the O(τ3) approximation to the Green’s func-tion as described in Refs. [17, 18], which provides ener-gies O(τ2). The time step adjustment is the following:from calculations at the relative large steps 0.001 and0.0005K−1, we obtain the Richardson extrapolated value13(4E(0.0005)− E(0.001)) , (7)based on the τ expansion E(τ) = E(0)+Cτ2+ · · · . Thisvalue turns out to be very close to the calculations withmuch smaller time steps, as it may be checked in the lastthree rows of Table I. This checks that the algorithmbehaves as O(τ2), as expected, and suggests to use thevalue τ = 0.0001K−1 for massive calculations with a neg-ligible bias.TABLE I: Determination of optimal time step τ from differ-ent evaluations of the binding energy B(N) of several clusters.The row labelled R.E. is the Richardson extrapolated valueobtained from the previous two rows. The statistical standarddeviation is indicated in parenthesis (error in the final digitshown). Energies and statistical errors are in K.τ B(10) B(20) B(30)0.001 183.47(5) 559.28(17) 1006.4(3)0.0005 185.91(6) 566.56(17) 1020.0(4)R.E. 186.72(9) 568.99(28) 1024.5(5)0.0001 186.93(6) 569.16(12) 1025.2(2)0.00002 186.72(3) 569.48(7) 1024.8(1)A further improvement of the calculation regards theestimate of the statistical error. Because of the sequentialMarkov chain nature of Monte Carlo algorithms, succes-sive samples are strongly correlated, and the typical wayof estimating the variance, σ2 = 〈H2〉 − 〈H〉2, may betoo optimistic. To avoid these correlations we computeda number of times (10, typically), the binding energies,with independent and randomized runs, and estimate thevariance from these results. This requires a considerableincreasing in computational time, but the obtained stan-dard deviations are very precisely computed. Specifically,we have used 1000 walkers with 105 steps plus 20000 sta-bilization steps in each walker.Hydrogen molecules interact through weak van derWaals forces that, nevertheless, are sufficiently strong tobound clusters with any number of molecules. Severalforms have been derived to describe the p-H2–p-H2 inter-action. Two of them are of particular interest becausethey combine ab initio properties with properties of thegas (or solid) as well as experimental information fromcollisions, one due to Silvera and Goldman [15], and theother to Buck et. al. [16], hereafter referred to as SGand BHKOS, respectively. The main difference amongthem is that the former contains a repulsive long-rangeterm (c9/r9) with the objective of providing an approxi-mation for the effective potential in a solid. Recent cal-culations have employed BHKOS potential [10], SG po-tential [12, 13, 14] or both [11]. We present here resultswith both interactions.TABLE II: DMC ground state (p-H2)N binding energies (inK) obtained with BHKOS interaction [16].N N N2 4.3114 15 371.17(3) 28 931.46(18)3 14.66(1) 16 408.56(5) 29 978.50(10)4 30.50(1) 17 446.68(3) 30 1025.66(14)5 50.33(1) 18 486.34(8) 31 1074.03(30)6 73.40(1) 19 527.46(8) 32 1122.19(23)7 98.76(1) 20 569.72(6) 33 1170.42(38)8 126.16(1) 21 612.27(11) 34 1219.26(25)9 155.51(2) 22 655.65(12) 35 1267.74(35)10 186.86(3) 23 700.50(13) 36 1317.87(15)11 220.98(2) 24 745.63(14) 37 1366.26(42)12 257.94(2) 25 791.57(16) 38 1416.12(22)13 297.80(9) 26 837.88(29) 39 1465.55(42)14 334.24(5) 27 884.71(15) 40 1516.00(37)TABLE III: DMC ground state (p-H2)N binding energies (inK) obtained with the SG interaction [15].N N N2 3.8456 15 341.12(5) 28 860.69(12)3 13.22(1) 16 375.81(6) 29 904.56(10)4 27.59(1) 17 411.56(4) 30 949.06(21)5 45.73(1) 18 448.24(5) 31 993.46(12)6 66.87(1) 19 486.31(5) 32 1038.04(25)7 90.14(1) 20 525.35(7) 33 1083.36(22)8 115.42(1) 21 564.92(10) 34 1128.38(32)9 142.46(2) 22 605.40(8) 35 1174.29(23)10 171.41(3) 23 646.60(8) 36 1220.32(22)11 202.74(1) 24 688.50(11) 37 1266.42(20)12 236.56(4) 25 731.02(11) 38 1312.22(21)13 272.53(5) 26 774.08(7) 39 1358.46(29)14 306.63(4) 27 817.34(14) 40 1405.45(27)The calculated DMC ground state binding energies aredisplayed in Tables II and III, for BHKOS and SG poten-tials, respectively. As usual, the numbers in parenthesisare the errors in the final digit shown, and correspondto the standard deviation. The binding energies for thedimer have been obtained by numerical integration of theSchro¨dinger equation.The inclusion of triplet correlations in the importancesampling function leads to a noticeably improvement ofthe variational energies, as showed in Fig. 1 for BHKOSpotential. The DMC energies are basically the same asthose of Ref. [10], with slightly more binding in the heav-ier clusters [20]. The main difference lies in the reductionof the standard deviation by typically a factor of 2.The total binding energies grow monotonically withthe number of constituents. In order to determine anenhanced stability related to magic sizes it is convenient3FIG. 1: (Color online) Comparison of the binding energiesper molecule obtained with diffusion Monte Carlo (DMC) anda variational Monte Carlo (VMC) calculation based on Eq. (1)with two- and three-body correlations. The interaction isBHKOS. The error bars are within the symbol size.0 5 10 15 20 25 3001020304050Number of moleculesBinding energy per particle (K)VMC two-bodyVMC two- and three-bodyDMCto analyze the variation with N of the dissociation en-ergy or chemical potential, defined from the ground stateenergies E(N) asµN = E(N − 1)− E(N) . (8)This quantity is plotted in Fig. 2 as a function of thenumber of molecules N . The main physical result of thisfigure is the presence of a neat peak at N = 13, indicat-ing the magical character of this cluster. Although thetwo used interactions give different total energies, theBHKOS potential providing more binding than the SGone, this peak is present for both interactions.FIG. 2: (Color online) DMC chemical potential (in K) of(p-H2)N clusters as a function of the number N of molecules.The error bars are within the symbol size.0 5 10 15 20 25 30 35 400102030405060Number of moleculesChemical potential (K)BHKOSSGBeyond N = 13 our calculations do not show any clearsignal of local enhancement of the chemical potential. Itshould be mentioned that even if the relative error of thetotal energies is around 10−4, the relative error in theseparation energies can be as high as 10−2, as a conse-quence of the strong cancelations appearing when com-puting µ. So, even after our formidable numerical effort,the absolute error of µ for N ≃ 40 may be as high as0.5 K. Having this fact in mind, the only possible struc-ture, apart from N=13, is N=36 for BHKOS potential.However, µ36 is ≃ 1.70± 0.6 K higher than its neighborsµ35,37 and one cannot exclude that it could simply bea statistical fluctuation. Consequently the only magicalcluster firmly established here is N=13, independent ofthe interaction.FIG. 3: (Color online) DMC and PIMC chemical potentials of(p-H2)N clusters as a function of the number N of molecules,calculated with SG interaction [15]. Filled circles and squaresare PIMC results of Ref. [14], the solid line corresponds toPIMC results of Ref. [13]. Error bars have not been drawnfor the clarity of presentation.0 5 10 15 20 25 30 35 400102030405060Number of moleculesChemical potential (K)T=0.5KT=1.5KPresent (T=0)T=1KIn contrast, PIMC calculations show very prominentpeaks at N = 26, 29, 32, 34 and 39, as show in Fig. 3.The PIMC results at T = 0.5 K and 1.5 K are fromRef. [14], and those at T = 1 K from Ref. [13]. Upto N ≈ 25, our DMC results are indistinguishable fromthe the PIMC ones of Ref. 14 at T = 0.5 K or those ofRef. 13 at T = 1 K. It is worth noticing that these PIMCresults at T = 0.5 and 1 K are essentially identical inthe calculated range of cluster sizes, with the noticeableexceptions of N = 23, 35 and 36. However, it should bekept in mind that PIMC error bars have not been drawnin Fig. 3 for the sake of clarity.The existence of peaks in the chemical potential seemsto be related to thermal effects. These could manifest inenhanced stability thresholds at finite temperature, sim-ilarly to what has been observed in 4He droplets [21].But according to Ref. [13] such thermal effects shouldbe associated to a coexistence of solid-like and liquid-likephases, with a dominance of the latter at low T, as aresult of both the zero-point motion and quantum per-mutation exchanges. To this respect it is worth recalling4that while DMC may be affected by the constraint im-posed by the importance sampling function, PIMC hasno such constraint. Actually, our importance samplingfunction is of the type used to describe liquid-like clus-ters. We have checked that the DMC ground state energyfor the magical (p-H2)13 cluster does not change employ-ing instead an importance sampling function where themolecules are localized at the vertex of the correspondingtruncated MacKay icosahedra. This cluster is definitelyliquid-like, in agreement with PIMC results. It wouldbe interesting to perform a similar solid-like DMC calcu-lation for N ≥ 26. Although very computationally de-manding, such a calculation could be useful to ascertainthe phase of p-H2 clusters in this size region.The authors acknowledge stimulating conversationswith J.P. Toennies and S. Montero and correspon-dence with J.E. Cuervo. This work is supportedby grants FIS2004-00912 (MCyT/FEDER, Spain), andACOMP07-003 (Generalitat Valenciana, Spain).[1] G. Tejeda, J.M. Ferna´ndez, S. Montero, D. Blume, andJ. P. Toennies, Phys. Rev. Lett. 92, 223401 (2004).[2] F. Baletto and R. Ferrando, Rev. Mod. Phys. 77, 319(2005).[3] P. Sindzingre, D.M.Ceperley and M.L.Klein, Phys. Rev.Lett. 67, 1871 (1991).[4] S. Grebenev, B. Sartakov, J.P. Toennies, and A.F.Vilesov, Science 289, 1532 (2000).[5] Y. Kwon and K.B. Whaley, Phys. Rev. Lett. 89, 273401(2002).[6] F. Paesani, R.E. Zillich and K.B. Whaley, J. Chem. Phys.119, 11682 (2003).[7] J. Tang and A.R.W. McKellar, J. Chem. Phys. 121, 3087(2004).[8] F. Paesani, R.E. Zillich, Y. Kwon, and K.B. Whaley, J.Chem. Phys. 122, 181106 (2005).[9] S. Baroni and S. Moroni, ChemPhysChem 6, 1884 (2005].[10] R. Guardiola and J. Navarro, Phys. Rev. A 74, 025201(2006).[11] J.E. Cuervo and P.N. Roy, J. Chem. Phys. 125, 124314(2006).[12] F. Mezzacapo and M. Boninsegni, Phys. Rev. Lett. 97,045301 (2006).[13] F. Mezzacapo and M. Boninsegni, Phys. Rev. A 75,033201 (2007).[14] S.A. Khairallah, M.B. Sevryuk, D.M.Ceperley, and J. P.Toennies, Phys. Rev. Lett. 98, 183401 (2007).[15] I.F. Silvera, V.V. Goldman, J. Chem. Phys. 69, 4209(1978)[16] U. Buck, F. Huisken, A. Kohlhase, D. Otten, and J. Scha-effer, J. Chem. Phys. 78, 4439 (1983)[17] J. Vrbik and S.M. Rothstein, J. Comput. Phys. 63, 130(1986).[18] S.A. Chin, Phys. Rev. A 42, 6991 (1990).[19] K.E. Schmidt, M.A. Lee, M.H. Kalos, and G.V. Chester,Phys. Rev. Lett. 47, 807 (1981).[20] There is a misprint in Table I of Ref. [10]. The entryfor B/N at N = 23 should read 30.43(2), instead of thequoted 29.94(2). We acknowledge J.E. Cuervo for havingdetected it. See also Ref. [13].[21] R. Bru¨hl, R. Guardiola, A. Kalinin, O. Kornilov, J.Navarro, T. Savas, and J.P. Toennies, Phys. Rev. Lett.92, 185301 (2004).

A diffusion Monte Carlo study of small para-Hydrogen clusters

Abstract

Similar works

Full text

Available Versions

Repositori d'Objectes Digitals per a l'Ensenyament la Recerca i la Cultura

Crossref

Directory of Open Access Journals