h?P Finite Element Approximation for Full-Potential Electronic Structure Calculations

Yvon MADAY

(In honor of the scienti fic heritage of Jacques-Louis Lions)

1 Introduction

The basic problem in quantum chemistry starts from the postulate of the existence of a time dependent complex function of the coordinatesxcalled the wave function Ψ that contains all possible information about the system we want to consider.The evolution of this wave function depends on its current state through the following equation proposed by Schr¨odinger:It involves a potential-energy functionVthat takes into account internal or external interactions as,for instance,those of electrostatic nature;for a single particle,it takes the form

The understanding of what the wave function represents was provided by Born who postulated,after Schr¨odinger,that|Ψ(x,t)|2dxrepresents the probability density offinding at timetthe particle at positionx.The wave function Ψ is thus normalized in such a way that the spatialL2norm of Ψ is 1.The strength of the concept comes from the fact that it applies to any system,in particular to molecules;the coordinatesxare then the positions of each particle(electrons and nuclei)of the system:hencexbelongs to3(N+M),whereNis the number of electrons andMis the number of nuclei.The Schr¨odinger’s equations contains all the physical information on the system it is applied to,it does not involve any empirical parameter except some fundamental constants of physics like the Planck constant,the mass and charge of the electrons and nuclei···.It is thus a fantastic tool to better understand,predict and control the properties of matter from the fundamental background.The very simple Schr¨odinger equation in appearance is however set in a much too high dimensional framework:1+3(N+M),so that it is not tractable for most problems of interest,except that a Quantum Monte Carlo(or QMC for short)approachis used to model and approximate the solutions.These QMC methods allow now to have access to properties other than the energy,including dipole and quadrupole moments,as well as matrix elements between different electronic states.Development and implementation of linear scaling QMC,analytical forces,wave function optimization,and embedding techniques are being pursued(see,e.g.,[35–36]).

For direct methods,though,simplifications need to be proposed to make this much too high dimensional problem accessible to numerical discretizations and simulations.Taking into account the time is quite easy from the principle of the separation of variables in case where the potentialVdoes not depend on time.As it is classical in this approach,the problem becomes time independent and takes the form of an eigenvalue problem:

whereEhas the dimension of an energy.

Through the variation principle,the various solutions to this(linear)eigenproblem,starting by the one associated with the smallest eigenvalue,are associated with a Hamiltonian energy,and,the ground state energy of the molecule corresponds to the smallest eigenvalue in(1.1).This interpretation through a variation principle does not simplify the matter but leads to tractable simplified models.The first one—known as the Born Oppenheimer approximation(see[5])—allows to separate the behavior of nuclei and electrons taking into account their large difference of masses.By considering the nuclei as fixed(or moving very slowly),the problem focuses on the behavior of the electrons—in the so-called electronic structure calculation—and is thus related to the wave function Ψ that depends onNvariables in3(the position of the electrons)and is parametrized by the position of theMnuclei in the associated Hamiltonian.

In order to comply with the Pauli principle of exclusion,the electronic wave function has to be antisymmetric with respect to the electron positions.The electronic problem thus consists in the minimization of the Schr¨odinger’s Hamiltonian over allL2normalized,antisymmetric wave functions.By minimizing instead on a smaller set offunctions provides a tractable problem at the price of yielding to a larger ground state energy.This is the matter of the Hartree Fock problem that consists in minimizing the actual Schr¨odinger’s energy over all wave functions that are written as a so-called Slater determinant,i.e.,a determinant:det[?i(xj)],where the one electron orbitals?i(i=1,···,N)are unknown functions over3.The minimization problem over such Slater determinants leads to a minimization problem involving a new energy.

Let us describe this model associated to a so-called closed-shell system with an even numberN=2Nof electrons,the electronic state is described byNorbitals Φ =(?1,···,?N)T∈(H1(R3))Nsatisfying the orthonormality conditions

and the associated electronic density

The factor 2 in the above expression accounts for the spin.In closed-shell systems,each orbital is indeed occupied by two electrons,one with spin up and one with spin down.

We then introduce the admissible space for molecular orbitals

In the case where the molecular system we consider is in vacuo and consists ofMnuclei of charges(z1,···,zM)∈({0})Mlocated at the positions(R1,···,RM)∈(3)Mof the physical space,and ofNpairs of electrons,the so-called Hartree Fock problem reads:Find Φ0such that

with

An alternative formalism,different in nature but that ends up to a similar mathematical problem,is based on the key result of Hohenberg and Kohn[30]that shows that ground state properties of a system is fully described by the electronic density.This led to the density functional theory,with the instrumental approach of Kohn and Sham[31].Indeed,from[30]the existence of an energy functional of the electronic density was established,this result is weakened however by the lack,even as of today,of knowledge of its proper functional form.It follows from the Hohenberg-Kohn theorem(see[30,37–38,56]),that there exists an exact functional,that is a functional of the electronic densityρthat provides the ground state electronic energy and density of theN-body electronic Schr¨odinger equation.The work of Kohn and Sham addressed this issue by providing approximations of the energy functional and laid the foundations for the practical application of DFT to materials systems.

The Kohn-Sham approach reduces the many-body problem of interacting electrons into an equivalent problem of non-interacting electrons in an effective mean field that is governed by the electron density.It is formulated in terms of an unknown exchange-correlation term that includes the quantum-mechanical interactions between electrons.Even though this exchangecorrelation term is approximated and takes the form of an explicit functional of electron density,these models were shown to predict a wide range of materials properties across various materials systems.The development of increasingly accurate and computationally tractable exchangecorrelation functionals is still an active research area in electronic structure calculations.

In the Kohn-Sham model,also described in the closed-shell configuration,the ground state is obtained by solving the minimization problem:Find Φ0such that

where the Kohn-Sham energy functional reads

The first term models the kinetic energy of Φ,the second term models the interactions between nuclei and electrons,and the third term models the interaction between electrons.The fourth term,called the exchange-correlation functional actually collects the errors made in the approximations of the kinetic energy and of the interactions between electrons by the first and third terms of the Kohn-Sham functional,respectively,as follows from the Hohenberg-Kohn theorem.The lack of precise knowledge for the Kohn-Sham functional is localized on this exchange-correlation term only.It therefore has to be approximated in practice.The local density approximation(or LDA for short)consists in approximating the exchange-correlation functional by

whereis an approximation of the exchange-correlation energy per unit volume in a uniform electron gas with density.The resulting Kohn-Sham LDA model is well understood from a mathematical viewpoint(see[1,33]).On the other hand,the existence of minimizers for Kohn-Sham models based on more re fined approximations of the exchange-correlation functional,such as generalized gradient approximations(see[1])or exact local exchange potentials(see[12])in the general case,is still an open problem.

Note that the Kohn-Sham problem can be split-up into two problems of minimization,with one among them being stated as a pure density problem.We first de fine the set of admissible densities:

then we propose the first problem

followed by the pure density functional problem

There are a lot of variations on the frame of the simulation of these equations.First we may be interested in simulating a molecule alone,small or big,the molecule may also have neighbors,and these can be taken into account exactly or in an average manner like for molecules in solvation(see[52]).When there are many molecules,these can be arranged in a periodic array that is exactly periodic or contains some local defects then the simulation will be done on a very large box composed of many cells,one of them containing the defect.In this case,the simulation domain,sometimes referred to as the supercell,is no longer the whole space R3,as in(1.5);it is the unit cell Γ of some periodic lattice of R3.In the periodic Kohn-Sham framework,the periodic boundary conditions are imposed to the Kohn-Sham orbitals(Born-von Karman PBC).Imposing PBC at the boundary of the simulation cell is the standard method to compute condensed phase properties with a limited number of atoms in the simulation cell,hence at a moderate computational cost.

Both minimization problems(1.2)–(1.5)lead to the resolution of a nonlinear eigenvalue problem,where the eigensolutions are atomic orbitals,function over R3,that thus become tractable to numerical simulations.In order to formulate these eigenproblems,we have to introduce the Hamiltonian for the Hartree-Fock or Kohn-Sham energies:

where

whereVCoulombandVExchangeare defined for anyψby

We notice thatand thus the Euler equations associated with the minimization problem(1.2)read

where theN×Nmatrix,whichis the Lagrange multiplier of the matrix constraint,is symmetric.

In fact,the problem(1.2)has an infinity of minimizers since any unitary transform of the Hartree-Fock orbitals Φ0is also a minimizer of the Hartree-Fock energy.This is a consequence of the following invariance property:

whereU(N)is the group of the real unitary matrices:

1Ndenoting the identity matrix of rankN.This invariance can be exploited to diagonalize the matrix of the Lagrange multipliers of the orthonormality constraints(see,e.g.,[18]),yielding the existence of a minimizer(still denoted by Φ0),such that

for some

Similarly,for the Kohn Sham problem,we introduce the associated Hamiltonian

wherehis the same as above and

The same analysis leads to an eigenvalue problem.By using again the invariance through unitary transforms(1.13),that still holds for the Kohn-Sham problem,we get the existence of a minimizer with a set of molecular orbitals still denoted as Φ0,such that

for some

2 About Numerical Methods

2.1 Generalities

For problems set in a periodic framework(analysis of crystals),the approximation by plane waves(Fourier)has traditionally been one of the popular approaches used for solving the Kohn-Sham problem since it allows for an efficient computation of the electrostatic interactions through Fourier transforms.In addition the plane waves(Fourier)approximation is a high order method that is fully deployed if the solutions to be approximated are very regular.Unfortunately,for full potential electronic structure calculations,the nuclear potential is not smeared out and induces singularities in the solutions(atomic orbitals and density)at the level of the nuclei,more precisely cusps in place of the nuclei and rapidly varying wave functions in their vicinity(see[24,29]).Another drawback of these methods lies in the nonlocality of the basis set that leads to a uniform spatial resolution which can be useless e.g.for materials systems with defects,where higher basis resolution is required in some spatial regions and a coarser resolution suffices elsewhere.In practice of such discretizations,the singular nuclear potentialVnucdefined by(1.4)is usually replaced with a smoother potentialVion;this amounts to replacing point nuclei with smeared nuclei.Not surprisingly,the smoother the potential,the faster the convergence of the planewave approximation to the exact solution of(1.2)or(1.5)(see[8]).The nuclear potentialVnucis replaced by a pseudopotential modeling the Coulomb interaction between the valence electrons on the one hand,and the nuclei and the core electrons on the other hand.The pseudopotential consists of two terms:a local componentVlocal(whose associated operator is the multiplication by the functionVlocal)and a nonlocal component.As a consequence,the second term in the Kohn-Sham energy functional(1.6)is replaced by

The pseudopotential approximation gives satisfactory results in most cases,but sometimes fails.Note that a mathematical analysis of the pseudopotential approximation is still lacking.Moreover,the core electrons need sometimes be considered since they are responsible for intricate properties.The full-potential/all-electron calculation is thus sometimes necessary.In order to overcome the convergence difficulties of the plane wave approximations,resulting from the cusp singularities one can augment the plane waves bases set as done in the augmented plane wave(or APW for short)method(see[42,50]),whichis among the most accurate methods for performing electronic structure calculations for crystals.We refer to[16]for the numerical analysis of the convergence based on the careful analysis of the properties of the cusp that we shall recall in Section 2.2.These APWs provide very good results,at the price however of two remaining drawbacks.The first one is that of the periodic framework that does not fit for single molecules or molecules in solvent.The second one is that the basis come from two different families and the locality of the plane waves(orthonormal basis)is lost.

For efficient computations in the case of all-electron calculations on a large materials system,approximation methods based on Gaussian basis are among the other most classical methods.An example is using the Gaussian package(see[25]).These approaches initially introduced on Hartree-Fock problem,have been developed both for reasons of accuracy and easiness of implementation due to the fact that product of these basis functions arising in nonlinear terms of the potential are easy to evaluate through analytical expressions.The basis functions are centered at each nuclei and are fitted so as to represent well the behavior of the atomic orbital at the level of the cusp and at infinity.There exist a large amount of know how in these methods,that benefit from highly optimized Gaussian basis functions on many molecules.When this expertise does not exist,the approximation properties of the Gaussian expansion are more questionable.We refer e.g.to[9–10]for the presentation and numerical analysis in this context.

Due to the large nonlinearities encountered in the energies involved in advanced Kohn-Sham models,the complexity of the computations,when it turns to implement the methods,scales asO(Nd),whereNis the number of degrees offreedom,anddcan be pretty large(d≥3).One way is to “squeeze” at most the numerical scheme,performing,at the mathematical level,what computational experts in simulations for electronic structure calculations design when they propose ad-hoc discrete basis(e.g.contracted Gaussian bases sets).The expertise here is based on the mathematical arguments involved in model reduction techniques(the reduced basis approximation),and we refer to[11,40]for a presentation of these techniques.They are based on adapted(not universal)discretizations and are shown to provide good approximations,but are still in their infancy.

There is thus room for the development of more robust approaches for electronic structure calculations,like for example finite element approximation of low or high order that,in other contexts(fluid mechanics,structure mechanics,wave,···),are of classical use.There has been already quite a lot of experiences in the domain of quantum chemistry even though the relative number of contributions is still small.We refer to[2,6,21,34,39,43,45–47,51,53–55,58–59]and the references therein for an overview of the contributions in this direction.

In order to be competitive with respect to plane-wave basis or Gaussian type basis,though,the full knowledge and expertise in the finite element machinery has to be requested,indeed,as appears in e.g.[6,28],the accuracy required for electronic structure calculation involves of the order of 100,000 basis functions per atom for P1finite elements,whichis far too expensive and the use of higher-order finite element methods is thus the only viable way.However,the use of high-order finite elements has some consequences on the complexity of the implementation,indeed these require the use of higher-order accurate numerical quadrature rules with larger stencils of points and leads also to increase in the bandwidth of the sti ffness matrices that grow cubically with the order of the finite-element,with a mass matrix that,contrarily to what happens for plane wave approximation,is not diagonal.In addition,using high order methods in regions where the solution presents singularities is a waste of resources.

The right question to raise is thus not accuracy with respect to number of degrees offreedom but accuracy with respect to run time.In this respect,the publication[57]analyses in full details on a variety of problems and regularity of solutions,the accuracy achieved by low to high order finite element approximations as a function of the number of degrees offreedom and of the run time.It appears,with respect to this second argument that the use of degrees between 5 and 8 is quite competitive.Of course,the answer depends on the implementation of the discretization method and the exact properties of the solution to be approximated but this indicates a tendency that is con firmed,both by the numerical analysis and by implementation on a large set of other applications.

A recent investigation in the context of orbital-free DFT indicates that the use of higherorder finite elements can signi ficantly improve the computational efficiency of the calculations(see[44,51]).For instance,a 100 to 1000 fold computational advantage was reported by using a fourth to sixth order finite element in comparison to a linear finite element.This involves a careful implementation of various numerical and computational techniques:(i)an a priori mesh adaption technique to construct a close to optimal finite element discretization of the problem;(ii)an efficient solution strategy for solving the discrete eigenvalue problem by using spectral finite elements in conjunction with Gauss-Lobatto quadrature,and a Chebyshev acceleration technique for computing the occupied eigenspace(see[44]).

As far as we are aware of,the current implementations of the finite element method involve uniform degree of the polynomial approximation.This results in an improved accuracy-per-node ratio that is still polynomial in the number of degrees offreedom.This is actually a bit disappointing since,as is explained in a series of papers by Fournais,S?rensen,Hoffmann-Ostenhof,and Hoffmann-Ostenhof[22–24,29]the solution is analytic(with exponential convergence to zero at infinity on unbounded domains)at least away from the position of the nuclei,where,if exact singular potential are used,the knowledge on singular behavior of the solution is rather well known(this one being of the shape,which results that the solution is not better thanaround the singularities.

In the finite element culture,such behavior—very regular except at some point where the behavior of the pointwise singularity is known—is know to allow for an exponential convergence with respect to the number of degrees offreedom.Indeed,in a series of papers written by Babu?ka and co-authors[3,26–27],a careful analysis is performed that leads to the conclusion that theh?Pversion of the finite element method allows for an exponential rate of convergence when solving problems with piecewise analytic data.In particular,in the three papers[3,26–27],the authors focus on the approximation of the functionover(0,1)forξ∈(0,1),this simple function is a prototype of pointwise singular behavior that can be present in the solution of regular problems in geometries with corners or edges,or for problems with nonregular coefficients.It is straightforward to check that when

We believe that it is interesting to summarize the conclusion of these papers as follows:

(1)For thePversion of the FEM(or spectral method(see[13])):Ifξ∈(0,1),the convergence of the best fit is of the order of,whereris related to the regularity of the function).Note that,ifξ=0(orξ=1),the convergence of the best fit is of the order of.This phenomenon is known as the doubling of convergence for singular functions(see,e.g.,[4]for more results in this direction).

(2)For theh?Pversion of the FEM(or spectral element method):The approximation is generally of the order of

(3)For theh?Pversion,with a graded mesh,i.e.,the size of the mesh diminishes as one gets closer to the singularities andPuniformly increasing:The approximation can be of exponential order with respect to the number of degrees offreedom.

(4)For the optimalh?Pversion of the finite element method:The approximation can be of exponential order with a better rate(with respect to the above rate)if the degreePthat is used in the largest elements increases while the graded meshis refined in the neighborhood of the singularity.Starting from a uniform mesh,the elements that contain the singularity are recursively refined;starting from the singularity,the degree of the approximation is equal to 1 and linearly increases with the distance to the singularities in the other elements;the error then scales like,whereNhis the number of degrees offreedom of the finite element approximation.

2.2 Regularity results

The natural question is then:What is the regularity of the density,the solution to the Hartree Fock or Kohn Sham problems?

It is proven(see,e.g.,the careful analysis of[26–28])that the solution to such systems is analytic(with exponential convergence to zero on unbounded domains)at least away from the position of the nuclei,where,if exact singular potential is used there,the solution is not better thanaround the singularities(this one being of the shape

For the same reasons as the doubling of convergence,if the finite element vertices are on the nuclei positions,then there is a doubling of convergence for thePandh?Pversion of the approximation leading to a convergence rate likeP?3for the solution obtained with polynomial degree 4,if the meshis uniform.The energy is approximated with another doubling of accuracy,i.e.,P?6···.This analysis deals only with the polynomial approximation without taking care of anyheffect and is consistent with the analysis of the paper[20].At this level,we want to emphasize on two points for which we refer to[16].The first one deals with a better understanding of the characteristics of the singularity of the solution of the Hartree Fock problem at the level of nuclei;indeed it appears that locally,when expressed in spherical coordinates around the nuclei,the solution is infinitely differentiable.The second is to indicate that,from this knowledge,it is actually possible to propose combined approximations that allows exponential convergence with respect to the number of degrees offreedom(see also the recent approach in[39]).

In order to be more precise on the regularity of the solutions to such systems,we are going to place ourselves in an adapted framework.We follow[53–54]to define,over any regular bounded domain ? that contains(far from the boundary)all nucleiRj,j=1,···,M,some weighted Sobolev spaces well suited for the numerical analysis of adapted finite element methods(as explained in[17],these weighted Sobolev spaces are well suited to characterize the singular behavior of solutions of general second-order elliptic boundary value problems in polyhedra).First,we define the subdomain ?0whichis the complementary in ? to the union of small enough ballsωjaround each nuclei positionRj,j=1,···,M.In addition,to each nuclei positionRj,j=1,···,M,we associate an exponentβj,and the following semi-norms for anym∈N:

whererjdenotes the distance toRjand norm

wheredenotes the derivative in the local coordinate directions corresponding to the multiindex.

The spaceis then the closure ofof all infinitely differentiable functions that vanish on the boundary of ?.

From the results stated above,we deduce that the solutions of the Hartree-Fock problemfor anyi(1≤i≤N)belong to such spaces(they are said asymptotically well behaved)and moreover

with.In[16],it is indicated that the same type of result can be assumed for the solution to the Kohn-Sham problem,at least for regular enough exchange correlation potential.In what follows we shall assume that the same regularity result holds for those solutions.

3 Galerkin Approximation

3.1 Generalities on the variational approximation

Let us consider a family offinite dimensional spacesXδ,with dimensionNδ.We assume that it is defined through the data of a finite basis set.Let us assume that these are subspacesof(for the time being this means that the discrete functions should be continuous).

The variational approximations of the Hartree-Fock or Kohn-Sham problems are

and

The solution to the Galerkin approximation procedure is determined by

Hence by the determination of the rectangular matrixC∈M(Nδ,N)contains theNδcoefficients of the molecular orbital?i,δin the basis{χμ}1≤μ≤Nδ.It is classical in this context to introduce the so-called overlap matrixSdefined as

so that the constraintson the discrete solutions read

or again in matrix form

Similarly,

wheredenotes the matrix of the operatorin the basis{χk}:

Finally,we can write the Coulomb and exchange terms by introducing first the notations

and,for any matrixXwith sizeNδ×Nδ

The Coulomb and exchange terms can respectively be expressed as

and

The discrete problem,is thus written equivalently as a minimization problem over the space

as

where for anyD∈MS(Nδ),

The energy can also be written in term of the so-called density matrix

leading to the problem

with

Similarly,the Kohn-Sham problem(3.2)reads

with

hereExc(D)denotes the exchange-correlation energy:

For general analysis of these discrete problems and the associated approximation results,we refer to[10–11,18–19,24,36,46,65]and the references therein.

3.2 The adapted h?P discrete spaces

Part 1Definition of theh?Pdiscrete spaces

The purpose of this section is to introduce the class ofh?Pfinite elements spaces for the approximation of the minimization problems(1.2)and(1.5)which we want to propose and analyze in this paper.They will be used to get an exponential convergence for the finite element approximation of the solution to the Hartree-Fock and Kohn-Sham problems.We start by truncating the domain R3in a regular bounded domain ? that,for the sake of simplicity,we shall consider to be a ball large enough to contain largely each nuclei,similarly as in papers where this class of approximations was proposed(see[3,26–27,48–49]).

The first step of the discretization consists in defining an initial triangulationT0of ?composed of hexahedral elementsKsubject to the following classical assumptions:

(1)

(2)Each elementKis the image of the reference cube(?1,1)3under a diffeomorphism,that,in most cases,is a homothetic transformation(except of course on the boundary of ?,but we shall not carefully analyze the approximation there since the solution is very regular at this level);

(3)The initial triangulation is conforming in the sense that the intersection of two different elementsKandK′is either empty,a common vertex,a common edge or a common face;

(4)The initial triangulation is regular and quasi-uniform in the sense that there exist two constantsκ1,κ2>0 with,wherehKdenotes the diameter ofK,ρKdenotes the diameter of the largest circle that can be inscribed intoKand

We assume in addition that the positions of the nucleiR1,R2,···,are the vertices of some element in the initial triangulationT0.Starting from this initial hexahedral triangulationT0of ?,for theh?Pprocedure,we define a family of so-calledσ-geometric triangulations.

First for the triangulation,we refine recursively(by partitioning into 4 hexahedra)each hexahedron that admits a nuclei at the vertex.This partitioning is based on a ratio(σ,1?σ)(with 0<σ<1)of each edge starting from the vertex that coincide at the nucleus,as explained in Figure 1.Note that the refinement only deals the elements that have a nuclei as a vertex,this refinement preserves the conformity with the non-refined elements.The first refinement allows to define the triangulation denoted asT1.Next,the same process applied toTiallows to defineTi+1.For any elementKof the new triangulationsTi(that of course are not quasi-uniform anymore),hKstill denotes the diameter ofK.

Figure 1 The partitioning of a hexahedron,that admits a nuclei at Rias a vertex,in a(σ,1 ? σ)-ratio(the left:before the refinement,the middle:after the refinement,the right:an element created in the refined hexahedron that maintains the conformity with adjacent elements).

Concerning the degree of the approximation that is used over each elements,we start at the levelT0by using polynomials inQ1,i.e.,tri-affine,over each element.Then each time a triangulation refinement is performed,the new elements,created at this stage,are provided withQ1polynomials,while,on the other elements that did not change,the degree of the polynomial is increased by 1 unit in each variable(or by a fixed factorσ>0 to be more general).In particular,the partial degree of the polynomial on the elements of a triangulationTithat have never been cut is≤1+iσand is uniform in all directions.

The discrete finite element spaceis composed of all continuous such piecewise polynomials(associated to the triangulationTi)that vanish on the boundary of ?.From the analysis in 1D explained in[26],exponential convergence for pointwise singular solutions can be obtained from the combination of theσ-geometric mesh refinement andσ-linear increase of the degree of the polynomials.

Part 2Analysis of theh?Papproximation

The difficulty in this analysis lies in the conjunction of three facts:

(1)The problem is set in three dimensions.

(2)Neither the mesh nor the degree from one element to the other is uniform.

(3)We are interested in approximation results for functions asymptotically well behaved.

For asymptotically well behaved functions,we modify the analysis proposed in[49]that dealt with discontinuous finite element methods.We thus start from the the existence of one dimensional operatorsdefined for any integerk≥0 and any integerp≥2k+1:Hk(?1,1)→p(?1,1)such that

These operators can actually be chosen such that the following proposition holds.

Proposition 3.1For every k∈N,there exists a constant Ck>0such that

For integers,p,k∈Nwith p>2k?1,κ=p?k+1and for u∈Hk+s(?1,1)with k≤s≤κ,there holds the error bound

for any j=0,1,···,k.

We then notice that,on the particular mesh of interestTi(except for those elements ofTi—that are actually also inT0—that are close to the boundary and for which there is no problem of regularity nor approximation),there are essentially two reference elements:A perfect cube bK=(?1,1)3and a truncated pyramidlike the one represented in Figure 2 below:

Figure 2 The reference truncated pyramidal element.

Over such a reference element,we want to propose a local reference quasi-interpolantand we follow,for this sake,the same construction as in[49]that is dedicated to the cube It uses the tensorization of the one dimensional operatorthat leads to the operator overdefined by

for which it is proved(see[49,Proposition 5.2]).

Theorem 3.1For any integer3≤s≤p,the operatorsatisfies

and

In order to get the same type of result,over,we modifying the operatordefined above overby using the affine transform from the cube to the truncated pyramid.This results in an operator with range equal to the set of all polynomials overwith partial degree≤p′with respect toxandybut≤3p′with respect to thezdirection.In order to be an operator of degree≤p,we choose,and denote bythis operator for which the following theorem holds.

Theorem 3.2For any integer3≤s≤p,the operatorsatisfies

where the constant cdTP≥1only depends on the shape of.

In order to construct a quasi-interpolant in the equivalentDGspace tobuilt overTi,the operatorwas then used in[49]by scaling it on each element of the meshTiwith the appropriate degree to propose a discontinuous approximation of the solution.Here we need to be a little bit more cautious since we want the approximation to be continuous sinceis a conforming approximation of

We nevertheless proceed as in[49].We first notice that every elementK∈Ti,is canonically associated to a reference elementthat is eitheror.The mapping that allows to go fromKtois denoted asχKand is composed of a rotation,i.e.,a homothetic transformation that thus preserve the polynomial degree.From these transformations,we first build from the operatorsanda totally discontinuous approximation of anyfunctionuwith the appropriate degree as inWe denote bythis operator.From the analysis performed in[49]based on the same regularity results as in(2.3),this first nonconforming approximation satisfies(see(5.21),(5.25)and(5.35)in[49]):

where,for anyK∈Ti,we denote by|Kthe pull back function associated withu|KthroughχK.The argument of the functionis thus points x∈.It can thus be projected with the appropriate reference operatorequal to either

We have now to make the approximation conforming(continuous)between two adjacent elements.This is done by lifting the discontinuities one after the other starting from the discontinuities at the vertex,then at the sides and then,finally at the faces.

Let us start with the vertices.We consider the set of all elements ofTithat share a common vertex a and denote them aswithj=1,···,J.The non-conforming approximationthus proposesJdistinct,but close,values.The rectification first consists in modifying this value so that the new approximation overwithj=2,···,Jis equal toAssume that for a givenwithj=2,···,J,the associatedis,and that the associated pull back transformation maps the vertex a onto(1,1,1).The rectification ofis obtained by adding a quantity

whereHp,1is the polynomialhereLpstands for the Legendre polynomial with degreep,andαis such thatHp,1(1)=1.

This modification,εa,j,is upper bounded by theL∞bound betweenandhence bounded by

From classical considerations we know that

Hence

The associated modification of the approximationofu,that we denote by recta,joveris thus upper bounded with an additional factorchK:

Let us continue on the rectification.We proceed similarly with the edge values,and then the faces values.Let us present the face rectification.We only have two elementsKandK′,that share a whole common face which we denote byFK,K′.The two approximations(already rectified at each vertex and edge)only differ from the internal values on this face.The difference is thus a functionε(x,y)(say)that vanishes on the boundary of the face,and that we are going to lift on the element that has the largest degree(sayK′).This lifting is again performed thanks to a functionHp,1():

the norm of which satisfies

By summing up these three type of corrections,we deduce that the new conforming approximation still satisfies

We finish up as in[49]where the term on the right-hand side above is bounded,with the regularity(2.3)by upper bounding this contribution by exp(?ci),whereiis the index of the triangulationTi,hence by expwhereNiis the total number of degrees offreedom of

Let us now introduce theorthogonal projection operatoroverthus defined as follows:

We can state as follows.

Theorem 3.3There exists a constant C0>0,such that,for all u that satisfies the regularity assumption(2.3),

4 A Priori Analysis

4.1 The Hartree-Fock problem

Part 1Prelimary analysis

Let us first start with the discretization of the Hartree-Fock problem.

Following(3.1),theh?Papproximation of the Hartree-Fock problem is

Remark 4.1The various integrals appearing in this energy should be—and are—generally computed through numerical quadrature.These affect(sometimes dramaticaly)the convergence of the discrete ground state to the exact one.We shall not investigate here this effect that is well described in e.g.[7]in more simple settings.

The lack of uniqueness for the minimization problem,as recalled in(1.13)is a difficulty for the error analysis that requires the understanding of the the geometry of the Grassmann manifoldM,this was first addressed in[41].For eachwe denote by

the tangent space toMat Φ,and we also define

Let us recall(see,e.g.,[41,Lemma 4])that

whereis the space of theN×Nantisymmetric real matrices.

The second order condition associated to the minimization problem(1.2)reads

where for allandin

It follows from the invariance property(1.13)that

This leads us,as in[41],to make the assumption thatis positive definite on,so that,as in[41,Proposition 1],it follows thatis actually coercive on(for thenorm).In all what follows,we thus assume that there exists a positive constantsuch that

Remark 4.2As noticed in[8],in the linear framework,the coercivity condition(4.3)is satisfied if and only if

(i)are the lowestNeigenvalues(including multiplicities)of the linear self-adjoint operator

(ii)There is a gapcΦ0>0 between the lowestNthand(N+1)steigenvalues ofh.

The topology of the Grassmann manifoldquotiented by the equivalence relation through unitary transformations(see e.g.[19])was analyzed in this context in[41]and in[8].In particular,

(1)Let Φ∈Mand Ψ∈M.IfMΨ,Φis invertible,thenis the unique minimizer to the problem

(2)The set

verifies

Hence,for any Φ∈Mand any Ψ∈MΦthere existsW∈Φ⊥⊥such that

whereis anN×Nsymmetric matrix,and the converse also holds.Similarly,at the discrete level,for anyand anythere existssuch that

whereis anN×Nsymmetric matrix,and the converse also holds.

In what follows,we shall compare,as in[8],the error between the solution to(4.1)and Φ0,the solution to(1.2),with its best approximation in

Lemma 4.1(1)If i∈Nis such that

then the unique minimizer of the problemis

In addition,

and for all i large enough,

(2)Let i∈such thatdimThen

The following Lemma 4.2(see[8])collects some properties of the function

Lemma 4.2Let

and S:(the space of the symmetric N×N real matrices)defined by

The function S is continuous on K and differentiable on the interior of K.In addition,

where‖·‖F(xiàn)denotes the Frobenius norm.For all(W1,W2,Z)∈K×K×(L2(?))Nsuch that

Now,we recall the following stability results that can be proved following the same lines as in[8].

Lemma 4.3There exists C≥0such that

(1)for all

(2)for all

In addition,for all(q,r,s)∈such thatand r≤min(q,s),and all0

(1)for allsuch that

Following the same lines as in[8,Lemma 4.7],we deduce from Lemma 4.3 that there existsC≥0 such that for all Ψ∈M,

with

Part 2Existence of a discrete solution

We now use the parametrization of the manifoldexplained in(4.9)to express a given minimizer of the problem(4.1)close to Φ0in terms of an elementWi,h?Pin a neighborhood of 0 expressed aswhere

Indeed,we can define

the minimizers of which are in one-to-one correspondence with those of(4.1).Then the same line as in the proof of Lemma 4.8 in[8]leads to the following lemma.

Lemma 4.4There exist r>0and i0such that for all i≥i0,the functionalEi,h?Phas aunique critical pointin the ball

Besides,is a local minimizer ofEi,h?Pover the above ball and we have the estimate

Since the solution Φ0is asymptotically well behaved,i.e.,satisfies(2.3),we obtain from the accuracy offered by theh?Pfinite elements spacesstated in Theorem 3.3 that there exists a constantCsuch that

Then we deduce the existence and uniqueness of a local minimizer to the problem(4.1)close to Φ0which satisfies that(see[8,(4.71)–(4.72)])forilarge enough,

Hence we have proven the following theorem.

Theorem 4.1LetΦ0be a local minimizer of the Hartree-Fock problem(1.2),and assume that it is asymptotically well behaved,i.e.,satisfies(2.3).Then there exist r0>0and i0suchthat,for anythe discrete problem(4.1)has a unique local minimizerin the set

and there exist constants C>0and c>0,independent ofisuch that

Finally,the discrete solutionsatisfies the Euler equations

whereand theN×Nmatrixis symmetric(but generally not diagonal).Of course,it follows from the invariance property(1.13)that(4.1)has a local minimizer of the formwithU∈U(N)for which the Lagrange multiplier of the orthonormality constraints is a diagonal matrix.

In order to get more precise results on the convergence rate of the eigenvalues,further analysis needs to be performed on the approximation properties in standard Sobolev spaces of the discrete space.As far as we know,these results concerning

(1)inverse inequalities,

(2)convergence properties of theL2(?)orthogonal projection operatorfor e.g.H1functions,

(3)convergence properties of theorthogonal projection operatorfor e.g.H2functions,

do not exist in optimal form,whichis what is required to get the doubling of convergence for the approximation of the eigenvalues with respect to the convergence ofThese results will be proven in a paper in preparation.

4.2 The Kohn-Sham problem

Following(3.2),theh?Papproximation of the Kohn-Sham problem is

Following the same lines as in the proof of the previous result,and the analysis of the plane wave approximation of the Kohn-Sham problem presented in[8],we can prove the following theorem.

Theorem 4.2LetΦ0be a local minimizer of the Kohn-Sham problem(1.5),and assume that it is asymptotically well behaved,i.e.satisfies(2.3),Then there exist r0>0and i0suchthat,for any i≥i0,the discrete problem(4.20)has a unique local minimizerin the set

and there exist constants C>0and c>0,independent ofisuch that

[1]Anantharaman,A.and Cancès,E.,Existence of minimizers for Kohn-Sham models in quantum chemistry,Ann.Inst.Henri Poincaré,26,2009,2425–2455.

[2]Bao,G.,Hu,G.and Liu,D.,An h-adaptive finite element solver for the calculation of the electronic structures,J.Comp.Phys.,231,2012,4967–4979.

[3]Babu?ka,I.and Suri,M.,ThehPandh?Pversions of the finite element method,an overview,Comput.Methods Appl.Mech.Engrg.,80(1),1990,5–26.

[4]Bernardi,C.and Maday,Y.,Polynomial approximation of some singular functions,Appl.Anal.,42(1–4),1991,1–32.

[5]Born,M.and Oppenheimer,J.R.,Zur Quantentheorie der Molekeln,Ann.Physik,84,1927,457–484.

[6]Bylaska,E.J.,Host,M.and Weare,J.H.,Adaptive finite element method for solving the exact Kohn-Sham equation of density functional theory,J.Chem.Theory Comput.,5,2009,937–948.

[7]Cancès,E.,Chakir,R.and Maday,Y.,Numerical analysis of nonlinear eigenvalue problems,J.Sci.Comp.,45(1–3),2010,90–117.

[8]Cancès,E.,Chakir,R.and Maday,Y.,Numerical analysis of the plane wave discretization of some orbitalfree and Kohn-Sham models,ESAIM:Mathematical Modelling and Numerical Analysis,46(2),2012,341–388.

[9]Cancès,E.,Defranceschi,M.,Kutzelnigg,W.,et al.,Computational quantum chemistry:a primer,Handbook of Numerical Analysis,Vol.X,North-Holland,Amsterdam,2003,3–270.

[10]Cancès,E.,Le Bris,C.and Maday,Y.,Méthodes Mathématiques en Chimie Quantique,Springer-Verlag,New York,2006.

[11]Cances,E.,Le Bris,C.,Nguyen,N.C.,et al.,Feasibility and competitiveness of a reduced basis approach for rapid electronic structure calculations in quantum chemistry,Proceedings of the Workshop for Highdimensional Partial Differential Equations in Science and Engineering,Montreal,2007.

[12]Cancès,E.,Stoltz,G.,Staroverov,V.N.,et al.,Local exchange potentials for electronic structure calculations,Mathematics in Action,2,2009 1–42.

[13]Canuto,C.,Hussaini,M.Y.,Quarteroni,A.and Zang,T.A.,Spectral Methods:Fundamentals in Single Domains,Springer-Verlag,New York,2006.

[14]Chen,H.,Gong,X.,He,L.and Zhou,A.,Convergence of adaptive finite element approximations for nonlinear eigenvalue problems,preprint.http://arxiv.org/pdf/1001.2344

[15]Chen,H.,Gong,X.and Zhou,A.,Numerical approximations of a nonlinear eigenvalue problem and applications to a density functional model,M2AS,33,2010,1723–1742.

[16]Chen,H.and Schneider,R.,Numerical Analysis of Augmented Plane Waves Methods for Full-Potential Electronic Structure Calculations,preprint,116.http://www.dfg-spp1324.de

[17]Costabel,M.,Dauge,M.and Nicaise,S.,Analytic regularity for linear elliptic systems in polygons and polyhedra,Math.Mod.Meth.Appl.Sci.,22(8),2012.

[18]Dreizler,R.M.and Gross,E.K.U.,Density Functional Theory,Springer-Verlag,Berlin,1990.

[19]Edelman,A.,Arias,T.A.and Smith,S.T.,The geometry of algorithms with orthogonality constraints,SIAM J.Matrix Anal.Appl.,20,1998,303–353.

[20]Gavini,V.,Knap,J.,Bhattacharya,K.and Ortiz,M.,Non-periodic finite element formulation of orbitalfree density functional theory,J.Mech.Phys.Solids,55,2007,669–696.

[21]Fang,J.,Gao,X.and Zhou,A.,A Kohn-Sham equation solver based on hexahedral finite elements,J.Comp.Phys.,231,2012,3166–3180.

[22]Fournais,S.,Hoffmann-Ostenhof,M.,Hoffmann-Ostenhof,T.and ?stergaard S?rensen,T.,The electron density is smooth away from the nuclei,Commun.Math.Phy.,228(3),2002,401–415.

[23]Fournais,S.,Hoffmann-Ostenhof,M.,Hoffmann-Ostenhof,T.and ?stergaard S?rensen,T.,Analyticity of the density of electronic wavefunctions,Arkiv fr Matematik,42(1),2004,87–106.

[24]Fournais,S.,Srensen,T.?.,Hoffmann-Ostenhof,M.and Hoffmann-Ostenhof,T.,Non-isotropic cusp conditions and regularity of the electron density of molecules at the nuclei,Annales Henri Poincare,8(4),2007,731–748.

[25]Gaussian web site.http://www.gaussian.com

[26]Gui,W.and Babu?ka,I.,Theh,Pandh?Pversions of the finite element method in 1 dimension,Numerische Mathematik,49(6),1986,577–683.

[27]Guo,B.and Babu?ka,I.,Theh?Pversion of the finite element method,Comp.Mech.,1(1),1986,21–41.

[28]Hermannson,B.and Yevick,D.,Finite-element approach to band-structure analysis,Phys.Rev.B,33,1986,7241–7242.

[29]Hoffmann-Ostenhof,M.,Hoffmann-Ostenhof,T.and S?rensen,T.?.,Electron wave functions and densities for atoms,Annales Henri Poincaré,2(1),2001,77–100.

[30]Hohenberg,P.and Kohn,W.,Inhomogeneous electron gas,Phys.Rev.,136,1964,B864–B871.

[31]Kohn,W.and Sham,L.J.,Self-consistent equations including exchange and correlation effects,Phys.Rev.,140,1965,A1133–A1138.

[32]Langwallner,B.,Ortner,C.and Süli,E.,Existence and convergence results for the Galerkin approximation of an electronic density functional,M3AS,20,2010,2237–2265.

[33]Le Bris,C.,PhD Thesis,Ecole Polytechnique,Paris,1993.

[34]Lehtovaara,L.,Havu,V.and Puska,M.,All-electron density functional theory and time-dependent density functional theory with high-order finite elements,J.Chem.Phys.,131,2009,054103.

[35]Lester Jr.W.A.,Recent Advances in Quantum Monte Carlo Methods,World Sientific,Singapore,1997.

[36]Jr.Lester,W.A.,Rothstein,S.M.and Tanaka,S.,Recent advances in Quantum Monte Carlo methods,Part II,World Sientific,Singapore,2002.

[37]Levy,M.,Universal variational functionals of electron densities,first order density matrices,and natural spin-orbitals and solution of the V-representability problem,Proc.Natl.Acad.Sci.USA,76,1979,6062–6065.

[38]Lieb,E.H.,Density functional for coulomb systems,Int.J.Quant.Chem.,24,1983,243–277.

[39]Lin,L.,Lu,J.F.,Ying,L.X.and E,W.N.,Adaptive local basis set for Kohn-Sham density functional theory in a discontinuous Galerkin framework I:Total energy calculation,J.Comp.Phys.,231(4),2012,2140–2154.

[40]Maday,Y.and Razafison,U.,A reduced basis method applied to the Restricted Hartree-Fock equations,Comptes Rendus Mathematique,346(3),2008,243–248.

[41]Maday,Y.and Turinici,G.,Error bars and quadratically convergent methods for the numerical simulation of the Hartree-Fock equations,Numer.Math.,94,2003,739–770.

[42]Martin,R.M.,Electronic Structure:Basic Theory and Practical Methods,Cambridge University Press,Cambridge,2004.

[43]Masud,A.and Kannan,R.,B-splines and NURBS based finite element methods for Kohn-Sham equations,Comput.Methods Appl.Mech.Engrg.,241,2012,112–127.

[44]Motamarri,P.,Nowak,M.R.,Leiter,K.,et al.,Higher-order adaptive finite element methods for Kohn-Sham density functional theory,2012,preprint.arXiv:1207.0167

[45]Pask,J.E.,Klein,B.M.,Fong,C.Y.and Sterne,P.A.,Real-space local polynomial basis for solid-state electronic-structure calculations:a finite element approach,Phys.Rev.B,59,1999,12352–12358.

[46]Pask,J.E.,Klein,B.M.,Sterne,P.A.and Fong,C.Y.,Finite element methods in electronic-structure theory,Comp.Phys.Comm.,135,2001,134.

[47]Pask,J.E.and Sterne,P.A.,Finite element methods in ab initio electronic structure calculations,Mod.Sim.Mat.Sci.Eng.,13,2005,R71–R96.

[48]Sch¨otzau,D.,Schwab,C.and Wihler,T.,hp-dGFEM for second-order elliptic problems in polyhedra.I:Stability and quasioptimality on geometric meshes,Technical Report 2009-28,SAM-ETH,Zrich,2009.

[49]Sch¨otzau,D.,Schwab,C.and Wihler,T.P.,hp-dGFEM for second-order elliptic problems in polyhedra.II:Exponential convergence.Technical report 2009-29,SAM-ETH,Zrich,2009.

[50]Singh,D.J.and Nordstrom,L.,Plane waves,Pseudopotentials,and the LAPW Method,Springer-Verlag,New York,2005.

[51]Suryanarayana,P.,Gavini,V.,Blesgen,T.,et al.,Non-periodic finite element formulation of Kohn-Sham density functional theory,J.Mech.Phys.Solids,58,2010,256–280.

[52]Tomasi,J.and Persico,M.,Molecular interactions in solution:an overview of methods based on continuous distributions of the solvent,Chem.Rev.,94(7),1994,2027–2094.

[53]Tsuchida,E.and Tsukada,M.,Electronic-structure calculations based on the finite element method,Phys.Rev.B,52,1995,5573-5578.

[54]Tsuchida,E.and Tsukada,M.,Adaptive finite element method for electronic structure calculations,Phys.Rev.B,54,1996,7602-7605.

[55]Tsuchida,E.and Tsukada,M.,Large-scale electronic-structure calculations based on the adaptive finite element method,J.Phys.Soc.Jpn.,67,1998,3844–3858.

[56]Valone,S.,Consequences of extending 1 matrix energy functionals from purestate representable to all ensemble representable 1 matrices,J.Chem.Phys.,73,1980,1344–1349.

[57]Vos,P.E.J.,Spencer,S.and Kirby,R.M.,FromhtoPefficiently:Implementing finite and spectral/h?Pelement methods to achieve optimal performance for low-and high-order discretisations,J.Comput.Phys.,229(13),2010,5161–5181.

[58]White,S.R.,Wilkins,J.W.and Teter,M.P.,Finite element method for electronic structure,Phys.Rev.B,39,1989,5819–5830.

[59]Zhang,D.,Shen,L.,Zhou,A.and Gong,X.,Finite element method for solving Kohn-Sham equations based on self-adaptive tetrahedral mesh,Phys.Lett.A,372,2008,5071–5076.

[60]Zhou,A.,Finite dimensional approximations for the electronic ground state solution of a molecular system,Math.Meth.App.Sci.,30,2007,429–447.