Research interests

My research work focuses on the numerical and theoretical analysis of numerical models and methods in quantum chemistry, whether they come from solid state physics or molecular chemistry. More precisely, they aim to estimate the errors that have several possible origins

  1. modeling errors: the equations which govern quantum chemistry are known - the Schrödinger equation is a high-dimensional linear PDE - which can only be solved for systems reduced to a few electrons. The models (Hartree-Fock or Kohn-Sham) reduce the problem to solving a nonlinear system of PDEs in 3D.

  2. discretization errors: for a given model, this error depends on the numerical scheme considered (plane waves, finite elements or finite differences).

  3. algorithmic errors: at a fixed model and scheme, the chosen algorithm introduces an error which may be difficult to assess. The nonlinear PDEs in quantum chemistry are solved by fixed point methods. The underlying problem being in the vast majority of cases non-convex, the convergence with respect to a numerical tolerance is unclear.

Tensor Networks

Tensor Networks (TN) have seen rapid and emerging interests. They aim to give a sparse representation of high-dimensional functions. The first instance of TN was Matrix Product States (MPS) or Tensor Trains that were introduced in the late 90s/early 2000s in physics as DMRG (Density Matrix Renormalization Group) to find the ground-state of statistical models (one-dimensional Hubbard or Heisenberg models). MPS has been used for the NN-body electronic Schrödinger problem where the question of the optimal topology (or best ordering) is not clear.

In the contribution, we investigate the behaviour of the singular values of reshapes of a particular class of functions (Slater determinants). The decay of these singular values drives the approximability of the function by an MPS. For this class of functions, if we denote by (σj)(\sigma_j) the singular values of the reshaped tensor Ψμ1μL/2μL/2+1μL\Psi_{\mu_1 \dots \mu_{L/2}}^{\mu_{L/2+1}\dots \mu_L}, then we show that there are 2N2^N positive singular values, where NN is the number of particles and they satisfy

σ1σ2N=σ2σ2N1=σ3σ2N2==Const. \sigma_1 \sigma_{2^N} = \sigma_2 \sigma_{2^N-1} = \sigma_3 \sigma_{2^N-2} = \cdots = \mathrm{Const}.

The constant is moreover explicit and depends on the ordering of the modes of the tensor. This provides a natural criterion to optimze a TT approximation of the state.

 Singular values of the reshaped tensor for different ordering schemes
Singular values of the reshaped tensor for different ordering schemes


Linear response in quantum chemistry

Properties of a molecule in its ground-state perturbed by a time-dependent potential can be approximated in the linear regime (assuming that the perturbation is small) using the Kubo formula. To be more precise, if ψ0\psi_0 is the ground-state of H=Δ+VH = -\Delta + V an operator acting on L2(Rd)L^2(\mathbb{R}^d), for some observable VOV_\mathcal{O} and ψ(t)\psi(t) solution to

itψ=Hψ+εf(t)VPψ,ψ(0)=ψ0, i\partial_{t} \psi = H \psi + \varepsilon f(t) V_{\mathcal P} \psi , \quad \psi(0) = \psi_{0},

the linear response function KK defined by

ψ(t),VOψ(t)=ψ0,VOψ0+ε(Kf)(t)+O(ε2), \langle \psi(t), V_{\mathcal O} \psi(t) \rangle = \Big\langle \psi_{0}, V_{\mathcal O} \psi_{0}\Big\rangle + \varepsilon (K \ast f)(t) + \mathcal{O}(\varepsilon^2),

is given by the Kubo formula

K(τ)=iθ(τ)VOψ0,ei(HE0)τVPψ0+c.c.. K(\tau) = -i \theta(\tau) \Big\langle V_{\mathcal O} \psi_{0}, e^{-i (H-E_{0}) \tau} V_{\mathcal P} \psi_{0}\Big\rangle + {\rm c.c.}.

The Fourier transform of KK is then

Kundefined(ω)=limη0+ψ0,VO(ω+iη(HE0))1VPψ0ψ0,VP(ω+iη+(HE0))1VOψ0. \widehat K(\omega) = \lim_{\eta \to 0^{+}} \Big\langle \psi_{0}, V_{\mathcal O} \Big(\omega +i\eta - (H-E_0)\Big)^{-1} V_{\mathcal P} \psi_{0}\Big\rangle - \Big\langle \psi_{0},V_{\mathcal P} \Big(\omega +i\eta + (H-E_0)\Big)^{-1} V_{\mathcal O} \psi_{0}\Big\rangle.

The above formula is in practice approximated by taking a positive η\eta and by truncating the space to a finite region (L,L)d(-L,L)^d. This necessarily yields a singular approximation of the linear response function Kundefined\widehat{K}. In the paper, the regularity of the response function Kundefined\widehat{K} is determined, depending on the decay of the potential VV. Moreover, the relationship between η\eta and LL guaranteeing the convergence of the approximation is shown, supported by numerical experiments.


The particle-hole random phase approximation in density functional theory

The ground-state of a molecule is given by the lowest eigenvalue of the operator HNH_N acting on i=1NL2(R3×Z2)\bigwedge_{i=1}^N L^2(\mathbb{R}^3\times \mathbb{Z}_2) with domain H2(R3N×Z2)H^2(\mathbb{R}^{3N}\times \mathbb{Z}_2)

HN(vext,w)=i=1N(12Δri+vext(ri))+1i<jNw(rirj), H_N(v_\mathrm{ext},w) = \sum\limits_{i=1}^N \Big(-\frac{1}{2}\Delta_{r_i} + v_{\mathrm{ext}}(r_i) \Big) + \sum\limits_{1 \leq i < j \leq N} w(r_i-r_j),

where NN is the number of electrons, vextv_\mathrm{ext} accounts for the interaction with the nuclei and ww is the electron-electron interaction. Solving the high-dimensional eigenvalue problem is quickly intractable for systems with more than a few electrons, but thankfully, the ground state energy can be computed using density functional theory (DFT). DFT states the existence of a universal functional FF of the electronic density ρ(r)=NR3(N1)Ψ0(r,r)2dr\rho(r) = N\int_{\mathbb{R}^{3(N-1)}} |\Psi_0(r,\overline{r})|^2 \, \mathrm{d}\overline{r} such that

E0=minΨ=1Ψ,HNΨ=minρINF(ρ)+vext,ρ, E_0 = \min_{\|\Psi\|=1} \langle \Psi, H_N \Psi \rangle = \min_{\rho \in \mathcal{I}_N} F(\rho) + \langle v_{\mathrm{ext}}, \rho\rangle,

where IN\mathcal{I}_N is the admissible set of electronic densities.

It is a tremendous reduction of dimensionality, at the expense that the functional FF is unknown... Over the past 60 years, great efforts have been put in the design of more and more accurate density functionals, that have been used for systems with tens of thousands electrons, but it appeared that the simple dissociation of the H2_2 molecule, i.e. a two-electron molecule remains a challenging system for most density functionals. At the dissociation limit, when the hydrogen atoms are pulled away to infinity, we expect that the energy of the system at the limit is twice the energy of a single hydrogen atom. This fails for example for restricted Hartree-Fock as the spin restriction introduces a spurious self-interaction energy that can only be lifted by breaking the spin symmetry. In other words, the correct depiction of the dissociation requires a superposition of two determinantal states, hence it cannot be captured at the level of restricted Hartree-Fock.

phRPA, also simply called RPA, is one of the few ansatz for the correlation energy that is able to dissociate H2_2 correctly. In the paper with Kyle Thicke, we show that by taking the RPA correlation energy, which can be derived from an adiabatic connection argument that is explained, on top of a restricted Hartree-Fock functional, we have the exact dissociation of the H2_2 molecule.


The random phase approximation of TDDFT

TDDFT or Time Dependent Density Functional Theory is the time-dependent equivalent of DFT. It states that the evolution of the electronic density ρΨ\rho^\Psi of the many-body Schrödinger equation

{iΨt(t)=(12Δ+i=1Nv(t,ri)+1i<jNw(rirj))Ψ(t)Ψ(0)=Ψ0, \left\lbrace \begin{aligned} \mathrm{i} \frac{\partial \Psi}{\partial t}(t) &= \Big( -\frac{1}{2} \Delta + \sum_{i=1}^N v(t,r_i) + \sum_{1 \leq i < j \leq N} w(r_i-r_j) \Big) \Psi(t) \\ \Psi(0) &= \Psi_0, \end{aligned} \right.

where Ψ0\Psi_0 is the ground-state of the rest Hamiltonian H=12Δ+i=1Nv(0,ri)+1i<jNw(rirj)H = -\frac{1}{2} \Delta + \sum_{i=1}^N v(0,r_i) + \sum_{1 \leq i < j \leq N} w(r_i-r_j), can be reproduced by the noninteracting evolution

{iΦt(t)=(12Δ+i=1Nvrg[{ρΦ(s)}s>0](t,ri))Φ(t)Φ(0)=Φ0, \left\lbrace \begin{aligned} \mathrm{i} \frac{\partial \Phi}{\partial t}(t) &= \Big( -\frac{1}{2} \Delta + \sum_{i=1}^N v_{\mathrm{rg}}[\lbrace \rho^{\Phi(s)} \rbrace_{s>0}](t,r_i) \Big) \Phi(t) \\ \Phi(0) &= \Phi_0, \end{aligned} \right.

with a potential vrgv_\mathrm{rg} called the Runge-Gross potential. This is the Runge-Gross theorem which has been proved for one-body potentials that are regular enough and analytical in time. In particular, it is not known if it holds for Coulomb potentials. Nevertheless, assuming that the Runge-Gross theorem is valid, it gives a way to approximate the excited states of the rest Hamiltonian HH via the density-density linear response function (DDLRF).

For a one-body observable vOv_\mathcal{O} and a one-body perturbation v(t,r)=vext(r)+εf(t)vP(r)v(t,r) = v_\mathrm{ext}(r) + \varepsilon f(t) v_\mathcal{P}(r), the DDLRF is defined by

Ψ(t),i=1NvO(ri)Ψ(t)=vO,ρΨ(t)=vO,ρΨ(0)+εvO,χvPf(t)+O(ε2). \Big\langle \Psi(t), \sum_{i=1}^N v_\mathcal{O}(r_i) \Psi(t) \Big\rangle = \langle v_\mathcal{O}, \rho^{\Psi}(t) \rangle = \langle v_\mathcal{O}, \rho^{\Psi}(0) \rangle + \varepsilon \langle v_\mathcal{O}, \chi v_\mathcal{P} \star f(t) \rangle + \mathcal{O}(\varepsilon^2).

The Fourier transform of the DDLRF χundefined\widehat{\chi} has poles that are located at an eigenvalue of HE0H-E_0

vO,χundefined(ω)vP=limη0+Ψ0,i=1NvO(ri)Ψ0(ω+iη(HE0))1i=1NvO(ri)Ψ0Ψ0,i=1NvO(ri)Ψ0(ω+iη+(HE0))1i=1NvO(ri)Ψ0. \begin{aligned} \langle v_\mathcal{O}, \widehat{\chi}(\omega) v_\mathcal{P} \rangle = \lim_{\eta \to 0_+} \Big\langle \Psi_0, & \sum_{i=1}^N v_\mathcal{O}(r_i) \Psi_0 \big( \omega + \mathrm{i} \eta -(H-E_0) \big)^{-1} \sum_{i=1}^N v_\mathcal{O}(r_i) \Psi_0 \Big\rangle \\ & \Big\langle \Psi_0, \sum_{i=1}^N v_\mathcal{O}(r_i) \Psi_0 \big( \omega + \mathrm{i} \eta +(H-E_0) \big)^{-1} \sum_{i=1}^N v_\mathcal{O}(r_i) \Psi_0 \Big\rangle. \end{aligned}

Assuming that the Runge-Gross theorem holds true, the DDLRF χ\chi is the solution to the TDDFT Dyson equation

χundefined(ω)=χundefined0(ω)+χundefined0(ω)fundefinedxc(ω)χundefined(ω), \widehat{\chi}(\omega) = \widehat{\chi}_0(\omega) + \widehat{\chi}_0(\omega) \widehat{f}_{\mathrm{xc}}(\omega) \widehat{\chi}(\omega),

where χ0\chi_0 is the noninteracting DDLRF of the frozen Hamiltonian

H0=12Δ+i=1Nvext(ri)+ρΦ(0)1(ri)+vxc[ρΦ(0)](ri). H_0 = -\frac{1}{2} \Delta + \sum_{i=1}^N v_\mathrm{ext}(r_i) + \rho^{\Phi}(0) * \frac{1}{|\cdot|}(r_i) + v_{\mathrm{xc}}[\rho^{\Phi}(0)](r_i).

fxcf_\mathrm{xc} is the exchange-correlation kernel formally defined as fxc(rt,rt)=δvrg(rt)δρ(rt).f_\mathrm{xc}(rt,r't') = \frac{\delta v_\mathrm{rg}(rt)}{\delta \rho(r't')}. The poles of the exact DDLRF can be directly inferred from those of the noninteracting DDLRF by solving the Dyson equation. Poles of χundefined(ω)\widehat{\chi}(\omega) are displaced and this is called pole shifting.

Since the Runge-Gross potential has no explicit expression, the exchange-correlation kernel has to be approximated. In the random phase approximation (RPA), it is simply given as

fxc(RPA)(rt,rt)=δttrr. f_{\mathrm{xc}}^{(\mathrm{RPA})}(rt,r't') = \frac{\delta_{tt'}}{|r-r'|}.

In the paper, we show that the number of poles of the RPA DDLRF is the same as those of noninteracting DDLRF. Moreover, because of the positivity of the Hartree kernel 1rr\frac{1}{|r-r'|}, the poles of the noninteracting DDLRF are "forward" shifted, i.e. for all k1k\geq 1, ωk(RPA)ωk\omega_k^{(\mathrm{RPA})} \geq \omega_k where (ωk(RPA))(\omega_k^{(\mathrm{RPA})}) and (ωk)(\omega_k) are the poles of the RPA and the noninteracting DDLRF. This is a mathematical explanation of the well known empirical fact that the noninteracting poles (i.e. the spectral gaps of the Kohn-Sham equations) underestimate the true transition frequencies.


Anderson acceleration

Anderson acceleration, also known as DIIS (Direct Inversion of the Iterative Space) in quantum chemistry, aims to accelerate the convergence of fixed-point iteration xk+1=g(xk)x_{k+1} = g(x_k). Such fixed-point problems are ubiquitous, and in quantum chemistry come from the Hartree-Fock or Kohn-Sham equations. DIIS is an extrapolation method where the next iterate is built using the knowledge of previous iterates. The DIIS sequence is given by

(ci(k))0im=arg minci=1i=0mcif(xki)2andxk+1=g(i=0mci(k)xki), (c^{(k)}_i)_{0 \leq i \leq m} = \argmin_{\sum c_i = 1} \left\|\sum\limits_{i=0}^m c_i f(x_{k-i})\right\|_2 \quad \text{and} \quad x_{k+1} = g\left(\sum_{i=0}^m c^{(k)}_i x_{k-i}\right),

where ff is an error function that locally vanishes only at the solution of the fixed-point problem. A usual choice is f(x)=xg(x)f(x) = x-g(x).

An important parameter in the efficiency of the method is the width mm of the history/number of iterates kept. Two variants are considered in our publication:

For both variants, we prove local linear convergence and superlinear convergence by tuning the DIIS numerical parameters.

 Comparison between the DIIS variants with fixed depth (FD-DIIS), with restarts (R-DIIS) and adaptive depth (AD-DIIS)
Comparison between the DIIS variants with fixed depth (FD-DIIS), with restarts (R-DIIS) and adaptive depth (AD-DIIS)


Solid state physics

In solid-state physics, the periodic boundary conditions motivate a plane-wave discretisation. However because of the Coulomb singularities at the nuclei, the eigenfunctions of the Hamiltonian Δ+V-\Delta + V have cusps which impedes the convergence. Pseudopotential methods are widely used, where the Coulomb singularities are smoothened. In the PAW method (Projector Augmented-Wave), introduced by Blöchl, the pseudopotential is introduced using an invertible operator acting locally around the nuclei. The PAW equations involve infinitely many terms that are truncated in the numerical treatment which causes a modelling error. This modelling error has been analysed for the one-dimensional model d2dx2Z0kZδkZakZδa+k-\frac{\mathrm{d}^2}{\mathrm{d}x^2} - Z_0 \sum_{k \in \mathbb{Z}} \delta_k - Z_a\sum_{k \in \mathbb{Z}} \delta_{a+k} for which explicit formulas of the eigenpairs are available.

Another approach consists in introducing an invertible operator TT acting locally around the nuclei. This operator is defined such that it maps smooth functions to functions with cusps satisfying the Kato cusp condition. This operator acts as a preconditonner for a plane-wave discretisation. This method is called the Variational Projector Augmented-Wave and has been analysed for one-dimensional and three-dimensional Hamiltonians.

 Comparison between a direct discretization and the VPAW method for $-\Delta - \frac{1}{|x-R|}-\frac{1}{|x+R|}$ with periodic boundary conditions
Comparison between a direct discretization and the VPAW method for $-\Delta - \frac{1}{|x-R|}-\frac{1}{|x+R|}$ with periodic boundary conditions