An ab initio path integral Monte Carlo simulation method for molecules and clusters: application to Li_4 and Li_5^+

A novel method for simulating the statistical mechanics of molecular systems in which both nuclear and electronic degrees of freedom are treated quantum mechanically is presented. The scheme combines a path integral description of the nuclear variables with a first-principles adiabatic description of the electronic structure. The electronic problem is solved for the ground state within a density functional approach, with the electronic orbitals expanded in a localized (Gaussian) basis set. The discretized path integral is computed by a Metropolis Monte Carlo sampling technique on the normal modes of the isomorphic ring-polymer. An effective short-time action correct to order $\tau^4$ is used. The validity and performance of the method are tested in two small Lithium clusters, namely Li$_4$ and Li$_5^+$. Structural and electronic properties computed within this fully quantum-mechanical scheme are presented and compared to those obtained within the classical nuclei approximation. Quantum delocalization effects are significant but tunneling turns out to be irrelevant at low temperatures.

A novel method for simulating the statistical mechanics of molecular systems in which both nuclear and electronic degrees of freedom are treated quantum mechanically is presented. The scheme combines a path integral description of the nuclear variables with a first-principles adiabatic description of the electronic structure. The electronic problem is solved for the ground state within a density functional approach, with the electronic orbitals expanded in a localized ͑Gaussian͒ basis set. The discretized path integral is computed by a METROPOLIS Monte Carlo sampling technique on the normal modes of the isomorphic ring polymer. An effective short-time action correct to order 4 is used. The validity and performance of the method are tested by studying two small lithium clusters, namely Li 4 and Li 5 ϩ . Structural and electronic properties computed within this fully quantum-mechanical scheme are presented and compared to those obtained within the classical nuclei approximation. Quantum delocalization effects turn out to be significant as shown by the fact that quantum simulation results at 50 K approximately correspond to those of classical simulations carried out at 150 K. The scaling factor depends, however, on the specific physical property, thus evidencing the different character of quantum and thermal correlations. Tunneling turns out to be irrelevant in the temperature range investigated ͑50-200 K͒. © 1998 American Institute of Physics. ͓S0021-9606͑98͒51421-3͔

I. INTRODUCTION
Light atoms, such as H, He, Li, or Be, cannot very often be treated as classical particles, particularly at low temperatures. As temperature decreases and the thermal de Broglie wavelength increases, the quantum character emerges, and a description in terms of classical coordinates and momenta breaks down. The most obvious manifestation of the quantum character of light atoms is a large zero-point energy ͑ZPE͒. A particle of mass m in a harmonic potential with characteristic frequency will have a ZPE of ប/2 and an associated spatial delocalization of ⌬xϭͱប/m. For instance, a proton in a typical bonding environment such as a H-O bond or a H 2 molecule, will have a ZPE of 0.15-0.25 eV and ⌬x between 0.2 and 0.3 Å. This represents a sizable effect which could be decisive in stabilizing a particular crystalline structure for a solid, or the ground state configuration of a molecule or a cluster. An even more interesting manifestation of quantum effects is the possibility that these light nuclei can tunnel across potential energy barriers, thus exploring classically forbidden regions of configuration space and giving rise to a variety of interesting quantum effects such as temperature-independent diffusion, exotic ground states, resonances in ion-surface scattering, and fluxional molecules. Signatures of quantum effects can also be observed in low-energy atomic collisions, or in protontransfer reactions in the gas and condensed phases.
To date, most studies that consider the quantum character of atomic nuclei are based on an empirical description of the interatomic interactions, or otherwise consist of extending and/or correcting a posteriori the results obtained within a classical nuclei approximation. Classical potentials are frequently not transferable from one environment to another, and are ill-suited to modeling phenomena involving significant electronic density redistribution, as in the making and breaking of chemical bonds. The natural route to overcome these limitations is to describe the interactions at a firstprinciples level, i.e., by including explicitly the electronic component in the description of the system. The recent de-velopment of such schemes, which address the question of the interplay between electronic structure and the quantum nature of light nuclei, has only very recently began to be realized, thus opening a fascinating field with important implications in many branches of physics, chemistry, and biology.
Since small clusters usually exhibit rich landscapes of isomeric forms within narrow energy bands, 1 they constitute good systems for studying the effects of quantum delocalization and associated tunneling behavior. Lithium clusters are particularly interesting because in addition to the small atomic mass, they are bound by metallic many-body interactions which cannot be adequately represented by means of classical interatomic potentials. 2 In the remainder of Sec. I we outline the methodology that we have developed to study this class of problems, and review the present understanding of Li clusters which is the test system for our method. In Section II we introduce the theoretical framework of our ab initio path integral approach and discuss the approximations involved. Section III is devoted to the details of the practical implementation of the path integral Monte Carlo ͑PIMC͒ and the electronic structure methods. In Sec. IV we present the validation of the electronic structure calculations, zero-temperature geometries, and electronic properties of Li 4 and Li 5 ϩ clusters. The results of our simulations for the classical and quantum Li 4 and Li 5 ϩ clusters at finite temperatures are presented in Sec. V. Section VI contains our conclusions and an assessment of the potential of this novel simulation tool.

A. Methodological aspects
The goal of the present work is to introduce a novel computational technique for studying the statistical mechanics of isolated systems like clusters and molecules containing light atoms. Our approach combines an imaginary-time path integral description of the nuclear degrees of freedom 3 with a first-principles density functional ͑DFT͒ description of the electronic structure. 4 Since the natural choice for investigating isolated systems is to use a localized basis set for the electronic orbitals, we adopt a Gaussian basis set. 5 In the present implementation the electronic structure is computed at the all-electron level, i.e., explicitly including core electrons. The sampling of the path integral is implemented using Monte Carlo ͑MC͒ techniques. The electronic energy is minimized for each nuclear configuration, and the MC rejection step is also performed using the energy calculated at the same level of sophistication.
Other schemes along these lines have been recently proposed by Marx and Parrinello 6 and Cheng et al. 7 At variance with our approach, these two methods use molecular dynamics ͑MD͒ for sampling the path integral, a plane-wave ͑PW͒ expansion for the electronic orbitals, and treat the electronion interaction at the pseudopotential level. PW expansions with periodic boundary conditions are more appropriate for extended systems such as solids or liquids though they can be adequately modified to deal with isolated systems. 6 The evaluation of real-time path integrals, which would include the full dynamical information, is an extremely difficult numerical task because the integrand is a rapid oscilla-tory function of the path. In the imaginary-time framework the statistical weight is real and positive definite, so that the integrals are well-conditioned, but the price is that the dynamics is not directly accessible. In the absence of real-time dynamical information, it is not particularly advantageous to sample the integral using MD in place of MC techniques. In particular, the MD technique requires elaborate thermostatting mechanisms 8 in order to overcome ergodicity problems in the sampling of the quasiharmonic degrees of freedom that appear in the path integral formulation ͑see Sec. II͒. In the present method we propose a METROPOLIS MC sampling technique on the normal modes of the polymer, which has less severe ergodicity problems, and is more convenient from the point of view of efficient evaluation of the path action ͑see Sec. III͒. 9 Moreover, the MC strategy is easy to interface with any electronic structure code.

B. Small lithium clusters
Structural and electronic properties of Li n and Li n ϩ clusters have been systematically investigated for nϭ2 -9 by means of ab initio configuration-interaction calculations. 10 A key observation was the existence of several isomers of comparable energy. Recently, ab initio classical MD simulations at the Hartree-Fock ͑HF͒ and nonlocal density functional ͑NLDF͒ level were carried out in order to explore the cluster dynamics as a function of temperature, to analyze isomerization reactions, and to study the melting transition. [11][12][13][14] The extent to which different levels of first-principles calculations are reliable for describing small Li clusters has been very recently discussed by Rousseau and Marx, 15 who concluded that either ab initio MP2 calculations or gradient corrected NLDF provide a reasonable potential energy landscape, while HF and LDF are inadequate.
This aspect is very important because the energy landscape determines the probability with which the cluster visits different possible geometries. In general, the ability to jump from one minimum to another will depend on the extent of both thermal and quantum fluctuations. Figure 1 shows three typical situations: At low temperatures quantum delocalization is the dominant mechanism. If the ZPE is higher than the barrier between the two minima ͑a͒, or smaller than the barrier but larger than the energy difference between them ͑b͒ ͑tunneling regime͒, the ground state wave function will sample both configurational minima. If, instead, it is smaller than the energy difference between the two minima ͑c͒, then the ground state will be basically unchanged from the classical one, though all structural distributions will be broadened by quantum delocalization. Thermal excitations can promote a classical-like situation like ͑c͒ to a thermally assisted tunneling regime.
The above picture is valid as long as the electronic ground state is nondegenerate, and this is realized for clusters with closed electronic shells. Open-shell clusters undergo Jahn-Teller distortions, and are characterized by degenerate ground states and pseudorotations. An adequate treatment of this situation would require the introduction of the concepts of conical intersections and geometric phases. 16 Though these phenomena are of great interest, and have indeed been detected in Li 3 ͑Ref. 17͒ and Li 5 , 18 we have preferred to avoid this extra complication at this stage in the development of our method.
Quantum nuclear effects in Li clusters were first addressed by Ballone and Milani, 19 who studied magic-number clusters ͑20, 40, and 92 atoms͒ by means of PIMC simulations using a simple, jellium model description of the electronic component. Their most interesting observation was the existence of a large number of isomers differing only in the location of the outermost atoms, such that quantum tunneling among these isomers led to a fluid like behavior. A very recent ab initio path integral MD study of Li 8 and Li 20 clusters by Rousseau and Marx 20 shows that this picture does not hold when the description of the electronic component is improved, thus demonstrating the necessity of going beyond the jellium-type description.
For our study we have chosen the following two examples: Li 4 has a single isomer and is well described by a quasiclassical picture. Li 5 ϩ has two isomers at an energy difference of about 150 K, and two other isomers at higher energies ͑see Sec. V͒. The zero-point energy is of the order of 100 K, so that Li 5 ϩ could constitute a good candidate for exhibiting significant quantum effects in terms of the sampling of configurational minima. The results presented in Sec. V will show that, in spite of this, Li 5 ϩ actually behaves in a very similar way to Li 4 .

II. THE AB INITIO PATH INTEGRAL PARTITION FUNCTION
The statistical mechanics of quantum many-body systems can be formulated in terms of the two-point density matrix, or imaginary-time propagator: ͑R,RЈ,␤ ͒ϭ͗R͉exp͑ Ϫ␤H͉͒RЈ͘, ͑1͒ whose trace is the partition function Z. ␤ϭ1/k B T is the inverse temperature. The path integral representation of the density matrix is given by where R(u) represents the configuration of an N-body system as a function of imaginary time u. The range of u is from 0 to ␤ប and the paths considered are restricted to those beginning at R(0)ϭR and ending at R(␤ប)ϭRЈ. The partition function can be similarly expressed as a path integral with contributions from all possible cyclic paths for which R(0)ϭR(␤ប). D͓R(u)͔ represents the differential element for all paths. The Euclidean action, S͓R(u)͔, associated with a path is defined as The first term corresponds to the kinetic energy contribution to the action, with m the mass of the particles. The generalization to heterogeneous systems, i.e., composed of two or more species of different masses, is straightforward.
In order to devise a feasible computational scheme, the path integral is typically discretized by representing the cyclic paths as a finite set of 3N-dimensional configurations, R i , at equispaced points in imaginary time between 0 and ␤ប. The degree of discretization is referred to as the Trotter number, P. The short-time or high temperature propagator, (R i ,R iϩ1 ;␤/ P), can be evaluated semiclassically at different levels of approximation. The contribution of the kinetic energy term to the short-time action is written in terms of a first-order finite difference between configurations on adjacent time slices, while the short-time integral of the potential energy together with higher-order corrections to the kinetic energy is replaced with an effective, quantum-corrected potential V eff (R). The sum of the two terms is referred to as the effective action. Therefore, the expression for the partition function of N interacting, distinguishable quantum particles with Trotter number P is given by

͑4͒
According to the level of approximation of the effective action, the number of slices P needed to achieve convergence in the partition function can be small enough that the problem is tractable, or large enough that the evaluation of the multidimensional integral becomes a hopeless task. It is therefore important to use the best possible effective action compatible with the computational complexity involved in its calculation. The simplest one, or primitive approximation, replaces the effective potential by the bare potential, which is equivalent to an end-point approximation for the short-time integral: and is correct only to order 2 . At the other extreme, the pair-action approximation provides a very accurate technique when the full, many-body potential can be reasonably approximated with a sum of pair potentials. This scheme has been very effectively exploited to investigate the properties of liquid and superfluid He down to temperatures of about 1 K, 21 a task that would not have been possible using the primitive action.
As mentioned in Sec. I, classical interaction potentials are computationally fast, but very often unreliable. Realistic interaction potentials can instead be obtained from more expensive first-principles techniques. In the latter, the electronic degrees of freedom are explicitly included in the Hamiltonian description of the system: Ĥ ͑ R,r͒ϭT n ϩT e ϩV ee ͑r͒ϩV en ͑r,R͒ϩV nn ͑R͒, ͑6͒ where r and R are the electronic and nuclear coordinates, T and V stand for kinetic and potential operators, while the subscripts e and n indicate electronic and nuclear components, respectively. The path integral representation for the partition function could then be developed using the coordinate basis for both the electrons and the nuclei. 22 However, standard electronic structure calculations are carried out in a wave function representation by resorting to the adiabatic separation of nuclear and electronic motion. It is therefore more convenient to expand the electronic component in the adiabatic basis set where electronic wave functions ͉ ␣ ͘ and total energies E ␣ (R) are obtained by diagonalizing the electronic Hamiltonian T e ϩV ee (r)ϩV en (r,R) ϩV nn (R). If ϭ␤/P, then the discretized partition function reads: Then, in the spirit of expression ͑4͒, the short-time propagator can be written as where E ␣ eff (R) is an effective potential which derives directly from the electronic structure. In the primitive approximation, which is what has been used in other ab initio path integral methods, 6,7 the effective potential is simply the total energy corresponding to adiabatic state ␣ of the electronic Hamiltonian, i.e., E ␣ (R).
The nuclear kinetic energy operator can give rise to nonadiabatic coupling matrix elements between adiabatic eigenstates. If these are negligible, but more than one Born-Oppenheimer ͑BO͒ surface is occupied due to thermal excitations, then the partition function will split into independent manifolds indexed by the BO electronic eigenstate. In the absence of degeneracies in the ground electronic state, the energy differences between electronic eigenstates are typically orders of magnitude larger than reasonable thermal kinetic energies. Consequently, only the ground electronic state contributes to the partition function and the electrons only enter at the level of replacing the total potential energy with the ground state first-principles effective potential E 0 eff (R) in Eq. ͑4͒:

͑10͒
This is the partition function that will be evaluated using Monte Carlo techniques. The expression for Z NP can be interpreted as the partition function of N classical polymers, each of P monomeric units or beads, with adjacent beads linked by harmonic springs with force constant mP/␤ប 2 . Beads on two separate cyclic polymers are coupled by the interaction potential only if they lie on the same time slice. Let us remark that the computation of excited electronic states is an open issue in density functional theory, and it is known that the usual approximations to exchange and correlation, like the LDA and even NLDF do not provide reliable excitation energies. Excited states could be calculated properly using very high-level ab initio methods. However, since we have developed a methodology in which the electronic degrees of freedom are dealt within DFT, we are for the time being not in a position to incorporate nonadiabatic couplings and excited electronic manifolds.

Effective action
We use a discretized time representation in which a path is described as a set of configurations, ͕R i ͖, iϭ1,...,P, at P-equispaced points in imaginary time. The effective shorttime propagator for two adjacent points along the path has been evaluated to fourth-order accuracy in ϭ␤/ P, so that the first-principles effective potential reads 23 where x i j is the three-dimensional vector of the coordinates of particle j in slice i, such that R i ϭ(x i1 ,x i2 ,...,x iN ). This form of the effective action leads to an error of the order of P(␤/ P) 5 in the partition function and allows us to significantly reduce the Trotter number required for convergence as compared to the primitive action ͑see Sec. V͒. The quantum correction to the potential requires the evaluation of the first-principles forces F j ϭϪ‫ץ‬E 0 /‫ץ‬x j in the ground electronic state. The cost of this operation is of the same order of magnitude as the rest of the electronic structure calculation, at least in most ab initio or density func-tional schemes. In contrast, second-and higher-order derivatives are sufficiently expensive to evaluate that the computational advantages of a more accurate short-time action are lost.

Normal modes sampling
The above expression for the partition function can be directly used to set up a METROPOLIS MC simulation scheme by assigning the appropriate Boltzmann weight to each configuration of the Nϫ P-dimensional isomorphic classical system. However, as quantum effects increase, the degree of discretization must be increased to maintain accuracy. Since the harmonic force constant between adjacent beads on the quantum polymer is mP/␤ប 2 , increasing the Trotter index results in increasingly stiff harmonic links and the computational problem of ensuring the ergodicity of the METROPOLIS walk becomes intractable. An intuitively appealing and computationally simple way for circumventing this difficulty comes from considering the normal modes of the quantum polymer. 24 In the absence of an interaction potential, all Cartesian degrees of freedom of the system are decoupled. For a single degree of freedom the harmonic intrapolymer potential is given by V p ϭ(mP/2␤ប 2 ) ͚ lϭ1 P (x l Ϫx lϩ1 ) 2 . Diagonalization of the second derivative matrix of this potential leads to the normal coordinates, ͕Q k ͖, kϭ1,...,P, In the normal mode representation, the kinetic energy contribution to the path action is The zero-frequency mode (kϭ P) corresponds to motion of the center of mass of the polymer and makes no contribution to the kinetic energy. All other normal modes would be Gaussian distributed with variance k 2 ϭ␤ប 2 /4mP sin 2 (k/P) if they corresponded to free particles. The potential energy term couples these normal modes and causes distortions from the free-particle distribution. The low-frequency modes correspond to large, collective motions of all beads of the polymer, while the high-frequency modes cause small, local path fluctuations. The normal modes are then used as METROPOLIS variables and the displacements scaled according to the Gaussian dispersions associated with each normal mode.

Observables
The canonical ensemble average of an observable O is given by where Ô is the corresponding quantum mechanical operator. If the operator Ô is diagonal in the coordinate representation, then ͗O͘ϭ ͵ dRO͑R͒͑R,R;␤ ͒.

͑15͒
Our MC strategy samples configurations R with probability proportional to ͑R,R;␤͒, such that equilibrium averages can be readily estimated via discrete summations. Dynamical variables that can be related to the partition function are also straightforward to obtain using thermodynamic estimators. With our choice of the short-time action, thermodynamic estimators of the total energy ͗E͘, the kinetic energy ͗K͘, and the potential energy ͗V͘ are given by

B. Electronic structure calculations
The calculation of the electronic energies is carried out within the framework of DFT. For a given nuclear configuration, Kohn-Sham single-particle equations 25 are solved self-consistently for the electronic density, and the total energy and forces are computed accordingly. Kohn-Sham orbitals are expanded in a Gaussian basis set.
The electronic density is also expanded in an additional Gaussian basis set. 26 The coefficients for the fit of the electronic density are computed by minimizing the error in the Coulomb repulsion energy. The use of this procedure results in an important speedup, since the cost of evaluating matrix elements reduces from O(N 4 ) to O(N 2 M) ͑N is the number of functions in the orbital set, and M the number of functions in the auxiliary set, typically of size comparable to N͒.
Matrix elements of the exchange-correlation potential are calculated by a numerical integration scheme based on grids and quadratures proposed by Becke. 27 During the selfconsistency cycle, the integration is performed on a set of coarse atom-centered, spherical grids. At the end of the selfconsistent procedure, the exchange-correlation energy is evaluated using an augmented, finer grid. This strategy of combining coarse and fine grids results in a considerable gain in computational efficiency, which is very important because this part is one of the main bottlenecks of the calculation.
The exchange-correlation term is described at a gradient corrected NLDF level. Correlation is given by the parametrization of the homogeneous electron gas of Vosko et al. 28 supplemented with the gradient corrections proposed by Perdew. 29 Gradient corrections to the exchange term are taken from Becke. 30 The first derivatives of the energy with respect to the nuclear coordinates, required by the fourth-order effective action, are evaluated by taking analytical derivatives of oneelectron and Coulomb terms, while the exchange-correlation contribution is obtained by numerical integration. 5

A. Validation of the basis set and optimized geometries
We have analyzed five different basis sets for Li 4 and Li 5 ϩ . The first one ͑labeled 1͒ is the standard 3-21 G basis. 31 The second set ͑labeled 2͒ is a double zeta plus polarization basis set, optimized for DFT calculations. 32 The third one ͑labeled 3͒ is the standard 6-311G basis, 33 the fourth set ͑labeled 4͒ is the 6-311G set augmented with a polarization function (6-311G*), and finally the fifth set ͑labeled 5͒ consists of a large uncontracted basis set (13s/9p/1d), proposed by Dunning. 34 The calculations performed with basis sets optimized for standard ab initio calculations ͑labeled 1, 3, 4, and 5͒ have been carried out using an uncontracted auxiliary basis set with a scheme (7s/3p/3d), as proposed in Refs. 32 and 10. The calculations performed with basis set 2 have been carried out using an auxiliary basis set proposed in Ref. 32, with a scheme (7s/2p/1d). Full geometry optimizations without symmetry constraints have been performed in all cases.
In agreement with previous work, 10,15 only one stable minimum with a rhombus geometry has been found for Li 4 ͑see Fig. 2͒, while for Li 5 ϩ we found four stable local minima ͑also shown in Fig. 2͒. The highest energy isomer, i.e., isomer I, consists of two triangles which lie on the same plane, joined by a shared central atom. The second isomer, i.e., isomer II, is similar to isomer I, but now the triangles lie on perpendicular planes. The third isomer, i.e., isomer III, has C 2 symmetry, and can be described as an isosceles triangle plus a dimer. The dimer is located perpendicular to the plane of the triangle, close to its shortest side. The fourth isomer, i.e., isomer IV, has a trigonal bipyramidal structure. It should be pointed out that neither isomer II nor isomer III have been reported in earlier works, probably because of symmetry constraints used during the geometry optimization procedure. The basis set dependence of the binding energies is shown in Fig. 3. Relevant structural parameters for all basis sets considered here are given in Table I The data in Fig. 3 and Table I show that calculations performed using basis set 2 yielded results that deviate considerably from the ones obtained using the very large basis set 5, which can be considered as almost converged. Even if the errors in computed bond lengths and binding energies are not too large ͑about 5% for bond lengths and 5%-10% for binding energies͒, calculations carried out using basis set 2 fail in reproducing the energy sequence of isomers of Li 5 ϩ . This can be ascribed to the fact that the description of the p shell using this set is rather poor, since it contains only one function. On the other hand, calculations performed with basis sets 1, 3, and 4 yield the same energy ordering for isomers of Li 5 ϩ and structural results within 2%-3% from the almost converged, basis set 5, values.  Basis set superposition error ͑BSSE͒ 35 calculations for Li 4 yield 14.36, 1.00, 0.96, 2.34, and 0.13 kJ/mol, for basis sets 1, 2, 3, 4, and 5, respectively. It is clear that in order to reduce BSSE, and obtain meaningful interaction energies, a better description of the p shell than the one provided by the small basis set 1 is required.
In view of these results, the intermediate basis set 3 ͑6-311G͒, which yields results close to the ones obtained with the larger sets is chosen to perform the electronic structure calculations required in the MC simulations. The relative energies with respect to the most stable isomer, isomer IV, for isomers I, II, and III, are 17.58, 15.16, and 9.67 kJ/mol, respectively, for calculations using basis set 3, compared with 19.38, 17.08, and 11.01 kJ/mol, respectively, using the large set 5.

B. Electronic properties: Dipole moments, Mulliken population charges, and eigenvalues
The Li 4 cluster is nonpolar ͑vanishing dipole moment͒ for symmetry reasons. Li 5 ϩ is a charged system, so its dipole moment depends on the choice of the origin. However, it is customary to evaluate the dipole moment using the center of charge as the origin, and in that case it provides a useful indicator of the asymmetry of the charge distribution, and could also be experimentally relevant. Isomers I and II are nonpolar, isomer IV is only slightly polar, but isomer III is considerably polar. The dipole moments of isomers III and IV, computed using basis set 3, are 1.163D and 0.017D, respectively. Mulliken population charges 36 are also useful indicators of the charge distribution. Results obtained with basis set 3 for Li 4 and the four isomers of Li 5 ϩ are shown in Table II. Significant differences are observed between different isomers-even between isomers I and II which are very similar both geometrically and energetically. Both quantities provide useful indicators of isomerization during the MC simulations in the Li 5 ϩ case. In addition, we present the two highest occupied and two lowest unoccupied molecular orbital energies for Li 4 and the four isomers of Li 5 ϩ in Table II. The HOMO-LUMO gap is quite large in all cases, with values around 1 eV. Therefore, it is unlikely that either thermal or quantum fluctuations will contribute to its closure.

A. Sampling strategy and convergence of the path integral with the degree of discretization
We used a simple METROPOLIS algorithm for the PIMC simulations with each trial move consisting of an attempt to move all the normal modes associated with all the particles.
Two different step sizes were used: ␦ c and ␦ s . The maximum displacement of the center of mass was set by ␦ c , and that of the normal modes of order k-associated with a length scale k -by k ␦ s . We have analyzed the possibility of introducing an additional convergence parameter k*, such that modes with kϽk* or kϾ PϪk* are moved with a relatively small step size, ␦ s k , while those associated with small length scale fluctuations are moved by amounts proportional to ␦ l k (␦ l Ͼ␦ s ). However, it turned out that, in this particular case, a single step size ␦ s for all values of k was efficient enough. This is often not the case when a large number of Trotter slices is used. The various parameters were adjusted to keep the overall acceptance ratio around 50% though occasional runs were used with acceptance ratios between 40% and 60%. The same displacement parameter for the center of mass ␦ c was used for both classical and quantum simulations.
Simulations on the Li 2 dimer using a classical potential fitted to first-principles calculations were used to check the convergence of various properties with the degree of discretization. Table III gives the PIMC results for the expectation value of the potential and kinetic energies using the primitive action and the fourth-order corrected form of the action. It can be observed that convergence to within the statistical error bars occurs with a Trotter number of just 4 when using the fourth-order correction, in contrast to 16 when using the primitive action. Errors in ͗V͘ are an order of magnitude less than those in ͗K͘. This is typical of the relatively facile convergence of expectation values of operators diagonal in the coordinates as opposed to those than must be evaluated using thermodynamic estimators. In our experience, structural quantities such as pair distribution functions converge even faster than the potential energy. The increase in the error bars-at constant number of MC configurations-as the Trotter number is increased is also typical of PIMC simulations. Based on our tests with Li 2 we used a Trotter number of 4 at 100 K and 8 at 50 K for the larger clusters.
We have performed a series of classical and quantum TABLE II. Mulliken populations and orbital energies of Li 4 and Li 5 ϩ computed with basis set 3 ͑in Å and eV͒. Atoms are labeled as in Fig. 2. The two lowest energy unoccupied and two highest energy occupied orbital energies are given.  III. Convergence data for Li 2 at 100 K. The number of MC configurations is 5.12 mil with acceptance ratios between 0.4 and 0.7. Error bars are given in brackets. The first two columns correspond to results using the primitive approximation and the last two to results obtained using the fourth-order correction to the effective action. Energies are expressed in K. Monte Carlo simulations for Li 4 and Li 5 ϩ clusters, with temperatures ranging from 50 to 200 K. Each of them consisted of 10 000-15 000 Monte Carlo steps, preceeded by around 1000 steps of thermalization. We stored atomic coordinates, energies, eigenvalues, Mulliken population charges, and dipole moment for each MC step, for later analysis.
In the following we will concentrate only on structural parameters, one-electron eigenvalues, and dipole moments, leaving aside other thermodynamical properties which would need longer simulations to reduce the statistical error bars to useful values.

B. Results for Li 4 and Li 5
؉ As mentioned in Sec. IV, Li 4 has a single, deep minimum at the 1 A singlet state, in the form of a planar rhombus. Due to this fact the cluster is very rigid and thermal and/or quantum effects basically sample the potential energy surface ͑PES͒ around the minimum. The Li 5 ϩ cluster constitutes a somewhat richer example due to the existence of several isomers. In particular, it is interesting to analyze the possibility of thermal activation and quantum tunneling between different regions of configuration space. In order to explore different regions of the PES, classical simulations were started from three different isomers, namely the ground state ͑isomer IV͒ and two higher energy configurations ͑isomers I and III͒. After a few thousand MC steps it was observed that the second and third simulations were attracted toward the ground state basin, showing that our MC strategy is quite effective in equilibrating and exploring configuration space, possibly because interconversion barriers were low. Based on this, quantum simulations were started from the ground state ͑isomer IV͒ and, in order to facilitate the detection of possible tunneling behavior, also from the first excited isomer ͑isomer III͒. Again, the isomerization toward the ground state was observed in the latter; the polymer moved as a whole, without showing any signature of tunneling. Let us mention that none of the higher-energy isomers appeared again during the simulation, although structures slightly reminiscent of isomer III ͑the closest in energy to the ground state͒ were observed. In other words, Li 5 ϩ appears to be unable to sample metastable regions of configurational space out from the ground state. Figure 4͑a͒ shows the pair distribution function g(r) of Li 4 in the classical case, and for different temperatures. It can be observed that the peaks are approximately centered at the optimized zero-temperature distances. Temperature effects consist basically of broadening the peaks; the first of them, corresponding to the first neighbors shell, almost disappears above 200 K. In Fig. 4͑b͒ we show the effects of the quantum nature of the nuclei by comparing simulations performed at 50 K using the classical and quantum schemes. It can be observed that quantum effects generate a pronounced broadening of the peaks, thus demonstrating the importance of their inclusion.
In Fig. 5 we show the pair distribution functions g(r) for Li 5 ϩ in the classical ͑a͒ and quantum ͑b͒ cases. The main panels contain the distribution averaged over all the five particles. Two groups of atoms can be identified: one composed of three atoms more strongly bound, which form the central triangle of the bipyramidal structure IV, and a second group composed of the two external atoms. The upper inset shows the partial g(r) corresponding to these two groups in the classical case at 100 K, and the lower inset in the quantum case at 50 K. No qualitative difference with Li 4 can be observed.
The trends discussed above also hold for the electronic properties, i.e., the distribution of one-particle eigenvalues. As can be observed in Fig. 6͑a͒, the HOMO and LUMO eigenvalue distributions for Li 5 ϩ exhibit significant broadening upon temperature increase. It is to be remarked that the widths are different for different eigenvalues, a fact that could be reflected in the temperature dependence of the optical photoabsorption spectrum. 37 Similar broadening can be observed in the lower panel, corresponding to the quantum and classical simulations for Li 5 ϩ performed at 50 K. As advanced above, the minimum HOMO-LUMO distance never becomes smaller than about 0.5 eV, so that quantum effects cannot promote an eigenvalue crossing which would result in a major modification of the electronic properties.
Another important quantity is the electric dipole moment, because of its experimental relevance. For Li 4 the dipole moment vanishes at zero temperature due to symmetry considerations. However, at finite temperature the cluster samples regions of the PES characterized by a finite dipole moment, such that the mean value is nonvanishing. A more pronounced effect of the same type can be observed in the quantum simulations. The averages computed using the classical and quantum schemes at 50 K are (0.17Ϯ0.07) and  It is interesting to note that where structural properties are concerned, the overall effect of considering quantum nuclei is qualitatively similar to the effect of increasing temperature in the classical simulations. For the closed-shell Li clusters considered in this work the classical temperature equivalent to the quantum system at 50 K is around 150 K. However, the dipole moment and eigenvalue distributions obtained using the quantum mechanical scheme at 50 K are qualitatively similar to those obtained within the classical scheme at about 100 K. This points out an important difference between the two types of correlations involved in the statistical mechanics of quantum systems, namely coherent quantum fluctuations as opposed to incoherent thermal fluctuations. These appear to behave in different ways according to the physical properties under consideration.
Further characterization of the wave packet behavior of the nuclei is provided by the imaginary-time correlation function, R 2 (tϪtЈ)ϭ͉͗r(t)Ϫr(tЈ)͉ 2 ͘, where 0ϽϭtϪtЈ Ͻ␤ប. The value at ϭ␤ប/2 is particularly important because it gives an estimation of the quantum delocalization of the nuclei. The delocalization length, or gyration radius, is about 0.15 Å for both Li 4 and Li 5 ϩ at 50 K. The imaginary time correlation function for Li 4 is shown in Fig. 7, averaged for atoms 1 and 2 ͑lower curve͒ and for 3 and 4 ͑upper curve͒. It can be noticed that the two atoms which are more strongly bound ͑1 and 2 in Fig. 2͒ are less delocalized than the other two, as may be expected.

VI. CONCLUSIONS AND OUTLOOK
We have introduced a novel method for simulating the statistical mechanics of quantum nuclei interacting through first-principles potentials, i.e., that derive directly from the electronic structure. The scheme presented and discussed here combines a path integral description of the nuclear variables with an adiabatic, ground state, density functional description of the electronic degrees of freedom. In the present scheme we have choosen a specific ͑NLDF-Gaussian͒ formulation to solve the electronic structure problem, but it is important to stress that any other implementation is perfectly valid and compatible with the present scheme, e.g., ab initio quantum chemical approaches like Hartree-Fock or MP2, and/or different localized ͑linear combination of atomic orbitals, linear muffin-tin orbitals͒ or extended ͑pseudopotential or augmented PW͒ basis sets. Moreover, the present scheme is extremely simple to interface with any electronic structure code, since the only input needed to compute the statistical weights are the ͑self-consistently͒ converged energy and forces. MD sampling schemes are more involved in this respect.
We have shown the adequacy and the performance of this methodology by simulating quantum nuclear effects in the clusters Li 4 and Li 5 ϩ . The number of imaginary-time slices needed to achieve convergence to a relative error of 0.5% in the nuclear kinetic energy is 4 at a temperature of 100 K and 8 at 50 K. The same level of accuracy is obtained only with 16 slices ͑at 100 K͒ if the primitive approximation to the action is used. This represents a gain of a factor of 4, which is very important due to the high computational cost of these simulations. The level of gain depends, however, on the shape of the potential energy surface sampled by the nuclei.
The results presented here for the above clusters show that, at temperatures below 50 K, quantum nuclear effects are crucial to account for their structural and electronic properties. Pair correlation functions are quite broadened with respect to the classical counterparts, to a level that similar distributions would be obtained for an effective temperature of about 150 K if only thermal, and not also quantum, fluctuations were considered. The results at Tϭ50 K can be considered to be representative of the ground state, as quantum nuclear effects largely overcome thermal motion. Electronic properties like one-electron eigenvalues and dipole moments show the same type of broadened distributions, although the effective classical temperature appears to be slightly lower, around 100 K, instead of 150 K. This is to emphasize that quantum effects cannot be readily mimicked by adding extra thermal fluctuations, because the correlations involved are of a completely different character: Quantum motion is coherent while thermal motion is incoherent.
The simulation technology presented here opens up the possibility of studying the role of light nuclei, especially protons, in biological and chemical systems. The advantage of describing the electronic component using a localized basis set ͑as Gaussians͒, as opposed to a PW basis set, for isolated clusters and molecules as well as for gas-phase reactions, is obvious because the vacuum around is easily taken into account. In addition, very often reactive condensed-phase systems and large biological molecules do not need a full quantum description because the relevant chemical processes occur in a circumscribed region of space. For example, enzymatic reactions require a first-principles electronic descrip-tion only in the vicinity of the active site, while the rest of the system can be treated by means of classical force fields. Therefore, hybrid schemes that combine quantum and classical mechanical descriptions in different spatial regions 38 are appropriate for these situations and, in conjunction with the present methodology, constitute a general computational approach appropriate for studying the effects of nuclear delocalization in the above situations. In fact, an insight into the problem of quantum hydrogen bonding in water has already appeared in the literature, 39 thus signaling the time for a new and exciting area of multidisciplinary research.