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

# Kinetic Formulation of the Kohn-Sham Equations for ab initioElectronic Structure Calculations

M. Mendoza ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials, Schafmattstrasse 6, HIF, CH-8093 Zürich (Switzerland)    S. Succi Istituto per le Applicazioni del Calcolo C.N.R., Via dei Taurini, 19 00185, Rome (Italy),
and Freiburg Institute for Advanced Studies, Albertstrasse, 19, D-79104, Freiburg, (Germany)
H. J. Herrmann ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials, Schafmattstrasse 6, HIF, CH-8093 Zürich (Switzerland) Departamento de Física, Universidade Federal do Ceará, Campus do Pici, 60455-760 Fortaleza, Ceará, (Brazil)
July 17, 2019
###### Abstract

We introduce a new approach to density functional theory based on kinetic theory, showing that the Kohn-Sham equations can be derived as a macroscopic limit of a suitable Boltzmann kinetic equation in the limit of small mean free path versus the typical scale of density gradients (Chapman-Enskog expansion). To derive the approach, we first write the Schrödinger equation as a special case of a Boltzmann equation for a gas of quasi-particles, with the potential playing the role of an external source that generates and destroys particles, so as to drive the system towards the ground state. The ions are treated as classical particles, using 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) model and compute, the exchange and correlation energies of simple atoms, and the geometrical configuration of the methane molecule. Excellent agreement with values in the literature is found.

###### pacs:
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 spatial coordinates and spin variables for electrons, where is the number of particles in the system. Since this is computationally very intensive, 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) DFTbook, developed by Hohenberg and Kohn hohen1 and Kohn and Sham kohn1. The Kohn-Sham approach to density functional theory allows an exact description of the interacting many-particle systems in terms of a series of effective one-particle systems, 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, due to a tremendous amount of intensive work, many physical and practical approximations of increasing accuracy, have continued to appear becke1; lee1; kaxiras1; exact_XC; carpari3; review1; hutter; carpari2; carpari.

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 -dimensional one-particle phase-space. This distribution function lives at the interface between the microscopic dynamics and the macroscopic description in terms of continuum fields representing average quantities over microscopic degrees of freedom ottinger2005.

In this Letter, we introduce a new approach to DFT, based on the kinetic theory, whereby the Kohn-Sham orbitals are regarded as density of quasi-particles, moving in an imaginary time, with the potential playing the role of a source or sink of quasi-particles. The dynamics of these quasi-particles 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). Excellent agreement with the literature is reported. Furthermore, due to the simple structure of the kinetic equation, it is argued that this approach might prove valuable also for developing an evolution equation for the total electronic density.

Let us begin by considering a many-body quantum system consisting of electrons. In the Kohn-Sham approach to DFT, this system is reduced to independent Schrödinger equations (Kohn-Sham equations), subject to a mean field potential that depends only on the total electronic density (Kohn theorem hohen1; kohn1), namely:

 iℏ∂ψi∂t=−ℏ22m∇2ψi+Vψi, (1)

where is the potential, the electron mass, and the index denotes the -th electronic orbital. The total electronic density is approximated as , where is the occupation number (equal to for an open-shell and for a closed-shell) and denotes or for even and odd numbers of electrons, respectively. Since the eigenfunctions of the Hamiltonian form a complete basis in Hilbert space, any wave function that describes the state of a given orbital can be expanded onto this basis,

 ψ′(→r,t)=∑nanexp(iϵntℏ)ϕn(→r), (2)

where are projection coefficients defined via the inner product, , and is the spatially dependent wave function, such that . We have assumed that the potential does not depend explicitly on time. It is well-known that by using the Wick rotation, which consists in replacing (with an fictitious time), one obtains a diffusion equation with source term,

 ∂ψi∂τ=ℏ2m∇2ψi−Vℏψi. (3)

Therefore, given an initial condition and using Eq. (2), we obtain , where we see that for confined systems (negative energy eigenvalues), every term in the sum grows proportional to , while for unconfined systems (positive energy eigenvalues) each term decreases as . In both cases, after some time, the only noticeable term in the sum will be the ground state, namely the one that grows fastest or decreases slowest, depending on the case. Therefore, the wave function will converge to , which, upon normalization, leads to the ground state. This time projection technique has been used to solve the Kohn-Sham system and obtain the ground state for different electronic configurations wick1; wick2.

However, since the diffusion equation can be derived from kinetic theory, it must be possible to recast the Kohn-Sham equations in the form of a kinetic equation in imaginary time. More precisely, the diffusion equation can be obtained from the Boltzmann equation,

 ∂f∂t+→v⋅∇f=Ω(f)+S(→r,→v), (4)

where is a single distribution function defined in the phase space. Here, represents a source term, and is the collision operator. For simplicity, we approximate it with the Bhatnagar-Gross-Krook (BGK) relaxation operator BGK: , where is the kinetic relaxation time, and the equilibrium distribution function, typically a local Maxwell distribution. Thus, by a Chapman-Enskog expansion chapman, for small values of the Knudsen number , defined as being the mean free path and the typical system size, we recover the diffusion equation for the density field ,

 ∂ρ∂t=τKD∇2ρ+S, (5)

where , is the zeroth order moment of the distribution function, and is defined according to the second order moment of the equilibrium distribution, . By comparing Eqs. (11) and (5), we observe that the Kohn-Sham equations emerge as the macroscopic limit of a distribution function of a gas of quasi-particles in phase space, with the identifications , , and .

Since we only require the zeroth, first, and second order moments of the equilibrium distribution (they are sufficient to recover the diffusion equation, and therefore the Kohn-Sham equations as well), we do not need to know the exact analytical expression of the equilibrium function (nor the intrinsic properties of the quasi-particles). Therefore, we can expand the distribution function onto an orthogonal basis of polynomials in velocity space up to second order. As a result, the distribution function, , and the source term, , related to each -th orbital can be written as

 fi(→r,→v,t)=w(→v)∞∑n=0a(n)n,i(→r,t)H(n)n(→v), (6)

and , where is the weight function, a tensor of order , namely the projection of the distribution function upon the tensorial polynomial of degree ,

 a(n)n,i(→r,t)=∫fi(→r,→v,t)H(n)n(→v)d3v. (7)

As mentioned before, we could choose any kind of polynomials and weight functions that satisfy the first three order moment constraints. However, for simplicity we use the Hermite polynomials with the weight , being a normalized temperature. By using the first three normalized Hermite polynomials, , , and , we can readily show that the equilibrium distribution function for the -th wave function becomes,

 feqi(→r,→v,t)=w(→v)ψi(→r)[1+D−θ2θ2(→v2−3θ)], (8)

where . Note that by taking , implying that behaves as the normalized temperature of the gas of quasi-particles, we obtain a simpler expression. However, we shall keep this general form for later applications, because, by increasing , i.e. the diffusivity , the particle dynamics can be accelerated and reach the equilibrium configuration faster.

Up to this point, there is no direct interaction between wave functions, and therefore they all reach the same ground state. In order to calculate excited states, hence impose the exclusion principle for fermions, we add an interaction potential, , to the Boltzmann equation for each wave function . This yields

 ∂fi∂t+→v⋅∇fi=−1τK[(1+VτKℏ)fi−feqi]−Vi, (9)

with , where . By introducing this potential, we dynamically and sequentially remove the contribution of excited levels , with , that do not belong to the respective orbital . This is equivalent to a dynamic Gram-Schmidt orthogonalization procedure. Note that the interaction potential vanishes once the system reaches the ground state, since all orbitals are then orthogonal. The role of this potential is to generate local quasi-particles in order to satisfy the orthogonality conditions between the orbitals.

To summarize, we have converted the system of original Kohn-Sham equations into kinetic equations. In this approach, we need to solve Eq. (9), using Eq. (8). Note that Eq. (9) has no second order spatial derivatives and therefore space and time go on the same footing. This offers a number of computational advantages, as we shall detail shortly.

The energy of each orbital can be calculated through the relation, .

In order to check the validity of our approach, we develop a numerical scheme and implement different simulations for simple atoms and molecules. As it stands, Eq. (9) looks computationally over-demanding, since it lives in a six-dimensional phase space. However, for the purpose of solving hydrodynamic problems, it is known that velocity space can be constrained to a handful of properly chosen discrete velocities , where runs over a small neighborhood of any given lattice site. This strategy has spawned a powerful technique, known as Lattice Boltzmann method, which has proven very successful for the simulation of a broad variety of complex flows Benzi1992145; chen1992recovery; succi2001lattice; electroLB; MUPHY. In the present work we employ the discrete velocities shown in Fig. 2. The result is the following set of Lattice Kinetic Kohn-Sham equations (LKKS):

 fij(→r +→vjδt,t+δt)−fij(→r,t)= (10) −δtτK[(1+VτKℏ)fij(→r,t)−feqij(→r,t)]−Vij,

where . Further details on the discretization and the values of the numerical parameters for the following simulations are presented in the Supplementary Material supp.

We first calculate the exchange and correlation energies for four different atoms, H, Be, B, and C. For this purpose, we have used a lattice size of cells. The measured values for the exchange and correlation energies are given in Table 1, together with the computational time and number of iterations taken by the discretized model. In Fig. 2, we show the orbitals obtained for the carbon atom, finding very good qualitative results.

The differences between the obtained and expected exchange and correlation energies are of the order of except for the case of the correlation energy of the carbon atom. The large error for this case is due to the larger spatial extension of orbital (see Fig. 2) which makes it sensitive to the boundary conditions. The computational time of these simulations is measured when the energy of the orbitals presents changes of less than between two subsequent steps.

Our kinetic equations can also be solved in dynamic fashion, i.e. by evolving the Kohn-Sham orbitals concurrently with the ionic motion (Concurrent Dynamics, CD for short). Although a rigorous theoretical support for such a procedure remains to be developed, we performed a study of how small the time step in the CD implementation of our algorithm should be, in order for the system to remain sufficiently close to the Born-Oppenhemier (BO) surface. To this purpose, we compare both versions of our kinetic scheme, BO and CD, against each other. For this case, we excite the H molecule and let it vibrate in its first mode. We have used a lattice size of cells. Note that by choosing the same time step, fs, (see Fig. 3), CD does not conserve the energy of the system, and the amplitude of the oscillations decays in time. However, by choosing an eight times smaller time step, the CD shows excellent agreement with BO, still performing three times faster than BO (CD took minutes and BO minutes on the same machine). The oscillation frequency of the vibrational mode obtained by the simulation is cm, which presents an error of in comparison with the experimental value of cm H2freq. This discrepancy is not inherent to our approach, but is rather due to discretization errors. Improvements of the numerical integration of the kinetic equations will be a subject of future work.

Finally, we build the methane molecule from scratch, by using the bare Coulomb potential. For this simulation, we use a lattice size of cells, and place the carbon atom in the center of the simulation zone. The hydrogen atoms are located randomly in space, and we let the system evolve to the ground configuration. After a few hours, we achieved the configuration shown in Fig. 1. Note that the angles are reproduced precisely, degrees, but the bond distance of pm is slightly shorter than the experimental value of pm ch4bond. This bond distance discrepancy is due to finite-size boundary effects, which compress the molecule. This statement is supported by further simulations with larger system sizes. The fact that we can reproduce the angles exactly with the bare Coulomb potential, is encouraging, since the use of pseudo-potentials requires orbital dependent parameters, which may significantly increase the complexity of the simulation.

Summarizing, we have introduced a new approach to DFT based on Boltzmann’s kinetic theory. We have shown that the Kohn-Sham equations can be regarded as a macroscopic limit of a first-order Boltzmann kinetic equation. In the kinetic Kohn-Sham equations, we assume that the ground state of a quantum system can be regarded as a gas of quasi-particles, which are generated or destroyed by the potential, the different wave functions interacting with each other in such a way as to match the orthogonality constraints. This opens a new perspective in the interpretation of multi-electron quantum systems.

The kinetic approach offers a number of computational advantages. First, since diffusion emerges from the underlying first-order propagation-relaxation microscopic dynamics, space and time always appear on the same footing (no second order spatial derivatives). This permits to advance the system in time with a time-step scaling only linearly with the mesh resolution, rather than quadratically. Second, since the information always travels along straight lines, defined by the discrete velocities , the streaming is computationally exact (zero round-off error). Third, owing to its excellent amenability to parallel computing, the present lattice kinetic formulation is expected to prove particularly valuable for large-scale simulations.

The lattice kinetic approach has shown excellent performance and accuracy. It can calculate satisfactorily molecular structures by using the bare Coulomb potential, obtaining the right geometric configuration, as it is the case for the methane molecule, without dealing with complicated orbital dependent pseudo-potentials. However, for large molecular structures, the inclusion of pseudo-potentials might well become necessary. We have also performed computational benchmarks and compared the performance with Car-Parrinello molecular dynamics carpari (CPMD), finding satisfactory results for the computational time to ground state, at least for simple molecules. The details of these benchmarks are provided in the Supplementary Material supp.

Finally, the present kinetic approach appears well positioned to solve the evolution equation for the total electronic density, by implementing the analogous of the Chapman-Enskog expansion for classical fluids. This, together with improvements of the discretization, such as local grid refinement, should permit to extend the present lattice kinetic approach to many challenging ab initio calculations, possibly including those currently handled by time-dependent density functional theory EGROSS.

###### Acknowledgements.
Prof. E. Kaxiras is kindly acknowledged for many invaluable discussions and practical suggestions, including the one of performing the methane molecule calculation. We thank Joost Vandevondele for discussions. We acknowledge financial support from the European Research Council (ERC) Advanced Grant 319968-FlowCCS. MM would also like to acknowledge many fruitful discussions with Tobias Kesselring.

## Appendix B Theoretical Background

The aim of this supplementary material is to show the discretized form of the kinetic equations, which in the macroscopic limit reproduce the Kohn-Sham equations with imaginary time,

 ∂ψi∂t=ℏ2m∇2ψi−Vℏψi. (11)

The potential energy of the system, , according to the theorem of Kohn, should only depend on the total electronic density , and, by using mean field approximations, contains four parts , where is the external potential due to ions and/or external interactions, is the interaction between electrons, and is the exchange-correlation interaction. For a multi-electronic system interacting with ions, we can write as, , where is the atomic number of the -th ion, is the charge of the electron, and the position of each ion. The electron-electron interaction, , is obtained by solving the Poisson equation for the electric potential , .

For the exchange-correlation potential , we use the functional derivative of the exchange-correlation energy . Here we use the approach by Becke becke1 and Lee et al. lee1. The Laplacians and gradients of the total density needed to calculated the exchange and correlation potentials, are computed by using an elegant fourth-order method recently proposed in Ref. fourthorder.

The ions move following classical molecular dynamics where the Hellman-Feynman forces due to the electron density are introduced via the electric potential , dictated by the spatial distribution of the electronic charge. In our case, the Hamiltonian of the -th ion is given by , where is the atomic number, the ion mass and the ion momentum. Due to the Born-Oppenheimer approximation, for each movement of the ions, the electronic density must be updated at the ground state. Here we will not consider photon or phonon emissions, i.e. the electrons adapt instantaneously to the ground state after the motion of the ions. Furthermore, for the purpose of this work, we will only take into account the electron-ion and ion-ion electrostatic interations, and therefore consider small time steps in order to conserve the total energy of the system kaxiras1. However, to be precise, one would need to use instead of just the electrostatic interaction, the total Hamiltonian of the electrons. Extensions in this direction will be subject of future research.

## Appendix C Lattice Kinetic Description

Our model is based on the lattice Boltzmann method. Therefore, we define the distribution function that describes the dynamics of the orbital and is associated with the velocity vector . For this purpose, we will use a special case of the Boltzmann equation () that lets the system evolve using just the equilibrium function,

 fij(→r+→vj,t+1)=wjψi(→r,t)[1+θ−D2θ2(3θ−→v2j)], (12)

where are discrete weights and is the lattice temperature, which are defined depending on the cell configuration, and (note that for the analytical equilibrium distribution in the main text, , here the factor will appear as a consequence of the discretization). Standard LB practice shows that in the continuum limit, the above equation converges to a diffusion-reaction equation (such is the Kohn-Sham with imaginary time, Eq. (11)) for the “density” . This is more memory intensive than a standard finite-difference scheme, but offers significant advantages in return. First, since diffusion is emergent from the relaxation to local equilibrium, i.e. the rhs of (12), no Laplacian is required, hence the time-step scales only linearly with the mesh size, rather than quadratically. Second, since the local equilibrium at the rhs is integrally shifted to the neighbor locations , the spatial transport component of the algorithm is virtually exact, i.e. zero round-off error. Third, the local equilibrium conserves the local density to machine roundoff. All of the above configures LB as a very fast, second-order accurate scheme for diffusion-reaction equation. Second order accuracy might appear poor as compared to, say, spectral methods, but the three aforementioned properties make the coefficients of second order errors so small that LB has indeed been shown to provide similar accuracy as spectral methods for, say, turbulence simulations on grids of the order .

The functions contain a forcing term in order to satisfy the second term in the rhs of Eq. (11) and the orthogonality conditions between electronic orbitals. Thus, we can write the wave function as,

 ψj=(ψ′j−Var,j)/(1+V). (13)

Note that the second term is the interaction potential that ensures the orthogonality condition of the wave functions and it vanishes once the system reaches the ground state. This potential is kind of a time-dependent Gram-Schmidt procedure. The auxiliary fields are obtained from the distribution functions , , where is the number of velocity vectors.

Increasing the lattice temperature, , implies to increase the number of velocity vectors, allowing a faster convergence but also introducing errors of higher order derivatives. After a systematic study of different configuration for the velocity vectors, and following the procedure in Ref. karli, we have found a new cell that gives the best performance and accuracy, D3Q25, which can be seen in Fig. 2 in the main text. The weights and velocity vectors, can be derived from the equations in Ref. karli, for the case of .

Thus, at each time step, we know the wave functions , and can calculate the electronic density using . For the computation of the exchange and correlation potentials, we cannot use , since it depends itself of the potentials via Eq. (13), leading to highly complicated system of non-linear algebraic equations. Therefore, as a reasonable guess, we take the density calculated with , .

For the Poisson equation, we implement an additional diffusion equation using as the distribution functions that model the electric potential . For this distribution, we have as evolution equation,

 φi(→r+→vi,t+1)=wiΦ[1+θ−ν2θ2(3θ−→v2i)], (14)

where is a given diffusion parameter that tunes how fast the potential converges to the solution of the Poisson equation. Due to the presence of charge density, we have to modify the potential with the expression, .

The energy level of each orbital can be calculated as, , which converges to a constant value once the system reaches the ground state. For the dynamics of the ions we use the Verlet method. In order to have an accurate value of the forcing term, we use a cubic interpolation to calculate the gradients of the potentials at the location of each ion.

## Appendix D Details of the Different Simulations

In order to validate our model, we have performed several simulations, whose details are presented here. First, we show the details of the simulation for calculating the exchange and correlation energies for four different atoms, H, Be, B, and C. Next, we implement a comparison with Car-Parrinello molecular dynamics carpari (CPMD) for the case of the hydrogen molecule, and finally, we present the numerical details of the simulation of the first vibrational mode of the hydrogen molecule and the respective details of the calculation of the geometry of the methane molecule.

### d.1 Exchange and Correlation energies for several atoms

In all cases, we take a lattice size of cells, and fix the Bohr radius to , , , and cells, for H, Be, B, and C, respectively. We have also set . This implies that we are fixing the relation between the Planck constant and the electronic mass via the parameter . The consequences of this assumption are only numerical, since we can always make a conversion of units where the natural constants, e.g. speed of light, Boltzmann constant, Plank constant, electronic charge and mass, etc., have arbitrary numerical values.

### d.2 Comparison with CPMD

Although a fair comparison with CPMD is not possible because it uses pseudo-potentials, we will perform simulations for small molecules where the valence electrons are also core electrons, i.e. hydrogen and helium molecules. In order to model the hydrogen molecule, we consider two hydrogen atoms separated by a distance (equilibrium configuration), with the Bohr radius, and we let the system optimize the wave function of the ground state, using both, CPMD and our model. Since our model starts with a constant wave function, we will use for CPMD, as initial pseudopotential and wave function, LDA exchange-correlation term, adjusted in the CPMD input file to optimize the wave function taking into account the exchange-correlation model proposed by Becke-Lee-Yang-Parr (BLYP) becke1; lee1. The reason for this choice is to make the comparison between both methods as fair as possible. In our model, we have set .

In Table 2, we show the results of the time consumption for our model for different system sizes. For , and , where we have chosen , and , respectively. We have used different lattice sizes to check convergence and boundary effects. The lattice is quite small and nevertheless can optimize properly the wave function. This is one of the main advantages of our model, the possibility to obtain reasonable accuracy by using few lattice sites. Also, in Table 2, we can see the performance of CPMD for the same configuration. Here, we have used different cutoffs for the plane wave expansion. Note that the times are of the same order as the ones of LB. In both cases, the tolerance of error has been adjusted to .

Let us now take two atoms of hydrogen and locate them at twice their equilibrium distance, , and let the system optimize the position of the atoms until reaching the ground state. For our model, we use again , and we have chosen a time step of ( ps), and the tolerance of both, the electronic density and ion dynamics, of . The lattice size is cells.

In the case of CPMD, we have implemented several simulations for different cutoffs and using for the spatial optimization a tolerance of and for the wave function optimization . In Table 3 are shown the results of the time consumption for both methods and also the final distance between the atoms. For this test, CPMD has shown to be almost one order of magnitude slower than our model, and even without achieving a very good accuracy for the equilibrium configuration which is around . These results show that our model can be very competitive for ab initio simulations of the ground state. All simulations for the comparison with CPMD have been run on a single core of an Intel Xeon E 5430 processor at 2.66 GHz.

### d.3 Concurrence and Born-Oppenheimer Dynamics

The numerical parameters are the same as in the previous simulations in Sec. D.2.

### d.4 Methane molecule

For building the methane molecule from scratch, we use a lattice size of cells, and we place the carbon atom in the center of the simulation zone. The hydrogen atoms are located randomly in space, and we let the system evolve towards the ground configuration. We have chosen the same parameters as before and a Bohr radius of cells. After hours, we have achieved the final configuration.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters