A variational framework for spectral discretization of the density matrix in Kohn Sham density functional theory

Kohn-Sham density functional theory (KSDFT) is currently the main work-horse of quantum mechanical calculations in physics, chemistry, and materials science. From a mechanical engineering perspective, we are interested in studying the role of defects in the mechanical properties in materials. In real materials, defects are typically found at very small concentrations e.g., vacancies occur at parts per million, dislocation density in metals ranges from $10^{10} m^{-2}$ to $10^{15} m^{-2}$, and grain sizes vary from nanometers to micrometers in polycrystalline materials, etc. In order to model materials at realistic defect concentrations using DFT, we would need to work with system sizes beyond millions of atoms. Due to the cubic-scaling computational cost with respect to the number of atoms in conventional DFT implementations, such system sizes are unreachable. Since the early 1990s, there has been a huge interest in developing DFT implementations that have linear-scaling computational cost. A promising approach to achieving linear-scaling cost is to approximate the density matrix in KSDFT. The focus of this thesis is to provide a firm mathematical framework to study the convergence of these approximations. We reformulate the Kohn-Sham density functional theory as a nested variational problem in the density matrix, the electrostatic potential, and a field dual to the electron density. The corresponding functional is linear in the density matrix and thus amenable to spectral representation. Based on this reformulation, we introduce a new approximation scheme, called spectral binning, which does not require smoothing of the occupancy function and thus applies at arbitrarily low temperatures. We proof convergence of the approximate solutions with respect to spectral binning and with respect to an additional spatial discretization of the domain. For a standard one-dimensional benchmark problem, we present numerical experiments for which spectral binning exhibits excellent convergence characteristics and outperforms other linear-scaling methods.

New variational principles for systems of partial differential equations

The question of finding variational principles for coupled systems of first order partial differential equations is considered. Using a potential representation for solutions of the first order system a higher order system is obtained. Existence of a variational principle follows if the original system can be transformed to a self-adjoint higher order system. Existence of variational principles for all linear wave equations with constant coefficients having real dispersion relations is established. The method of adjoining some of the equations of the original system to a suitable Lagrangian function by the method of Lagrange multipliers is used to construct new variational principles for a class of linear systems. The equations used as side conditions must satisfy highly-restrictive integrability conditions. In the more difficult nonlinear case the system of two equations in two independent variables can be analyzed completely. For systems determined by two conservation laws the side condition must be a conservation law in addition to satisfying the integrability conditions.

Part I. On the breaking of nonlinear dispersive waves. Part II. Variational principles in continuum mechanics

A model equation for water waves has been suggested by Whitham to study, qualitatively at least, the different kinds of breaking. This is an integro-differential equation which combines a typical nonlinear convection term with an integral for the dispersive effects and is of independent mathematical interest. For an approximate kernel of the form e^(-b|x|) it is shown first that solitary waves have a maximum height with sharp crests and secondly that waves which are sufficiently asymmetric break into "bores." The second part applies to a wide class of bounded kernels, but the kernel giving the correct dispersion effects of water waves has a square root singularity and the present argument does not go through. Nevertheless the possibility of the two kinds of breaking in such integro-differential equations is demonstrated.

Difficulties arise in finding variational principles for continuum mechanics problems in the Eulerian (field) description. The reason is found to be that continuum equations in the original field variables lack a mathematical "self-adjointness" property which is necessary for Euler equations. This is a feature of the Eulerian description and occurs in non-dissipative problems which have variational principles for their Lagrangian description. To overcome this difficulty a "potential representation" approach is used which consists of transforming to new (Eulerian) variables whose equations are self-adjoint. The transformations to the velocity potential or stream function in fluids or the scaler and vector potentials in electromagnetism often lead to variational principles in this way. As yet no general procedure is available for finding suitable transformations. Existing variational principles for the inviscid fluid equations in the Eulerian description are reviewed and some ideas on the form of the appropriate transformations and Lagrangians for fluid problems are obtained. These ideas are developed in a series of examples which include finding variational principles for Rossby waves and for the internal waves of a stratified fluid.

I. Electronic processes in α-sulfur. II. Polaron motion in a D.C. electric field

Part I: The mobilities of photo-generated electrons and holes in orthorhombic sulfur are determined by drift mobility techniques. At room temperature electron mobilities between 0.4 cm²/V-sec and 4.8 cm²/V-sec and hole mobilities of about 5.0 cm²/V-sec are reported. The temperature dependence of the electron mobility is attributed to a level of traps whose effective depth is about 0.12 eV. This value is further supported by both the voltage dependence of the space-charge-limited, D.C. photocurrents and the photocurrent versus photon energy measurements.

As the field is increased from 10 kV/cm to 30 kV/cm a second mechanism for electron transport becomes appreciable and eventually dominates. Evidence that this is due to impurity band conduction at an appreciably lower mobility (4.10^-4 cm²/V-sec) is presented. No low mobility hole current could be detected. When fields exceeding 30 kV/cm for electron transport and 35 kV/cm for hole transport are applied, avalanche phenomena are observed. The results obtained are consistent with recent energy gap studies in sulfur.

The theory of the transport of photo-generated carriers is modified to include the case of appreciable thermos-regeneration from the traps in one transit time.

Part II: An explicit formula for the electric field E necessary to accelerate an electron to a steady-state velocity v in a polarizable crystal at arbitrary temperature is determined via two methods utilizing Feynman Path Integrals. No approximation is made regarding the magnitude of the velocity or the strength of the field. However, the actual electron-lattice Coulombic interaction is approximated by a distribution of harmonic oscillator potentials. One may be able to find the “best possible” distribution of oscillators using a variational principle, but we have not been able to find the expected criterion. However, our result is relatively insensitive to the actual distribution of oscillators used, and our E-v relationship exhibits the physical behavior expected for the polaron. Threshold fields for ejecting the electron for the polaron state are calculated for several substances using numerical results for a simple oscillator distribution.

I. Quantum mechanics of one-dimensional two-particle models. II. A quantum mechanical treatment of inelastic collisions

Part I

Solutions of Schrödinger’s equation for system of two particles bound in various stationary one-dimensional potential wells and repelling each other with a Coulomb force are obtained by the method of finite differences. The general properties of such systems are worked out in detail for the case of two electrons in an infinite square well. For small well widths (1-10 a.u.) the energy levels lie above those of the noninteresting particle model by as much as a factor of 4, although excitation energies are only half again as great. The analytical form of the solutions is obtained and it is shown that every eigenstate is doubly degenerate due to the “pathological” nature of the one-dimensional Coulomb potential. This degeneracy is verified numerically by the finite-difference method. The properties of the square-well system are compared with those of the free-electron and hard-sphere models; perturbation and variational treatments are also carried out using the hard-sphere Hamiltonian as a zeroth-order approximation. The lowest several finite-difference eigenvalues converge from below with decreasing mesh size to energies below those of the “best” linear variational function consisting of hard-sphere eigenfunctions. The finite-difference solutions in general yield expectation values and matrix elements as accurate as those obtained using the “best” variational function.

The system of two electrons in a parabolic well is also treated by finite differences. In this system it is possible to separate the center-of-mass motion and hence to effect a considerable numerical simplification. It is shown that the pathological one-dimensional Coulomb potential gives rise to doubly degenerate eigenstates for the parabolic well in exactly the same manner as for the infinite square well.

Part II

A general method of treating inelastic collisions quantum mechanically is developed and applied to several one-dimensional models. The formalism is first developed for nonreactive “vibrational” excitations of a bound system by an incident free particle. It is then extended to treat simple exchange reactions of the form A + BC →AB + C. The method consists essentially of finding a set of linearly independent solutions of the Schrödinger equation such that each solution of the set satisfies a distinct, yet arbitrary boundary condition specified in the asymptotic region. These linearly independent solutions are then combined to form a total scattering wavefunction having the correct asymptotic form. The method of finite differences is used to determine the linearly independent functions.

The theory is applied to the impulsive collision of a free particle with a particle bound in (1) an infinite square well and (2) a parabolic well. Calculated transition probabilities agree well with previously obtained values.

Several models for the exchange reaction involving three identical particles are also treated: (1) infinite-square-well potential surface, in which all three particles interact as hard spheres and each two-particle subsystem (i.e. BC and AB) is bound by an attractive infinite-square-well potential; (2) truncated parabolic potential surface, in which the two-particle subsystems are bound by a harmonic oscillator potential which becomes infinite for interparticle separations greater than a certain value; (3) parabolic (untruncated) surface. Although there are no published values with which to compare our reaction probabilities, several independent checks on internal consistency indicate that the results are reliable.

Part I. Approximate Hartree-Fock wavefunctions one-electron properties and electronic structure of the water molecule. Part II. Perturbation-variational calculation of the nuclear spin-spin isotope coupling constant in HD

Part I

Several approximate Hartree-Fock SCF wavefunctions for the ground electronic state of the water molecule have been obtained using an increasing number of multicenter s, p, and d Slater-type atomic orbitals as basis sets. The predicted charge distribution has been extensively tested at each stage by calculating the electric dipole moment, molecular quadrupole moment, diamagnetic shielding, Hellmann-Feynman forces, and electric field gradients at both the hydrogen and the oxygen nuclei. It was found that a carefully optimized minimal basis set suffices to describe the electronic charge distribution adequately except in the vicinity of the oxygen nucleus. Our calculations indicate, for example, that the correct prediction of the field gradient at this nucleus requires a more flexible linear combination of p-orbitals centered on this nucleus than that in the minimal basis set. Theoretical values for the molecular octopole moment components are also reported.

Part II

The perturbation-variational theory of R. M. Pitzer for nuclear spin-spin coupling constants is applied to the HD molecule. The zero-order molecular orbital is described in terms of a single 1s Slater-type basis function centered on each nucleus. The first-order molecular orbital is expressed in terms of these two functions plus one singular basis function each of the types e^-r/r and e^-r ln r centered on one of the nuclei. The new kinds of molecular integrals were evaluated to high accuracy using numerical and analytical means. The value of the HD spin-spin coupling constant calculated with this near-minimal set of basis functions is J_HD = +96.6 cps. This represents an improvement over the previous calculated value of +120 cps obtained without using the logarithmic basis function but is still considerably off in magnitude compared with the experimental measurement of J_HD = +43 0 ± 0.5 cps.