I. Acknowledgments

• Achi Brandt
• Ruth Pachter
• Jian Wang, Karthik Iyer, Michael Merrick, Rob Coalson
• AFOSR and NSF
• Ohio Supercomputer Center

II. Electronic Structure Methods

A. Basis-Set Methods (Quantum Chemistry)

1. Methods

   • Typically localized basis functions: gaussians (GTO) or atomic orbitals (STO).

2. Advantages

   • Analytical forms for many integrals for electron-electron interactions.
   • Input of atomic information when using atomic orbitals.
   • Controlled convergence with basis set size.
   • Variational theorem.

3. Disadvantages

   • Less-structured Hamiltonian matrices, and scaling can be severe.

1. Methods

   • Expand eigenfunctions in plane-wave basis.
   • Utilize FFTs for electrostatics problems and updating the orbitals.

2. Advantages

   • Efficiency of FFTs.
   • Lack of dependence of basis on atomic positions.
   • Control of convergence with cutoff energy of shortest wavelength mode.

3. Disadvantages

   • Difficulty if widely varying length scales.
   • Charged systems.
   • Completely nonlocal basis.
   • Implicitly periodic systems (supercells).

1. Methods

   • Finite differences.
   • Finite elements.
   • Wavelets.

2. Advantages

   • Convergence controlled only by grid resolution and order of approximation.
   • Handle multiple length scales with adaptive refinements.
   • Finite or periodic systems.
   • Locality of each function update and linear scaling/parallel algorithms.
   • Multigrid algorithms and efficiency.

3. Disadvantages

   • Domain size required?
   • Accuracy of representation?

A. Poisson problems and biophysical electrostatics

• Finite difference and finite element methods.
• Mainly SOR and conjugate gradient solvers on single level (Honig, McCammon).
• More recently application of multigrid techniques (Holst and Saied, Coalson and Beck, Tomac and Gräslund).
• Mesh refinements: Bai and Brandt, Beck, Gygi and Galli, Kaxiras, Goedecker.

B. Electronic Structure

• Goedecker (linear scaling algorithms): Rev. Mod. Phys. 71, 1085 (1999).
• Arias (wavelet methods): Rev. Mod. Phys. 71, 267 (1999).
• Beck [finite differences (FD) and finite elements (FE)]: Rev. Mod. Phys. (submitted).

C. Eigenvalue solvers with fixed potential

• Rabitz (1983): FD multigrid method for single eigenfunctions, fixed eigenvalue.
• Seitsonen (1995): high-order FD, conjugate-gradient method.
• Kolb (1993): 2-d FE method, high accuracy for molecules.
• Ackermann (1994): adaptive multilevel FE method.
• Pask (1999): FE method for solids.

D. Self-consistent problems

• Laaksonen (1985), Becke (1989): FD methods with radial grids for molecules.
• Chelikowsky (1994): high-order FD methods for Kohn-Sham equations. Examined accuracy vs. plane waves.
• Hoshi and Fujiwara (1995): FD techniques incorporating linear-scaling algorithms.

E. Multigrid techniques

• White, Wilkins, and Teter (1989): FE method using MG acceleration.
• Bernholc (1991, 1995+): MG for electronic structure, FD representation.
• Ancilotto (1999): algorithm similar to Bernholc's.
• Davstad (1992): 2d MG solver for Hartree-Fock. FD representation.
• Beck (1995+): MG for Poisson equation, nonlinear MG solver for Kohn-Sham equations. All particles on grid. FD representation.
• Martin (1999): one-way MG for Kohn-Sham.

F. Mesh refinements

• Adaptive coordinates: Gygi and Galli (1995), Kaxiras (1997).
• Mesh refinements: Fattebert (1999), uses MG acceleration in solution; Ono and Hirose (1999), double-grid method; Beck (1999).

G. Finite elements for electronic structure

• Heinemann (1987): FE for diatomics, high accuracy.
• Kopylow (1998): MG solver with 2d FE method, diatomics.
• White, Wilkins, Teter (1989): FE for atoms and diatomics, 3d, MG acceleration.
• Gillan (1995+): FE-type method with linear-scaling algorithm. Large systems.
• Tsuchida and Tsukada (1995+): FE method with linear scaling. Large systems.

Typical of real-space solvers without MG acceleration (or linearized MG) is requirement of 20-50 or more self-consistency iterations to reach ground state.

H. Time-dependent DFT calculations

• TDDFT in real time for optical response: Yabana and Bertsch (1996+), high-order FD methods.
• TDDFT eigenvalue methods for excitation energies: Chelikowsky (1999), high-order FD methods.

IV. Density Functional Theory for Electronic Structure

• Consider the KS equations in the Local Density Approximation (LDA):

  [- 1 2 Ñ2 + veff(r) ] yi(r) = ei yi(r).

  (1)

• The one-electron effective potential is:

  veff(r) = ves(r)+ vxc(r(r)),

  (2)

• The electrostatic portion of the potential (nuclear potential plus Hartree potential) given by:

  ves(r) = ó õ r(r¢) | r - r¢ | d r¢ - å i  Zi | r - Ri | .

  (3)

• This potential can be obtained by real space numerical solution of the Poisson equation:

  Ñ2 ves(r) = - 4 prtot (r),

  (4)

• r_tot 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.

• We have also developed MG methods for solving the nonlinear Poisson-Boltzmann equation of electrostatics:

  Ñ·(e(r) Ñf(r)) = - 4p[rs (r) + q _ n   +  e-bqf(r)-v(r)- q _ n   -  ebqf(r)-v(r) ]

  (5)

• All of the above are nonlinear problems.

• 2nd order matrix representation:

  1 h2 é ê ê ê ê ê ê ë -2 1 0 0 ¼ 1 -2 1 0 ¼ 0 1 -2 1 ¼ 0 0 1 -2 ¼ : : : : ··· ù ú ú ú ú ú ú û é ê ê ê ê ê ê ë f(x1) · · · f(xN) ù ú ú ú ú ú ú û = -4p é ê ê ê ê ê ê ë r(x1) · · · r(xN) ù ú ú ú ú ú ú û .

  (6)

• The Laplacian operator is discretized with a high order finite difference representation (Hamming). In the present work, a 12th order form is mainly used. Table I shows the coefficients up through 12th order.
• The impact of the order is shown in Figs. 1, 2, and 3, for the eigenvalues, first moments, and virial ratios of the 1s, 2s, and 2p states of the hydrogen atom. The calculations were performed on a 65³ grid domain where the fine grid spacing was h=0.5au.

VI. Finite-element discretization

• Localized polynomial basis set.
• Therefore, variational theorem obeyed.
• For example, if piecewise linear functions in 1d, problem looks very much like (2nd order) FD discretization.
• Adaptable.
• Number of terms along row of matrix proportional to p³, vs. 3p + 1 for FD.

VII. FAS-FMG Algorithm

• 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. (1983).

• The discrete equations can be represented on the current finest grid as:

  Lh Uh = fh.

  (7)

• For the Poisson problem L^h is the finite difference Laplacian on the fine grid with spacing h, U^h is the electrostatic potential v_es(r), and f^h is -4pr_tot(r).

• Upper case for the potential denotes the exact numerical solution. Below, lower case implies the current approximation.

• In the eigenvalue equations U^h becomes U_i^h, the eigenfunctions, where the index i indicates the orbital number.

• The eigenvalue operator L^h is the Hamiltonian minus l_i, where l_i is the eigenvalue for orbital i, and there is no source term f^h. The eigenvalue has the same value on all levels at convergence and therefore is not labeled with a grid size.

• The discrete equations on the coarse level are modified by the inclusion of a defect correction term:

  LHuH = IHh fh + tH ,

  (8)

• I^H_h is the restriction operator which takes a local average of the fine grid function.

• The defect correction is:

  tH = LH IHh uh - IHh Lh uh .

  (9)

• The fine grid function is then corrected via:

  uh ¬ uh + IhH (uH - IHh uh) ,

  (10)

• I_H^h is the interpolation operator.

• For the eigenvalue problem, constraints must be imposed on the coarsest level to maintain eigenfunction separation:

  áuiH, IHh ujh ñ = áIhH uih , IhH ujh ñ

  (11)

• On coarsest scale, require very little computational overhead.
• The eigenvalue is also updated on the coarsest level by computing the Rayleigh quotient:

  li = < LHuiH-tiH,uiH > < uiH,uiH >

  (12)

• Subspace orthogonalization via Ritz projection.
• At end of MG correction cycle, Gram-Schmidt orthogonalization is first performed, followed by diagonalization of the Hamiltonian in the basis of the occupied orbitals:

  wT L wzi -li zi = 0

  (13)

• w symbolizes the q ×N_g matrix of eigenfunctions on the grid, and the z_i are the coefficients used to improve the occupied subspace.

• The FMG approach is initiated by first iterating the self consistent problem on the coarsest level used for the eigenvalue problem.

• The problem is then interpolated to the next finer level, where the MG correction cycles begin.

• The FMG cycle is shown in Figure 4.

VIII. H Atom Convergence Behavior

• 12th order, 65³ domain, h = 0.5au.
• Potential from Poisson solution for single charge in center of domain.
• Eigenvalues: -0.50050 (1s), -0.12504 (2s), -0.12496 (2p). Three 2p states degenerate to 10 decimal places.
• Residuals can be driven to machine precision zero.
• Eigenvalues converged to 5 decimal places after a single FAS-FMG V-cycle. Total of 6 applications of Hamiltonian on the eigenfunctions on fine scale.
• Major cost is relaxation of orbitals on fine level.
• Total of about 90 sec on 450 MHz Pentium II machine.

IX. Numerical Results

• Calculations on atoms, ions, and small molecules to test the convergence and accuracy of the method.

• The convergence behavior is illustrated in Fig. 5 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 (Ancilotto, et al., 1999) on the 8 valence electron C₂ molecule. Wang and Beck, J. Chem. Phys. (in press).
  
  Figure 5: Convergence behavior. The top curve is the Car-Parrinello result of Ancilotto et al.. The second curve is the MG result from that work. 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.

• We present the ionization potentials for three atoms and the eigenvalues of the CO molecule computed at the Xa level for comparison with previous numerical work (Table III).
  Table 2: Calculations on Atomic Ionization Potentials and Eigenvalues for the CO molecule. The LDA calculations of Tong and Sham used a form for the correlation energy which interpolates between the Wigner expression at low density and the Gell-Mann-Brueckner form at high density. For the CO calculation, the grid spacing was h = 0.1335 au and the bond length was taken as 1.13Å. All energies are in hartrees.

  Ionization Potentials
  Atom	h	Tong and Sham	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
  Orbital	h	TeVelde	FAS-FMG
  1s(O)	0.1335	-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

• Large 129³ 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 (Laaksonen), and 0.10 D in finite difference pseudopotential calculations (Chelikowsky).

X. Scaling

• q orbitals and N_g fine grid points in the domain
• Relaxation of orbitals (qN_g), relaxation of potential (N_g)
• Gram-Schmidt process (q² N_g)
• Ritz projection (q² N_g to construct the matrix and q³ to solve)
• Computation of the eigenvalues on the coarse grid (qN_g^H)
• Solution of the constraint equations on the coarse grid (q²N_g^H to construct the constraint matrix and q³ to solve).
• The most costly portion of the algorithm for the small 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.

XI. High Order Mesh Refinement

• We have developed a generalization of Bai and Brandt's 2nd order multigrid mesh refinement method to arbitrary orders. [TLB, J. Comp. Chem. 20, 1731 (1999)]
• It will allow for a high resolution treatment in the regions around the nuclei and will result in a large decrease in the required total number of grid points over the composite domain.
• The method maintains the favorable scaling properties of the MG method.
• Figure 6 displays a schematic of the refined grid, and Figure 7 shows numerical results for a three dimensional singular source problem.

XII. Conclusions

• The discretization scheme is relatively simple.
• So long as the condition of zero correction at convergence is satisfied in the FAS-FMG algorithm, the solvers converge rapidly to the solution.
• The time per self-consistency update is comparable to the plane-wave method. The nonlinear solver requires less time per update than the linearized version due to handling of constraints on the coarsest level.
• Nonlinear multiscale solvers require at most a few self-consistency iterations to obtain the ground state solution.
• The Poisson solvers scale linearly with system size. For the eigenvalue problem, the method scales as q² N_g (q is the number of eigenfunctions and N_g is the total number of grid points on the fine scale) if the orbitals span the entire domain. If a localized orbital representation is possible, the method scales linearly.
• Even on uniform meshes of reasonable size, good accuracies for the valence states can be obtained in high order finite difference all electron computations. Future work will focus on mesh refinements in the region of the core and the incorporation of real space pseudopotential methods.
• The real space methods are readily modified for adaptive discretization schemes without losing the favorable convergence and scaling properties of the MG method.
• Due to the locality of the iterations, parallel algorithms can be developed in a straightforward way.
• Only one parameter controls the convergence in terms of numerical accuracy, the fine grid spacing h.

References

[]

R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).

[]

M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias and J.D. Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992).

[]

S. R. White, J. W. Wilkins, and M. P. Teter, Phys. Rev. B 39, 5819 (1989).

[]

J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72, 1240 (1994).

[]

K. Davstad, J. Comp. Phys. 99, 33 (1992).

[]

F. Gygi and G. Galli, Phys. Rev. B 52, R2229 (1995).

[]

E. L. Briggs, D. J. Sullivan, and J. Bernholc, Phys. Rev. B 54, 14362 (1996).

[]

F. Ancilotto, P. Blandin, and F. Toigo, Phys. Rev. B 59, 7868 (1999).

[]

N. A. Modine, G. Zumbach, and E. Kaxiras, Phys. Rev. B 55, 10289 (1997).

[]

K. A. Iyer, M. P. Merrick, and T. L. Beck, J. Chem. Phys. 103, 227 (1995); T. L. Beck, Intl. J. Quant. Chem. 65, 477 (1997).

[]

T. L. Beck, J. Comp. Chem. (accepted). Preprint located at xxx.lanl.gov/abs/physics/9805025.

[]

A. Brandt, S. F. McCormick, and J. Ruge, SIAM J. Sci. Stat. Comput. 4, 244 (1983).

[]

W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).

[]

S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).

[]

R. W. Hamming, Numerical Methods for Scientists and Engineers, (Dover, New York, 1962).

[]

B. Y. Tong and L. J. Sham, Phys. Rev. 144, 1 (1966).

[]

G. te Velde and E. J. Baerends, Phys. Rev. B 44, 7888 (1991).

[]

L. Laaksonen, D. Sundholm, and P. Pyykko, Intl. J. Quant. Chem. 27, 601 (1985).

[]

D. Bai and A. Brandt, SIAM J. Sci. Stat. Comput. 8, 109 (1987).