Examination of the Feynman-Hibbs Approach in the Study of Ne N -Coronene Clusters at Low Temperatures
R. Rodríguez-Cantano, R. Pérez de Tudela, M. Bartolomei, M. I. Hernández, J. Campos-Martínez, T. González-Lezana, P. Villarreal, J. Hernández-Rojas, J. Bretón
EExamination of the Feynman-Hibbs Approach in the Study ofNe N -Coronene Clusters at Low Temperatures Roc´ıo Rodr´ıguez-Cantano, Massimiliano Bartolomei, Marta I. Hern´andez, ∗ Jos´e Campos-Mart´ınez, Tom´as Gonz´alez-Lezana, and Pablo Villarreal
Instituto de F´ısica Fundamental (IFF-CSIC), Serrano 123, 28006 Madrid, Spain
Ricardo P´erez de Tudela
Lehrstuhl f¨ur Theoretische Chemie, Ruhr-Universit¨at Bochum, 44780 Bochum, Germany
Javier Hern´andez-Rojas and Jos´e Bret´on
Departamento de F´ısica and IUdEA,Universidad de La Laguna, 38203 Tenerife, Spain (Dated: August 28, 2018) a r X i v : . [ phy s i c s . a t m - c l u s ] A p r bstract Feynman-Hibbs (FH) effective potentials constitute an appealing approach for investigations ofmany-body systems at thermal equilibrium since they allow us to easily include quantum correctionswithin standard classical simulations. In this work we apply the FH formulation to the study ofNe N -coronene clusters ( N = 1-4, 14) in the 2-14 K temperature range. Quadratic (FH2) and quartic(FH4) contributions to the effective potentials are built upon Ne-Ne and Ne-coronene analyticalpotentials. In particular, a new corrected expression for the FH4 effective potential is reported.FH2 and FH4 cluster energies and structures -obtained from energy optimization through a basin-hoping algorithm as well as classical Monte Carlo simulations- are reported and compared withreference path integral Monte Carlo calculations. For temperatures T > N isotropic harmonic oscillators allows us to propose a means ofestimating the cut-off temperature for the validity of the method, which is found to increase withthe number of atoms adsorbed on the coronene molecule.KEYWORDS: quantum effects, Feynman-Hibbs potentials, van der Waals clusters, polycyclicaromatic hydrocarbons ∗ Electronic address: marta@iff.csic.es . INTRODUCTION Quantum effects are very important in a variety of many-body systems at thermal equi-librium, especially for light molecules and/or low temperatures. Among the problems ofinterest, we mention the behavior of water[1], the storage of small molecules in nanoporousmaterials[2–5] and the study of molecules trapped inside low temperature matrices, eithersolid[6, 7] or superfluid as in the case of He nanodroplets[8, 9]. The path-integral formula-tion of statistical mechanics[10] has provided the framework for the development of accuratepath-integral Monte Carlo (PIMC) methods[11, 12] to study these systems. However, thesemethods are computationally very demanding so various approximate approaches have beendeveloped over the years to overcome this drawback.One of the simplest approximations is the use of Feynman-Hibbs (FH) effectivepotentials[13] in classical Monte Carlo (CMC) or molecular dynamics simulations. Thesepotentials are given as a temperature- and mass-dependent expansion of the intermolecularpotentials in powers of ¯ h . In this way, the approach provides an easy and appealing meansto include quantum corrections in a purely classical simulation. This formulation has beenapplied to the study of both homogeneous[14–17] and heterogeneous systems as, for instance,the sieving of H and D in microporous materials[5, 18–22]. In the case of homogeneousmedia, the validity of this method has been investigated in detail by Ses´e[23, 24] from com-parisons with “exact” PIMC calculations of Lennard-Jones systems (He, Ne, Ar, D , CH ).Similarly, Calvo et al [15] found that the quadratic FH effective potential reproduces quitewell the thermodynamic properties (melting temperatures, heat capacities, caloric curves,etc.) of Ne, Ar, and Xe rare gas clusters. More recently, Kowalczyk et al [25] assessedthe FH approach for supercritical He at 10 K and concluded that the FH potentials areonly suitable at low fluid densities, suggesting that previous applications to dense para-H in nanoporous materials at low temperatures should be revised. It is thus interesting toinvestigate the validity of the FH approach for heterogeneous systems such as fluid-solidmixtures, for which specific studies are scarcer.3n this work we study the performance of the FH approach for Ne N -coronene ( N =1 − ,
14) clusters at temperatures ranging from 2 to 14 K, by comparing CMC calculationsusing FH potentials with the accurate PIMC method. This system can be considered as aprototype for studies of van der Waals interactions between small molecules and carbona-ceous substrates[26] and, in this way, its findings may serve as a guide for future simulationsof the storage of gases by new porous carbon materials[27]. In addition, the interactionbetween coronene and other polycyclic aromatic hydrocarbons (PAHs) with Ne atoms isinteresting in connection with the spectroscopy of these molecules in Neon matrices at lowtemperatures ( ≈ h ) and quartic (¯ h )order, hereafter referred as FH2 and FH4 potentials, respectively. As in our recent study ofHe N -coronene clusters[31], minimum energies of the bare interaction potential are obtainedby means of a basin-hopping (BH) approach[32], and the optimized structures are thenused as seeds for the CMC and PIMC calculations. In this work, we additionally run BHand CMC calculations with the FH2 and FH4 potentials in order to assess the extent ofimprovement with respect to the use of the bare potentials.The paper is organized as follows. Section 2 presents the FH effective potentials ap-plied to the Ne N -coronene interaction. BH, CMC and PIMC computational methods arebriefly reviewed in Section 3. Results (energies and structures at different temperatures) arereported and discussed in Section 4. Finally, concluding remarks are given in Section 5. II. NE N -CORONENE INTERACTION POTENTIALSA. Bare potentials The coronene molecule is assumed to be rigid and fixed to the reference frame. Theorigin of the coordinate system is placed at the coronene center of mass, with the z axis4eing perpendicular to the molecular plane and the x axis being overimposed to two of theC-C bonds of this molecule. The position of i − th Ne atom is given by the Cartesian vector r i . The total potential of the Ne N -coronene system is given as a sum of pairwise interactions, V ( r , ..., r N ) = N (cid:88) i =1 V Ne − Cor ( r i ) + N (cid:88) i The FH effective potentials for an A-B interacting system are obtained from the followinggaussian average of the V AB potential [13, 23] V FHAB ( r ) = (cid:18) µπ ¯ h β (cid:19) / (cid:90) d u V AB ( r + u ) exp (cid:18) − µβ ¯ h u (cid:19) , (6)where r is the vector joining A and B, µ is the A-B reduced mass, and β = 1 /k B T , k B being the Boltzmann constant. Eq.6 is based on a variational treatment of the path-integraland involves the neglect of exchange effects[23]. Closed expressions of the FH potential6re obtained upon a Taylor expansion of V AB ( r + u ) around the vector r in powers of theCartesian components of u and solving the corresponding integrals of Eq.6. We have[24] V FH(2 p )AB ( r ) = p (cid:88) n =0 n ! (cid:18) β ¯ h µ (cid:19) n ∇ n [ V AB ( r )] , (7)where ∇ = 1 and p = 0 , , u . It can be seen that at high temperatures theeffective potentials tend to the bare potential and the classical regime is correctly recovered.However, as temperature tends to zero these potentials become unphysically high, whichwill inevitably impose limits to the model.In the case of the Ne-Ne interaction, V AB ( r ) depends on | r | ≡ ρ and the FH2 potentialbecomes V FH2Ne − Ne ( ρ ) = V Ne − Ne ( ρ ) + ¯ h β µ (cid:18) V (cid:48)(cid:48) Ne − Ne ( ρ ) + 2 ρ V (cid:48) Ne − Ne ( ρ ) (cid:19) , (8)where V Ne − Ne is the potential of Eq.2 and V (cid:48) Ne − Ne and V (cid:48)(cid:48) Ne − Ne are their first and secondderivatives, respectively, with respect to ρ .The quartic Ne-Ne effective potential (FH4) writes V FH4Ne − Ne ( ρ ) = V FH2Ne − Ne + 12 (cid:18) ¯ h β µ (cid:19) (cid:18) V (cid:48)(cid:48)(cid:48)(cid:48) Ne − Ne ( ρ ) + 4 ρ V (cid:48)(cid:48)(cid:48) Ne − Ne ( ρ ) (cid:19) , (9)where, analogously, V (cid:48)(cid:48)(cid:48) Ne − Ne and V (cid:48)(cid:48)(cid:48)(cid:48) Ne − Ne are the third and fourth derivatives of V Ne − Ne withrespect to ρ .It is worth pointing our that Eq. 9 differs from previously published expressions[18, 19,21, 25], which we believe are incorrect. More details are provided in the Appendix.The form of the Ne-coronene FH2 potential is more complicated because the atom-bondpotentials U k of Eq.4 depend both on the distance to the bond center ρ k and the cosine ofthe Jacobian angle, c k . The result is 7 FH2Ne − Cor ( r ) = V Ne − Cor ( r )+ ¯ h β m Ne (cid:88) k (cid:26) ∂ U k ∂ρ k + 2 ρ k ∂U k ∂ρ k + (1 − c k ) ρ k ∂ U k ∂c k − c k ρ k ∂U k ∂c k (cid:27) . (10)Note that in Eq.10 the Ne mass is written instead of the reduced mass of Ne-coronene. Wehave made this choice because the coronene molecule is fixed to the reference frame and itis thus assumed that it has an infinite mass.We have found that the contribution to V FH2Ne − Cor due to the derivatives of U k with respectto c k are very small. To simplify the calculations, we have assumed that the contributionsfrom these derivatives to the terms of ¯ h order are negligible. In this way, the FH4 potentialis written as V FH4Ne − Cor ( r ) = V FH2Ne − Cor ( r ) + 12 (cid:18) ¯ h β m Ne (cid:19) (cid:88) k (cid:26) ∂ U k ∂ρ k + 4 ρ k ∂ U k ∂ρ k (cid:27) . (11)Nevertheless, the exact expression for the correction of ¯ h order is given in Eq. 25, and ithas been additionally checked that the contribution to that equation of the derivatives ofthe atom-bond potential with respect to c k is negligible.More details about the derivation of Eqs. 8-11 are provided in the Appendix. All therequired derivatives of the ILJ functions were obtained analytically as well as the gradientsof the potentials of Eqs.8-11 (needed for the application of the BH algorithm referred be-low). In any case, recent advances in techniques of automatic differentiation[36] might proveadvantageous for this kind of simulations.The Ne-Ne bare potential as well as the FH2 and FH4 ones at 6 K are presented in Fig. 1.The ILJ form of the bare potential is very realistic, as we have compared it with the Tang-Toennies potential[37] and found that both potentials would appear as indistinguishable inFig. 1. Indeed, the well depth and equilibrium distance of the Tang-Toennies potential are3.646 meV and 3.090 ˚A, very close to the values of Table I. On the other hand, it can beseen in Fig. 1 that the FH corrections significantly modify the bare potential, the effective8otentials being more “repulsive”: the equilibrium distance changes from 3.09 to 3.24 (3.28)˚A and the well depth, from 3.66 to 3.03 (2.84) meV, as one goes from the bare to the FH2(FH4) potentials, respectively.The corresponding potentials for Ne-coronene are shown in Fig.2 as functions of the y coordinate, while z and x are fixed at the absolute minimum of the classical potential (at3.21 and 0 ˚A, respectively). Again, it can be seen that the effect of the quadratic correctionis non negligible, for instance, the minimum energy moves from -27.83 (bare potential) to-25.68 meV (FH2 potential) but, on the other hand, the FH4 potential is very close to theFH2 one. III. CLASSICAL AND QUANTUM-MECHANICAL CALCULATIONSA. General procedure and notation In this work we compare the performance of the Ne N -coronene FH2 and FH4 effectivepotentials with that of the bare interaction potentials. First, we will study the equilibriumgeometries of these clusters by means of the BH approach. The corresponding calculationsare denoted by BH, BH-FH2 and BH-FH4 for the bare, FH2 and FH4 potentials, respectively.Note that, while the BH calculations are independent of the temperature, BH-FH2 and BH-FH4 must be repeated for the different temperatures of the study. The BH equilibriumgeometries are used as initial configurations of the CMC calculations at each temperature,which will be denoted as CMC, CMC-FH2 and CMC-FH4, for the bare, FH2 and FH4potentials, respectively. The resulting energies and configurations are compared with thePIMC calculations, where just the bare interaction potential is employed.For the BH calculations using the bare interaction potential, the zero point energy (ZPE)was also computed within the harmonic approximation and added to the BH equilibriumenergies (details can be found in Ref.[31]). This is a convenient estimation of the quantumeffects of the system when the thermal effects become small, as will be discussed in SectionIV. These calculations are denoted by BH+ZPE.9 . BH minimization Likely candidates for the global potential energy minima of Ne N -coronene clusters werelocated using the BH scheme [32], which is also known as the “Monte Carlo plus energyminimization” approach of Li and Scheraga [38]. This method transforms the potentialenergy surface into a collection of basins and explore them by hopping between local minima.This technique has been used successfully for both neutral [32, 39, 40] and charged atomicand molecular clusters [41–46], along with many other applications.[47] In the size rangeconsidered here the global optimization problem is feasible at a reasonable computationalcost. A total of 5 runs of 5 × BH steps each were performed for all clusters sizes.The global minimum was generally found in fewer than 10 BH steps. The optimizationtemperature was chosen between 8 and 10 K. C. PIMC and CMC calculations Details of the PIMC method employed here can be found in our study on He N -coroneneclusters [31] and in previous literature[11, 48–51]. The basic assumption involves express-ing the density matrix at a temperature T as a product of M density matrices at highertemperatures M T : ρ ( R , R M ; β ) = (cid:90) d R . . . d R M − M − (cid:89) α =0 ρ ( R α , R α +1 ; η ) , (12)where η = β/M . R α is the vector which collects the 3 N positions of the N Ne atoms: R α ≡ { r α , . . . , r αN } , being r αi the position vector of the i − th Ne atom at the time slice or imaginary time α . The total Hamiltonian ˆ H of the system with the coronene molecule fixedto the origin of coordinates can be written as:ˆ H = − ¯ h m Ne N (cid:88) i =1 ∇ i + V ( R ) . (13)where V is the interaction potential of Eq.1.10he internal energy is obtained by means of the virial estimator [52, 53] as: (cid:104) E ( T ) (cid:105) = 3 N β − (cid:42) M M − (cid:88) α =0 N (cid:88) i =1 ( r αi − r C i ) · F αi − M M − (cid:88) α =0 V ( R α ) (cid:43) . (14)where r C i = M − (cid:80) M − α =0 r αi is the centroid of the i th particle and F αi is the force experiencedby the i particle on the α slice. The integration is carried out via a Metropolis Monte Carloalgorithm, as an average over a number of paths {R , R , . . . , R M , R M +1 } sampled accordingto a probability density proportional to the factorized product of M density matrices ofEq.12. Exchange effects are neglected. The number of beads vary between M = 1 (forthe CMC, CMC-FH2 and CMC-FH4 calculations) to a maximum of 150 for the lowesttemperature PIMC simulations. Depending on M , the number of steps varies between 10 to 10 . The staging sampling method has been employed[54] involving a number of eightbeads in each movement for the PIMC simulations. The final average energy is obtained byextrapolation to the M → ∞ following a parabolic law[55, 56]. IV. RESULTS AND DISCUSSIONA. Cluster energies and structures at 6 K Before tackling the study of the FH approximation as a function of the temperature, westart presenting results at 6 K, an intermediate value in the range addressed here and incoincidence with the temperature of various experiments on PAHs isolated in Ne matrices[7,29, 30].In Table II, energies of various Ne N -coronene clusters, N =1-4 and 14, at 6 K as obtainedfrom the BH and CMC approaches and using the bare, FH2 and FH4 effective potentialsare reported and compared with the PIMC results. Various arrangements (or “isomers”)( n a , n b ) are examined for each number of Ne atoms, where n a and n b refer to the numberof atoms placed above and below the coronene plane, respectively. It can be seen that boththe BH and the CMC energies tend to the PIMC energies as one goes from using the bare to11he FH2 and FH4 potentials. However, given the difference between the BH and the CMCenergies, it is clear that thermal effects are non negligible at this temperature so, among allthe calculations, the CMC-FH4 energies are the ones giving the closest agreement with thereference PIMC results. In addition, notice that CMC-FH4 agrees with PIMC as to whichis the most stable isomer ( n a , n b ) for a given cluster size N . For example, for N = 4, themost stable arrangement within the FH4 potential is the (3,1) one, in agreement with thePIMC calculation, whereas (4,0) gives the absolute minimum when using the bare and FH2potentials. This result can be explained by a more repulsive character of Ne-Ne interactionwhen it is considered at the FH4 level (Fig.1), thus making an arrangement with a largerdensity of Ne atoms relatively less stable (as the (4,0) one).Also from Table II, it is worth mentioning that the addition of the ZPE to the BH energiesgives a fair agreement with the PIMC energies, and that the BH+ZPE calculations correctlypredict the relative stability of the different isomers for a given N .Computation of the energies per atom helps us to quantify the performance of the FHapproach. Results are depicted in Fig.3, where it is clear the improvement of the CMC-FH energies with respect to the CMC ones. In more detail, note that, for a given N ,the differences between the CMC and PIMC energies are larger for those clusters havinga larger number of atoms on a given side of the molecule. In this case, quantum effectsincrease because the number of effective Ne-Ne interactions increases as well (interactionsbetween atoms sitting on different sides of the molecule are negligible). For example, for(3,0), with three effective Ne-Ne interactions, the relative error of the CMC calculation is 14%, whereas for the (2,1) cluster, with just one Ne-Ne interaction, the error reduces to 11 %.The FH2 and FH4 potentials certainly amend the classical result, although the errors are alsolarger when the number of interacting atoms increase. Indeed, the errors of the CMC-FH2calculation are 5 % and 2 % for the (3,0) and (2,1) clusters, respectively. It is importantto mention that both Ne-Ne and Ne-coronene interactions contribute significantly to thequantum corrections: for the (3,0) cluster, 24 % and 76 % of the FH2 correction (computedat the corresponding optimal geometry) are due to Ne-Ne and Ne-coronene, respectively,12hereas these numbers become 12 % and 88% in the case of the (2,1) arrangement.Besides increasing the overall energies, there are some changes in the geometries of theNe N -coronene clusters when the interaction potentials are modified. Probability densities asfunctions of the coordinates parallel to the coronene plane, D ( x, y ), as obtained from PIMC,CMC-FH4 and CMC, are shown in Fig.4 for (7,7) Ne -coronene at 6 K (results for the FH2potential are not shown as they are quite close to the FH4 ones). These plots have beenobtained by means of a histogramming procedure on the x and y coordinates, accumulatingthe probability density along the z coordinate. For each side of the molecule, one Ne atomis placed above the central hollow, while the other six atoms are located near the bordersof the outer hexagons. The distributions from the different calculations are quite similar.However, it can be noticed that the peaks of the PIMC and CMC-FH4 probability densitiesare somewhat wider and that the average distance between those peaks is slightly largerthan in the case of the CMC calculation. These features can be further examined by meansof the one-dimensional distribution D (cid:0) [ x + y ] / (cid:1) as shown in Fig.5 (it is computed byaccumulating the probability density D ( x, y ) along the angle ϕ given by tan ϕ = y/x ). TheCMC-FH4 distribution is in better agreement with the PIMC calculations than the purelyclassical one. Fig.5 also depicts the distance of the outer Ne atoms with respect to thecoronene symmetry axis ( z ) as obtained from the BH and BH-FH4 optimized geometries,and it can be seen that these positions coincide with the maxima of the correspondingCMC distributions. The FH4 potential involves larger equilibrium distances than the barepotential (for instance, see Fig. 1), a feature that can explain the shift in the peak positionof the CMC-FH4 with respect to the CMC one. The peak of the PIMC distribution is inthe middle, suggesting that the FH4 potential is overestimating this effect. B. Temperature dependence of the cluster energies Further insight is gained into the FH approach by studying the cluster energies as func-tions of the temperature for the different cluster sizes. In Fig.6, the behavior of the (1-4,13) and (7,7) cluster energies as obtained from the different potentials and methods is testedagainst PIMC in the temperature range 2-14 K. To make a more coherent comparison, theenergies have been shifted by the minimum energy (BH) and divided by the number ofatoms, (cid:101) E ( T ) = E ( T ) − E BH N . (15)Although not all of the ( n a =1-4,0) arrangements correspond to the absolute minimum energy(see Table II), we have chosen this sequence in order to study the FH approach as a functionof the number of Ne atoms over a given face of the molecule.In the higher temperature range of Fig. 6 ( T > T = 10 K. Note that “un-shifted” energies were taken for the calculation of the relative errors (i.e., E instead of (cid:101) E ).Note that the CMC errors rise with the number of Ne atoms, as already discussed above.It can also be seen from the Table that the CMC-FH2 and CMC-FH4 calculations roughlyhalve the CMC errors independently of the number of Ne atoms. Therefore, although theFH2 and FH4 potentials do not provide a perfect agreement with the PIMC energies, theydo introduce quantum corrections in a steady way.As temperature decreases ( T < ∼ < ∼ k B ,in agreement with a model of N classical harmonic oscillators in a three-dimensional (3D)space. In addition, the PIMC energies at low temperatures are in a fairly good agreementwith the ZPE computed here within the harmonic approximation (dotted lines in Fig. 6).In this way, we have tested a model of N 3D isotropic harmonic oscillators against presentcalculations. It is assumed that each of the N Ne atoms moves under an effective potential V har ( r ) = m Ne ω N r / 2, where ω N is a characteristic frequency which varies with N . Itsvalue is obtained from equating the computed ZPE per atom to 3¯ hω N / 2. In this way, the“average interaction” undergone by each atom -which depends on the interaction with boththe molecular substrate and the remaining atoms in the cluster- is approximated by anisotropic harmonic potential. Within this model, the quantum mechanical average energyper atom is [57] (cid:101) E qhar ( T ) = 32 ¯ hω N coth ¯ hω N k B T , (16)to be compared with the PIMC energy. Applying Eq. 7 to the harmonic potential V har , theFH2 potential is (cid:101) E BH − FH2har ( T ) = 18 (¯ hω N ) k B T , (17)which will be associated with the BH-FH2 energy. Note that the FH4 potential is identicalto the FH2 one within the present harmonic approximation. The FH2 thermal energy per15tom is obtained by adding 3 k B T , (cid:101) E CMC − FH2har ( T ) = (cid:101) E BH − FH2har ( T ) + 3 k B T, (18)which will be related to the CMC-FH2 energies. Results are shown in Fig. 7 for the (2,0)and (4,0) clusters. Although anisotropy and anharmonicity effects should probably be addedto attain a more quantitative agreement, the model energies compare fairly well with thecomputed ones and, in this way, this simple model does provide an adequate zero-orderdescription of the behavior of these clusters. C. Applicability of the FH2 and FH4 effective potentials The results of Fig. 6 indicate that both FH2 and FH4 potentials are useful for T > He Lennard-Jones systems, Ses´e[24] found that the FH2 potential generally performsbetter than the FH4 one. In this work, we have found that the FH4 potential is useful justin a narrow range of temperatures (around 6 K), since at larger temperatures its correctionsbecome negligible whereas at lower temperatures it worsens the FH2 estimations.It is thus important to determine the temperature range of validity of these effectivepotentials as they can lead to erroneous results when applied below a critical temperature T (cid:63) . Ses´e has found that, for relatively low densities, this temperature can be deduced fromthe condition λ B /σ < . 5, where λ B = ¯ h (cid:112) π/mk B T (cid:63) is the thermal de Broglie wavelengthand σ is the Lennard-Jones collision diameter. Here, taking σ = 2.76 ˚A (the diameter of theNe-Ne ILJ potential), we obtain T (cid:63) ≈ T (cid:63) that might be extrapolated to relatedsystems where PIMC calculations would become too time-consuming. This temperatureis defined as the crossing between the approximate ( (cid:101) E CMC − FH2har ) and the quantum ( (cid:101) E qhar )energies. Given that at low temperatures (cid:101) E qhar ( T ) ≈ ¯ hω N , the ZPE per atom of thesystem, a simple quadratic equation follows, whose lowest root is T (cid:63) = 32 ¯ hω N (cid:32) − (cid:112) / k B (cid:33) . (19)The resulting temperatures range from T (cid:63) = 2.8 to 4.8 K as we go from the (1,0) to the(7,7) cluster, respectively. It is worth noting that this model correctly takes into account theshift of T (cid:63) with the number of rare gas atoms (density), as found for other systems[23, 25].Indeed, T (cid:63) is proportional to the ZPE per atom, which increases with the number of Neatoms, as can be seen in Fig. 6. This trend is possibly due to a rise of the frequencies of theaverage interaction undergone by each Ne atom as the number of rare gas atoms surroundingit increases. V. CONCLUDING REMARKS The Feynman-Hibbs (FH) approach has been applied to the study of Ne N -coronene clus-ters ( N = 1-4, 14) at low temperatures ( T = 2-14 K) and using realistic analytical potentials.The suitability of the quadratic (FH2) and quartic (FH4) effective potentials has been in-vestigated by comparing basin-hoping (BH) optimizations and classical Monte Carlo (CMC)calculations of cluster energies and structures with benchmark path-integral Monte Carlo(PIMC) calculations.For T > T < ∼ VI. APPENDIX: DERIVATION OF THE FEYNMAN-HIBBS EFFECTIVE PO-TENTIALS In this paragraph we give a more detailed account of the derivation of the quadratic (FH2)and quartic (FH4) effective potentials for the interaction potentials of this study (Eqs. 8-11).Since these potentials are of the type V ( ρ ) or V ( ρ, cos θ ), it is convenient to consider theoperators of Eq.7 in spherical coordinates. In particular, the Laplacian is given, in terms of r ≡ ( r, cos θ, φ ), as ∇ = ∂ ∂r + 2 r ∂∂r + 1 r (cid:18) (1 − c ) ∂ ∂c − c ∂∂c + 1(1 − c ) ∂ ∂φ (cid:19) , (20)18here c = cos θ .To apply Eq.7, we use the property of invariance of the Laplacian operator under anytranslation or rotation in the three-dimensional space. In the case of the Ne-Ne interaction,the origin of the coordinate system is displaced to coincide with one of the atoms, and, since V Ne − Ne only depends on the radial distance ρ , ∇ V Ne − Ne = d V Ne − Ne dρ + 2 ρ dV Ne − Ne dρ , (21)leading to the FH2 correction of Eq.8.For the Ne-coronene interaction, taking into account Eq. 4 ∇ V Ne − Cor ( r ) = (cid:88) k ∇ U k ( ρ k , c k ) . (22)For every bond k , ∇ U k can be easily performed after a transformation from the originalCartesian system to a new one where the origin is at the bond center and the z axis isaligned with the bond axis. In this way, ρ k and c k coincide with the radial distance and thecosine of the polar angle of the new reference system and, taking into account Eq.20 andthe independence of the potential with respect to the azimuthal angle, ∇ U k = ∂ U k ∂ρ k + 2 ρ k ∂U k ∂ρ k + 1 ρ k (cid:18) (1 − c k ) ∂ U k ∂c k − c k ∂U k ∂c k (cid:19) , (23)which can be readily related to Eq. 10.To obtain the FH4 potential of the Ne-Ne interaction, the ∇ operator from Eq. 7 isapplied to V Ne − Ne as ∇ V Ne − Ne = (cid:18) d dρ + 2 ρ ddρ (cid:19) (cid:18) d V Ne − Ne dρ + 2 ρ dV Ne − Ne dρ (cid:19) = d V Ne − Ne dρ + 4 ρ d V Ne − Ne dρ , (24)19n expression which can be immediately identified with the FH4 corrections of Eq.9.At this point, it is worthwhile noting that other authors[18, 19, 21, 25] have reported anextra term, especifically ρ dVdρ , in the formal expressions of the FH4 potentials (Eq. 2 of Ref.[18] for example). We have checked that the contributions depending on other derivativesof the potential out of the third and fourth derivatives are canceled out in the calculationof ∇ V , and therefore, that the correct FH4 potential expression is that given by Eq. 9 ofthis work.Finally, the ∇ operator applied to the atom-bond potential U k ( ρ k , c k ) is ∇ U k = ∂ U k ∂ρ k + 4 ρ k ∂ V k ∂ρ k + 2 ρ k (cid:20) (1 − c k ) ∂ U k ∂ρ k ∂c k − c k ∂ U k ∂ρ k ∂c k (cid:21) + 1 ρ k (cid:20) (1 − c k ) ∂ U k ∂c k − c k (1 − c k ) ∂ U k ∂c k − − c k ) ∂ U k ∂c k (cid:21) . (25)Eq. 11 is readily obtained once all the derivatives with respect c k are neglected, as alreadydiscussed elsewhere. Acknowledgments The work has been funded by Spanish MINECO grants FIS2013-48275-C2-1-P, FIS2014-51993-P, and FIS2013-41532-P. R. R.-C. thanks support from FIS2013-48275-C2-1-P. Al-location of computing time by CESGA (Spain) and support by the COST-CMTS ActionCM1405 “Molecules in Motion (MOLIM)” are also acknowledged. [1] F. Paesani and G. A. Voth, J. Phys. Chem. B , 5702 (2009).[2] D. P. Broom and K. M. Thomas, MRS Bulletin , 412 (2013).[3] J. Cai, Y. Xing, and X. Zhao, RSC Advances , 8579 (2012).[4] A. Mart´ınez-Mesa, L. Zhechkov, S. N. Yurchenko, T. Heine, G. Seifert, and J. Rubayo-Soneira,J. Phys. Chem. C , 19543 (2012).[5] T. X. Nguyen, H. Jobic, and S. K. Bhatia, Phys. Rev. Lett. , 085901 (2010). 6] M. Bahou, Y. J. Wu, and Y. P. Lee, Angew. Chem. Int. Ed. , 1021 (2015).[7] I. Garkusha, A. Nagy, J. Fulara, M. F. Rode, A. L. Sobolewski, and J. P. Maier, J. Phys.Chem. A , 351 (2013).[8] S. Yang and A. M. Ellis, Chem. Soc. Rev. , 472 (2013).[9] S. Szalewicz, Int. Rev. Phys. Chem. , 273 (2008).[10] R. P. Feynman, Statistical Mechanics (Benjamin, New York, 1972).[11] D. M. Ceperley, Rev. Mod. Phys. , 279 (1995).[12] C. Chakravarty, Int. Rev. Phys. , 421 (1997).[13] R. P. Feynman and A. Hibbs, Quantum Mechanics and Path-Integrals (McGraw-Hill, NewYork, 1965).[14] B. Guillot and Y. Guissani, J. Chem. Phys. , 10162 (1998).[15] F. Calvo, J. P. K. Doye, and D. J. Wales, J. Chem. Phys. , 7312 (2001).[16] N. Tchouar, F. Ould-Kaddour, and D. Levesque, J. Chem. Phys. , 7326 (2004).[17] M. Abbaspour and Z. Borzouie, Fluid Phase Equilibria , 167 (2014).[18] A. V. A. Kumar and S. K. Bhatia, Phys. Rev. Lett. , 245901 (2005).[19] A. V. A. Kumar, H. Jobic, and S. K. Bhatia, J. Phys. Chem. B , 16666 (2006).[20] D. Noguchi, H. Tanaka, A. Kondo, H. Kajiro, H. Noguchi, T. Ohba, H. Kanoh, and K. Kaneko,J. Am. Chem. Soc. , 6367 (2008).[21] D. Liu, W. Wang, J. Mi, C. Zhong, Q. Yang, and D. Wu, Ind. Eng. Chem. Res. , 434 (2012).[22] C. I. Contescu, H. Zhang, R. J. Olsen, E. Mamontov, J. R. Morris, and N. C. Gallego, Phys.Rev. Lett. , 236102 (2013).[23] L. M. Ses´e, Mol. Phys. , 1297 (1994).[24] L. M. Ses´e, Mol. Phys. , 931 (1995).[25] P. Kowalczyk, L. Brualla, P. A. Gauden, and A. P. Terzyk, Phys. Chem. Chem. Phys. ,9182 (2009).[26] Y. Jiao, A. Du, M. Hankel, and S. C. Smith, Phys. Chem. Chem. Phys. , 4832 (2013).[27] M. Bartolomei, E. Carmona-Novillo, and G. Giorgi, Carbon , 1076 (2015).[28] C. Joblin, L. d’Hendecourt, A. L´eger, and D. D´efourneau, Astron. Astrophys. , 923 (1994).[29] M. Steglich, C. J¨ager, G. Rouill´e, F. Huisken, H. Mutschke, and T. Henning, ApJL , L16(2010).[30] I. Garkusha, J. Fulara, P. J. Sarre, and J. P. Maier, J. Phys. Chem. A , 10972 (2011). 31] R. Rodr´ıguez-Cantano, R. P. de Tudela, M. Bartolomei, M. I. Hern´andez, J. Campos-Mart´ınez,T. Gonz´alez-Lezana, P. Villarreal, J. Hern´andez-Rojas, and J. Bret´on, J. Chem. Phys. ,224306 (2015).[32] D. J. Wales and J. P. K. Doye, J. Phys. Chem. A , 5111 (1997).[33] F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti, and F. Vecchiocattivi, Phys.Chem. Chem. Phys , 5489 (2008).[34] F. Pirani, M. Albert´ı, A. Castro, M. M. Teixidor, and D. Cappelletti, Chem. Phys. Lett. ,37 (2004).[35] M. Bartolomei, E. Carmona-Novillo, M. I. Hern´andez, J. Campos-Mart´ınez, and F. Pirani, J.Phys. Chem. C , 10512 (2013).[36] A. Griewank and A. Walter, Evaluating Derivatives: Principles and Techniques of AlgorithmicDifferenciation (Society of Industrial and Applied Mathematics, 2008).[37] K. T. Tang and J. P. Toennies, J. Chem. Phys. , 4976 (2003).[38] Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. U.S.A , 6611 (1987).[39] J. Hern´andez-Rojas, F. Calvo, J. Bret´on, and J. Gomez Llorente, J. Phys. Chem. C ,17019 (2012).[40] S. Acosta-Guti´errez, J. Bret´on, J. M. G. Llorente, and J. Hern´andez-Rojas, J. Chem. Phys. , 074306 (2012).[41] J. Hern´andez-Rojas, J. Bret´on, J. M. Gomez Llorente, and D. J. Wales, J. Phys. Chem. B , 13357 (2006).[42] J. P. K. Doye and D. J. Wales, Phys. Rev. B , 2292 (1999).[43] J. Hern´andez-Rojas and D. J. Wales, J. Chem. Phys. , 7800 (2003).[44] J. Hern´andez-Rojas, J. Bret´on, J. M. Gomez Llorente, and D. J. Wales, J. Chem. Phys. ,12315 (2004).[45] J. Hern´andez-Rojas, J. Bret´on, J. G. Llorente, and D. Wales, Chem. Phys. Lett. , 404(2005).[46] J. Hern´andez-Rojas, F. Calvo, F. Rabilloud, J. Bret´on, and J. M. Gomez Llorente, J. Phys.Chem. A , 7267 (2010).[47] D. J. Wales, Energy Landscapes (Cambridge University Press, Cambridge, 2003).[48] Y. Kwon, D. M. Ceperley, and K. B. Whaley, J. Chem. Phys. , 2341 (1996).[49] R. Rodr´ıguez-Cantano, D. L´opez-Dur´an, R. P´erez de Tudela, T. Gonz´alez-Lezana, G. Delgado- arrio, P. Villarreal, and F. A. Gianturco, Comp. Theor. Chem. , 106 (2012).[50] R. Rodr´ıguez-Cantano, R. P´erez de Tudela, D. L´opez-Dur´an, T. Gonz´alez-Lezana, F. A.Gianturco, G. Delgado-Barrio, and P. Villarreal, Eur. Phys. J. D , 119 (2013).[51] R. Rodr´ıguez-Cantano, T. Gonz´alez-Lezana, P. Villarreal, and F. A. Gianturco, J. Chem.Phys. , 104303 (2015).[52] M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. Phys. , 5150 (1982).[53] K. R. Glaesemann and L. E. Fried, J. Chem. Phys. , 5951 (2002).[54] M. Sprik, M. Klein, and D. Chandler, Phys. Rev. B , 4234 (1985).[55] E. Cuervo and P.-N. Roy, J. Chem. Phys. , 124314 (2006).[56] R. P´erez de Tudela, D. L´opez-Dur´an, T. Gonz´alez-Lezana, G. Delgado-Barrio, P. Villarreal,F. A. Gianturco, and E. Yurtsever, J. Phys. Chem. A , 6892 (2011).[57] M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford UniversityPress, 2010).[58] M. I. Hern´andez, M. Bartolomei, and J. Campos-Mart´ınez, J. Phys. Chem. A , 10743(2015). igures E ne r g y ( m e V ) Ne−Ne distance (Å) Ne−Ne potential Bare potentialFH2 (6K)FH4 (6K) FIG. 1: Bare Ne-Ne interaction potential (meV) as a function of the interatomic distance (in ˚A),compared with the quadratic and quartic Feynman-Hibbs effective potentials at 6 K. E ne r g y ( m e V ) y (Å) yBare potentialFH2 (6K)FH4 (6K) FIG. 2: Ne-coronene (Bare, FH2 and FH4 at 6K) interaction potentials (in meV), as functions ofthe y Cartesian coordinate (depicted in the inset), where the x and z coordinates are fixed at theabsolute minimum of the bare potential (0. and 3.21 ˚A, respectively). 28 -26 -24 Energy per atom (meV) N u m be r o f N e a t o m s , N CMCCMC-FH2CMC-FH4PIMC (3,1)(4,0)(2,2)(3,0)(2,1)(1,1)(2,0)(n a ,n b ) {{{ FIG. 3: Energies per atom (in abscissas) of Ne N -coronene at 6 K as obtained from CMC, CMC-FH2, CMC-FH4 and PIMC, for the different sizes N and isomers ( n a , n b ) (in ordinates). IG. 4: Probability densities D ( x, y ) of (7,7) Ne -coronene at 6 K. Left, middle and right panels:PIMC, CMC-FH4 and CMC calculations, respectively. +y ) [Å]T = 6 K BH BH−FH4PIMCCMCCMC−FH4 FIG. 5: Normalized probability density D (cid:0) [ x + y ] / (cid:1) of (7,7) Ne -coronene at 6 K, for CMC,CMC-FH4 and PIMC calculations. In addition, the distance of the outer Ne atoms to the coronenesymmetry axis ( z ), as obtained from the BH and BH-FH4 calculations, are depicted by arrows. E ne r g y ( m e V ) BH-FH2BH-FH4CMCCMC-FH2CMC-FH4PIMCzpe E ne r g y ( m e V ) FIG. 6: Energies per atom, shifted as in Eq.15, of Ne N -coronene clusters as functions of temperaturefor the different methods used in this work. Insets depict the classical optimal geometries of theclusters. (1,0), (2,0), (3,0), (4,0) and (7,7) are shown in the left upper, middle upper, right upper,left lower and middle lower panels, respectively. Calculations have been performed at 2, 4, 6, 10,and 14 K (lines are guides to the eye). See text for discussion. E ne r g y ( m e V ) BH-FH2CMCCMC-FH2PIMC E ne r g y ( m e V ) FIG. 7: Comparison between the computed BH-FH2, CMC, CMC-FH2 and PIMC energies (sym-bols) and the corresponding energies obtained with a model of N harmonic oscillators (solid lines).Upper and lower panels: (2,0) and (4,0) clusters, respectively. See text for details. ables ABLE I: Parameters of the Ne-Ne, Ne-CC and Ne-CH pair potentials[33, 35]. Distances are in˚A, energies in meV, and γ is dimensionless.Pair ρ e ε γ Ne-Ne 3.094 3.660 9.0 ρ ⊥ e ρ (cid:107) e ε ⊥ ε (cid:107) γ Ne-CC 1.297 1.809 3.643 3.995 8.5Ne-CH 2.544 1.782 3.316 3.655 9.0 ABLE II: Ne N -coronene energies (in meV) at 6 K for the different calculations of this work(see text for details). Various isomers ( n a , n b ) are studied, where n a ( n b ) is the number of Neatoms above (below) the coronene plane. PIMC error bars (in meV), associated to the M → ∞ extrapolation procedure, are given in parenthesis. Standard deviation of the CMC energies (notshown) is about 0.01 meV. For a given N , the absolute minimum energy within each method isshown in boldface. N ( n a , n b ) BH BH+ZPE BH-FH2 BH-FH4 CMC CMC-FH2 CMC-FH4 PIMC1 (1,0) -27.83 -24.35 -26.04 -25.80 -25.99 -24.19 -23.95 -23.89 (0.01)2 (1,1) -55.72 -48.75 -52.12 -51.65 -52.04 -48.44 -47.96 -47.85 (0.10)(2,0) -55.36 -47.26 -51.13 -50.44 -52.20 -48.10 -47.42 -46.72 (0.03)3 (3,0) -87.42 -73.85 -80.51 -79.23 -82.44 -75.49 -74.20 -72.10 (0.01)(2,1) -83.28 -71.69 -77.25 -76.32 -78.28 -72.36 -71.44 -70.62 (0.01)4 (4,0) -117.29 -97.03 -106.68 -104.58 -110.50 -99.97 -97.88 -94.81 (0.10)(3,1) -115.38 -98.32 -106.66 -105.14 -108.60 -99.80 -98.28 -96.04 (0.20)(2,2) -110.90 -94.66 -102.41 -101.01 -104.54 -96.34 -94.96 -93.50 (0.20)14 (7,7) -415.15 -332.53 -369.99 -360.34 -392.46 -346.83 -336.94 -322.44 (0.02) ABLE III: Relative errors (in %) of the CMC energies with respect to the PIMC ones ( | E CMC − E PIMC | /E PIMC ) as obtained from the different potentials and clusters ( n a , n b ). T (K) Potential (1,0) (2,0) (3,0) (4,0) (7,7)Bare 7 10 12 14 2010 FH2 2 4 6 7 10FH4 2 4 5 6 8(K) Potential (1,0) (2,0) (3,0) (4,0) (7,7)Bare 7 10 12 14 2010 FH2 2 4 6 7 10FH4 2 4 5 6 8