Many if not most computational physics and chemistry problems require consideration of a large range of length scales. In protein folding, both short and long ranged interactions and competitions between them lead to the final configuration of the molecule.[1,2] For the interaction between a protein and a nucleic acid or a charged interface, specific ionized groups may contribute significantly to the binding, while distant regions of the molecule have less importance.[3] If one studies the electronic structure of large molecules, there is a concentration of electron density around the atomic nuclei and between atoms in the chemical bonding regions, while often large portions of space exhibit very low and smoothly varying density.[4] In numerical simulations of fluid dynamics, there may exist specified regions which require a higher resolution treatment locally.[5]

Multiscale methods provide one approach for tackling computational problems exhibiting a range of length scales.[6,7] These methods were developed in order to overcome convergence problems in iterative solutions to partial differential equations. By utilizing approximations from coarser grids, components of the error on a wide range of length scales can be decimated, typically leading to linear scaling computing time with system size. The underlying differential equations can be represented in various ways including finite differences and finite elements. For the present method, high order finite difference representations are employed.

Based on the physical examples given above and many others requiring variable resolution, it is appropriate to develop mesh refinement strategies in which fine gridding is focused only in those spatial regions which require it.[8] The goal is to maintain the linear scaling property of the multigrid method while minimizing the prefactor in the scaling relation. The computational philosophy adopted here is to generate a sequence of nested regular grids; this strategy allows one to use the existing multigrid routines over the mesh refinement patches with no significant changes. It also allows for the implementation of accurate high order difference equations on a composite mesh. This approach is to be contrasted with the many curved grid generation techniques[10] which have been widely applied in engineering applications.

Bai and Brandt[8] addressed several pivotal issues concerning extension of the FAS multigrid technique for locally refined meshes. First, they developed a l-FMG altorithm in order to restore linear scaling behavior which can be lost when many levels of refinements are used and thus the coarser global grids are themselves visited in a way which scales linearly with the number of levels. The work to accuracy exchange rate l is the Lagrange multiplier for the grid optimization equations. Second, they determined that the same interpolation order can be used at the mesh boundaries as is used over the rest of the domain. Third, they found that local relaxation sweeps near structural singularities (to be differentiated from source singularities) can restore convergence rates to those observed away from the singularities. Finally, a second order conservative differencing method was developed for interior source singularities, when it was found that the only important factor for obtaining accurate solutions far from the singularity was to correctly reproduce the source strength around that singularity. These mesh refinement techniques were tested on several model two dimensional problems.

In the present work, the Bai and Brandt conservative differencing method is extended to high orders and three dimensions. The method is then successfully tested on a three dimensional source singularity, the Coulomb potential. First, a review of the FAS multigrid method is given. Next, a simple approach for generating the high order formulas for the Laplacian and the interpolation near the boundaries is reviewed. Then, the conservative FAS differencing forms are derived from general considerations of balancing the fluxes at the boundaries. The three dimensional conservative forms are obtained simply by locally averaging the one dimensional fluxes over surfaces at the boundaries. With the inclusion of the high order flux corrections, the summation of the defect correction over the mesh refinement is zero to machine precision. The boundary corrections lead to accuracy within the patch which is the same as that for a high order uniform fine mesh covering the whole domain; with no corrections, serious errors occur over the whole composite domain. The accuracy outside the refinement zone is improved over that for the high order uniform coarse mesh over the whole domain. The correct high order behavior is thus obtained over the entire composite domain. The only exceptions are the points one grid spacing away from the singularity, where high order is not obtained on uniform nor composite meshes.

High orders are often required in order to obtain numerically accurate results on three dimensional grids of reasonable size. While the methods derived here are general for elliptic problems, part of our motivation stems from the recent development, along with other groups, of ab initio multigrid methods for quantum chemistry.[11,12,13] In our all electron approach, all particles are represented numerically on the grid, including the nuclei. In previous work, we have carried out density functional calculations on atoms and small molecules on uniform grids, and it is clear that the majority of the numerical errors originate from the regions around the nuclei. These errors are due both to the finite size of the nucleus on the grid and to the poor numerical representation of the core electron orbitals. The methods developed in this paper will be incorporated in our electronic structure codes both for the Poisson solver and for the eigenvalue solver.

2  Full Approximation Scheme Algorithm

The FAS multigrid technique allows for solution of nonlinear problems and is ideally suited for the mesh refinement methods presented below.[6,8,9] A two dimensional schematic of the composite mesh examined here is presented in Figure 1, with two nested patches within a full domain. In the test calculations here, four or five levels were employed with the coarsest three covering the entire domain (see Figure 2). The Poisson equation to be solved on the finest patch of the four level problem is written as:

For this case, L^h4 is the finite difference Laplacian on the finest scale, U^h4 is the exact grid solution for the potential on that scale, and f^h4 is -4p times the charge density for a three dimensional problem. The current approximation to the exact solution will be written in lower case: for example the approximate h4 solution is u^h4. On the next coarser level h3, the level h4 patch covers only a portion of the level h3 domain. The equation to be solved on the h3 level is:

where I_h4^h3 is the full restriction operator which performs a weighted local average of the fine scale function, and t^h3 is the level h3 defect correction given by:

The fine scale function can then be corrected as follows:

and further iterations are subsequently performed on the h4 level. Here I_h3^h4 is the interpolation operator. One way to understand the function of the defect correction is to observe that, if the exact solution U^h4 were passed to level h3, no correction would be made. That is, the defect correction causes the level h3 equation to `optimally mimic' the level h4 problem. Note that t^h3 is only defined over the coarse grid points within the interior region of the h4 level patch, with zero values outside.

When using multiple scales, the defect correction includes an additional contribution from the previous scale. Here the example is given of the level h2 defect correction computed during the final level h4 V cycle:

By performing these coarse grid correction cycles recursively to coarser and coarser scales, errors of all wavelengths can be effectively removed with only several iterations necessary on the fine scale.

Before deriving the high order conservative forms for the composite mesh computations, a simple general procedure for obtaining the required high order Laplacian and interpolation operators is summarized.

3  Direct Method for High Order Formulas

For completeness, a summary is given here of a general method for developing any of the formulas required in a high order multigrid method. Hamming[14] outlines a direct method for obtaining numerical formulas of arbitrary order. Sample a function at N points. Define the Lagrange sample polynomials:

which are of order N-1 in x. Then

and

A polynomial which passes through the N sample points y_i is:

Expand the sample polynomial as follows:

where c_i,N = 1. Then it is easy to see that the product of the coefficient matrix and the Vandermonde matrix of the sample points (second matrix below) is the matrix [p_i(x_j)]:

Thus, except for a normalization factor of p_i(x_i) for each row of the coefficient matrix, the left hand matrix is the inverse of the Vandermonde matrix.

In the direct method, Hamming shows that there is a simple connection between the Vandermonde matrix, the desired weights for the approximation, and a vector of `moments':

These moments result from allowing the operation of interest to act on the sequence of functions: 1, x, x², x³, ¼x^N-1. For example, for the second derivative operator centered at x = 0, the first four elements of the moment vector are 0, 0, 2, 0. Similarly, moments can also be obtained for the operations of integration and interpolation. For integration, the moments are the integrals of 1, x, x², ¼ over the sampling domain, and for interpolation, they are these elementary polynomials themselves.

Since the normalized version of the coefficient matrix in Eqn. 11 is the inverse of X, the weights for each of the approximations in the multigrid process can be calculated to any desired order by one matrix-vector multiply. The normalized coefficient matrices for N sampling points are obtained by expanding the sampling polynomial (Eqn. 10) and dividing each row by the normalizing factor p_i(x_i). They are termed `universal matrices' by Hamming due to their generality, i.e. they depend only on the sampling points, not on the formula to be approximated. The matrices are tabulated in Ref. [14] up to seven sampling points, which allows for computation up to the 6th order Laplacian. In the present work, simple C codes have been written which generate the eight, nine, etc. point matrices as well. The weight vectors for the Laplacians through 8th order are given in Table 1. The three dimensional versions are generated from the sum of the three orthogonal x, y, z axes. We have utilized Laplacians up to 8th order in previous multigrid work on uniform domains.[13] The interpolation weight vectors for even numbers of sampling points are listed in Table 2 through 8th order. Similary, formulas for any higher orders can be obtained from the universal matrices. The high order interpolation formulas only need be used when setting the function values on and outside the refinement boundaries which are fixed during iterations over the patch. It is essential that the order of the interpolation match the order of the Laplacian at these boundaries. Lower order interpolations are adequate during the rest of the multigrid processing.

4  High Order Conservative Forms

When solving for the potential on coarser scales which contain a mesh refinement patch at the next finer level, it is clear from Eqn. 2 that, if the sum of t over the interior domain is not zero, additional sources have been introduced. This is in fact the case, which can be shown by examination of the t terms in a one dimensional example; most of the interior terms do cancel, but nonzero contributions remain at the patch boundaries. The terms which remain are of the form of one dimensional flux operators. Without correcting for these new sources, the solution will be polluted over the whole domain. The method of Bai and Brandt corrects for these sources by introducing local opposing fluxes at the boundary. In this section, their second order method is extended to high orders.

First, the problem is illustrated schematically by using a continuous notation (in the grid notation, all integrals go over to sums). The coarse scale is labelled by H and the fine by h. It is desired to satisfy:

where the integration is over the whole patch domain D, including the boundaries. However, it is true that

This integration is only over the interior region of the patch. Therefore:

Here t_b^H is a boundary term designed to oppose locally the additional terms due to nonconservation of source and the S integration is over a narrow strip at the surface. This implies:

The form of t_int^H is the difference of the Laplacian acting on the coarse scale function minus a local average of the Laplacian acting on the fine scale function, Eqn. 3. Therefore, converting a volume integral into a surface integral:

The brackets áñ signify a local average (restriction) of the fine scale Laplacian acting on the function, and the gradient operators Ñ_b are obtained by noncancellation of terms near the boundary of the volume integral. The final expression shows that the boundary t_b^H generates a flux which locally opposes the flux from the additional sources in the interior. Therefore, after collecting the correct units from the two scales, it is apparent that the form for t_b^H is:

where H is the coarse grid spacing, a is the numerical prefactor to the Laplacian (see Table 3), and now the gradients are one dimensional operators directed outward from the surface (determined below). Here the gradients simply represent the unitless coefficients since H² and a have been moved to the other side. For example, on a one dimensional domain (corresponding to the 2nd order Laplacian):

on the left boundary, and

on the right.

Since the process of full restriction is a weighted local average over the fine scale function (27 points in three dimensions), the averaging in Eqn. 17 can be viewed as follows. First, average over the direction normal to the boundary surface, compute the two gradient terms on the rhs of Eqn. 18, and then average over the other two directions. The full restriction weights for this process along one fine scale dimension are w = [1/4 1/2 1/4]. There is no requirement for high order restriction operators, as long as the same restriction method is used consistently. Therefore, the coefficients for the two gradient terms in Eqn. 18 can be determined by solving the one dimensional problem. In two dimensions, the fine scale gradient operator Ñ_b^h is averaged over three points along the boundary line with weights [1/4 1/2 1/4]. This yields a local average of the fine scale flux through the boundary. A similar procedure was followed in the work of Bai and Brandt.[8] In three dimensions, the local flux average is over a square centered on the location of the coarse scale gradient. The weights are 1/4 for the center, 1/8 for the edges, and 1/16 for the corners.

The one dimensional version of the flux difference in Eqn. 18 was solved for high orders by examination of the cancellation of terms near the boundary. The result for the left hand side of the coarse scale gradient on a left boundary is given by:

where n_L is the number of points in the Laplacian to the left of the center and the c_-nL+j are the Laplacian coefficients from Table 1. The rhs side of the gradient is antisymmetric with respect to these coefficients. For a right side boundary, all the signs are reversed.

The locally averaged (in one dimension) fine scale gradient coefficients are:

For the fine scale coefficients, the central term always cancels completely, so both gradient operators Ñ_b^H and Ñ_b^h are centered about the fine grid location one point inside the patch boundary. All of the coefficients up through 8th order are listed in Table 3, and the terms for a 6th order left boundary are shown in Fig. 3 to illustrate the locations. Similarly, conservative forms can be derived for higher orders if desired.

5  Computational Details and Numerical Results

The computational test case presented here is for the 4th order form. The three coarsest scales covered the whole domain, while the finest one or two were nested patches. On the three coarsest (full domain) scales, the boundaries were set by fixing the potential at the analytical value for a singular source in three dimensions f(r) = 1/r. The boundary was fixed with one additional term outside the physical boundary since the Laplacian has two terms beyond the center in one dimension. Iterations were performed over all the interior points of the full domain or patch. The FAS-MG technique was used in the form of the series of nested V cycles as shown in Fig. 2. SOR iterations were employed for all relaxation steps, with w = 1.2. The optimal relaxation parameter was determined empirically for the high order case. Full weighting restriction and linear interpolation were used, except 4th order interpolation was performed over the patch regions, including the required points beyond the boundaries. These points were set such that the Laplacian and defect correction were defined over the entire interior of the patch. The boundary potential terms for the patches were reset during the correction step after each visit to coarser scales.

The code was written in C with double precision arithmetic, utilizing the prescription of Ref. [15] for dynamic memory allocation. The test calculations were run on a Pentium 133 MHz processor laptop with 40Mb of RAM, requiring a total of roughly 15 relaxation sweeps on the fine scale and 3 seconds of total processor time for convergence. The `exact' grid results were obtained by repeated loops around the final V cycle of the FAS procedure, until the residuals were on the order of machine precision zero. The coarsest (full domain) scale had 5 points on a side, the next two finer scales 9 and 17, and the two nested patches both had 9 points on a side. To examine the order of the method, computations were performed on a full domain coarse grid corresponding to the that of the composite mesh (level h3), followed by one finer full domain grid with the spacing halved. The accuracy of the composite mesh method was then determined by comparison with the high order coarse and fine uniform grid calculations used in the determination of the order.

The increased accuracy obtained using fourth order equations vs. second order is displayed for uniform domain computations in Figure 4, in which the absolute errors of the solution are presented away from the singularity. That the fourth order Laplacian leads to fourth order behavior was confirmed in the uniform grid computations described above. Except for the set of points one grid spacing away from the singularity, the correct order is obtained over the entire domain.

Then computations were performed on the four level composite domain with a single refinement patch centered at the origin. To test the effect of the boundary correction on conservation, the integral of the defect correction over the refinement patch was computed with and without the boundary terms. Without the boundary correction, the integral was 0.8 in magnitude, while with the boundary terms, the integral was zero to double precision accuracy. The impact of the conservative boundary correction on the accuracy of the solution is apparent in Figure 5; serious errors are incurred over the whole domain in the absence of the boundary corrections. The accuracy of the method can be determined by comparison with the separate fine and coarse uniform domain results. Figure 6 shows that the accuracy within the patch is virtually identical to that for the uniform fine domain results. In Figure 7, the errors outside the refinement patch are displayed. The refinement mesh leads to increased accuracy outside the refinement on the coarser level in comparison with the uniform coarse level results. The numerical results are presented in Table 4. Finally, test computations were also successfully performed on a five level problem with two nested refinement patches. The resulting potential is plotted in Figure 8 to illustrate the accuracy of the method in relation to the numerical errors. These numerical results thus confirm that the conservative mesh refinement technique developed here leads to results of the correct high order within the refinement region, while increasing the accuracy on the coarse domain outside the refinement zone.

6  Summary

A general technique has been presented for carrying out high order mesh refinement multigrid calculations. The FAS method for composite domains was first summarized. Then, Hamming's direct method for generating high order formulas was outlined. Both the high order Laplacian and interpolation coefficients were obtained from the universal matrices of the direct method. Since the sum of the defect correction over the interior of the patch is nonzero, high order conservative forms were derived by analysis of the one dimensional problem. The two and three dimensional forms can be obtained by averaging locally over three points on a line or nine points on a square, respectively. The new method was successfully tested for the fourth order case on a Poisson problem, the source singularity in three dimensions.

The high order mesh refinement methods should allow for accurate computations on three dimensional domains which require a range of length scales. We are developing a quantum chemical Density Functional Theory (DFT) multigrid method for ab initio calculations.[13] So far, our fully numerical three dimensional calculations have been performed on uniform domains, treating both the electrons and nuclei with the high order grid approximations. As a test computation, we examined the CO molecule, and obtained good results in all electron computations. However, it is apparent from those results that the crude treatment of the nuclei and the core electrons limits the accuracy of the method. We plan to incorporate the high order composite mesh techniques into the quantum chemistry method to obtain more accurate results in the region of the nuclei, and to investigate the impact of those improvements on the eigenfunctions, eigenvalues, and total molecular energies. Beyond the quantum chemical applications, these methods should prove helpful in other large scale electrostatics calculations in biophysics and in multigrid fluid dynamics computations, especially for cases where increased resolution is required only over small local regions of space.

Acknowledgments

I thank Achi Brandt, Dov Bai, and Michael Merrick for many helpful discussions. I would like to acknowledge the support of NSF grant CHE-9632309. I also thank Daan Frenkel and Bela Mulder for support during a sabbatical leave at the FOM Institute in Amsterdam during the fall of 1996.

References

[1]

Z. Luthey-Schulten, B. E. Ramirez, and P. G. Wolynes, J. Phys. Chem. 99, 2177 (1995).

[2]

K. A. Dill, S. Bromberg, K. Yue, K. M. Fiebig, D. P. Yee, P. D. Thomas, and H. S. Chan, Protein Sci.  4, 561 (1995).

[3]

N. Ben-Tal, B. Honig, C. Miller, and S. McLaughlin, Biophys. J. 73, 1717 (1997).

[4]

R. A. Friesner, Ann. Rev. Phys. Chem. 42, 341 (1991).

[5]

J. A. Michelsen, in Multigrid Methods III, eds. W. Hackbusch and U. Trottenberg (Birkhäuser, Berlin, 1991), p. 301. See other articles in this collection as well.

[6]

A. Brandt, Math. Comp. 31, 333 (1977).

[7]

W. L. Briggs, A Multigrid Tutorial, (SIAM, Philadelphia, 1987).

[8]

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

[9]

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

[10]

Computational Fluid Dynamics, ed. H. Deconinck (Von Karman Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium, 1995); Mathematical Aspects of Numerical Grid Generation, ed. J. E. Castillo (SIAM, Philadelphia, 1991); Modeling, Mesh Generation, and Adaptive Numerical Methods for Partial Differential Equations, eds. I. Babusa, et al. (Springer-Verlag, New York, 1995).

[11]

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

[12]

E. L. Briggs, D. J. Sullivan, and J. Bernholc, Phys. Rev. B 52, R5471 (1995); ibid. 54, 14362 (1996); J. Bernholc, E. L. Briggs, D. J. Sullivan, C. J. Brabec, M. Buongiorno Nardelli, K. Rapcewicz, C. Roland, and M. Wensell, Intl. J. Quant. Chem. 65, 531 (1997).

[13]

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

[14]

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

[15]

W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T Vetterling, Numerical Recipes in C: The Art of Scientific Computing (Cambridge Univ. Press, Cambridge, 1992).

Points	Order	Prefactor	Coefficients
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.

Points	Order	Prefactor	Coefficients
N=2	2nd	2				1
N=4	4th	16			-1	9
N=6	6th	256		3	-25	150
N=8	8th	2048	-5	49	-245	1225

Table 2: Coefficients for interpolation. One side of the symmetric weight vector is shown. Each coefficient term should be divided by the prefactor.

Level	Points	Order	Prefactor	Coefficients
H	N=2	2nd	1				1
h							1
H	N=4	4th	12			-1	15
h						-1	14
H	N=6	6th	180		2	-25	245
h					2	-23	220
H	N=8	8th	5040	-9	119	-889	7175
h				-9	110	-770	6286

Table 3: Coefficients for conservative forms. One side is shown. Each term on the other half of the gradient operator has the opposite sign. The set of coefficients is for a left boundary. All the signs are reversed for a right boundary. The locations of the terms for a 6th order example are shown in Fig. 3.

r	Exact	Grid Exact	Trunc. Err.	Fine Trunc. Err.	Coarse Trunc. Err.	MG Err.
2	0.5	0.520842737	0.020842737	0.020888618		0.020875723
4	0.25	0.257687505	0.007687505	0.007710007	0.010443407	0.007726450
6	0.1[`6]	0.167993560	0.001326893	0.001183652		0.001389060
8	0.125	0.126091267	0.001091267	0.000232278	0.003854094	0.001232428
12	0.08[`3]	0.0834252199	0.0000918866	0.0000295553	0.0005908833	0.0001265184
16	0.0625	0.0625022779	0.0000022779	0.0000072359	0.0001151041	0.0000094832
20	0.05	0.0499998476	-0.0000001524	0.0000023526	0.0000347405	0.0000016978
24	0.041[`6]	0.0416667925	0.0000001258	0.0000008825	0.0000131217	0.0000007387

Table 4: Numerical results for the FAS-MG composite mesh computations at several distances from the singular source. The edge of the patch is at r = 8. The last column is for a single V cycle MG computation with a total of 16 relaxation sweeps on the fine scale. Notice that for the points in the outer regions of the domain, the single cycle MG errors are not strictly less than the anomalously small truncation errors for the composite domain; however, the errors are considerably smaller than the truncation errors on the uniform coarse scale domain (they are of the same magnitude as the uniform fine grid results at those points).

Figure 1: Schematic two dimensional cut through the three dimensional composite mesh.

Figure 2: Four level FAS cycle.

Figure 3: Locations and values for the coefficients used to generate the 6th order conservative form on a left boundary. The boundary is located on the coarse scale H at the point i and on the fine scale h at the point j = 2i.

Picture Not Created.

Figure 4: Uniform domain calculations illustrating the improved accuracy of the fourth order approximation over the second order case.

Picture Not Created.

Figure 5: Impact of nonconservation on the accuracy. The refinement patch edge is at r = 8. The diamonds are for the nonconservative calculation and the crosses for the conservative case.

Picture Not Created.

Figure 6: Comparison of the accuracy for the composite grid with that on the uniform fine and coarse scales. The refinement patch edge is at r = 8. The diamonds are for the composite grid, the crosses for the uniform fine grid, and the squares for the uniform coarse grid. Note the potential value at the r = 8 point on the composite grid is set by the coarse grid, which leads to the larger error there. The interior patch points are at r = 2,4,6.

Picture Not Created.

Figure 7: Comparison of accuracy outside the refinement patch with the uniform coarse grid results. The crosses are the uniform coarse grid results and the diamonds are for the composite grid computation.

Figure 8: Plotted are the analytical 1/r potential and the numerical results from the conservative mesh refinement multigrid computation. The two fine patches span the ranges -8. to 8. and -4. to 4. The lower curve gives the magnitude of the difference between the exact and numerical results, illustrating the larger errors near the source singularity.