Research

Electronic structure calculations of large systems using divide-and-conquer (DC) method

Introduction

The electronic structure calculation is the computer simulation clarifying the quantum behavior of electrons and nuclei that rules any and all properties of molecules. The fundamental equations determining the quantum behavior of molecular systems (i.e., Schrödinger and Dirac equations) were established in 1920s. Paul Dirac (Nobel laureate in Physics, 1933) stated in 1929 that "The fundamental laws necessary for the mathematical treatment of the whole chemistry are fully known." But he simultaneously said that "application of these laws (to an atomic or molecular system) leads to equations that are too complex to be solved." The development of the electronic structure theory oriented toward the construction of models approximating these laws for applying them to atomic and molecular systems.

Even if one uses these approximated theories, one comes up to give up the practical treatment of large systems because the computational time rapidly increases as the system size increases. In addition, the computational time becomes longer as the accuracy of the approximation becomes higher (HF → MP2 → CCSD, see Fig. 1). To accurately treat the phenomena of realistic systems such as biomolecules and chemical reactions on surfaces, we demand a trick to drastically reduce the computational cost of highly accurate approximation. We have developed the divide-and-conquer (DC) method to solve this challenging issue.

Figure 1. CPU time for HF, MP2, and CCSD calculations of polyenes C_nH_n₊₂.

DC method based on one-particle approximation

René Descartes, the Father of Modern Philosophy, stated in Discourse on Method that "Divide each of the difficulties under examination into as many parts as possible, and as might be necessary for its adequate solution^[1]." Namely, the difficulty in calculating a large system can be solved by dividing the issue into several small-system calculations. Weitao Yang^[2,3] firstly imported this philosophy to the electronic structure calculations and proposed the DC method.

The procedure of the DC method in the framework of one-particle approximation is depicted in Fig. 2. First, the entire system is divided into disjoint subsystems. To improve the description of the subsystem, atoms adjacent to the central region (called the buffer region) are taken into consideration when constructing subsystem molecular orbitals (MOs). Hartree-Fock (HF) or Kohn-Sham (KS) equation is solved in each union of the central and buffer regions (called the localization region) to obtain subsystem MOs. At this step, the electrostatic potential derived from the outside of the localization region is also taken into account. The density matrix and energy of the entire system are obtained by defining the unique Fermi level to keep the number of electrons in the entire system and averaging the subsystem density matrices. This procedure is iterated until the density matrix and energy converge. In summary, the DC method provides the subsystem MOs and the density matrix and SCF energy of the entire system.

Figure 2. Schematic of the DC SCF procedure.

Our systematic assessments of DC-HF method revealed that the computational time is drastically reduced from the conventional HF method, and that the energetic error introduced by the DC approximation can be controled to be less than a few kcal/mol by varying the size of buffer region even if an atom is adopted as the central region^[4,5]. Furthermore, we have improved the SCF convergence of the DC method by means of DC-DIIS method^[4] and DC-FON method^[6]. However, this procedure can only provides the solution of one-particle problems such as HF and DFT. We further extended to the electron correlation method, which is indispensable for accurate electronic structure calculations, by means of subsystem MOs.

DC method for high-accuracy calculations including electron correlation

Generally speaking, we require MOs for the electron correlation calculations. It is, in principle, possible to evaluate electron correlation by means of subsystem MOs. However, the correlation energies of the central regions (not the localization regions) are required for avoiding overlapped counting (Fig. 3). Then, we utilized the idea of the energy density analysis (EDA)^[7], which has been developed in Nakai group, for extracting the correlation energy of the central region.

Figure 3. Schematic of the DC-CCSD procedure of polyene C₆₀H₆₂.

We have applied this method to MP2^[8,9], CCSD^[10], and CCSD(T)^[11] theories and have demonstrated its effectiveness. As an example, Fig. 4 shows the computational time for DC-CCSD calculation of polyene chains C_nH_n₊₂. In the conventional CCSD method (black line), the computational time increases in proportion as sixth power of the number of atoms. By applying the DC method (red line), we accomplished reducing the computational cost drastically, where the time scales almost linearly with respect to the number of atoms.

Figure 4. CPU time for DC and conventional CCSD calculations of polyenes C_nH_n₊₂.

Table 1 compares the energy obtained from the DC and conventional CCSD calculations of polyene chains. The energy error introduced by DC approximation is less than 0.1 mhartree in this system. Although the magnitude of the error depends on the system and the basis set, its systematic improvement is possible by changing the size of the buffer region.

Table 1. Comparison of DC and conventional CCSD energies [hartree] of polyenes C_nH_n₊₂ (6-31G).

Concluding remarks

By utilizing the DC method, we have accomplished reducing the computational cost of highly-accurate electronic structure calculations of large systems. Although we did not mention in this review, we have proposed the DC-DM MP2 method^[12] where the MP2 correlation energy is evaluated as the functional of DC-HF density matrix^[13]. We emphasize here that these efforts of reducing the computational cost makes the impossible possible.

The evolution of the electronic structure calculation requires the advance in the computer technology in addition to the theoretical and algorithmic improvements of electronic structure methods. Even a leading-edge computer architecture often desires novel theories and algorithms. Recently, the efficient parallelization of the computation is becoming important for fast calculations using cluster computers. The DC method is suitable to the parallel cluster computations because each subsystem can be treated independently. We beleve that the solid cooperation of highly accurate electron correlation theories, fast parallel computers, and the DC method opens up a new paradigm for the theoretical chemistry.

- DC method is implemented in GAMESS program package developed by Gordon group at Iowa State University.

References

[1] R. Descartes, Discourse on Method, Project Gutenberg.

[2] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).

[3] W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).

[4] T. Akama, M. Kobayashi, and H. Nakai, J. Comput. Chem. 28, 2003 (2007).

[5] T. Akama, A. Fujii, M. Kobayashi, and H. Nakai, Mol. Phys. 105, 2799 (2007).

[6] T. Akama, M. Kobayashi, and H. Nakai, Int. J. Quantum Chem. 109, 2706 (2009).

[7] H. Nakai, Chem. Phys. Lett. 363, 73 (2002).

[8] M. Kobayashi, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007).

[9] M. Kobayashi and H. Nakai, Int. J. Quantum Chem. 109, 2227 (2009).

[10] M. Kobayashi and H. Nakai, J. Chem. Phys. 129, 044103 (2008).

[11] M. Kobayashi and H. Nakai, J. Chem. Phys. 131, 114108 (2009).

[12] M. Kobayashi and H. Nakai, Chem. Phys. Lett. 420, 250 (2006).

[13] M. Kobayashi, T. Akama, and H. Nakai, J. Chem. Phys. 125, 204106 (2006).