My interest

My current research interest is to develop novel, efficient and reliable numerical algorithms and mathematical software tools for ground state and excited state electronic structure theories. I find this field particularly exciting because
  1. Scientifically, it is amazing that many (electronic, mechanical, optical) properties calculated from such first principle theories (i.e. the only input is atomic specicies and atomic positions. There is in principle no dependence on empirical or fitting parameters) can be directly compared to experimental results.
  2. Mathematically, the theories are relatively well defined. They also provide mathematical and computational challenges not commonly seen in other branches of applied mathematics. The complexity of electronic structure theories can range from high (like density functional theory) to very high (like configuration interaction), and hence numerical methods with reduced complexity are urgently needed in order to simulate ever larger relastic systems.
  3. This field had been active in physics and chemistry before computers even existed, and yet the field is still at its very early stage from the perspective of applied mathematics. In my view, "language barriers" are sometimes the biggest obstacle preventing mathematicians from understanding the problem. There is great potential for numerical analysis and numerical linear algebra tools to deeply impact this field.

Introductory material

With Jianfeng Lu, we developed a summer school course, aimed at introducing electronic structure theory specifically to students and researchers with mathematical backgrounds. In particular, it does not assume that you already know quantum mechanics (we start from scratch and then move very quickly).

The video of the lectures is available on the MSRI website (2.5 hours each day for 10 days).

Additional reading, matlab codes and projects can be found on the LBNL site.


Recent Research Directions


Please see the "Publication" and "Presentation" tabs for more details.



Adaptive compression schemes

The adaptively compressed exchange operator (ACE) formulation reduces the computational cost associated with the Fock exchange operator in Hartree-Fock and hybrid functional Kohn-Sham calculations. The ACE formulation does not depend on the size of the band gap, and thus can be applied to insulating, semiconducting as well as metallic systems, and can reduce the cost of hybrid functional calculations by 5-10 fold without loss of accuracy. In an iterative framework for solving Hartree-Fock-like systems, the ACE formulation only requires moderate modification of the code.

The ACE formulation for Fock exchange calculation is adopted in a number of community electronic structure software packages such as Quantum ESPRESSO and PETot.

The concept of adaptive compression also finds its use in compressing the polarizability operator in phonon calculations. The adaptively compressed polarizability (ACP) formulation reduces the cost of phonon calculations from \(\mathcal{O}(N^4)\) to \(\mathcal{O}(N^3)\).


QRCP based compression schemes

QR factorization with column pivoting (QRCP) is an overlooked and yet powerful technique for electronic structure calculation. We have developed the selected columns of the density matrix (SCDM) formulation for finding localized representation of the subspace spanned by occupied Kohn-Sham orbitals for insulating systems using QRCP. SCDM explicitly uses the fact that density matrices associated with insulating systems decay exponentially along the off-diagonal direction in the real space representation, and avoid the usage of an optimization procedure. SCDM is simple to code (two lines of MATLAB code for demonstration purpose!), and is efficient, robust and massively parallelizable.

The interface of SCDM with Wannier90 package is under progress.


Multiscale methods (embedding)

For mesoscopic systems, the global computational domain is very large. Due to the limit of computational power, only certain degrees of freedoms can be treated accurately with high fidelity quantum mechanical theories. The rest of the domain must be modeled using some lower fidelity theories. In physics, this is referred to as the ``embedding'' problem. We developed the PEXSI-\(\Sigma\) method (see introduction to the PEXSI method below) for coupling quantum mechanical calculations to classical calculations. The PEXSI-\(\Sigma\) method is a discrete analog of the Dirichlet-to-Neumann mapping (DtN) method, where the DtN operator is constructed using a reference Green's function from a crystal configuration.


Spectral density estimation

Methods for estimating the spectral density (a.k.a. density of states) only use matrix-vector multiplication, and can be powerful techniques when diagonalization of the full matrix is prohibitively expensive. They can be useful for computing excited state properties without fully diagonalizing the matrix, or even without having the matrix explicitly constructed. This is the case e.g. when the absorption spectrum is to be computed in linear response time-dependent density functional theory (TDDFT), or when the Bethe-Salpeter equation (BSE) is to be solved.

The convergence of randomized methods can also be accelerated using the spectrum sweeping method, which sweep through the entire spectrum by building such low rank decomposition at different parts of the spectrum.


Fast Green's function methods

The pole expansion and selected inversion (PEXSI) method is an alternative way for solving KSDFT without diagonalizing the Hamiltonian operator. The pole expansion converts the problem of evaluating the electron density to the problem of evaluating the selected elements of Green's functions (i.e. inverse matrices). The PEXSI method can be used for both metallic and insulating systems. PEXSI scales as \(O(N)\) for quasi-1D systems, \(O(N^{1.5})\) for quasi-2D systems, \(O(N^2)\) for 3D bulk systems, and therefore overcomes the \(\mathcal{O}(N^3)\) barrier in standard Kohn-Sham dFT calculations.

The PEXSI software package has been integrated into community electronic structure packages such as SIESTA, CP2K, DGDFT, BigDFT, and Atomistix ToolKit, and can harness more than 100,000 cores on the NERSC system at DOE for systems over 20,000 atoms. It has been used for studying large scale 2D materials such as graphene nanoflake, phospherene nanoflake, and for predicting phospherene nanoflake heterojunction as highly efficient solar cells materials.


Adaptive basis set

The Discontinuous Galerkin Density Functional Theory (DGDFT) uses adaptive local basis functions to discretize the Kohn-Sham Hamiltonian. These basis functions are localized in the real space, and are discontinuous in the global domain. DGDFT can yield both systematically improvable results while keeping a relatively small basis set, with rigorous a posteriori error analysis.

The DGDFT package can scale beyond 100,000 processors and solves electronic structure problem with more than 20,000 atoms. The recent wall clock time per SCF for a graphene system with more than 10,000 atoms is around one minute.