Solving the Two-Electron Integral Problem in Hartree-Fock-Theory

Theoretical methods for studying fundamental biochemical processes at the electronic level hold tremendous promise for enabling rational drug design, the interpretation of biophysical measurements, and understanding the complex mechanisms of life. A primary difficulty in such theories is that, as the size of the systems become larger and more complex, the empirical methods become less reliable. One promising approach to avoid empiricism is to employ first principles quantum mechanical methods. Although the quantum mechanical equations governing chemistry are well known, their mathematical solution is extremely demanding of computer time-so demanding that it is impractical to apply these methods to most biochemical problems.

University of Minnesota Supercomputer Institute researchers Matt Challacombe and Eric Schwegler, who had been working with Chemistry Professor Jan Almlöf before his death in January, 1996, have developed novel methods for reducing the cost of first principle calculations and have made an initial application to the p53 tumor suppressor tetramerization monomer shown in Figure 1. Mutations in the p53 tumor suppressor are the most frequently observed genetic alterations in human cancer. A central model in first principles quantum chemistry is the Hartree-Fock theory. In this model, each electron moves in the average potential from all other electrons, and the exchange interaction (a manifestation of the Pauli exclusion principle) is exactly accounted for. Missing from the Hartree-Fock approximation is the dynamical correlation resulting from their pairwise Coulomb repulsion. However, the Hartree-Fock approximation is often the starting point for theories that do account for electron correlation. Within the density functional theory, Hartree-Fock exchange appears to be a necessary ingredient in the best functionals, which are able to achieve chemical accuracy (1 kcal/mol).

Figure 1 Electrostatic potential isosurfaces of the p53 tumor suppressor tetramerization monomer. The isosurfaces correspond to ±0.06 atomic units, with the red surface corresponding to positive values. This calculation involved 3836 basis functions and was carried out on an IBM RS6000 model 590 workstation. Timing for one Fock build is 3.3 central processor unit (cpu) hours.

The Two-Electron Problem

The Hartree-Fock self-consistent field equations are solved iteratively, requiring assembly of a Fock matrix with each cycle. Conventional implementations of Hartree-Fock theory are dominated by computation of exchange and Coulomb contributions to the Fock matrix. Calculation of the exchange and Coulomb matrices are expensive because they involve the computation and handling of two-electron integrals, which are four index quantities. Calculation of the two-electron integrals (and therefore the Fock matrix) is formally an N^4 process, where N is the number of basis functions.

In practice, the use of Gaussian basis sets and the thresholding of negligible integrals employed by state-of-the-art quantum chemistry programs yield computational complexities that are approximately N^3, and which approach a quadratic limit as N becomes large. Nevertheless, applications of Hartree-Fock theory to systems with more than 2000 basis functions are rare; one of the largest calculation to date (N = 2700) was possible only through the exploitation of high (icosahedral) symmetry. For macromolecules, conventional ab initio methods are inapplicable.

Linear Scaling

These researchers have recently introduced linear scaling methods for computation of the Fock matrix that are valid for insulating systems (not conductors). Linear scaling computation of the Fock matrix is achieved with specialized methods that exploit the unique characteristics of the exchange and Coulomb interactions.

The exchange interaction is a quantum phenomenon. For insulators, this interaction is short range, a feature that has been exploited to obtain linear scaling methods for computing the exchange matrix. These methods are described in more detail in the Journal of Chemical Physics, 105, p. 2726, (UMSIReport 96/17, March 1996).

Figure 2 Average timings for construction of the Fock matrix two-electron components corresponding to a sequence of water clusters. Average times for one MONDO SCF cycle are also shown. Linear fits to the Fock matrix two-electron component times are shown as an aid to gauging linearity. Correlation coefficients of both fits are r = 0.999.

Linear scaling computation of the Coulomb matrix is accomplished with hierarchical clustering technologies that rely on rapidly convergent series expansions to reduce complexity. This method is equivalent to computing the gravitational interaction between the Earth and the Moon by using the product of their masses divided by the separation of their center of masses, rather than adding up the individual interactions between all particles that compose the Earth and the Moon. This method is a tree-code implementation of the multipole series and provides an alternative to fast multipole re-expansion methods used earlier by other researchers. This method is described in UMSI Report 96/168, September 1996. These linear scaling methods have been interfaced with a modified version of the ab initio quantum chemistry code HONDO 95.3, developed by Michel Dupuis of Pacific Northwest Laboratory. This modified version is called MONDO. MONDO Hartree-Fock calculations using the 3-21G basis set have been performed on a series of approximately spherical water clusters taken from an equilibrated molecular dynamics simulation. Timings for computation of the exchange and the Coulomb matrices are shown in Figure 2, together with average times for one MONDO SCF cycle. These timings were obtained on an IBM RS/6000 model 590 workstation.

Ab Initio Electronic Structure of Proteins

With linear scaling methods, large systems of biophysical relevance are accessible to first principles quantum chemical studies. One important molecular property readily obtained from Hartree-Fock theory is the electrostatic potential, which is widely implicated in molecular recognition and binding. Electrostatic potential isosurfaces of the p53 tumor suppressor tetramerization monomer are shown in Figure 1, and were derived from MONDO Hartree-Fock calculations using the 3-21G basis set.

Linear scaling methods already exist for computation of the exchange-correlation matrix in density functional theory. In combination with the linear scaling methods for Hartree-Fock calculations discussed here, these techniques open the door to accurate quantum chemical methods that are appropriate for the study of large biochemical systems.

For a more detailed discussion of this research, please see the two University of Minnesota Supercomputer Institute Research Reports mentioned in this article (UMSI 96/17 and UMSI 96/168). They are available by contacting our Research Report Coordinator, Susan Kilber Kalenze, kilber@msi.umn.edu.

In This Issue:

Electronic Structure Workshop

Two-Electron Problem in Hartree-Fock Theory

Upcoming Symposia

People & Visitors

Research Scholars

Research Reports