SciDAC Review - Speed and Precision in Quantum Chemistry

COMPUTATIONAL CHEMISTRY

SPEED and PRECISION in QUANTUM CHEMISTRY

Figure 1. A molecular model for amorphous polyisobutylene. Accurate calculation of Raman spectra would provide additional validation for this molecular model against experimental data. However, this demands high-quality basis sets which are not practical in conventional schemes for molecules of this size.

MADNESS is a simulation framework powered by innovative mathematical tools.With a combination of guaranteed speed and precision unprecedented in computational chemistry, MADNESS is capable of solving quantum chemical problems that are too complicated for conventional methods. It will become even more powerful as versions are developed for use on the petascale computers.

The field of quantum chemistry—like so many areas of science—has made extraordinary progress due to advanced computing and simulation technology. “The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble.” This statement was made by Paul Dirac, who later became the 1933 Nobel Laureate in Physics. Dirac’s statement rings true today, as equations that represent chemical and physical phenomena are still very complicated and difficult to solve. Yet, tremendous gains in theory, algorithm, and computational power over the course of eight decades have provided a number of simulation tools capable of answering some chemical questions. When applied to small molecules these methods are completely quantitative, but with larger molecules (figure 1) the answers become less precise. Such simulation tools are constructive with careful application, such as for use in prediction of trends, rather than determination of absolute values. Despite impressive progress to date, more powerful simulation tools are needed to usher in the next era of scientific discovery.

Further advances require the coupling of new theories and advanced algorithms with the very largest, most powerful computers—a goal ubiquitous across the SciDAC program.These challenges have motivated a collaborative effort led by Dr. George I. Fann and Dr. Robert J. Harrison of ORNL, and Dr. Gregory Beylkin of University of Colorado–Boulder. Dr. Fann is the Principal Investigator on the recently awarded (see “SciDAC-2: The Next Phase of Discovery,” p16) SciDAC Scientific Application Partnership (SAP), “Advanced Mathematics for Electronic Structure,” a project tightly coupled to “Advanced Methods for Electronic Structure,” a SciDAC project that focused on advancing quantum chemical methods. Dr. Harrison served as PI on the latter and is currently a co-investigator on the former.Dr. Fann and Dr. Harrison are working with a diverse group of chemists, mathematicians, and computer scientists to lead computational chemistry to the petascale. From this collaborative network of scientists emerges Multiresolution ADap- tive NumErical Scientific Simulation (MADNESS), an advanced computational quantum chemistry software package. With the new mathematical tools comprising MADNESS, complex quantum chemical problems in many dimensions can be computed with guaranteed speed and precision.

MADNESS serves as a powerful scientific simulation framework, applicable to a number of important chemical problems. Potential applications include clean energy innovations for automobiles and industry, improvements to the efficiency of drug development in the pharmaceutical field, and analysis of the dynamics of molecules and electrons in powerful laser fields. Although the development of MADNESS is currently motivated primarily by questions in chemistry, it is a general purpose framework with expected extensions into many disciplines and interdisciplinary applications. Furthermore, the mathematical tools underlying the software may eventually be applied broadly, paving the way to discovery in many fields of science.

Experiment, Simulation, and Theory
Quantum chemistry is concerned with calculating the electronic structures and energies of atoms and molecules. Chemical behavior—how atoms and molecules interact—can be predicted with great confidence from this information (sidebar "Quantum Chemistry Basics"). To gather the information critical to understanding such systems, quantum chemists rely on a combination of experiment, simulation, and theory.

Modern experimental techniques such as femto-second laser spectroscopy, nuclear magnetic resonance, and atomic force microscopy are very powerful and can reveal some of the necessary information. Despite these tremendous capabilities, many important chemical questions remain beyond the reach of experimental study. For example, the rational design of new catalysts requires a detailed understanding of the reaction mechanism and, in particular, the rate limiting step, but these ephemeral structures are often too short-lived to be observed. Moreover, many chemical transformations occur at interfaces or inside nanostructured pores and surfaces, and are thus hidden from most experimental probes. “Experimental chemistry is left making indirect inferences about what is actually happening,” explains Dr. Harrison. However, missing information is often available from computer simulations. Indeed many insights—such as the detailed nature of a given chemical bond or the interpretation of X-ray spectra—are only accessible from theory. “Thus modern chemistry has emerged as a close partnership between experiment and computer simulation, underpinned by rigorous theory.”

Costs of Simulation
Much interesting and economically important chemistry actually occurs in very poorly characterized systems, and such systems require advanced methods. Unfortunately, the capabilities of simulation are severely constrained by both the construction of molecular structures with full atomic detail and the computational cost of running fully quantitative models.

Quantum chemists generally use mathematical methods executed by computer programs to determine the energies and structures of molecules. These programs employ a wide variety of methods to approximate wave functions, which are the solutions to complex wave equations. Nearly all of these methods start from effective one-electron models in which electrons “move” independently, experiencing only an average potential due to the influence of other electrons, nuclei, and external fields. These one-electron wave functions are termed molecular orbitals, and since molecules are comprised of atoms, a very effective technique has been to expand the molecular orbitals in terms of atom-centered atomic orbitals (sidebar "Atomic Orbital Methods"). Indeed, this linear combination of atomic orbitals (LCAO) approach has been responsible for many of the successes of the first fifty years of quantum chemistry. However, applying this to large molecular systems requires calculating the contributions of the many associated functions. Unfortunately, these calculations are both computationally intensive and numerically erratic.

Any question that can be asked about a stationary molecule can be answered by the molecular wave function, which depends upon information about all the electrons, such as spatial and spin coordinates. The wave function is found by solving the wave equations fundamental to quantum chemistry—the non-relativistic Schrödinger equation, and the Dirac equation for relativistic problems. Yet, when applied to real molecules, these equations are too mathematically complicated for exact solution, and so must instead be approximated. Much of computational quantum chemistry involves the development of methods that provide accurate approximations of the wave function at a practical cost, meaning they must be computationally tractable. But accuracy is very expensive in terms of computational cost, so nearly all studies require a careful tradeoff, balancing computational feasibility against the precision required to answer specific chemical questions.

Scaling Behavior
Ideally, the cost of calculations is expected to increase only linearly, O(N), as the size of the system increases. This assumption is based upon the observation that electron correlation is local, meaning that long-range interaction between electrons is effectively screened. However, as larger molecular systems or more accurate chemical simulations are pursued, the complexity of the involved equations grows in a nonlinear manner, such as quadratically O(N²), cubically O(N³), quartically O(N⁴), and greater (see figure 10, p63). This rise in complexity that accompanies increases in either system size or required numerical accuracy renders the concept of scaling behavior very important to computational chemistry researchers.

Scaling behavior is a measure of how the required computational effort climbs with increases in the size or complexity of the system being simulated. Problems in quantum chemistry are notorious for being extremely computationally intensive. For example, the CCSD(T) method—a “gold standard” of quantum chemistry—is known from both theory and extensive benchmark calculations to predict with chemical precision the properties of most molecules. However, the computational cost of this method scales as the seventh power of the system size O(N7) and the inverse fourth power of the requested precision (e-4). Thus, doubling the size of the system increases the computational cost by a factor of 128, and requesting just another decimal place of precision multiplies this by a factor of 10,000. This kind of scaling behavior renders studies impractical unless limited to small systems. Improving scaling behavior of such methods is a priority, and efforts to do so are always ongoing.

Methods of Approximation
More widely applicable methods, such as Hartree–Fock (HF) and density functional theory (DFT), are not capable of yielding the exact result, but are nevertheless capable of making quantitative predictions (sidebar “Hartree–Fock and Density Functional Theory”). Formally, the computational cost of these methods scale as the fourth power of the system size [quartically, O(N4)] or greater, but reduced scaling methods have been available for some time. Achieving reduced scaling behavior often means accepting a low-level of precision. Practically this allows fully predictive simulations for small systems, but computational feasibility forces conventional quantum chemistry studies to make increasingly severe and uncontrolled approximations in order to simulate larger systems. The result is an undesirable lack of precision.

Figure 5. The illustrations above (V₀, V₁, and V₂) show grids of varying degrees of detail. The mathematical notation describes how the grids are related in a telescoping series of spaces. The space called V0 is the coarsest grid, and would be used in smooth areas where information can be ignored without affecting overall accuracy very much. Finer grids, such as V₁, V₂….V_n, are necessary to represent more detailed areas, where information cannot be neglected without sacrificing accuracy.

Improving the scaling behavior of quantum chemistry calculations without sacrificing speed and accuracy is a primary motivation for the development of MADNESS. “It seems counterintuitive that we can improve both precision and speed at the same time,” remarks Dr. Harrison, “but it should be remembered that the methods of mainstream chemistry originated in the mid-1950s when computers were still literally graduate students chained to their hand-cranked calculators for the first few hours of each day.” At the time, a compact representation of molecular wave functions was crucial; expanding molecular wave functions in terms of atomic wave functions is both compact and physically advantageous (sidebar "Atomic Orbital Methods," p57). The problems arise when higher precision is pursued, especially in large systems that require the addition of more functions for each atom. Because the functions on each atom overlap with each other, the atom-centered expansions become numerically ill-conditioned as more functions are added. Moreover, the overlap of diffuse functions with atoms a long distance away increases the computational cost very substantially. In contrast, fully numerical approaches have been found to be more robust at high precision. “Application of such approaches has so far been limited by speed, the ability to reliably trade speed for precision, and a path to computing in beyond three dimensions,” explains Dr. Harrison. “This is where MADNESS comes into the picture.”

The Math Behind MADNESS
The MADNESS computational scheme combines a number of specialized mathematical methods in order to achieve its impressive performance.

The “M” in MADNESS stands for “Multiresolution,” and the “ADNE” is for “ADaptive NumErical” methods. In Scientific Simulations—the “SS” in the MADNESS acronym—physical objects of interest must be represented mathematically.

Mathematical representations can be constructed with varying degrees of detail. A representation will be more accurate when it contains more information about detail, but the added information also demands greater computational times. The complex geometries of physical objects—molecules, in the case of quantum chemistry—means that the optimal mathematical representation might be more detailed for some areas of the object and less detailed in other areas.

Figure 6. In the first panel is a visualization of a molecular orbital of the benzene dimer computed with the local density approximation using MADNESS. In the next two panels, the same molecular orbital of the benzene dimer is shown from two different perspectives and with the adaptive mesh displayed.

Multiresolution analysis is used in numerous applications, including image processing and solving differential equations. The central idea of multiresolution analysis is perhaps most succinctly conveyed by what mathematicians call a telescoping series of representations. In figure 5, the space named Vn is the most accurate and informative basis, or finest grid, of the telescoping series and V0 is the coarsest, least accurate grid of the series. When the physical object is finely detailed everywhere, its representation cannot be sparse anywhere, meaning no information can be neglected without sacrificing accuracy. In this case, the finest grid must be used by the method. However, if some regions are smooth with little variation, these can be represented accurately on a coarser mesh at some level of refinement because differences at finer scale lengths will be too small to matter very much. Adaptive methods discard differences smaller than some threshold and generate a sparse representation which fits, or adapts to, the object being represented. Essentially, this is done by retaining fine scale information only where it is necessary. Eliminating unnecessary information lowers complexity, which in turn lowers the computational cost.

Setting the threshold value controls the error introduced by technique. A smaller threshold will keep more information, meaning a more accurate calculation, but also a more computationally demanding simulation. A larger threshold setting reduces the computational cost by lowering the information content, and thus the precision. Thus, these methods allow scientists to trade speed for precision, and vice versa.

Multiresolution adaptive numerical methods allow for variable representation of areas that vary at different scales, with more information included for fine scale details and less information at areas of coarser scale. Non-adaptive numerical methods are not as computationally efficient because the entire representation would include detail at the finest scale, even in areas where so much information was not necessary. Hierarchical representations from multiresolution analysis provide MADNESS with capabilities for fast computation, and the separated representations of operators (sidebar "Separated Representations") make these methods efficient in three, six, and higher dimensions. The algorithm utilized in MADNESS is called a discontinuous spectral element method, and is implemented to provide a sampling scheme for adaptive grid refinement methods. The mathematical workings of these methods are what allow MADNESS to complete powerful calculations with speed and guaranteed precision.

Code Logistics
Another parameter affecting the performance of conventional quantum chemistry calculation methods is the large quantity of computer code required. Mathematical models of physical phenomena, such as the Schrödinger equation (figure 3), can often be written out in relatively few lines of math notation. However, computer programs for the solution of such equations can contain many thousands of lines of code. This disparity between mathematical notation and computer code is what another SciDAC PI, Dr. Phil Colella, articulately calls the “semantic gap.” One example is the widely available computational chemistry software package NWChem, which contains over one million human-generated lines and several million machine-generated lines of code.

Developing, debugging, optimizing, and maintaining such unwieldy code is a tremendous logistical exercise. Calculations must be validated numerically and portability of the code between different platforms and operating systems must be ensured. These tasks are quite formidable with very large codes. The MADNESS collaborators have kept this in mind and made these logistical concerns a high priority.

Parallel Computing
In today’s era of supercomputing, another particularly important aspect of code development is optimization of the division of computational effort between each of the microprocessors in large parallel computing environments. These tasks are usually beyond the capabilities of a single person and require specialized teams of computer programmers, further adding to the logistical cost of scientific simulations.

Figure 8. A substituted cobalt-porphyrin has been shown to catalyze the insertion of nitrogen and oxygen into carbon–carbon double bonds to form three-membered rings. The most striking feature of this process is that it preserves with high yield the relative position of the substituents enabling the synthesis of chiral species with high efficiency. Chemical catalysis is of immense economic importance to the US chemical and petrochemical industries, and is of direct interest to the DOE due to its role in efficient green or low-emission energy systems such as hydrogen fuel cells.

Internally, MADNESS has adopted the message passing interface (MPI) standard for parallel programming because it is the most portable, especially with the very latest computers. “But MPI is a very low-level programming tool and very cumbersome to use for the dynamic computation within MADNESS,” says Dr. Harrison.“ Thus, we have layered on top of MPI some software so that each processor has an independent queue of tasks to perform, and it can add tasks to either its own queue or that of another processor.” In addition to introducing additional parallelism—to utilize multicore chips, for example—users now have the freedom to move work to the data in cases where it is more efficient than moving the data to the work. Furthermore, because the task scheduling looks after dependencies between the tasks, the application code has become much easier to write. “We simply generate the tasks, expressing how they relate to each other, and then compute.”

User Interface
As the size and complexity of codes have increased, so has the demand for a user-friendly interface that requires less programming. A big advance within MADNESS is the ability to translate nearly directly the calculus of mathematics—the operations that combine and transform functions such as addition, multiplication, differentiation, and integration—into a numerical calculus necessary for running scientific simulations. Popular programs like Mathematica® and Matlab® perform high-level symbolic manipulations of algebra and calculus directly, thereby dispensing with the time-consuming and typically error-prone process of translating math formulae into computer code.

Figure 9. A molecular orbital of the benzene molecule projected onto a surface with the adaptive mesh also displayed.

While MADNESS does not perform symbolic manipulations, it does eliminate the tedious and very error prone step of re-expressing a few lines of math as potentially tens of thousands of lines of computer code for manipulating complicated data structures. In MADNESS, scientific applications are composed in terms of operations on functions, rather than on the raw numbers and matrices of coefficients. This is a result of the guaranteed precision associated with each operation.

MADNESS hides much of the complexity of computing on a massively parallel computer by providing several levels at which applications can be programmed. At the highest level the user is nearly totally unaware he/she is using a parallel computer. However, this will only be efficient for very large problems. At the next level down, the user composes the application as operations on many functions at the same time. This provides necessary additional parallelism and more opportunity to overlap communication and computation. At the lowest level, the user can directly access the coefficients and data structures for new functionality or maximum possible performance.

Testing and Validation
In the process of validating and testing the prototype Python code (sidebar, "Coding Languages," p61), MADNESS had been successfully applied to complex problems. Figure 9 displays the charge density profile in a single benzene molecule, upon which is superimposed the three-dimensional grid implied by the use of wavelet eigenfunctions of different widths. The machines on which the program was developed were the 32-processor IBM Power-4 nodes and 256-processor SGI Altix at ORNL. The scaling behavior of the computational effort of this prototype program is found to be between linear (proportional to the system size) and quadratic (proportional to the square of the system size), which represents a tremendous improvement over traditional HF calculations that scale quartically O(N⁴) for high precision (figure 10). The production code is targeting near-linear scaling for all systems by solving for localized molecular orbitals.

One initial test of fundamental importance was whether the results of the electronic structure calculations depend on the linear dimension of the simulation volume. To test this effect, HF calculations were performed on a series of two-electron ions isoelectronic with helium: He, Be²⁺, Mg¹⁰⁺, and Ca²⁰⁺. This test revealed that a threshold in the program needed to scale with the size of the system. With this change, the resulting ground-state energies were found not only to be independent of the cell size, but matched published literature values. HF calculations on the simplest molecules, H₂⁺ and H2, were also in full agreement with literature results.

Figure 10. This graph illustrates how quickly demand for computational resources increases as rates of scaling increase. Examples are given for linear O(N), quadratic O(N²), cubic O(N³), and quartic O(N⁴) scaling as the system being simulated doubles in size. The linear O(N) method is most desirable, requiring twice the computational resources to simulate a system twice as large To simulate the very same system, a method that scales quartically O(N⁴) demands sixteen times the computational power.

The multiresolution adaptive approach was also calibrated for DFT calculations for closed-shell systems, such as the neutral atoms He, Be, Mg, Ca, and Sr; results agreed exactly with the values tabulated by National Institute of Standards and Technology (NIST). The MADNESS team also performed DFT calculations on the water and benzene molecules, at geometries close to their equilibrium configurations. The calculations on water were performed to investigate the translational invariance of the energy—to see whether the energies changed when the nuclear positions were displaced from the so-called dyadic points, the boundaries of the mathematical functions involved. They found not only that the energy was translationally invariant, but also that the size of the orbitals—as measured by the number of functions required for a converged calculation of the wave function—was also essentially constant.

Results and Achievements
More recently, MADNESS has prototyped a complete set of fully functional calculations. These include electron energies, gradients, linear-response for polarizability, optical rotation, and excitation energies. Additionally, Dr. Harrison and his colleagues have begun to develop tools for the study of electron correlation. For what is believed to be the first time ever, “we have solved a two-electron Schrödinger equation directly in six dimensions without use of any symmetry to simplify the problem,” says Dr. Harrison. Madness has also been used for determining stacking energies of benzene molecules and DNA base pairs (figure 13). The MADNESS team is also collaborating with ORNL physicists to examine the time-evolution of electronic systems in intense laser fields motivated by experiments planned at the Linac Coherent Light Source (LCLS).

Towards Petascale
In addition to these major achievements in quantum chemistry applications, the project scientists are also working towards practical, computational improvements to MADNESS. Pushing to adapt MADNESS to massively parallel systems requires a shift in programming language, and the distributed parallel programming paradigm used by MADNESS has been described above. The prototype code used Python (sidebar "Coding Languages," p61), but for exploiting both finer grain parallelism and massive parallelism the C++ language must be used, as “Python is just not suitable for that,” explains Dr. Harrison. The advantages of Python can be retained by layering Python on top of C++, even in the highest-level C++ interface. An early version of the parallel code has been demonstrated on the 4096-processor Cray XT3 at ORNL. “But, our ultimate goal,” says Dr. Harrison, “are the anticipated petascale systems.” Expected to include multicore chips and at least 100,000 processors, the petascale machines planned for installation by 2010 present a profound challenge. Yet, the challenge was anticipated, and in preparation for such systems the MADNESS code was architected from the outset with capabilities to handle MPI sub-communicators, which enable extra levels of coarse-grain parallelism.

Figure 13. H. Sekino of Toyohashi, Japan used MADNESS to eliminate basis set incompleteness errors in conventional atomic orbital basis studies of the stacking energies of DNA base pairs. Shown here, from multiple points of view, is a cytosine dimer at a representative geometry.

Future Goals
A goal in 2007 for the MADNESS collaboration is the implementation of a full suite of chemistry tools on top of the production version of MADNESS. The new tools will be distributed as part of the high-performance computational chemistry software package NWChem. Adding capabilities for six and higher dimensions to MADNESS is also in the works. In terms of cutting-edge quantum chemistry, Dr. Harrison and his colleagues plan to complete several detailed benchmark calculations on electron correlation.

For quantum chemists, MADNESS already represents a powerful simulation system that will lead to advancement of the field. The innovative computational tools developed for MADNESS hold even greater promise for a broad range of disciplines. With these advances, an exciting new era of computational physics is emerging.

Contributors:
Dr. Robert Harrison, ORNL; Dr. George I. Fann, ORNL

Further Reading: http://www.csm.ornl.gov/ccsg/html/projects/madness.html