12478
J. Phys. Chem. 1995,99, 12478-12482
Multigrid Method for Electrostatic Computations in Numerical Density Functional Theory Michael P. Merrick, Karthik A. Iyer, and Thomas L. Beck* Department of Chemistry, University of Cincinnati, Cincinnati, Ohio 45221-01 72 Received: March 20, 1995@
A multigrid method is presented for the real space numerical solution of the Poisson equation arising in ab-initio quantum computations. The specific algorithm presented is designed to compute the total electrostatic potential in real space density functional theory electronic structure calculations. Using a high-order finite difference approximation for the Laplacian and representing the nuclei as distributed charges, accurate electrostatic potentials due to both nuclei and electrons are obtained simultaneously in contrast to separate calculations for each component. Computations are presented for finite and periodic model problems in order to illustrate the accuracy of the approximation in addition to the speed and linear scaling properties of the algorithm.
I. Introduction Quantum mechanical calculations have played a crucial role in the development of chemistry and physics. They are becoming more generally useful for routine chemical calculations with the advent of efficient electronic structure algorithms and fully quantum many-body simulation methods. These abinitio calculations impact several scientific fields, ranging from fundamental physics to materials research and biophysics. Until recently, one trait shared by all quantum methods was that calculations could not be performed on very large systems. Even with increases in computing capabilities due to hardware advances and more efficient methods, few of these algorithms have favorable scaling properties, and this places severe restrictions on the accessible system size. Our group is currently developing a linear scaling numerical real space density functional theory (DFT) algorithm. An important part of any ab-initio quantum chemical method is the efficient calculation of the electrostatic potential due to nuclear and electronic distributions. In DFT, an external potential due to the nuclei yields a unique ground state electron density.’ This ground state density is obtained when the total energy functional is minimized. The total energy of the ground state consists of a sum of kinetic, total electrostatic, and exchange-correlation energies. The total electrostatic potential is comprised of three components: the nuclear-nuclear (nn), nuclear-electron (ne), and electron-electron (ee) interactions. The ee and ne components are usually obtained using Fourier transform methods2 or by analytic solution of integrals representing the interactions in some bask3 The scaling behavior of these types of calculations is on the order N log N for the FIT-based approaches and of substantially higher order for basis set methods. Becke et aL4 have developed a numerical integration scheme to solve Poisson’s equation. The electronic distribution is decomposed into several single center problems on angular meshes, one for each nucleus. A numerical solution is obtained in spherical coordinates for each of the individual components and then summed, reconstructing the total potential. For a static system of nuclei, the energy due to the nn term is typically calculated once using an Ewald summation method. The original Ewald method is an order Mi* procedure but has been improved to order N log N in the fast Fourier Poisson formulation of York and YangeS The fast multipole method is @Abstractpublished in Advance ACS Abstracts, July 15, 1995.
0022-365419512099- 12478$09.00/0
an order N alternative approach for potential energy calculations in a system of point charges6 In contrast to the above methods which solve separately for each component of the electrostatic potential, here we develop a linear scaling multigrid algorithm to solve for the total electrostatic potential due to the nuclear and electronic distributions simultaneously. A distributed nuclear approximation is used in contrast to the usual 6 function representation of that portion for the charge distribution. A rationale for using a distributed nucleus can be given in light of a multilevel calculation itself. When solving a set of equations on a grid, the problem must first be represented discretely. On a coarse grid, the representation of the 6 function leads to a substantially smeared charge. As we pass to ever finer scales, the charge correspondingly shrinks in spatial extent. Eventually, the length scale of the nucleus must become insignificant on the scale of the electron distributions. It is the purpose of the present work and a second paper7 using the Poisson solver for atomic and molecular electronic structure computations to determine whether that limit can be obtained with a physically reasonable fine grid spacing. Similar theoretical ideas for smearing the charge have been found useful in electrostatic problems in molecular fluids and Once the nuclei are represented in this manner, the electrostatic problem can be solved in one fast linear scaling computation, since the electrons and nuclei are now on an equal footing. The discrete and local nature of the updates leads to algorithms which can be easily adapted for use on highly parallel architectures. The algorithm would be especially useful in large scale molecular dynamics or Monte Carlo simulations where the nuclear positions are altered at each step. This requires repeated calculations of the total potential. Much has been written about the scaling and accuracy of the multigrid (MG) method when applied to matrix equations of the form Au = f,e.g., Poisson’s equation:
Our algorithm delivers the electrostatic potential at every grid point for finite or periodic systems. Using a high-order finite difference approximation for and representing the nuclei as distributed charges, we have successfully implemented an MG algorithm to solve for the total electrostatic potential in one step as part of an ab-initio Kohn-Sham DFT calculation. The energies of several atomic systems as well as the binding energy 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 33, 1995 12479
Multigrid Method for Electrostatic Computations of a two-atom system have been calculated to high accuracy. That work appears in a separate papere7 In this paper, solutions to three model problems are presented in order to illustrate that accurate physics ean be obtained quickly using a distributed nuclear approximation with an MG Poisson solver. First, a singular source is placed within a nonuniform neutralizing distribution in order to mimic the electrostatics of an atomic system. Second, the Madelung constants for two cubic salts (NaC1 and CsC1) are computed. These results show the method's ability to accurately and rapidly calculate quantities associated with long range forces. Finally, electrostatic maps are generated for CsCl.and sodalite, a zeolite of high symmetry.'O These electrostatic maps may be of interest in many areas of study to explain the energetic interactions that are localized in a region of importance. In the case of the simple zeolite, the catalytic activity may be viewed as resulting from the electrostatic topography of the crystal lattice. In the next section, an overview of the method is given; readers are referred to the works of Brandt, Briggs, and Hackbusch for further details concerning MG techniques.' In sections I11 and IV, the specific examples and results are presented. Finally, the results are summarized and future plans discussed in the concluding section.
and r is a quantity called the residual. The residual is a measure of how well the solution satisfies the original difference equation: r=f-Av
(3)
The error can be obtained on coarse grids and used to improve the fine grid solution in a process called coarse grid correction. On the coarse grids the long wavelength modes of an error appear more oscillatory and are quickly damped by the smoothing algorithm. The necessary work is also substantially less on a coarse grid. The order N scaling, accuracy, and speed of the MG method are achieved when the errors on the various grids are themselves improved in an iterative fashion by successive coarse grid corrections. Grid transfer is performed by interpolation and restriction operators. Interpolation is defined as the process of transferring a function to a fine grid from a coarse grid and restriction the reverse of that. Our propagation method uses a weighted Gauss-Seidel (SOR) smoothing algorithm which is represented here for illustration using a three-point finite difference approximation to the VZ operator in one dimension:
11. Multigrid Method for Solution of Poisson Equation Our MG algorithm is designed to compute the electrostatic potential of atomic systems for DFT calculations using a distributed nucleus approximation. In the model problems presented we use either a single cube or a Gaussian to represent the nuclear charges. Both functional representations are normalized to the total charge. In three dimensions, the value of the charge density in the cube approximation is simply a h 3 at the grid point where the nucleus is centered; h is the grid spacing. In this manner, the charge distribution is confined to exactly one grid point. The total distribution is the sum of the electronic and nucleur contributions. The Gaussian functions are represented in the usual manner. As the grid spacing or the width of the Gaussian tends toward zero, increasingly accurate representations of the 6 function are obtained, i.e., both forms can be considered prelimit Q functions. The MG method developed by Brandt" and extended by others is a highly efficient technique used to solve the algebraic equations resulting from a discretized partial differential equation. Originally applied to boundary value problems, it is an alternative to direct methods used to solve sparse linear systems of equations. Direct methods are often difficult to apply to many interesting and important problems. As will be seen, the MG method is fast, accurate, and easy to apply. MG is an extension of various iterative techniques and wits developed to overcome a shortcoming common to all real space iterative methods, namely, critical slowing down (CSD). This phenomenon results from an iterative method's inability to effectively remove long wavelength components of an error. However, an advantage of the iterative methods is that they rapidly damp the short wavelength components of the error. Once the short wavelength parts of the error are removed, the technique stalls and improvements to the solution are made only very slowly. The MG technique overcomes this shortcoming in a unique way. For a set of equations represented by Au = 5 it can be shown that the error, e = u - v (u is the exact solution and v is its current approximation), associated with a given approximate solution satisfies Ae = r
(2)
This equation is usually called the residual equation, where A is the discrete operator of the original problem, e is the error,
The weighting parameter is denoted by w , and m is the iteration number. As shown in the equation, the method uses the updated values of the solution as soon as they are calculated in addition to the weighting to improve the convergence properties of the algorithm. In this paper, calculations were performed using fourth- and sixth-order finite difference approximations to the operator. We found that propagators of an order higher than 2 were required to obtain physically accurate results. The numerical coefficients for the approximations were obtained from Vichnevet~ky.'~ An accommodative scheme was used to recognize when the smoothing algorithm stalls, which initiates grid transfer. Intergrid transfers were performed by linear interpolation and full weighting restriction processes.
111. Model Problems Two types of problems were examined to exhibit the method's ability to handle the three parts of the total electrostatic potential in an ab-initio calculation. A singular source in an exponentially decaying neutralizing background mimics a finite atomic problem with an analytic solution for compari~on.'~ The charge distribution (atomic units),
(5) leads to a potential
The continuous problem was modified such that the singularity at r = 0 in the second term was removed and the nucleus was represented as a cube of charge. The value at the origin was chosen such that the grid computed integral of e ( r ) over the interval was zero. The electronic charge density in a real atom is of course finite at the nucleus. The parameters Z and a were both given values of 1. The computed solution leads to a large but finite value of the potential at the origin, in contrast to the infinite analytical result there. The lattice calculations used both types of distributed nuclear representations mentioned above. The accuracy of the Madelung constants,I6 a,calculated from
12480 J. Phys. Chem., Vol. 99,No. 33, 1995
-ae2
Epr= r
Menick et al.
(7)
-
0.8
analytic Q-0 calc
clearly indicates the method's ability to handle long range force calculations. The energy per ion pair, Epr,used to compute the Madelung constant was obtained from
where Ecalcis the total energy calculated using the solution @(I-) returned by the MG calculation in the equation for the electrostatic energy:
w dr
= '/zje(r)
(9)
is the nuclear self energy, which is dependent on the representation, and Npris the number of pairs in the periodic box. The self energy can be calculated analytically for the Gaussian representation or by a one-time boundary value calculation using the same Poisson solver. The CsCl structure had two atoms represented by Gaussians, one in the center of the box and eight on each comer, each of which is shared by 8 other unit cells. The CsCl model represented the true unit cell with two atoms per cell in the body centered cubic lattice form. The NaCl type structure was constructed with eight full Gaussians or eight cubes in each box centered at I/,& increments in three dimensions. The sodalite unit cell has space group I", and consists of 36 atoms. We used the cubic charge representation for the 12 Si and 24 0 nuclei in the unit cell.'" The 36 cubes of charge were placed on the nearest grid point. Periodic boundary conditions were enforced in the usual manner. There is of course no requirement of regular lattices, but they offer a good means of testing the accuracy of the present approach.
IV. Results We fist examine the numerical solution of eq 5, a singular source in a neutralizing background, to begin to assess the accuracy of the method. This problem contains interactions of the ee and en types. A fourth-order representation was employed for the V operator. The finest scale had M3grid points. In Figure 1, the computed result is compared with the analytical solution for grid points away from the origin. It is clear the numerical solution is close to the analytical result, except for grid points adjacent to the origin, where small deviations occur. These results are encouraging; however, the true test of the accuracy of the method is obtained from use of the Poisson solver in an electronic structure computation. We have obtained accurate grid solutions to the Kohn-Sham equations for multiorbital atomic problems, and these results are presented in another paper.7 The grid solutions for various cube widths are compared in Figure 2 to show that the errors spread in a zone around the cube of charge as it is distributed over a larger volume. This is to be expected. Confining the cube of charge to increasingly smaller volumes provides a better approximation to the 6 function, which is regained in the limit of h -,0. The important conclusion here is that numerically accurate solutions to the Poisson equation involving a singular source can be obtained using the distributed nucleus approximation, a high-order finite difference representation, and our multigrid solver. This model problem illustrates where the errors are to be expected and how they vary with the distributed nucleus approximation. The grid solution can be obtained in a matter of seconds (=20 for 323 problem) on a SPARCstation ELC, starting from an initial guess for #(r) of zero everywhere.
Figure 1. Comparison of the analytic and computed solutions of the Poisson equation resulting from the charge distribution of eq 5.
0.8 -h
M2h Q--i3 3h
0.8
Figure 2. Variations in the calculated solution as the cube of charge is distributed over two and three grid spaces in each direction (x, y, and z). TABLE 1: Energies (hartree) and Madelung Constants Obtained Using a Cubic Charge Distribution with the NaCl Lattice Model for Various Gridso 15.220 236 2 31.753 488 9 64.819 178 7
16.532 660 4 33.065 658 9 66.131 242 4
1.747 938 1.747 599 1.747 566
16 32 64
the energies were used in eqs 7 and 8 to compute the Madelung constants (a). The Laplacian was represented with a fourth-order approximation. The interval was held constant at 10.654 713, and the grid spacing varied as 10.654713" where N is the number of grid points. The number of pairs (Npr)was 4, e2 is 1, and r = 5.327 356 5. Atomic units were used in all lattice calculations. The literature value of a is 1.747 558.
We now discuss the solutions for periodic problems possessing long range forces. Here we have used both the cubic and Gaussian representations to test the relative merits of each. Problems involving sets of bare charges (nn contribution) without the distributed background of the previous problem provide a stringent test for a real space solver. We first used the cubic representation of the charges to compute the Madelung constant for NaCl, with eight charges in the periodically replicated box. The results in Table 1 show that, even on a 163grid, the Madelung constants are obtained to high accuracy. As expected, the accuracy of the calculation is increased as the fine grid spacing is reduced. Increases in accuracy are also seen when the representation of V is improved, as shown in Table 2. The accuracy is noticeably higher for the sixth-order
J. Phys. Chem., Vol. 99,No. 33, 1995 12481
Multigrid Method for Electrostatic Computations
TABLE 2: Data in Table 1 When a Sixth-Order Approximation to v2 Is Used Ecslc
Eselt
a
14.489 926 7 30.292 051 6 61.896 248 2 125.104 640 30
15.802 102 9 31.604 196 1 63.208 392 1 126.416 784 2
1.747 607 1.747 565 1.747 565 1.747 565
N 16 32
64 128
Mlnlnum -0.17
TABLE 3: Data Obtained Using a Gaussian Charge Distribution and a Sixth-Order Approximationa Eca~c Exit a N a9 2.454 912 3 2.449 228 5 2.449 124 0 2.449 122 0
3.767 036 8 3.761 361 6 3.761 263 2 3.761 263 2
1.747 530 1.747 550 1.747 558 1.747 561
16 32
64 128
1.739 850 1.747 420 1.747 559 1.747 561
a In addition to the data in Table 1, the Madelung constants calculated with the analytic self energy (aa)are given.
TABLE 4: Data in Table 3 Using a Gaussian Cha e Distribution and a Sixth-Order Approximation to g f o r a CsCl Lafflcd
0.00
2.46
4.92
7.38
9.84
Figure 3. Electrostaticpotential in a plane midway between the cesium _____
0.737 730 4 0.736 797 4 0.736 780 1
0.941 258 4 0.940 332 9 0.940 315 7
1.762 600 1.762 667 1.762 668
16 32
64
1.754 293 1.762 609 1.762 671
the interval was held constant at 10.0 and the grid spacing varied as lO.O/N, where N is the number of grid points. The number of pairs (Npr)was 1, ez is 1, and r = 8.660 245. The literature value of a is 1.762 670. a
approximation on the coarse grids but becomes less apparent as the grid spacing is reduced. In all further lattice calculations, a sixth-order approximation was used. The average computing time for the 163grid was roughly 2 s on a SPARCstation 5. To test the scaling properties, the number of ions in a NaCl type lattice was increased from 8 to 64 as the number of grid points was changed from 323 to M3. The average ratio of computing times on a Cray-YMP was 8.07 for complete multigrid solution starting from zero potential everywhere. The ratio of the two times definitely points toward linear scaling with system size, as expected for a full multigrid solver. Defining the work-unitI2 as the amount of work necessary to perform a fine grid iteration, the average number of work-units for a calculation of the type described above was 13 with an average of 8 fine grid iterations. The computations were stopped on the basis of a residual criterion of ( ~ k , r ; ) ' ' * / nI10-6. Distributions which do not depend on the grid spacing were also used to test the generality of the algorithm. We chose to represent the nuclei here using Gaussian functions, which in principle allows the use of an analytically calculated self energy. The self energy for a Gaussian with (3 = 0.6 au of length is 0.470 158 hartree. The results obtained using a Gaussian representation for the NaCl lattice model are given in Table 3. Also given is the Madelung constant calculated using the analytic self energy in eq 8. In general, the Madelung constants calculated using the numerical self energies are more accurate. This is of course because the numerically generated solution for the entire system contains within it the self energy corresponding to the grid representation. The two values calculated in different fashions approach each other as the number of grid points is increased, and the grid representation of the charge distribution becomes more accurate. Mixing of analytical and numerical results is not reasonable until the discrete representation of a problem approaches the continuous case. Table 4 gives the same data for the CsCl lattice. Here the calculation box was the unit cell itself, so each comer contained just one eighth of a Gaussian for the C1 atoms, and the Cs atom was located at the center of the cell. Again, the method yields highly accurate
ion and four chlorine ions. There is a negative maximum at the comers and a positive maximum in the center.
Madelung constants.I7 The Poisson solver delivers the potential at every grid point (which is required for minimizing total energies in electronic structure calculations). As an illustration, Figure 3 displays a contour map of the electrostatic potential in a plane midway between the positive ion and the four negative ions. These results show that high accuracy can be obtained for the long range potentials even though the representation of the charge is rather rough. However, with a crude representation, the self energy must be obtained numerically in a onetime calculation. A slightly more complicated test case concerns the mineral, sodalite which does not possess the simple CsCl or NaCl structures. An Ewald calculation for the case where the ions were placed on the nearest grid points of a 643 grid of unit length gave 1120.38 au of energy compared to 1120.21 au from our multigrid algorithm using 36 cubes of charge (0.015% difference). The energy obtained using the continuous space positions was 1124.17 au (0.34% difference from the Ewald result for grid coordinates). The difference between the continuous space total and that for the ions placed on the grid shows that a substantially greater error is incurred by placing the problem on a grid than by solving the problem with our numerical Poisson grid method (see below). The Ewald program was obtained from a standard FORTRAN recipe.I8 A mapping of the potential in the plane of the Si atoms and another shifted slightly toward the 0 atoms is presented in Figure 4. This is the potential a unit charge would feel in the plane.
V. Discussion The results of the model problems presented show promise for use in a wide range of applications. The algorithm is relatively simple and fast, and scales linearly in computer time with the number of charged particles. The discrete nature of the method leads to an algorithm which is highly parallelizable. Although crude, the representation used for the nuclei leads to results that we consider surprising in their accuracy. The simple cube of charge representation yields results that are comparable in accuracy to those obtained using a well represented Gaussian function. Desired accuracies can be obtained by choosing appropriate finite difference representations of the Laplacian and fine grid sizes. Our examples show that, with finite difference representations of fourth or sixth order, realistic results can be obtained with a moderate number of grid points. It should be
12482 J. Phys. Chem., Vol. 99, No. 33, 1995
SODALITE
Figure 4. (a) Electrostatic potential in the plane containing the silicon atoms. The Si atoms are connected by bridging 0 atoms, two of which are above and the other two are below the shown plane. The electrostatics of the represented plane is dominated by the Si atoms, yet one can discern the presence of the bridging oxygens in the electrostatic maps. (b) Electrostatic potential in the plane between the plane of the Si atoms and the plane containing the bridging oxygens. The presence of the oxygen atoms is more clearly visible in this map.
noted that the higher order representations will have larger bandwidths and thus increase the number of computations necessary for each update. For example, the sixth-order finite difference representation of lP has a total of 19 terms in three dimensions in contrast to 7 for the second-order formula. In addition to the speed and scaling properties of the method, the algorithm yields the electrostatic potential at every grid point. This is a distinct advantage over the methods which give the total energy or the potential at a point. It is ideally suited for real space Kohn-Sham DFT calculations using a minimization strategy, since the potential at every point is required for each propagation step.* When the Kohn-Sham equations are also represented directly on the grid at the same level of accuracy as the Poisson equation, the entire solution process can be carried out consistently and efficiently in real space. Having addressed all three terms (ee, en, nn) of the electrostatic potential in our model calculations and having carried out a series of computations on atomic and simple molecular species at the LDA level: we believe this combined approach holds promise in numerical real space DFT. The clear advantage of a real space approach is that one can utilize localized orbital methods to overcome the M bottleneck due to the required orbital orthonormalization. If this is done and if a multigrid solver is used for the KohnSham equations as well, a fully linear scaling algorithm results for the total energy calculation. For large scale Monte Carlo and molecular dynamics simulations, special methods may be required to properly assign molecular positions. For example, our results on the zeolite problem show that greater errors result from placing the problem on the grid than from the numerical approximations in our grid solver. Methods have been developed in the simulation of plasmas for distributing charges on neighboring grid points to accurately compute particle forces.19 In addition, one of the most important features of the multigrid method is its high degree of adaptability. In a large scale simulation, the movement of a particle in one region of the box may not require that
Merrick et al. the total potential be recomputed over the entire box due to charge relaxation and resulting screening effects. Using an algorithm such as the one we've described, it is straightforward to perform a boundary value calculation on the portion of the grid where the potential changes. The fixed boundary for the new calculation will be chosen some distance from the center of the region of interest. This volume may be small in comparison to the size of the entire problem, leading to large increases in computational speed. In addition, composite grids will give the ability to focus on regions of interest where the potential changes rapidly. Problems requiring periodic boundary conditions in one, two, and three dimensions can be performed easily. Therefore, this method should provide a viable alternative to recently developed approaches for carrying out Ewald sums for three-dimensional problems with two-dimensional periodic boundaries, namely surface problems involving charged particles or dipoles.20 Electrostatic maps of large systems can be generated relatively quickly with modest grid sizes; the present algorithm and related methods should thus prove useful in biological studies of electrostatic driving forces in active or recognition sites.
Acknowledgment. We thank Achi Brandt, Rob Coalson, and J. Thomas King for helpful discussions. This research was supported by NSF Grant CHE-9225123. We would like to thank the Ohio Supercomputer Center for a grant of time on the Cray-YMP, on which some of these computations were done. T.L.B. thanks AFOSR and Dr.Ruth Pachter for a summer fellowship at Wright Patterson Air Force Base. References and Notes (1) Parr, R. G.; Yang, W. Density-Functional Theory Of Atoms And Molecules; Oxford University Press: Oxford, 1989. (2) Remler, D. K.; Madden, P. A. Molecular Physics 1990, 70, 921. (3) Szabo, A.; Ostlund, N. S. Modem Quantum Chemistry-Introduction to Advanced Electronic Structure Theory; McGraw-Hill: New York, 1982. (4) Becke, A. D.; Dickson, R. M. J. Chem. Phys. 1988, 89, 2993. (5) York, D.; Yang, W. J . Chem. Phys. 1994, 101, 3298. (6) Greengard, L.; Rokhlin, V. J. Comput. Phys. 1987, 73, 325. (7) Iyer, K. A.; Menick, M. P.; Beck, T. L. J. Chem. Phys., in press. (8) Rosenfeld, Y. In Strongly Coupled Plasma Physics; Rogers, F. J., Dewitt, H. E., Eds.; NATO AS1 Series B; Plenum Press: New York, 1987; Vol. 154, p 543. (9) Onsager, L. J . Phys. Chem. 1939, 43, 189. (10) Meier, W. M.; Olson, D. H. Atlas of Zeolite Structure Types, 2nd ed.; Butterworth: London, 1987. (1 1) Brandt, A. Math. Comput. 1977, 31, 333. (12) Briggs, W. L. A Multigrid Tutorial; SIAM: Philadelphia, PA, 1987. (13) Hackbusch, W. Multi-Grid Methods and Applications; SpnngerVerlag: Berlin, 1985. (14) Vichnevetsky, R. Computer Methods for Partial DifSerential Equations; Prentice-Hall: Englewood Cliffs, NJ, 1981; Vol. 1. (15) Arfken, G. Mathematical Methods for Physicists; Academic Press: San Diego, CA, 1985; p 921. (16) Ashcroft, N. W.; Mermin, N. D. Solid State Physics; Saunders College: Philadelphia, PA, 1976; Chapter 20. (17) The box lengths were 10.65 and 10.0 atomic units for the NaCl
and CsCl models, respectively. The values of (J were chosen such that the Gaussians would not overlap. Various values were tested to ensure that the calculations were independent 0. The presented data were obtained using a value of (J = 0.6 au for both models. (18) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (19) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Adam Hilger: Bristol, 1988. (20) Hautman, J.; Klein, M. L. Mol. Phys. 1992, 75, 379.
JP950789D
