PRL 113, 096402 (2014)

week ending 29 AUGUST 2014

PHYSICAL REVIEW LETTERS

Kinetic Formulation of the Kohn-Sham Equations for ab initio Electronic Structure Calculations 1

M. Mendoza,1,* S. Succi,2,† and H. J. Herrmann1,3,‡

Computational Physics for Engineering Materials, Institute for Building Materials, ETH Zürich, Schafmattstrasse 6, HIF, CH-8093 Zürich, Switzerland 2 Istituto per le Applicazioni del Calcolo C.N.R., Via dei Taurini, 19 00185 Rome, Italy and Institute for Advanced Computational Science, Harvard University, 29 Oxford Street, Cambridge, Massachusetts 02138, USA 3 Departamento de Física, Universidade Federal do Ceará, Campus do Pici, 60455-760 Fortaleza, Ceará, Brazil (Received 10 May 2013; published 27 August 2014) We introduce a new connection between density functional theory and kinetic theory. In particular, we show that the Kohn-Sham equations can be reformulated as a macroscopic limit of the steady-state solution of a suitable single-particle kinetic equation. We derive a Boltzmann-like equation for a gas of quasiparticles, where the potential plays the role of an external source that generates and destroys particles, so as to drive the system towards its ground state. The ions are treated as classical particles by using either the Born-Oppenheimer dynamics or by imposing concurrent evolution with the electronic orbitals. In order to provide quantitative support to our approach, we implement a discrete (lattice) kinetic model and compute the exchange and correlation energies of simple atoms and the geometrical configuration of the methane molecule. Moreover, we also compute the first vibrational mode of the hydrogen molecule, with both Born-Oppenheimer and concurrent dynamics. Excellent agreement with values in the literature is found in all cases. DOI: 10.1103/PhysRevLett.113.096402

PACS numbers: 71.10.-w, 31.15.A-, 51.10.+y

The calculation of physical properties of interacting many-body quantum systems is one of the major challenges in chemistry and condensed matter. In principle, this task requires the solution of the Schrödinger or Dirac equations for 3N spatial coordinates and N spin variables for electrons, where N is the number of electrons in the system. Since this is computationally prohibitive for all but smallest N’s, the development of approximate models to describe these systems is in continued demand. A very successful formalism in this context is provided by density functional theory (DFT) [1] developed by Hohenberg and Kohn [2] and Kohn and Sham [3]. The Kohn-Sham (KS) approach to density functional theory allows an exact description of the interacting many-particle systems in terms of a series of effective one-particle equations coupled through an effective potential which depends only on the total electronic density. In particular, the ground-state energy of the system is a functional of the electron density, whose exact expression is, however, not known due to the complicated nature of the many-body problem. However, over the years, thanks to a tremendous amount of intensive work, many physical and practical approximations of increasing accuracy have continued to appear [4–12]. Kinetic theory, on the other hand, is the tool of choice for the study of transport phenomena in dilute media. It takes a mesoscopic point of view by defining the probability distribution function of finding a particle at a given position with a given momentum in the six-dimensional one-particle phase space. 0031-9007=14=113(9)=096402(5)

In this Letter, we introduce a new connection between DFT and kinetic theory, whereby the KS orbitals are regarded as the density of quasiparticles moving in imaginary time, with the potential playing the role of a source or sink of quasiparticles. The dynamics of these quasiparticles is governed by a kinetic equation, which, upon Wick rotation to imaginary time, recovers the Kohn-Sham equations in the macroscopic limit of small mean-free path. Using this kinetic approach, we compute the exchange and correlation energies of simple atoms and molecules, particularly the methane molecule (see Fig. 1).

FIG. 1 (color online). Methane molecule, CH4 . The blue (outer) and red (inner) isosurfaces denote low and high electronic density, respectively. Using our model, we have obtained for the angles between bonds, 109.47 deg, and a C-H bond distance of 108.6 pm.

096402-1

© 2014 American Physical Society

PRL 113, 096402 (2014)

PHYSICAL REVIEW LETTERS

Excellent agreement with the literature is reported. In order to perform ab initio molecular dynamics calculations in an efficient way, we chose a convenient collision operator in the kinetic equation deriving a time-dependent equation for each orbital that can be evolved concurrently with the ionic motion (concurrent dynamics, CD for short). The dynamics of the ions by using CD is compared with the respective Born-Oppenheimer (BO) dynamics, finding very good agreement at less computational time. In the KS approach to DFT, the full quantum many-body problem is reduced to a system of N electrons obeying N time-independent Schrödinger equations (KS equations) subject to a mean field potential that depends only on the total electronic density [2,3], namely, ϵi ψ i ¼ −ðℏ2 =2mÞ∇2 ψ i þ Vψ i , where V is the KS potential, m the electron mass, ψ i the ith KS orbital, and ϵi the corresponding energy.PThe total electronic density is approximated as ρ ≃ Ni jψ i j2 , and the KS potential is given by V ¼ V ext þ V ee þ V xc, where V ext is the external potential due to ions and/or external interactions, V ee is the interaction between electrons, and V xc is the exchangecorrelation interaction.PFor a system with N I ions, we write ~ I jÞ, where ZI is V ext ð~rÞ ¼ −ð1=4πε0 Þ NI I ðZI e2 Þ=ðj~r − R the atomic number of the Ith ion, e the charge of the electron, and RI the position of each ion. The electronelectron interaction V ee ¼ eΦ is obtained by solving the Poisson equation for the electrostatic potential Φ, ∇2 Φ ¼ eρ=ε0 . For the exchange-correlation potential V xc , we use the approach by Becke [4] and Lee et al. [5]. Note that the KS potential energy is the connection between our approach and the density functional theory. We start our approach by writing the following kinetic equation, ∂f=∂t þ v~ × ∇f ¼ ΩðfÞ þ Sð~r; v~ Þ, where f ¼ fð~r; v~ ; tÞ is a single distribution function defined in phase space. Here, Sð~r; v~ Þ represents a source term, and Ω is the collision operator. For simplicity, we approximate it with the Bhatnagar-Gross-Krook (BGK) relaxation operator [13]: ΩBGK ðfÞ ¼ −ðf − f eq Þ=τK , where τK is the kinetic relaxation time, and f eq the equilibrium distribution function, typically a local Maxwell distribution. Thus, by a Chapman-Enskog expansion [14], for small values of the Knudsen number, Kn Kn ¼ λ=L ≪ 1, λ ∝ τK being the mean-free path and L the typical system size, we recover the diffusion equation for the Rdensity R field Ξ, ∂Ξ=∂t ¼ τK D∇2 Ξ þ S, where Ξ ¼ fd~v ¼ f eq d~v is the of the distribution function, S ¼ R zeroth-order moment R Sd~v and DΞ ¼ f eq v~ 2 d~v (corresponding to the average kinetic energy of the quasiparticles). By identifying ϕ ≡ Ξ, ðℏ=2mÞϕ ≡ τK DΞ, and −Vϕ=ℏ ≡ S, one recovers the time-dependent Schrödinger equation in imaginary time, for which the KS orbitals ψ i are the stationary-state solutions, such that any wave function ϕ describing the state of a given system can be expanded in this basis. Within the framework Pof the Wick rotation, this can be written as ϕð~r; tÞ ¼ n an exp ð−ϵn t=ℏÞψ n ð~rÞ, where an

week ending 29 AUGUST 2014

are projection coefficients defined via the inner product, an ¼ hψ n jϕi. Therefore, given an initial condition, we see that for negative (positive) energy eigenvalues, every term in the sum grows (decays) proportional to expð−ϵi t=ℏÞ. After some time, the only noticeable term in the sum will be the ground state. Therefore, the wave function ϕ will converge to ϕ ≃ a0 expð−ϵ0 t=ℏÞψ 0 , which, upon normalization, leads to the ground state. This time-projection technique has been used to solve the KS system and obtain the ground state for different electronic configurations [15,16], including diffusion Monte Carlo algorithms [17]. Thus, we observe that the KS orbitals emerge as the macroscopic limit of a distribution function of a gas of quasiparticles in phase space. The energy of each KS orbital can be calculated through the relation, ϵi ¼ −ðℏ=2Þ∂ logðhψ i jψ i iÞ=∂t, when the system has reached its ground state. The kinetic energy K can be computed with no the relation, R need of Laplacian operators, using R Ki ¼ ψ i ∇ × J~ i d~v=hψ i jψ i i, with J~ i ¼ f i v~ d~v. Since we require only the zeroth-, first-, and secondorder moments of the equilibrium distribution (they are sufficient to recover the diffusion equation), there is no need to know the exact analytical expression of the equilibrium distribution or the intrinsic properties of the quasiparticles. Therefore, one can expand the distribution function onto any orthogonal basis of polynomials in velocity space up to second order. As a result, the distribution function f i and the source term Si related to each ith orbital can be written as f i ð~r; v~ ; tÞ ¼ P ðnÞ ðnÞ wð~vÞ ∞ r; tÞHn ð~vÞ, and Si ¼ −ðV=ℏÞχð~vÞψ i , n¼0 an;i ð~ where wð~vÞ and χð~vÞ are suitable weight functions, and ðnÞ an;i a tensor of order n, namely, the projection of the ðnÞ

distribution function upon the tensorial polynomial Hn ð~vÞ R ðnÞ ðnÞ of degree n, an;i ð~r; tÞ ¼ f i ð~r; v~ ; tÞHn ð~vÞd3 v. For simplicity, we use the Hermite polynomials with the weight wð~vÞ ¼ expð−~v2 =2θÞ=ð2πθÞ3=2 , θ being a normalized temperature. By using the first three pffiffiffi normalized Hermite polynomials, H0 ¼ 1, Hk1 ¼ vk = θ, and Hkl 2 ¼ pffiffiffi k l kl ðv v − θδ Þ= 2θ, we can readily show that the equilibrium distribution function f eq i for the ith wave function ~ becomes f eq ð~ r ; v ;tÞ ¼ wð~ v Þψ ð~ v2 −3θÞ=2θ2 , i rÞ½1þðD−θÞð~ i where D ¼ ℏ=ð2mτK Þ. The weight function χð~vÞ is given by χð~vÞ ¼ wð~vÞ½1 − ð~v2 − 3θÞ=2θ. The expression for χð~vÞ ensures that the first and second moments of the source term are zero to prevent changes in the momentum and thermal energy of the system. In fact, one can include further terms in the expansion to impose that moments up to higher orders, N , of the distribution and source term also vanish i.e., R eq n (increasing the accuracy of our Rapproach), f v~ d~v ¼ 0, with 2 < n ≤ N , and χð~vÞ~vn d~v ¼ 0, with 0 < n ≤ N . In order to calculate excited states we, hence, impose the exclusion principle for fermions; we add an interaction

096402-2

∂f i f − f eq i þ Si − V i ; ð1Þ þ v~ × ∇f i ¼ − i ∂t τK P with V i ¼ χð~vÞ j

Kinetic formulation of the Kohn-Sham Equations for ab initio electronic structure calculations.

We introduce a new connection between density functional theory and kinetic theory. In particular, we show that the Kohn-Sham equations can be reformu...
326KB Sizes 2 Downloads 7 Views