The plane-wave pseudopotential method has proven to be an excellent computational strategy for solving large-scale electronic structure problems in condensed phases[1,2]. Notable strengths of the method are the ability to use fast Fourier transforms for updating the self-consistency equations, lack of dependence of the basis on atom positions, and the clear control of convergence with the cutoff energy determined by the shortest-wavelength mode. However, the method can encounter difficulties in treating widely varying length scales. This issue is relevant for surfaces, clusters, or systems where the pseudopotential varies rapidly near the nucleus such as first row elements or transition metals. Also, charged systems present significant complications in a plane-wave code. In recent years, considerable effort has been focused on alternative real-space approaches utilizing finite elements or finite difference representations[3,4,5,6,7,8,9,10,11,12]. Advantages of these approaches include the ability to handle finite or periodic systems with equal effort and the locality of each iteration step. Locality leads to simplicity in developing domain decomposition parallel algorithms. In addition, it is relatively straightforward to implement adaptive grid-refinement techniques in order to focus effort in spatial regions with large variations in the computed functions, for example near the nuclei in electronic structure computations[4,7,10,14]. Finally, representation directly in real space allows for the use of multigrid (MG) algorithms with their excellent convergence characteristics and scaling properties[15].
Several groups have employed MG algorithms for various portions of the solution process for self-consistent electronic structure computations. White et al. [3] developed a finite element method for electronic structure which utilized an MG solver for the Poisson equation; they also presented preliminary results on application of MG methods for one-orbital self-consistent eigenvalue problems. Davstad[6] discretized the Hartree-Fock equations for diatomic molecules and solved the resulting two-dimensional equations with a combination of the MG method and a Krylov subspace method for the coarsest-grid equations. Gygi and Galli[7] proposed an adaptive-coordinate approach which places increased resolution near the nuclei using curved grids. They solved the Poisson portion of the problem with MG techniques. Briggs et al.[8] have developed an MG solver for the Kohn-Sham (KS) equations which utilizes MG methods for both the eigenvalue and Poisson problems. Their MG eigenvalue solver uses a double-discretization approach and a linearized version of the MG method. Recently, Ancilotto et al.[9] have presented a similar method for solution of the KS equations, and have applied the method to calculations on charged lithium clusters. They show that the MG approach is as accurate and more efficient than the corresponding Car-Parrinello method. Modine et al.[11] developed an adaptive-grid real-space method which employs MG preconditioning for solution of the Poisson equation. Finally, we have developed high-order MG methods for solving the KS equations including all the electrons[12,13]. Typical of the MG solvers to date is the requirement of 20 or more self-consistency cycles to obtain convergence. In this paper, we present a high-order multigrid algorithm for solving the self-consistent Kohn-Sham equations which obtains convergence in an order-of-magnitude less numerical effort. The method utilizes the nonlinear Full Approximation Scheme (FAS), full multigrid (FMG) eigenvalue technique of Brandt et al.[15] with the further inclusion of new methods for the self-consistency portion of the problem.
The first step in development of the MG algorithm is the discretization of the problem in real space. Here we employ a high-order finite difference representation on a Cartesian grid. This approach has been shown to yield accurate results in pseudopotential and all-electron calculations[5,12,13]. Consider the KS equations[16] in the Local Density Approximation (LDA):
| (1) |
The one-electron effective potential is:
| (2) |
with the electrostatic portion of the potential (nuclear potential plus Hartree potential) given by:
| (3) |
This potential can be obtained by real-space numerical solution of the Poisson equation:
| (4) |
where rtot is the total charge density on the grid due to the electrons and nuclei. The exchange-correlation potential in our spin-restricted LDA computations is calculated with the VWN functional[17].
The Laplacian operator is discretized with a high-order finite difference representation[18]. In the present work, a 12th-order form is used. The high-order expression yields a large gain in accuracy, which reduces the number of required grid points significantly [13]. The same Laplacian is used for solution of the Poisson and KS equations, resulting in a consistent level of accuracy throughout the solver. Due to this consistent representation on all levels, the only parameter controlling the accuracy of the solution is the fine grid spacing h. The other quantities are diagonal in the coordinate representation. The real-space discretization scheme offers efficient calculation of electron densities and, as a result, the exchange-correlation potentials. The nuclear charge density is the discretized form of the Dirac delta function. The wavefunctions are vectors of length Ng, the total number of grid points on a given level. Integration is performed by simple trapezoid summation on the three-dimensional domains. The boundary conditions on the orbitals for the finite systems examined here are set to zero, while the electrostatic potentials on the boundary are obtained via multipole expansion of the total charge density to quadrupole order.
The KS eigenvalue problem is nonlinear since one simultaneously solves for both the eigenvalues and the eigenfunctions. Therefore, we have employed the FAS eigenvalue technique of Brandt et al.[15] This method allows for solution of nonlinear problems with efficiencies similar to linear ones. Here we outline the FAS method of solution for the eigenvalue and Poisson problems. Both solvers are required for the self-consistent problem.
The discrete equations can be represented on the current finest grid as:
| (5) |
For the Poisson problem Lh is the finite difference Laplacian on the fine grid with spacing h, Uh is the electrostatic potential ves(r), and fh is -4prtot(r). Upper case for the potential denotes the exact numerical solution. Below, lower case implies the current approximation. In the eigenvalue equations Uh becomes Uih, the eigenfunctions, where the index i labels the orbital number. The eigenvalue operator Lh is the grid Hamiltonian Hh minus li, where li is the eigenvalue for orbital i, and there is no source term fh.
In the FAS method, the desired functions are represented on each grid level. The discrete equations on coarser levels, however, are modified by the inclusion of a defect correction term which is required to yield zero correction at convergence. Suppose we have an approximate solution on the fine grid h, uh, which is obtained by relaxation sweeps on the grid:
| (6) |
We have used Gauss-Seidel updates for these smoothing steps on all levels. The coarse-grid problem with grid spacing H = 2h can then be constructed as:
| (7) |
where IHh is the restriction operator which takes a local average of the fine grid function. Full-weighting restriction is employed here, which weights the local points according to trapezoid-rule integration. The defect correction is:
| (8) |
The coarse-grid equation is iterated to obtain an improved solution on that level. In a two-grid method, the fine-grid function is then corrected via:
| (9) |
where IHh is the interpolation operator. These operations can be recursively extended to yet coarser grids. By inclusion of multiple levels, errors of all wavelengths are rapidly damped in the fine-grid function. The only change required on grids two or further from the finest grid is that an additional term IhH th must be included in the defect correction to incorporate information from the previous level. A multigrid V-cycle consists of 1) initial iterations on the finest level, 2) passage to the sequence of coarser levels where the coarse-grid equations are iterated, and 3) return to the fine scale via correction steps and further iterations on each new fine level (Fig. 1).
In solving the eigenvalue equations in KS theory, some additional components are required. On all grids except the coarsest, the orbital equations are updated just as outlined above. However, on the coarsest level, constraints must be imposed to maintain eigenvalue separation (otherwise all states begin to collapse to the ground state). If one had the exact solution on a fine grid Uih and then restricted the orbitals to the next coarser level, the functions would no longer be orthogonal. Hence, the constraints on the coarsest level can be implemented by solving a linear matrix equation:
| (10) |
These equations, since they are implemented only on the coarsest scale, require very little computational overhead, and can be solved either via exact matrix inversion or by iterative methods. The matrix size is of order q ×q, where q is the number of orbitals. The eigenvalue is also updated on the coarsest level by computing the Rayleigh quotient:
| (11) |
Since the eigenvalues are the same on all levels at convergence, they are not labeled with a grid level.
Finally, the convergence can be improved by including subspace orthogonalization via Ritz projection. At the end of an MG correction cycle, a Gram-Schmidt orthogonalization is first performed on the finest level, followed by diagonalization of the Hamiltonian in the basis of the occupied orbitals:
| (12) |
where w symbolizes the q ×Ng matrix of eigenfunctions on the grid, and the zi are the coefficients used to improve the occupied subspace. This Ritz projection matrix problem is thus also of size q ×q.
The FAS-FMG approach outlined above typically requires a total of only 10-20 relaxation steps on the finest level to obtain convergence to within the errors due to the discrete representation (for both Hartree and eigenvalue problems). Higher accuracy is achieved by increasing the order of the approximation and/or decreasing the grid spacing (the energy resolution increases with the inverse of h2). The convergence of the accuracy with finite difference order and fine mesh spacing has been addressed previously in Refs. [5,13]. Results of chemical accuracy can be obtained in pseudopotential calculations on grids of order 653 or smaller with high-order finite difference methods.
The FMG approach is initiated by first iterating the self-consistent problem on the coarsest level used for the eigenvalue problem. The initial approximation for the eigenfunctions is a set of random numbers of magnitude one. The coarsest level is chosen so that the eigenvalues have the same ordering as on the finest level. For the systems studied here, a coarse grid of 93 or 173 is required, depending on the problem. The Poisson solver extends the levels up to a 33 grid where only the one central point is iterated. In most of these calculations the finest grid employed is of size 653, so the eigenvalue solver comprises three or four levels, while a total of six levels is utilized for the Poisson solver on the finest scale. During the initial iterations on the coarest level, the eigenfunctions are first relaxed and orthogonalized via a Gram-Schmidt process, and then the potential is updated. Once an initial approximation is obtained on the coarse level, the solution is interpolated to the next finer level, where the MG correction cycles are initiated.
At the beginning of each MG V-cycle, the effective potential is updated once upon entry to the new finest level with the FAS method. Subsequently, the potential is updated at the end of each eigenfunction correction cycle; since the changes to the potential are relatively smooth, the potential corrections are performed with a simple V-cycle using the previous potential as input. Schematic diagrams of the FMG solver are presented in Figs. 2 and 3. The Poisson equation is typically solved with only a few relaxation sweeps on the finest scale on either side of the V-cycle. During each eigenfunction correction cycle, the eigenfunctions are relaxed 3 times on each side of the V-cycle. The solution of the entire electrostatic portion of the problem thus requires similar computational effort as for the update of a single eigenfunction. This algorithm differs from those of Refs. [8,9] by inclusion of updates of the eigenvalues and enforcement of the constraints on the coarsest level only without resorting to linearization of the equations; therefore, while the convergence behavior is more rapid, the total numerical effort in one V-cycle is less. One self consistency V-cycle for a 653 domain (14 electrons) requires roughly one minute of CPU time on a 350 MHz Pentium-II processor, including the Poisson and eigenfunction updates. The solution time increases linearly with Ng for fixed q, and as mentioned above, the energy accuracy is inversely proportional to h2.
We have explored two approaches for the eigenfunction correction cycles. In the first, each eigenfunction is carried through the V-cycle starting with the lowest energy function, and the effective potential is updated at the conclusion of the cycle. This sequential method leads to a rapidly convergent effective potential since the low-lying states are quickly stabilized. In the second approach, all eigenfunctions are corrected simultaneously, and the effective potential is updated upon conclusion of the V-cycle. For small systems the first approach is preferred since the Poisson updates are quite inexpensive. However, as we discuss below, the scaling of the two approaches differs, and the second method is preferable for larger systems. Its convergence behavior is slightly less dramatic, requiring three or four self-consistency updates as opposed to two for the first method. The algorithm can readily be adapted for either approach depending on the problem.
We have performed calculations on atoms, ions, and small molecules to test the convergence and accuracy of the method. The numerical results presented here were obtained with the second approach discussed above. The nuclei were placed on individual grid points; in future work, grid refinement techniques[14] will allow for placement of multiple nuclei at arbitrary locations. The convergence behavior is illustrated in Fig. 4 in calculations on the 4 electron Be atom and the 14 electron CO molecule. Comparison is made to recent MG and Car-Parrinello pseudopotential calculations[9] on the 8 valence electron C2 molecule. The FAS-FMG approach leads to dramatic acceleration; we emphasize that the present calculations include the nuclear singularity in the effective potential, making it a substantially more challenging scenario for convergence to the ground state. The slightly reduced acceleration for the CO molecule (relative to the Be atom) is linked to the convergence behavior of the two core states, which can be alleviated with mesh refinements.
We computed total energies for atoms and ions to obtain atomic ionization potentials at the LDA-VWN level (Table I). The total energy obtained for the He atom is -2.8345 au vs. the exact value of -2.8348 au. The computations on the ions are performed with equal effort as for the neutral species, and the IP results are of satisfactory accuracy[19]. Finally, we present the eigenvalues of the CO molecule computed at the Xa level for comparison with previous numerical work[20] (Table I). The computation was performed on a large 1293 domain since the dipole moment of CO is a sensitive function of the overall domain size. We obtained a dipole value of 0.25 D C-O+ compared to the previous result of 0.24 D in fully converged numerical calculations[21], and 0.10 D in finite difference pseudopotential calculations[5]. As discussed in Ref. [3], the most severe errors in real space all-electron calculations occur in the regions around the nuclei. The limitations on the accuracy of our uniform domain results is a consequence of these errors. Two improvements are suggested: inclusion of pseudopotential techniques to remove the core and/or local grid refinements in the neighborhood of the nucleus[7,10,14]. We have generalized a second-order multigrid local mesh refinement method[22] to arbitrary order[14] and will include these refinements in the KS solver in future work. Since the refinements are truly local, they require only modest computational overhead.
To conclude, we discuss the scaling properties of the various steps in the FAS-FMG algorithm for the case where the orbitals span the whole domain. With q orbitals and Ng fine grid points in the domain, the scalings of the important components of the algorithm are as follows: relaxation of orbitals (qNg), relaxation of potential (Ng), Gram-Schmidt process (q2 Ng), Ritz projection (q2 Ng to construct the matrix and q3 to solve), computation of the eigenvalues on the coarse grid (qNgH), and solution of the constraint equations on the coarse grid (q2NgH to construct the constraint matrix and q3 to solve)[23]. The most costly portion of the algorithm for the systems examined here is the relaxation step for the orbitals. For large systems where orbital localization is possible, each of the steps becomes linear scaling with the number of electrons since every operation is then local in space (except the generation of the Hartree potential which can be solved with linear scaling using MG methods)[24]. Other linear-scaling algorithms exist for the eigenvalue problem in the context of basis-set methods [25,26]. The two algorithms for implementing self consistency discussed above differ in that the first approach, while exhibiting stronger convergence to the ground state, leads to a qNg scaling which persists even with localized orbitals.
Highly efficient methods have also been developed to accelerate convergence in basis-set self-consistent field calculations. The direct inversion in the iterative subspace (DIIS) extrapolation minimizes the error vector of the off-diagonal blocks of the Fock matrix using information from previous steps of iteration[27]. This algorithm leads to tight convergence with roughly ten self-consistency steps. However, it requires the storage of the eigenvectors from several previous steps. Also, care must be taken to obtain a good initialization of the orbitals or the method can be unstable[28,29]; several steps with an alternative procedure (such as conjugate gradients) are typically taken to obtain a good guess before initiating the DIIS extrapolation. A direct efficiency comparison between the real-space and basis-set methods is more difficult than between real-space and plane-wave approaches due to the very different underlying structures of the methods.
We have presented a new method for solution of the KS equations in real space which utilizes nonlinear multigrid techniques to solve the self-consistent equations. The purpose has been to illustrate the rapid convergence behavior of the FAS-FMG algorithm in relation to related electronic structure methods. The efficiency is a result of the preconditioning on coarser levels and nonlinear treatment during the MG correction cycles. The combination of real-space discretization with FAS-FMG multiscale techniques of solution leads to order-of-magnitude improvement in convergence behavior relative to other MG methods and the Car-Parrinello approach, and it should thus prove a useful numerical stategy for solving large-scale electronic structure problems.
We would like to thank Prof. Achi Brandt for many helpful discussions concerning this work. The research was supported by NSF grant CHE-9632309.
| IONIZATION POTENTIALS | |||
| h | Ref. 19 | FAS-FMG | |
| He® He+ | 0.15 | 0.969 | 0.974 |
| Li® Li+ | 0.35 | 0.198 | 0.193 |
| O® O2+ | 0.12 | 1.843 | 1.821 |
| EIGENVALUES | |||
| Orbital | Ref. 20 | FAS-FMG | |
| 1s(O) | -18.745 | -18.832 | |
| 1s(C) | -9.912 | -9.940 | |
| s(2s) | -1.045 | -1.060 | |
| s*(2s) | -0.489 | -0.494 | |
| p(2p) | -0.413 | -0.407 | |
| s(2p) | -0.304 | -0.302 |
.
Figure 1. Multigrid V-cycle. The iterations begin on the left side on the finest level. The problem is restricted to the next coarser level (R) and iterations are performed there. This process is continued until the coarsest level is reached. After iterations on the coarsest level, the approximation on the next finest level is corrected (C) followed by relaxation steps. This is continued until the finest scale is reached. Only a few iterations are required on each level on each side of the V-cycle.
Figure 2. The FMG scheme. Iterations begin on the left side of the diagram on the coarsest level. The next finer-grid solution is initiated with the interpolated functions from the previous coarse grid approximation. The algorithm proceeds by moving through the sequence of finer resolution grids with a multigrid V-cycle on each level. The self-consistent solution is obtained with two or a few passages through the final V-cycle. The effective potential is updated (the Hartree and exchange-correlation potentials are generated) at each point labeled with a filled circle.
Figure 3. Flow chart of the FMG self-consistent eigenvalue solver. The process is initiated with iterations on the coarsest level. Here the eigenfunctions and the effective potential are updated in each iteration step (the cost of each iteration is 1/64 that on the finest level in a 3-level problem). The Hartree potential is computed with an FAS multigrid solver. The problem is then interpolated to the next-finer level to begin the multigrid procedure. The Hartree potential is generated on the current finest level with the Poisson multigrid solver, and the exchange-correlation potential is computed directly on the grid to yield the new effective potential (EP). A new set of eigenvalues is obtained from a Ritz projection. Then the fine scale eigenfunctions are relaxed (EF Update) for a few steps. The eigenfunctions are corrected and the eigenvalues updated during the eigenvalue problem coarse-grid correction (CGC) V-cycle. Afer returning to the fine level, the EP is updated and a final Ritz projection is performed. The variables k and n label the grid level and the number of self-consistency iterations, respectively; kmax is the finest grid level and ntot is the total number of self-consistency cycles on each level.
Figure 4. Convergence behavior. The top curve is the Car-Parrinello result of Ref. 9. The second curve is the MG result of Ref. 9. The next is our result for the CO molecule with the FAS-FMG solver. The bottom curve is the FAS-FMG result for the Be atom.
. .
. .
