Quantum chemical calculations involve solving the Schrödinger eigenvalue equation for many interacting electrons. The problem requires self consistent solution since the effective potential due to the nuclei and electrons depends on the eigenfunctions. A numerical solver therefore must include methods for solving both the Poisson and eigenvalue equations. In this talk, a method is discussed which discretizes the Kohn-Sham equations in real space with high order finite difference representations. The Poisson and eigenvalue problems are then both in sparse matrix form. The equations are solved with a Full Approximation Scheme (FAS), Full Multigrid (FMG) method which allows for solution of nonlinear problems and includes efficient preconditioning from the coarser levels. The resulting algorithm rapidly locates the ground state electron density, requiring only two or three self consistency cycles on the fine scale. Numerical results are presented which compare the nonlinear MG solver to previous MG and Car-Parrinello results; the present algorithm leads to a significant increase in efficiency.

Papers discussing this research can be found at http://xxx.lanl.gov/abs/cond-mat/9905422 and http://xxx.lanl.gov/abs/physics/9805025.

1  Acknowledgments

I would like to thank Prof. Achi Brandt and Dr. Jian Wang for many helpful discussions concerning this work. The research was supported by NSF grant CHE-9632309.

2  Electronic Structure Methods

2.1  Basis Set Methods

1. Methods

   • Typically localized basis functions: gaussians or atomic orbitals.

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

   • Very large matrices can result, and scaling can be severe.

2.2  Plane Wave Methods

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.
   • Periodic systems only.

2.3  Real Space Methods

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 iterations and parallel algorithms.
   • Multigrid algorithms and efficiency.

3. Disadvantages

   • Domain size required?
   • Accuracy of representation?

3  Density Functional Theory for Electronic Structure

• Consider the KS equations[13] 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) = å i  Zi | r - Ri | + ó õ r(r¢) | r - r¢ | d r¢.

  (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[14].

• 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)

• We have used this PB solver in simulations of flexible polyelectrolyte chains to study effects of counterion condensation on chain properties. The chain structure computed at more accurate PB level differs substantially from linearized Debye-Hückel calculations.

4  Discretization

• The Laplacian operator is discretized with a high order finite difference representation[15]. In the present work, a 12th order form is used. Table 1 shows the coefficients up through 8th order.
  
  

  Points	Order	Prefactor
  N=3	2nd	1				1	-2
  N=5	4th	12			-1	16	-30
  N=7	6th	180		2	-27	270	-490
  N=9	8th	5040	-9	128	-1008	8064	-14350

  Table 1: Coefficients for the Laplacian. One side plus the central point are shown. Each coefficient term should be divided by the prefactor. The Laplacian is symmetric about the central point.

• The 12th order form has coefficients -50, 864, -7425, 44000, -222750, 1425600, and -2480478 with prefactor of 831600.

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

5  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.[12]

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

  Lh Uh = fh.

  (6)

• 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 ,

  (7)

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

  (8)

• The fine grid function is then corrected via:

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

  (9)

• 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 ñ

  (10)

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

  (11)

• 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

  (12)

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

                  x      x \      x
     x    x \    x \    /   \    /
 x  x \  /   \  /   \  /     \  /
x \/   \/     \/     \/       \/
                                         

6  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.

7  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[8] on the 8 valence electron C₂ molecule.
  
  Figure 5: Convergence behavior. The top curve is the Car-Parrinello result of Ref. 8. The second curve is the MG result of Ref. 8. 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. This figure is rotated 90 degrees due to tth converter.

• 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[17](Table 2).
  
  

  Atom	h	Ref. 16	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	Ref. 17	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

  Table 2: Calculations on Atomic Ionization Potentials and Eigenvalues for the CO molecule. The LDA calculations of Ref. 16 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.

• 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[18], and 0.10 D in finite difference pseudopotential calculations[4].

8  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.

9  High Order Mesh Refinement

• We have developed a generalization of Bai and Brandt's[19] 2nd order multigrid mesh refinement method to arbitrary orders.[11]
• 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.

. .

10  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 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

[1]

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

[2]

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

[3]

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

[4]

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

[5]

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

[6]

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

[7]

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

[8]

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

[9]

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

[10]

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).

[11]

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

[12]

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

[13]

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

[14]

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

[15]

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

[16]

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

[17]

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

[18]

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

[19]

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