We describe an efficient implementation of an algorithm for computing selected elements of a general sparse symmetric matrix A that can be decomposed as A = LDL^T, where L is lower triangular and D is diagonal. Our implementation, which is called SelInv, is built on top of an efficient supernodal left-looking LDL^T factorization of A. We discuss how computational efficiency can be gained by making use of a relative index array to handle indirect addressing. We report the performance of SelInv on a collection of sparse matrices of various sizes and nonzero structures. We also demonstrate how SelInv can be used in electronic structure calculations.

We present an efficient parallel algorithm and its implementation for computing the diagonal of $H^-1$ where $H$ is a 2D Kohn-Sham Hamiltonian discretized on a rectangular domain using a standard second order finite difference scheme. This type of calculation can be used to obtain an accurate approximation to the diagonal of a Fermi-Dirac function of $H$ through a recently developed pole-expansion technique \cite{LinLuYingE2009}. The diagonal elements are needed in electronic structure calculations for quantum mechanical systems \citeHohenbergKohn1964, KohnSham 1965,DreizlerGross1990. We show how elimination tree is used to organize the parallel computation and how synchronization overhead is reduced by passing data level by level along this tree using the technique of local buffers and relative indices. We analyze the performance of our implementation by examining its load balance and communication overhead. We show that our implementation exhibits an excellent weak scaling on a large-scale high performance distributed parallel machine. When compared with standard approach for evaluating the diagonal a Fermi-Dirac function of a Kohn-Sham Hamiltonian associated a 2D electron quantum dot, the new pole-expansion technique that uses our algorithm to compute the diagonal of $(H-z_i I)^-1$ for a small number of poles $z_i$ is much faster, especially when the quantum dot contains many electrons.

Kohn-Sham density functional theory (KSDFT) is by far the most widely used electronic structure theory in condensed matter systems. Density functional perturbation theory (DFPT) studies the response of a quantum system under small perturbation, where the quantum system is described at the level of first principle electronic structure theories like KSDFT. One important application of DFPT is the calculation of vibration properties such as phonons, which can be further used to calculate many physical properties such as infrared spectroscopy, elastic neutron scattering, specific heat, heat conduction, and electron-phonon interaction related behaviors such as superconductivity . DFPT describes vibration proper- ties through a polarizability operator, which characterizes the linear response of the electron density with respect to the perturbation of the external potential. More specifically, in vibration calculations, the polarizability operator needs to be applied to d × NA ∼ O(Ne) perturbation vectors, where d is the spatial dimension (usually d = 3), NA is the number of atoms, and Ne is the number of electrons. In general the complexity for solving KSDFT is O(Ne3), while the complexity for solving DFPT is O(Ne4). It is possible to reduce the computational complexity of DFPT calculations by “linear scaling methods”. Such methods can be successful in reducing the computational cost for systems of large sizes with substantial band gaps, but this can be challenging for medium-sized systems with relatively small band gaps.

In the discussion below, we will slightly abuse the term “phonon calculation” to refer to calculation of vibration properties of condensed matter systems as well as isolated molecules. In order to apply the polarizability operator to O(Ne) vectors, we need to solve O(Ne2) coupled Sternheimer equations. On the other hand, when a constant number of degrees of freedom per electron is used, the size of the Hamiltonian matrix is only O(Ne). Hence asymptotically there is room to obtain a set of only O(Ne) “compressed perturbation vectors”, which encodes essentially all the information of the O(Ne2) Sternheimer equations. In this dissertation, we develop a new method called adaptively compressed polarizability operator (ACP) formulation, which successfully reduces the computational complexity of phonon

1

2

calculations to O(Ne3) for the first time. The ACP formulation does not rely on exponential decay properties of the density matrix as in linear scaling methods, and its accuracy depends weakly on the size of the band gap. Hence the method can be used for phonon calculations of both insulators and semiconductors with small gaps.

There are three key ingredients of the ACP formulation. 1) The Sternheimer equations are equations for shifted Hamiltonians, where each shift corresponds to an energy level of an occupied band. Hence for a general right hand side vector, there are Ne possible energies (shifts). We use a Chebyshev interpolation procedure to disentangle such energy dependence so that there are only constant number of shifts that is independent of Ne. 2) We disentangle the O(Ne2) right hand side vectors using the recently developed interpolative separable density fitting procedure, to compress the right-hand-side vectors. 3) We construct the polarizability by adaptive compression so that the operator remains low rank as well as accurate when applying to a certain set of vectors. This make it possible for fast computation of the matrix inversion using methods like Sherman-Morrison-Woodbury.

In particular, the new method does not employ the “nearsightedness” property of electrons for insulating systems with substantial band gaps as in linear scaling methods. Hence our method can be applied to insulators as well as semiconductors with small band gaps.

This dissertation also extend the method to deal with nonlocal pseudopotentials as well as systems in finite temperature. Combining all these components together, we obtain an accurate, efficient method to compute the vibrational properties for insulating and metallic systems.

This dissertation concerns the quantum many-body problem, which is the problem of predicting the properties of systems of several quantum particles from the first principles of quantum mechanics. Included under this umbrella are various problems of fundamental importance in quantum chemistry, condensed matter physics, and materials science. Of particular interest is the electronic structure problem, the problem of determining the state of the electrons in a system with fixed atomic nuclei. Since direct numerical solution of the many-body Schrödinger equation is intractable even for systems of moderate size, a diverse array of approximate methods has been developed. The broad goals of this dissertation are to improve the mathematical understanding of certain widely-used approximations, as well as to propose new methods. Roughly speaking, we consider three (overlapping) categories of methods: Green's function methods, embedding methods, and variational methods.

One can understand Green's function methods in terms of many-body perturbation theory, which computes series expansions of physical quantities about a non-interacting reference system. These expansions can be expressed graphically in terms of Feynman diagrams, which can in turn be reorganized, in some cases, into an expansion in terms of so-called bold diagrams. Green's function methods can be specified by choosing a subset of bold diagrams to approximate the sum. At the same time, such methods can be understood in terms of an object known as the Luttinger-Ward (LW) functional, which admits a representation in terms of the bold diagrams. Many aspects of these constructions are purely formal, and indeed the existence of the fermionic LW functional as a single-valued functional has recently been called into question. To contribute to the understanding of these issues, we provide rigorous proofs of the combinatorial construction and analytic interpretation of the bold diagrams in the simplified setting of a classical field theory. In this setting we also provide a rigorous non-perturbative construction of the LW functional via convex duality and prove several key properties, including continuity up to the boundary of its domain and asymptotics in the limit of large interaction.

Quantum embedding methods, meanwhile, view a large system as being composed of smaller fragments that are treated with high accuracy and embedded in the larger system in a mutually consistent way. Inspired by a connection between the boundary analysis of the LW functional and embedding, we perform similar analysis for the 1-RDM theory for fermionic systems, which is also developed via convex duality, illustrating a relation to fermionic embedding methods such as the density matrix embedding theory (DMET).

Another embedding method of note is the dynamical mean-field theory (DMFT), which is at the same time a Green's function method that can be understood in terms of the LW functional. DMFT relies on the solution of impurity problems, which specify the embedding of an interacting system into a non-interacting bath. Underlying DMFT is a result about the sparsity pattern of the self-energy matrix for impurity problems, which to our knowledge has not been proved in the literature. We provide a rigorous proof of this result in various classical and quantum settings. We go on to investigate the fermionic DMFT in depth, identifying the key mathematical structures that appear in the algorithmic loop for solving it and using these to prove the well-posedness of this loop, in a certain sense.

Finally, we introduce a suite of new approaches to the quantum many-body problem that provide variational lower bounds to the ground-state energy. These methods, which combine the themes of convexity and embedding, are based on novel convex relaxations of the variational principles for the ground-state energies of many-body systems. To begin, we recover a second-quantized version of the formalism of strictly correlated electrons (SCE), which yields an exact expression for the exchange-correlation functional in Kohn-Sham density functional theory in the limit of infinite Coulomb repulsion in terms of the solution of a multi-marginal optimal transport problem. We introduce a semidefinite relaxation method for approximately solving this problem and obtaining a lower bound for the ground-state energy. The ideas underlying this relaxation are generalized considerably, outside the context of SCE, to yield much tighter lower bounds, which we validate numerically for both quantum spin systems and fermionic systems. We also describe how these relaxation methods can be interpreted as embedding methods via convex duality.

This work is concerned with the accurate numerical simulation of the many-electron problem, which involves the modeling of the electron wavefunction, from which all of the properties of a chemical or condensed matter system can, in principle, be computed. This problem poses a number of challenges, including the effective parametrization of the wavefunction space. The combination of neural networks and quantum Monte Carlo methods has arisen as a promising path forward for highly accurate electronic structure calculations. Previous proposals have combined equivariant neural network layers with a final antisymmetric layer in order to satisfy the antisymmetry requirements of the electronic wavefunction. However, to date it is unclear if one can represent arbitrary antisymmetric functions of physical interest, and it is difficult to precisely measure the expressiveness of the antisymmetric layer.

In the first chapter of this dissertation, we begin by introducing the electronic structure problem and the variational nature of finding the lowest energy wavefunction, or ground state. We describe Metropolis Monte Carlo sampling techniques, as well as a simplifying reduction to the number of degrees of freedom. We conclude the first chapter with a brief discussion of some optimization techniques used to address the variational ground state problem once a trial parametrization has been established.

In the next chapter, we then introduce the form of some modern neural-network based trial wavefunctions. We attempt to investigate the expressiveness of the antisymmetric layers by proposing explicitly antisymmetrized universal neural network layers. This approach has a computational cost which increases factorially with respect to the system size, but we are nonetheless able to apply it to small systems to better understand how the structure of the antisymmetric layer affects its performance. We first introduce a generic antisymmetric (GA) neural network layer, which we use to replace the entire antisymmetric layer of the highly accurate ansatz known as the FermiNet. We also consider a factorized antisymmetric (FA) layer which more directly generalizes the FermiNet by replacing the products of determinants with products of antisymmetrized neural networks.

We next investigate the numerical performance of these explicitly antisymmetrized ansatzes. We demonstrate that the FermiNet-GA architecture can yield effectively the exact ground state energy for small atoms and molecules. We find, interestingly, that the resulting FermiNet-FA architecture does not outperform the FermiNet. This strongly suggests that the sum of products of antisymmetries is a key limiting aspect of the FermiNet architecture. To explore this further, we also investigate a slight modification of the FermiNet, called the full determinant mode, which replaces each product of determinants with a single combined determinant. We find that the full single-determinant FermiNet closes a large part of the gap between the standard single-determinant FermiNet and FermiNet-GA on small atomic and molecular problems. Surprisingly, on the nitrogen molecule at a dissociating bond length of 4.0 Bohr, the full single-determinant FermiNet can significantly outperform the largest standard FermiNet calculation with 64 determinants, yielding an energy within 0.4 kcal/mol of the best available computational benchmark.

In the final chapter, we introduce the VMCNet repository, which was used to implement the numerical experiments previously described. VMCNet is intended to be a flexible, general purpose VMC framework which interfaces natively with the JAX library for rapid prototyping, performance benefits due to just-in-time XLA compilation, and easy dispatch to multiple GPU systems. We describe both the Python API, intended for more complex use-cases that require customization of finer details of the VMC loop, and the command-line interface, which provides a more streamlined and encapsulated way to run variational Monte Carlo experiments, including those described in this dissertation.