Numerical Methods for Electronic Structure Calculations of Materials

This is the published version. Copyright 2010 Society for Industrial and Applied MathematicsThe goal of this article is to give an overview of numerical problems encountered when determining the electronic structure of materials and the rich variety of techniques used to solve these problems. The paper is intended for a diverse scientific computing audience. For this reason, we assume the reader does not have an extensive background in the related physics. Our overview focuses on the nature of the numerical problems to be solved, their origin, and the methods used to solve the resulting linear algebra or nonlinear optimization problems. It is common knowledge that the behavior of matter at the nanoscale is, in principle, entirely determined by the Schrödinger equation. In practice, this equation in its original form is not tractable. Successful but approximate versions of this equation, which allow one to study nontrivial systems, took about five or six decades to develop. In particular, the last two decades saw a flurry of activity in developing effective software. One of the main practical variants of the Schrödinger equation is based on what is referred to as density functional theory (DFT). The combination of DFT with pseudopotentials allows one to obtain in an efficient way the ground state configuration for many materials. This article will emphasize pseudopotential-density functional theory, but other techniques will be discussed as well

Novel Numerical Models of Electrostatic Interactions and Their Application to S-Nitrosothiol Simulations

Atom-centered point charge model of the molecular electrostatics remains a major workhorse in the atomistic biomolecular simulations. However, this approximation fails to reproduce anisotropic features of the molecular electrostatic potential (MEP), and the existing methods of the charge derivation are often associated with the numerical instabilities. This work provides an in-depth analysis of these limitations and offers a novel approach to describe electrostatic interactions that paves the way toward efficient next-generation force fields. By analyzing the charge fitting problem from first principles, as an example of the mathematical inverse problem, we show that the numerical instabilities of the charge-fitting problem arise due to the decreasing contribution from the higher multipole moments to the overall MEP. This insight suggests that if the point charges are arranged over the sphere using Lebedev quadrature, the resulting point charge model is able to exactly reproduce multipoles up to a given rank. At the same time, point charge values can be derived without fitting to the MEP, avoiding numerically unstable method of the charge derivation. This approach provides a systematic way to reproduce multipole moments up to any rank within the point charge approximation, which makes this model a computationally efficient analog of the multipolar expansion. Moreover, the proposed charged sphere model can be also used in the multi-site expansions with the expansion centers located at each atom in a molecule. This provides a natural approach to expand the traditional atom-centered point charge approximation to include higher-rank atomic multipoles and to account for the anisotropy of the MEP. We applied the proposed charged sphere model to S-nitrosothiols (RSNOs)—a class of biomolecules that serves to store and transmit nitric oxide, a biologically important signaling molecule. We showed that when the atom-centered charged spheres are optimized together with the Lennard-Jones parameters, the resulting force field can accurately reproduce the anisotropic features of the intermolecular interactions that play a crucial role in the biological regulation of RSNO chemistry. Overall, the developed charge model is a promising approach that can be used in the biomolecular simulations and beyond, e.g. in the multipolar force fields for atomistic and coarse-grained simulations

Numerical methods for phase retrieval

In this work we consider the problem of reconstruction of a signal from the magnitude of its Fourier transform, also known as phase retrieval. The problem arises in many areas of astronomy, crystallography, optics, and coherent diffraction imaging (CDI). Our main goal is to develop an efficient reconstruction method based on continuous optimization techniques. Unlike current reconstruction methods, which are based on alternating projections, our approach leads to a much faster and more robust method. However, all previous attempts to employ continuous optimization methods, such as Newton-type algorithms, to the phase retrieval problem failed. In this work we provide an explanation for this failure, and based on this explanation we devise a sufficient condition that allows development of new reconstruction methods---approximately known Fourier phase. We demonstrate that a rough (up to $\pi/2$ radians) Fourier phase estimate practically guarantees successful reconstruction by any reasonable method. We also present a new reconstruction method whose reconstruction time is orders of magnitude faster than that of the current method-of-choice in phase retrieval---Hybrid Input-Output (HIO). Moreover, our method is capable of successful reconstruction even in the situations where HIO is known to fail. We also extended our method to other applications: Fourier domain holography, and interferometry. Additionally we developed a new sparsity-based method for sub-wavelength CDI. Using this method we demonstrated experimental resolution exceeding several times the physical limit imposed by the diffraction light properties (so called diffraction limit).Comment: PhD. Thesi

Simulations of Chemical Catalysis

This dissertation contains simulations of chemical catalysis in both biological and heterogeneous contexts. A mixture of classical, quantum, and hybrid techniques are applied to explore the energy profiles and compare possible chemical mechanisms both within the context of human and bacterial enzymes, as well as exploring surface reactions on a metal catalyst. A brief summary of each project follows. Project 1 — Bacterial Enzyme SpvC The newly discovered SpvC effector protein from Salmonella typhimurium interferes with the host immune response by dephosphorylating mitogen-activated protein kinases (MAPKs) with a -elimination mechanism. The dynamics of the enzyme substrate complex of the SpvC effector is investigated with a 3.2 ns molecular dynamics simulation, which reveals that the phosphorylated peptide substrate is tightly held in the active site by a hydrogen bond network and the lysine general base is positioned for the abstraction of the alpha hydrogen. The catalysis is further modeled with density functional theory (DFT) in a truncated active-site model at the B3LYP/6-31 G(d,p) level of theory. The truncated model suggested the reaction proceeds via a single transition state. After including the enzyme environment in ab initio QM/MM studies, it was found to proceed via an E1cB-like pathway, in which the carbanion intermediate is stabilized by an enzyme oxyanion hole provided by Lys104 and Tyr158 of SpvC. Project 2 — Human Enzyme CDK2 Phosphorylation reactions catalyzed by kinases and phosphatases play an indispensable role in cellular signaling, and their malfunctioning is implicated in many diseases. Ab initio quantum mechanical/molecular mechanical studies are reported for the phosphoryl transfer reaction catalyzed by a cyclin-dependent kinase, CDK2. Our results suggest that an active-site Asp residue, rather than ATP as previously proposed, serves as the general base to activate the Ser nucleophile. The corresponding transition state features a dissociative, metaphosphate-like structure, stabilized by the Mg(II) ion and several hydrogen bonds. The calculated free-energy barrier is consistent with experimental values. Project 3 — Bacterial Enzyme Anthrax Lethal Factor In this dissertation, we report a hybrid quantum mechanical and molecular mechanical study of the catalysis of anthrax lethal factor, an important first step in designing inhibitors to help treat this powerful bacterial toxin. The calculations suggest that the zinc peptidase uses the same general base-general acid mechanism as in thermolysin and carboxypeptidase A, in which a zinc-bound water is activated by Glu687 to nucleophilically attack the scissile carbonyl carbon in the substrate. The catalysis is aided by an oxyanion hole formed by the zinc ion and the side chain of Tyr728, which provide stabilization for the fractionally charged carbonyl oxygen. Project 4 — Methanol Steam Reforming on PdZn alloy Recent experiments suggested that PdZn alloy on ZnO support is a very active and selective catalyst for methanol steam reforming (MSR). Plane-wave density functional theory calculations were carried out on the initial steps of MSR on both PdZn and ZnO surfaces. Our calculations indicate that the dissociation of both methanol and water is highly activated on \ufb02at surfaces of PdZn such as (111) and (100), while the dissociation barriers can be lowered significantly by surface defects, represented here by the (221), (110), and (321) faces of PdZn. The corresponding processes on the polar Zn-terminated ZnO(0001) surfaces are found to have low or null barriers. Implications of these results for both MSR and low temperature mechanisms are discussed

A Survey of Stochastic Simulation and Optimization Methods in Signal Processing

International audienceModern signal processing (SP) methods rely very heavily on probability and statistics to solve challenging SP problems. SP methods are now expected to deal with ever more complex models, requiring ever more sophisticated computational inference techniques. This has driven the development of statistical SP methods based on stochastic simulation and optimization. Stochastic simulation and optimization algorithms are computationally intensive tools for performing statistical inference in models that are anal ytically intractable and beyond the scope of deterministic inference methods. They have been recently successfully applied to many difficult problems involving complex statistical models and sophisticated (often Bayesian) statistical inference techniques. This survey paper offers an introduction to stochastic simulation and optimization methods in signal and image processing. The paper addresses a variety of high-dimensional Markov chain Monte Carlo (MCMC) methods as well as deterministic surrogate methods, such as variational Bayes, the Bethe approach, belief and expectation propagation and approximate message passing algorithms. It also discusses a range of optimization methods that have been adopted to solve stochastic problems, as well as stochastic methods for deterministic optimization. Subsequently, area as of overlap between simulation and optimization, in particular optimization-within-MCMC and MCMC-driven optimization are discussed

Towards the Understanding of Fluorescence Quenching Mechanisms : Molecular Dynamics Simulations of Dye-Quencher Interactions in Biomolecular Systems

The presence of antibodies directed against p53 in human blood serum is a specific and independent marker of cancer. The present thesis is a simulative study of two fluorescent dyes (rhodamine 6G and MR121) used in the development of fluorescence-based immunoassays for the detection of p53 antibodies. The "selective" fluorescence quenching property of tryptophan residues present in dye-conjugated peptide chains enables monitoring of conformational dynamics and antibody binding events by means of fluorescence spectroscopy. The molecular mechanisms of the quenching of fluorescent dyes are mostly unknown. Here MD simulations are used in combination with existing results from ensemble fluorescence experiments in order to obtain predictive theoretical insight into dye/quencher interactions. A new automated refinement method was developed for deriving reliable molecular mechanics force field parameters for small- to medium-sized molecules using reference data from high level ab initio quantum chemical calculations. Using this method CHARMM force field parameters for the two dyes were derived. The parameters were then used to perform MD simulations on two simplified, but realistic dye quencher systems: MR121/TRP and R6G/TRP. Results of these simulations have given important insight on the bimolecular interactions between the dyes and the quencher. A quantitative interpretation of the quenching mechanism based on it’s dependence on dye/quencher distance and orientation has emerged. These results were then applied to the interpretation of simulation data of an epitope from the tumor suppressor protein p53 which was labelled at the N-terminus first with MR121 and then with R6G. In all phases of the thesis care was taken to confront and/or combine the theoretical results with available experimental data

Point-set manifold processing for computational mechanics: thin shells, reduced order modeling, cell motility and molecular conformations

In many applications, one would like to perform calculations on smooth manifolds of dimension d embedded in a high-dimensional space of dimension D. Often, a continuous description of such manifold is not known, and instead it is sampled by a set of scattered points in high dimensions. This poses a serious challenge. In this thesis, we approximate the point-set manifold as an overlapping set of smooth parametric descriptions, whose geometric structure is revealed by statistical learning methods, and then parametrized by meshfree methods. This approach avoids any global parameterization, and hence is applicable to manifolds of any genus and complex geometry. It combines four ingredients: (1) partitioning of the point set into subregions of trivial topology, (2) the automatic detection of the local geometric structure of the manifold by nonlinear dimensionality reduction techniques, (3) the local parameterization of the manifold using smooth meshfree (here local maximum-entropy) approximants, and (4) patching together the local representations by means of a partition of unity. In this thesis we show the generality, flexibility, and accuracy of the method in four different problems. First, we exercise it in the context of Kirchhoff-Love thin shells, (d=2, D=3). We test our methodology against classical linear and non linear benchmarks in thin-shell analysis, and highlight its ability to handle point-set surfaces of complex topology and geometry. We then tackle problems of much higher dimensionality. We perform reduced order modeling in the context of finite deformation elastodynamics, considering a nonlinear reduced configuration space, in contrast with classical linear approaches based on Principal Component Analysis (d=2, D=10000's). We further quantitatively unveil the geometric structure of the motility strategy of a family of micro-organisms called Euglenids from experimental videos (d=1, D~30000's). Finally, in the context of enhanced sampling in molecular dynamics, we automatically construct collective variables for the molecular conformational dynamics (d=1...6, D~30,1000's)

Protein Structure Elastic Network Models and the Rank 3 Positive Semidefinite Matrix Manifold

This thesis is a contribution to the study of protein dynamics using elastic network models (ENMs). An ENM is an abstraction of a protein structure where inter-atomic interactions are assumed to be modelled by a Hookean potential energy which is a function if inter-atomic distances. This model has been studied by various authors, and despite being a very simple model, can nonetheless provide a realistic understanding of protein dynamics. For example, it was shown by Tirion, that the Hookean potential energy can reproduce the normal mode fluctuations of the more complicated semi-empirical potential. In addition, it was shown by Tekpinar and Zheng that an ENM can correctly model the order of local conformational changes before global conformational changes during ATP-driven conformational changes. The purpose of this thesis is to provide a second mathematical formulation for modelling ENMs. This thesis suggests removing the square-root in the Hookean potential which leads to a positive semidefinite (PSD) potential that is a function of quadrances rather than distances. There are many similarities between the two approaches, but also many differences. One main difference is PSD matrices are linearly related to quadrance, the square of distance, which opens the way to model the PSD potential using perceptrons whose weight matrix is a rank 3 PSD matrix. This interesting consequence is left as a topic of future research. The PSD potential is just as appropriate for modelling ENMs as observed by the following two agreements: The PSD potential produces normal mode fluctuations that agree with the Hookean potential introduced by Tirion. This agreement suggests both potentials provide the same information about a protein structure's flexibility. The generalization of the Hookean iENM potential (introduced by Tekpinar and Zheng) to the PSD iENM potential also interpolates the local conformational changes before the global conformational changes, in agreement with the original Hookean observations. Recall that the equations of motion in classical mechanics is formulated using an abstract Riemannian manifold. This abstraction gives modellers the flexibility to consider different Riemannian manifolds appropriate to the problem. After the introduction of the Hookean potential, the study of protein dynamics still uses the 3n dimensional Euclidean space as the Riemannian manifold, the same Riemannian manifold used by the semi-empirical potential. This is because both the semi-empirical potential and the Hookean potential assume the atomic coordinates of a protein structure are represented by a 3n by 1 vector. However, with the introduction of the PSD potential, the protein structure's atomic coordinates are represented as a point on the rank 3 n by n PSD matrix manifold. Consequently, a new Riemannian manifold for modelling protein dynamics has been proposed. In order to model protein dynamics on the rank 3 PSD matrix manifold, the equations of motion needs to be defined. This thesis presents the geometric objects: horizontal projection, gradient, Hessian, and retraction required for formulating the equations of motion for protein structures as an optimization problem on the rank 3 PSD matrix manifold. These formulas are a modification of the original formulas introduced by Journée et al. to allow constraints relevant to a protein structure to be described. Rosen's correction to the constraint manifold was already introduced in 1961, and was reintroduced by Goldenthal et al. in 2007 under the name of ``the fast projection algorithm''. Rosen's, Goldenthal et al.'s, and Journée et al.'s work are all closely related but were developed independently. This thesis makes their relationship more apparent