Monte Carlo implementation of density-functional theory

PHYSICAL REVIEW B 86, 085115 (2012) Monte Carlo implementation of density-functional theory K. Putteneers* and F. Brosens† Theory of Quantum and Complex Systems, Universiteit Antwerpen, Universiteitsplein 1, 2610 Wilrijk, Antwerpen, Belgium (Received 20 April 2012; published 13 August 2012) We propose a conceptually easy and relatively straigthforward numerical method for calculating the groundstate properties of many-particle systems based on the Hohenberg-Kohn theorems. In this “density-functional Monte Carlo” method a direct numerical minimization of the energy functional is performed by a Monte Carlo algorithm in which the density is simulated by a distribution of Bernoulli walkers. The total number of particles is conserved by construction, unlike for other implementations of density-functional theory. The feasibility of the method is illustrated by applying it to a nanoshell. DOI: 10.1103/PhysRevB.86.085115 PACS number(s): 71.15.Mb, 02.70.Ss, 31.15.E−, 71.10.−w I. INTRODUCTION B. Existing implementations and their practical problems Implementations of density-functional theory1 (hereafter denoted as DFT) are widely used in chemistry and physics2,3 to find ground-state properties of systems of interacting particles, including properties of metallic nanoshells4,5 discussed in this article because of their interesting applications. Although widely used, the existing implementations show some practical problems. These problems, and the proposed solution to them, can only be understood from the basics of DFT, which we therefore first briefly repeat. In current implementations of DFT, the constraint of a fixed number of particles is taken care of by mimimizing the functional   [n(r),λ] = Ev [n] − λ n(r) dr − N , (3) A. The Hohenberg-Kohn theorems The theory originates from two theorems proven by Hohenberg and Kohn1,2 concerning the ground state of interacting particles in an external potential v(r). The first theorem states that the full many-particle ground state is a unique functional of the particle density n(r), at least for nondegenerate ground states and up to an arbitrary constant for v(r). An important aspect of this theorem is the one-to-one correspondence between the particle density distribution n(r) and the value of the energy functional  (1) Ev [n(r)] ≡ v(r)n(r) dr + F [n(r)] in which F [n(r)] ≡ |T̂ + Ûee | = T [n] + Uee [n] is the expectation value of the kinetic and the interaction energy operators T̂ and Ûee in the many-particle state |. Obviously, if the functional (1) is evaluated in the ground-state density nGS (r), it yields the ground-state energy EGS ≡ Ev [nGS ]. Less obvious is the second Hohenberg-Kohn theorem. It states that the energy functional Ev [n] has as its minimum value the correct ground-state energy associated with v(r) if the number of particles N[n] ≡ n(r) dr is kept constant: min Ev [n(r)] n(r) N[n(r)]=N = Ev [nGS (r)]. (2) The Hohenberg-Kohn theorems formally justify that the particle density n(r) can be used as the basic ingredient in solving quantum-mechanical many-particle ground-state problems, as was assumed in the (extended) Thomas-Fermi theory.6,7 1098-0121/2012/86(8)/085115(5) which leads to the Euler equation δEv [n(r)] = λ, δn(r) (4) where λ is a Lagrange multiplier, to be determined such that n(r) dr = N. A widespread methodology to realize this minimization was developed by Kohn and Sham8 in the ansatz F [n] ≈ Ts [n] + Ecoul [n] + EX [n] + EC [n], (5) in which Ts [n] denotes the kinetic energy of noninteracting particles, Ecoul [n] the direct Coulomb interaction, EX [n] the exchange energy, and EC [n] the correlation energy which contains all quantum-mechanical many-particle effects not included in the other terms. Kohn and Sham noticed that if in this ansatz Eq. (4) is evaluated for noninteracting and interacting particles, the resulting equations are formally the same. Based on the fact that the former system can be described by the solutions of single-particle Schrödinger equations, Kohn and Sham developed a self-consistent scheme in which “effective” single-particle Schrödinger equations are solved,  1 2  − 2 ∇ + veff (r) ψi (r) = ǫi ψi (r) (6) with veff (r) = v(r) + vcoul (r) + δEX [n(r)] δEC [n(r)] + , δn(r) δn(r) where the particle density is constructed from the orbitals,  2 n(r) = N i=1 |ψi (r)| , and the direct Coulomb potential is calculated from the particle density using the Poisson equation ∇ · [ε(r)∇vcoul (r)] = 4π n(r). In this approach, Ts [n] is calculated as the sum of the single-particle kinetic energies. The energy eigenvalues ǫi are in fact Lagrange multipliers, resulting from the condition that the N effective single-particle wave functions are orthonormal. Because N single-particle Schrödinger equations have to be solved self-consistently, this orbital method becomes exceedingly demanding if N 085115-1 ©2012 American Physical Society PHYSICAL REVIEW B 86, 085115 (2012) K. PUTTENEERS AND F. BROSENS increases, except if a sufficiently large reduction of the number of orbitals is possible by symmetry arguments, e.g., in crystals. For systems that are too large for the Kohn-Sham approach, one nowadays aims at solving Eq. (4) as directly as possible. However, apart from rather well investigated expressions for EX [n] and EC [n], as needed in the Kohn-Sham scheme, one needs an expression for the kinetic energy functional. Hitherto, no expression has been found for which the results are as satisfactory as those of Kohn-Sham calculations and the research for more adequate functionals is ongoing.9–12 A second problem is that a correct inclusion of quantummechanical many-particle effects and realistic potentials can lead to nontrivial differential, integral, or integrodifferential equations which often can only be solved with additional approximations. Furthermore, somewhat artificial precautions are required to guarantee that the particle density remains non-negative everywhere. Finally, an input value for the Lagrange multiplier λ is needed. This constant is being related for atoms to the negative of the first ionization energy and for solids to the chemical potential.13 Even if one accepts this interpretation also in case approximate functionals are used, the need for such an external input parameter limits the applicability of the method to systems for which it is known from experiment or from other calculations. II. THE METHOD As an implementation of DFT that overcomes all but one of the above-mentioned problems, we propose a direct numerical minimization of the energy functional—without the detour over Eqs. (3) and (4)—by a Monte Carlo algorithm in which the density is garanteed to be non-negative in the entire domain and the number of particles is kept constant by construction. This “density-functional Monte Carlo” (DFMC) method is a pure implementation of the Hohenberg-Kohn theorems (1) and (2). Like more conventional orbital-free DFT methods, the results of DFMC are currently less accurate than those of Kohn-Sham calculations because of the lack of an adequate expression for the kinetic energy functional. This is especially true for small systems. Therefore, the DFMC method is at this moment mainly useful for systems that are too large for orbital-based calculations, e.g., nanoshells4,5,14 as discussed in this article, stellar matter,15 quantum dots,16 and other confined systems.17 When more accurate functionals become available, the quality of the results will increase because the DFMC method is not restricted to the density-functional forms that exist at this moment. Then the DFMC method can be used also for smaller systems. The DFMC method which we propose to implement the Hohenberg-Kohn theorems directly follows the standard Monte Carlo procedure: (1) Construct a trial distribution. (2) Modify the distribution. (3) Accept or reject the new distribution depending on the corresponding energy. (4) Repeat steps 2 and 3 until a minimum value of the energy is found. In DFMC calculations the density distribution of a system is simulated by entities (“walkers”) on a mesh, that can, for example, be an equidistant grid. For an interval characterized by the indices {jα ,jα + 1} in each direction α, there exists a one-to-one correspondence between the particle density nj and the number (Nw )j of walkers in the interval volume Vj , which is given by the relation nj = N (Nw )j Nw Vj (7) with Nw being the total number of walkers. The procedure is initiated by generating an initial walker distribution, drawn according to a trial density profile ntrial (r) of the particles. The number of mesh points as well as the number of walkers should be sufficiently high as to describe the underlying physics with a satisfactory accuracy. The density profile is changed by letting the walkers diffuse over the mesh: in a step k a walker is moved to a neighboring mesh point with equal probability in each direction. We have chosen to describe this diffusion process by Bernoulli walks rather than by Gaussian deviates because of the easy bookkeeping of the corresponding density changes. Unlike more conventional orbital-free DFT, the condition n(r)  0 is automatically satisfied since no walker is moved from an interval with zero concentration. After each step k the energy difference E (k) = (k) (k) (k−1) (k−1) E [n (r)] − E [n (r)] between the energy of the new density profile n(k) (r) and that of the previous one, n(k−1) (r), is calculated. To calculate the functional, any appropriate integration method can be used. If E (k)  0, the new density profile n(k) (r) is accepted. If E (k) > 0, the new density profile is accepted with probability Pa , which gives the possibility of escaping from local minima. This acceptance probability can be constructed by simulated annealing18 or as a threshold acceptance.19,20 The system configuration is allowed to stabilize by gradually lowering the acceptance probability Pa . In order to fulfill the condition of the second Hohenberg-Kohn theorem that the number of particles has to be constant, no branching or killing accelerators, as often used in the quantum Monte Carlo method, are applied. Despite the introduction of a probability Pa for accepting moves with increasing energy, the possibility of ending up in a local minimum cannot be excluded. Therefore several minimization runs have to be performed. III. APPLICATION TO A NANOSHELL We developed this density-functional Monte Carlo method in order to calculate the ground-state density and work function of a metallic nanoshell.14 This spherical symmetric core-shell nanosize particle can be used in, e.g., biomedicine21 and highpressure measurements.22 For the calculations we used the following model and system parameters. In the shell the conduction electrons are considered to move in a uniform neutralizing background with density ρb = 3/(4π rs3 ) with rs being the Wigner-Seitz radius. Realistic shell parameters were considered. The core radius of the nanoshell was taken to be RC = 40 nm, the shell radius RS = 55 nm, the relative permittivity of all layers ε = 1, and the shell was considered to be gold with Wigner-Seitz radius rs = 3. With these parameters the number of conduction electrons is of the order of N ∼ 26 × 106 , which is beyond the scope of Kohn-Sham calculations. Because of 085115-2 PHYSICAL REVIEW B 86, 085115 (2012) MONTE CARLO IMPLEMENTATION OF DENSITY- . . . the spherical symmetry of the system, the density is only radially dependent, n(r) = n(r), and the problem can be treated quasi-one-dimensionally. We used a radial mesh with intervals characterized by an index j which represent a threedimensional volume Vj = 4π (rj3+1 − rj3 )/3, with 20 mesh points per nm. All expressions below are given in atomic units. We used 100 000 walkers to construct a density profile and started with an almost homogeneous electron concentration in the shell. This initial density profile resembles the results obtained by Kohn-Sham calculations4 for small nanoshells, which we consider as a benchmark and which show an almost homogeneous distribution in the shell with a small spill-out into the core and the environment. At present there exist many approximate expressions for (the parts of) the universal functional F [n], Eq. (5); see, e.g., Refs. 3 and 11 and references therein. Because the discussion regarding which functional is most suitable for a given system is far beyond the scope of this proof-of-concept article, we limited the calculations to the basic local-density approximation for F [n]. For the kinetic energy per volume ETs [n] of the noninteracting electrons we considered the Thomas-Fermi expression ETs [n] = (3/10)(3π 2 )2/3 n5/3 , and the exchange energy density was taken to be the Hartree-Fock result23 EX [n] = −(3/4)(3/π )1/3 n4/3 . For the correlation contribution we used the combination proposed in Ref. 24 and used in Ref. 4 within the high-density limit [local Wigner-Seitz radius rs (r) < 1] EC = [A ln(rs ) + B + Crs ln(rs ) + Drs ]n with25 A = 0.0311, B = −0.048, C = 0.0020, and D = −0.0116, and for densities √ corresponding to rs (r)  1 the expression EC = γ n/[1 + β1 rs + β2 rs ] with26,27 γ = −0.1423, β1 = 1.0529, and β2 = 0.3334. The Hartree contribution EH [n] =  1 ρexc (r)φH (r) dr was calculated from the excess charge 2 density ρexc (r) = ρb − n(r) and from the closed form for the Hartree potential φH (r) of a spherical symmetric system  +∞ 4π  r ′ with ε = 1: φH (r) = r ρexc (r ′′ )r ′′2 dr ′′ dr ′ . Because r ′2 0 FIG. 1. (Color online) Minimization run with a fixed threshold probability Pa = 0.5 to illustrate the time scale of a minimization in order to derive adequate minimization parameters. Dashed black line: minimum energy found after the indicated number of iterations. Full red line: last accepted energy of the iteration. One iteration consists of a sweep over the mesh with a possible move of one walker in each interval. FIG. 2. Minimum energies found in 150 minimization runs. The outliers on the figure result from minimization runs ending up in a local minimum. of the rather dense mesh, a simple trapezoidal rule was used to integrate the energy density. The parameters for the minimization procedure were chosen based on a run with a fixed acceptance probability Pa = 0.5 which is shown in Fig. 1. We started the actual calculations with Pa = 0.5 and, based on Fig. 1, we decided to lower the acceptance probability with a power 1/0.6 every 2000 iterations. The minimization run was ended if the threshold reached a value less than 0.001 and for 1000 iterations no lower energy minimum was found. With these parameters one minimization run, written in FORTRAN, took about one minute on a 64-bit desktop (AMD dual core running at 3.0 GHz). We performed 150 runs of which the calculated energy minima are shown in Fig. 2. The outliers on this figure are attributed to local minima. In averaging for a final result, we excluded FIG. 3. (Color online) Calculated ground-state density of a nanoshell with (RC ,RS ) = (40,55) nm from averaging over 78 retained run results corresponding to the global minimum. The relative error is of the order 10−4 –10−3 , which results in error bars that are not visible in the figure. 085115-3 PHYSICAL REVIEW B 86, 085115 (2012) K. PUTTENEERS AND F. BROSENS FIG. 4. (Color online) Calculated ground-state density of a nanoshell with (RC ,RS ) = (55,95) a.u. calculated with the DFMC method (black solid line) and with Kohn-Sham calculations30 (blue dashed line). Error bars are not visible in the figure. all results with a minimum energy above the median so that no local minima are included. To determine the work function W , we calculated the ionization energies Eionization due to the removal of from one to six walkers. Observing a quadratic dependency of these ionization energies as a function of the number of removed electrons N and taking into account that Eionization ( N = 0) = 0, we fitted these six calculated points to the function Eionization ( N) = b × ( N) + c × ( N )2 with the leastsquares method. The high χ 2 probability Q = 0.9999 of the fit confirmed that Eionization ( N ) is quadratic. The work function was then calculated as W = √ Eionization ( N = 1) and its standard deviation as σW = σb2 + σc2 + 2 cov(b,c) with cov(b,c) being the covariance between the parameters b and c. Based on the median of the minimum energies from the 150 runs, 78 results were retained for averaging. For the work function of the nanoshell under consideration we found a value of W = (2.23 ± 0.10) eV. Given the relatively large dimensions of the nanoshell, its work function can be expected to be of the order of the bulk work function of gold, for which measured values lie in the range [5.31–5.47] eV depending on the crystal face being measured.28 The calculated work function of several eV thus has indeed the correct order of magnitude. The difference with the experimental value can be attributed entirely to the used model for the background and the approximate form of the kinetic energy functional. Kohn-Sham calculations for a flat surface, based on the uniform background model with rs = 3 and the same expression for the exchange and correlation energy, result in a work function29 W = 3.35 eV, which is also a few eV lower than the experimental bulk value. Furthermore, orbital-free calculations for a flat surface with the Thomas-Fermi kinetic energy and with the same parameters for rs and the exchange and correlation energy, lead to a work function value29 of W = 2.11 eV. Expecting a somewhat higher work function for the nanoshell because of confinement, we consider the DFMC value for the work function to be consistent with the calculations performed for a flat surface. The difference with experiment is thus not due to the used method but due to the used approximations. Note that the work function, a quantity of the order of a few eV, is obtained from the total energy of a system, a quantity with a magnitude of several 107 eV. This shows that the Monte Carlo calculations reach an acceptable level of accuracy for the energy. The resulting density is shown in Fig. 3. For further verification of the results, we have applied the DFMC method to a nanoshell with (RC ,RS ) = (55,95) a.u. for which KohnSham calculations were performed in Ref. 4 with the same approximations for the background and for the exchange and correlation energy as used in this article. Comparison of our results in Fig. 4 with the results of the Kohn-Sham calculations30 show that the resulting density profile is globally the same: almost homogeneous inside the shell with a spill-out into the core and the environment which is of comparable accuracy in both treatments. The observable differences can be entirely attributed to the different treatment of the kinetic energy in both approaches. IV. CONCLUSION We conclude that density-functional Monte Carlo (DFMC) is a relatively simple method for calculating the groundstate properties of rather large many-particle systems which overcomes several problems in existing implementations of density-functional theory. We show that the DFMC method allows one to treat spherical nanoshells with realistic dimensions, which are inaccessible to orbital-based DFT methods, by calculating the ground-state density and the work function of a nanoshell that contains about 3 × 107 conduction electrons. ACKNOWLEDGMENTS The authors gratefully acknowledge J. T. Devreese for critically commenting on the manuscript and P. Nordlander and V. Kulkarni for providing the numerical data from Kohn-Sham calculations. Research was funded by a grant of the Agency for Innovation by Science and Technology (Agentschap voor Innovatie door Wetenschap en Technologie) Flanders, by the Fund for Scientific Research (Fonds Wetenschappelijk Onderzoek) Flanders, Project No. G.0365.08, and the WOG Project No. WO.033.09N. * 2 † 3 Corresponding author: katrijn.putteneers@ua.ac.be fons.brosens@ua.ac.be 1 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). J. P. Perdew, A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. Scuseria, and G. I. Csonka, J. Chem. Phys. 123, 062201 (2005). 085115-4 PHYSICAL REVIEW B 86, 085115 (2012) MONTE CARLO IMPLEMENTATION OF DENSITY- . . . 4 E. Prodan and P. Nordlander, Chem. Phys. Lett. 349, 153 (2001); 352, 140 (2002); E. Prodan, A. Lee, and P. Nordlander, ibid. 360, 325 (2002). 5 E. Prodan, P. Nordlander, and N. J. Halas, Chem. Phys. Lett. 368, 94 (2003). 6 L. H. Thomas, Proc. Camb. Philos. Soc. 23, 542 (1927). 7 E. Fermi, Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Natur. 6, 602 (1927). 8 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 9 Y. A. Wang and E. A. Carter, in Theoretical Methods in Condensed Phase Chemistry, edited by S. D. Schwartz, Progress in Theoretical Chemistry and Physics Vol. 5 (Kluwer, New York, 2000). 10 J. P. Perdew and L. A. Constantin, Phys. Rev. B 75, 155109 (2007). 11 C. Huang and E. A. Carter, Phys. Rev. B 85, 045126 (2012). 12 J. Xia, C. Huang, I. Shin, and E. A. Carter, J. Chem. Phys. 136, 084102 (2012). 13 R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, New York, 1989). 14 R. D. Averitt, D. Sarkar, and N. J. Halas, Phys. Rev. Lett. 78, 4217 (1997). 15 M. Rotondo, J. A. Rueda, R. Ruffini, and S.-S. Xue, Phys. Rev. C 83, 045805 (2011). 16 R. Pino, A. J. Markvoort, and P. A. J. Hilbers, Physica B 325, 149 (2003). 17 G. S. Ho, C. Huang, and E. A. Carter, Curr. Opin. Solid State Mater. Sci. 11, 57 (2007). 18 S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, Science 220, 671 (1983). 19 G. Dueck and T. Scheuer, J. Comput. Phys. 90, 161 (1990). 20 P. Moscato and J. Fontanari, Phys. Lett. A 146, 204 (1990). 21 R. Bardhan, S. Lal, A. Joshi, and N. J. Halas, Acc. Chem. Res. 44, 936 (2011). 22 N. Van den Broeck, K. Putteneers, J. Tempere, and I. F. Silvera, J. Appl. Phys. 110, 114318 (2011). 23 P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930). 24 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 25 M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957). 26 D. Ceperley, Phys. Rev. B 18, 3126 (1978). 27 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 28 CRC Handbook of Chemistry and Physics, edited by W. M. Hayes, 92nd ed. (CRC Press, Boca Raton, FL, 2012). 29 J. P. Perdew and Y. Wang, Phys. Rev. B 38, 12228 (1988). 30 V. Kulkarni and P. Nordlander (private communication). 085115-5