8356
J Phys. Chem. 1990, 94, 8356-8359
The h = 0 Term in Coulomb Sums by the Ewald Transformation M. W. Deem, J. M. Newsam,* and S. K. Sinha Exxon Research and Engineering Company, Clinton Township, Route 22 East, Annandale, New Jersey 08801 (Received: January I I , 1990; In Final Form: April 26, 1990)
An explicit, surface integral formulation of the h = 0 term in the reciprocal space component of the Ewald transformation for Coulomb sums in the general case is derived, and the magnitude of this term is calculated for six representative noncentrosymmetric crystal structures.
Introduction
Atom-level simulations of ionic or partially ionic systems generally include a treatment of the Coulomb (charge-charge) interactions:
There is a substantial and growing literature on applications to extended crystal structures (refs 1-9 provide a small number of examples), for which the energy conventionally computed is that for a central unit cell (which is assumed representative of the crystal energy in the asymptotic limit of a large crystal). The sum performed on a unit cell by unit cell basis is then
where ri is the position of atom i (which has charge qi) in the unit cell, I are the lattice vectors inside the crystal volume V ( L ) for which L is the characteristic scaling length, and the prime indicates that terms with I = 0, i = j (the self-self terms) are skipped. This sum does not converge unless Cqi = 0. Even then the sum is only slowly and conditionally convergent (Le., the value of the sum depends on the order in which the terms are added). The convergence can be appreciated by representing the sum over ij terms as the interaction between two dipoles, one at 0 and one at I, which is 0 ( r - 3 ) ,whereas the space of the I vectors is 0 ( r 3 ) . Methods to facilitate convergence of lattice sumsl(t12have been explored by a number of authors, notably, Ewald,I3 Evjen,I4 B e r t a ~ t , Nijboer '~ and De Wette,I6 and Williams" (including ~~
~~
~~~~
(1) Catlow, C. R. A.; Mackrodt, W. C. Computer Simulation of Solids; Lecture Notes in Physics No. 166; Springer-Verlag: Berlin, 1982. (2) Catlow, C. R. A. Solid Stare tonics 1983, 8, 89-107. ( 3 ) Price, G. D.; Parker, S. C. Phys. Chem. Miner. 1984, 10, 209-216. (4) Matsui, M.; Matsumoto, T. Acta Crystallogr. 1985, 841, 377-382. (5) Fujino, T.; Morss, L. R. J . Solid State Chem. 1987, 67, 131-141. (6) Jackson. R. A.; Catlow, C . R. A. Mol. Simulation 1988, I , 207-224. (7) Ooms, G.; van Santen, R. A.; den Ouden, C. J . J.; Jackson. R. A,; Catlow, C. R. A. J . Phys. Chem. 1988, 92, 4462-4465. ( 8 ) Karasawa, N.; Goddard. W. A. J . Phys. Chem. 1989, 93, 7320-7327. (9) Dove, M. T. A m . Mineral. 1989, 74, 774-779. (10) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: London, 1954. ( 1 I ) Tosi, M. P. In Solid State Physics: Advances in Research and Applicarions Volume 16; Seitz. F., Turnbull, D., Eds.; Academic Press: New York, 1964; pp 1-120. (12) Glasser, M. L.; Zucker, 1. J. In Theoretical Chemistry: Aduances and Perspecfiues Yolunie 5; Academic Press: New York, 1980; pp 67-139. (13) Ewald, P. P. Ann. Phys. 1921, 64, 253-287. (14) Evjen, H. M . Phys. Reu. 1932, 39, 675. (15) Bertaut. F. J . f h y s . Radium 1952, 13, 499-505. (16) Nijboer, B. R. A.; De Wette, F. W. Physica 1957, 23. 309-321. (17) Williams. D.E. Acta Crysrallogr. 1971. A27. 452-455.
0022-3654/90/2094-8356$02.50/0
generalizations to cases where the dependence of contributions on distance scales as higher inverse powers of r ) . These most commonly used Ewald transformation formulas are correct for centrosymmetric crystal structures but give results that disagree with an explicit evaluation of eq 2 in cases in which the unit cell has a net dipole moment.'8-22 Earlier authors have derived expressions for the necessary shape-dependent correction term for simple c a s e ~ . I ~ -A~ lgeneral expression has also been derived,22 although in the form of a relatively cumbersome six-dimensional integral. We present here an alternative derivation of the correction to the classical Ewald sum (appropriate in the general case) that takes the form of a surface integral, which can be determined analytically in simple cases or evaluated numerically in the general case. The technique developed by Ewald and further exploited by subsequent authors involves the introduction of a convergence function, 4 ( r ) , and separation of the sum into two components
-
-
Conditions for a useful convergence function are that (1) d(r) should tend rapidly to zero as r a,( 2 ) the limit as r m of ( 1 - 4 ( r ) ) / rshould be finite, and (3) the three-dimensional Fourier transform of ( I - 4 ( r ) ) / rmust be computable. Satisfaction of these conditions permits evaluation of the second term in eq 3 in reciprocal (Fourier) space essentially by an application of a 3-D Poisson summation formula. For a Coulomb ( l / r ) sum the conventional choice of convergence function, + ( r ) , is the complementary error function 2 @ ( r )= erfc ( a r ) = -
(4)
With this convergence function23 ( a determines how much of the sum is performed in real space and should be optimized for ~
~~~
(18) Harris, F. E. Theor. Chem. ( N . Y . ) 1975, 1, 147-218. (19) Redlack, A.; Grindlay, J. J . Phys. Chem. Solids 1975, 36, 73-82. (20) De Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London 1980, ,4373, 27-56. (21) Perram, J. W.; De Leeuw, S. W. Physica 1981, 109A, 237-250. (22) Smith, E. R. Proc. R . SOC.London 1981, A375, 475-505. (23) The complementary error function meets the requirements on the convergence function as
erfc ( a r )
-?r'I2ar
as r
-m
A proof of relationship iii, using a convergence function, polar coordinates, interchange of the order of integration, and Cauchy's theorem is available from the authors; see also ref 16.
0 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94, No. 21, 1990 8357
h = 0 Term in Coulomb Sums run time for a particular structure, vide infra, but a = 0.3 is usually a good initial value), the restriction on I = 0, i = j , in the second component on the right-hand side of eq 3, S2,can be removed by subtracting the limiting values of the excluded terms
Development
Explicit indication of the order of summation allows for further manipulations on S2. To specify the order, we introduce a summation cutoff function, SCL(r):
The real exponential in the second term, S,, can be ignored, for FT,{SC,l(g) exists only in a small region about the origin in reciprocal space for L large. S4 = 2nU !- L-m lim
11
y F T 3 { S C L , ( g - ~ d7 ) b(g) dg (17)
Using the definition of b(g) above selects the value of the integrand at g = 0, yielding (7)
We further require that lim FT,ISC,(r)l(g) = b(g)
L--
(8)
which is satisfied for any reasonable V(L). The second summation in S2 (eq 6) can be rephrased as an integral by the use of delta functions
s, =
This expression is formally equivalent to eq 3.14 in the earlier work of Smith.22 One means of making use of this expression would be, where tractable, to specify explicitly the shape function, SCL, and then take its Fourier transform, FT3{SCL).22 We can, however, reduce eq 18 to a more insightful form. Using Parseval's second theorem again produces an integral in real space. We exploit the following Fourier transform pair:24
(9)
Use of Parseval's second theorem
s
SCL(r) dr
(20)
f(r) g(r) dr =
l F ( g ) G(-g) dg, F(g) =FT,tf(r)l,
G(g) = FTJg(r)l (10)
the convolution theorem F T 3 W g(r)l = FT,tf(r)l*
FT,lg(r)l
(11)
and the space translation property FT3tf(r+ro)} = e2*jrroFT,{(f(r))
(12)
yields (the Fourier transform of the real space delta functions, 111, is the set of reciprocal space delta functions, {h],divided by the unit cell volume, UI6
Recognizing the familiar structure factor, F(g), defined as F ( g ) = xqie-2*ig.r,
(14)
I
we can write S? =
At all reciprocal space vectors except h = 0 the term in the innermost brackets has a valid limit, and using the limiting value of FT3(SCL](g)we can reduce the integral to a sum plus the integral over the h = 0 term.
We demonstrate that the contributions to this integral are nonzero only at the surface of the summation volume and that the term S4therefore arises physically from interaction with a surface charge. Thus, consider grouping the values of r inside VL into sets of N points r = rr - rJ k, for which ri and k are constant. The value of the integrand in eq 21 for this fixed ri is then EqJ/lkl.As x q , = C q l = 0, this contribution is zero. This grouping cannot, however, be done at the surface, for certain of the vectors r would then protrude from the volume. (The fl singularities at r = r, - r, contribute nothing to the integral, for they are less singular than delta functions and exist only at a finite number of points; if the real exponential term in eq 16 had been retained, these singularities would, in any event, have been avoided.) To evaluate the integral in eq 20 for the surface contribution, then, consider V' E V" E V ( L ) ,for which the surfaces of V'and V"are defined such that the integral over V ( L )has no contribution from within V" and the integral over Y" has no contribution from within V'. The integration is taken as written over r E V'. Then, at the V'boundary, a switch is made to the asymptotic approximation of integrating the interaction between two dipoles d = x q l r l ,one at 0 and one at r. The error associated with this representation, neglect of dipolequadrupole and higher pole-pole terms, tends to 0 as L m, for such terms are on the surface O(L-"), n I3, whereas the surface area is O ( L 2 ) . The uniformly polarized volume V ( L ) - V' is analyzed as a surface charge on a[ V(L)-V'l. The surface charge on the inner surface of V(L)-V' must be canceled by an asymptotically equal and opposite surface charge on dV', for the contribution to the integral (20) over V ( L ) is zero within V". Thus, the only remaining term is the surface charge on a V ( L ) .
+
-
(24) This relationship is obtained by some algebraic manipulation in polar coordinates.
8358
Deem et al.
The Journal of Physical Chemistry, Vol. 94, No. 21, 1990
TABLE I" material
space group
unit cell contents
ref
quartz
P322 1
Si,O,
26
u-cristobalite
P4,2?
Si,Ox
26
zeolite $
P4122
si640128
27
BaTiO!
P4mm
BaTiO,
28
KTiOPO,
Pn2,a
K8Ti8P8040
29, 30
Ba M n F,
A2,am
Ba4Mn4F,6
31
explicit
Ewald (including h = 0)
h=O
-28.961 2 -28.751 8 -41.7227 -42.5 184 -696.497 -696.565 -10.5855 -10.6045 -71 . I 797 -73.3420 -1 6.5248 - I 6.5203
-28.9621 -28.7519 -41.7223 -42.5187 -696.496 -696.568 -10.5853 - 10.6045 -7 I ,1794 -73.3428 -16.5248 -16.5203
3.9533 4.1635 2.3489 1.5525 3.9181 3.8470 1.7961 1.7769 9.9952 7.8318 0.0560 0.0605
Sb PC
S
P S P
S P
S P
S P
'Equation 22 was checked for this table by comparing the results of explicit summations out to a large distance from the central unit cell (typically 250 A. column 5) with the result produced by eq 23 (column 6). The magnitude of the h = 0 term included in the Ewald sum is given in column 7. Formal valence charges were used, save for the phosphate group in KTiOPO, for which Pt and 0- were assumed. Results are given simply in units of (electronic charge)* ,&-I; the conversion factor 332.06 converts to the conventional units of kcal mol-'. The code used to compute the numbers given in this table will be supplied to interested readers. bSpherical SC for which the h = 0 term is calculated using eq 25 first derived by de Leeuw et al.20 'Parallelopipedal SC for which the h = 0 term given by eq 23 was evaluated to 15 significant figures by adaptive 16-point Gauss-Legendre quadrature
The general form of the h = 0 term, when the summation is performed on a unit cell by unit cell basis, is thus .SA =
- lim S - d d-r S d.n 2u
L-.m
av(L)
r3
This surface charge integral happens to be independent of L. The full expression of the Ewald transformation is then written
The historic answer for the Coulomb sum, S, is the sum of only the first three terms in eq 23.I3-I7 (Incidentally, the reciprocal space sum need only be performed over half of reciprocal space, for IF(g)l = lF(-g)l.) As noted earlier,'8-22the magnitude of the h = 0 correction term depends on the shape of the summation volume and on the unit cell dipole moment (and, hence, on how the unit cell is defined. For a real crystallite the choice is not arbitrary as the contents must be contiguous and, for termination in complete unit cells, the appropriate choice is dictated by the termination planes.) The term can be substantial for noncentrosymmetric structures (Table I). The proper Ewald transformation, eq 23, is the only nonvanishing term in the expansion of S, eq 2, in the asymptotic limit of a large crystal, L a. Consideration of this treatment in terms of its application to a finite (but still ideal) crystal indicates that terms that scale as various inverse powers of L become considerable. For example, S4will contribute a next higher term which is O(L-'). Choosing an analytically tractable form of the cutoff function, SC,, (see, for example, ref 22) will likely permit asymptotic analysis of the convolution integrals in eq 15, thus generating finite-crystal corrections to the h # 0 terms as well. Such a derivation has not been pursued here, for we anticipate that accelerated evaluations of the finite sum in eq 225will prove
-
(25) Greengad, L. F. The Rapid Evaluation ofPotenria1 Sums in Particle Systems; MIT: Harvard, MA, 1988. (26) Hyde, B. G.: Anderson, S. Inorganic Crystal Structures: Wiley: New York, 1988. (27) Newsam, J . M.; Treacy, M. M. J.; Koetsier, W. T.; deGruyter, C. B. Proc. R . Soc. London 1988, A420, 375-405. (28) Frazer, 8. C.:Danner. H. R.: Pepinsky, R . Phys. Rev. 1955, 100, 745-746. (29) Stucky, G .D.; Phillips, M. L. F.; Gier, T. E. Chem. Mater. 1989, I , 492-509
more useful for these special small crystal cases. Further, eq 2 may no longer be meaningful for very small crystals, for the energies of noncentral unit cells will then need be considered. Examples As an example of the application of this formula to a specific case, consider performing the Coulomb sum over all unit cells a. S, in polar coordinates within a sphere of radius R as R becomes
-
This expression is identical with that derived by de Leeuw et al.*O for the special (simplest) case of a summation of cubic unit cells within a sphere. Extensions to other summation cutoff volumes may require numerical integration of eq 22. We have, by way of further example, evaluated the term for some topical and representative materials summed with both spherical and parallelopipedal cutoff volumes (Table I). Efficient, explicit, direct space alternatives to the traditional Ewald approach might be found among Euler-like transformations. Numerical validation of such explicit summations should be based on results from application of the full Ewald result (23) in which the surface charge effects have been included. These results relate to an ideal, isolated crystal. Formal application of Born-von Karman periodic boundary conditions implies omission of the h = 0 surface charge term. In a real environment, substantial surface charges arising from the particular distribution of charge density within the individual (bulk) unit cells may not be maintained. Surface reconstruction, termination in incomplete unit cells, scavenging of appropriately charged species, and interactions with other crystallites will all serve, in general, to diminish the surface charge. For the aluminosilicate framework systems containing sorbates that we are currently simulating, we thus generally choose to treat the h = 0 term as zero, even though for an "ideal" crystal in isolation this term can be substantial (Table I). Conclusion
An explicit, surface integral formulation of the h = 0 term in the reciprocal space component of the Ewald transformation for (30) Voloshina, 1. V.; Gerr, R. G.; Antipin, M. Y.; Tsirel'son, V. G.; Pavlova, N. 1.; Struchkov, Y. T.; Ozerov, R. P.; Rez, I. S. Kristallografya 1985, 30, 668-676. ( 3 1 ) Keve, E. T.; Abrahams. S. C . ; Bernstein, J. L . J . Chem. Phys. 1969, 5 1 . 4928-4926.
J . Phys. Chem. 1990, 94, 8359-8362 Coulomb sums in the general case has been derived by using simple Fourier transformations. The magnitude of this term has been calculated for six representative noncentrosymmetric crystal structures. In each case the term is substantial, up to 17% of the total lattice energy sum. In applying the Ewald transformation
8359
to simulations of real systems, it appears likely that surface charges will largely be canceled in the real environment, although this assumption should be weighed in the context of the other simplifications made in such applications, such as a general neglect of finite crystal size and the effects of defects and disorder.
Isothermal Compressibility of SPC/E Water Kazi A. Motakabbir and M. Berkowitz* Department of Chemistry. Vniuersity of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27599 (Received: January 25, 1990; In Final Form: May 3, 1990)
Molecular dynamics computer simulations on rigid SPC/E water molecules were performed. The goal of the simulations is to study the behavior of isothermal compressibility, which was calculated in the simulations at T = 298, 273, and 248 K. The calculated isothermal compressibilitiesat these temperatures display a trend contrary to the experimental observations. The hydrogen-bonded network in water was also investigated. No correlation between the temperature dependence of isothermal compressibility and the number of hydrogen-bonded pentagons was observed.
Introduction
Among many unusual features displayed by liquid water, the anomalous behavior of thermodynamic response functions is one of the most interesting. The anomaly in the behavior of the response functions such as constant-pressure heat capacity ( Cp), thermal expansivity coefficient (aT),and isothermal compressibility ( K ~ ) especially . at temperatures below 0 O C , is believed to be due to the onset of some sort of cooperative phenomena, implying the possible presence of anomalous fluctuations in density.' Since the isothermal compressibility KT is a direct measure of density fluctuations, we shall concentrate our attention on this quantity. While for most liquids KT decreases when the temperature decreases, for liquid water KT displays a minimum at T = 50 "C. Below this temperature KT rises, especially very rapidly in the supercooled region. The described anomalous behavior of the compressibility is attributed to the superposition of "normal" fluctuations and anomalous fluctuations that diverge at temperature Ts.2 The nature of these anomalous fluctuations and their divergence at T, is not clear. Several theoretical models try to explain these anomalous density fluctuations. Stanley and Teixeira3 ascribe the cooperative nature of the anomalous fluctuations to the clustering of 4-fold hydrogen-bonded water molecules. The patches of such clusters are characterized by local order and reduced local density. The Stanley-Teixeira model is a percolation model with the percolation threshold for 4-folded molecules achieved at Ts. Another model is due to S t i l l i ~ ~ g ewho r , ~ proposed the existence of a significant concentration of bulky unstrained polyhedra in water. The cooperativity of the fluctuations is again due to the geometry; for example, a pair of dodecahedra can be stabilized by sharing the common pentagonal face. The model of Speedy5 is similar to Stillinger's. While Stillinger does not specify the particular polygonal structure, Speedy concentrates his attention on pentagonal hydrogen-bonded rings. These rings, according to Speedy, have a tendency to self-replicate, thus generating low-density patches. In Speedy's model the density fluctuations grow proportionally to the number of live-membered hydrogen-bonded rings until the limit of stability is approached at Ts. ( I ) Angell, C. A. In Wafer, A Comprehensive Treafise;Franks, F., Ed.; Plenum Press: New York, 1982; Vol. 7. (2) Speedy, R. J.; Angell. C. A. J . Chem. Pfiys. 1976, 65, 851. (3) Stanley, H. E.; Teixeira, J . J . Chem. f h y s . 1980, 73, 3404. (4) Stillinger, F. H . Science 1980, 209, 451. ( 5 ) Speedy, R. J. J . Pfiys. Chem. 1984, 88, 3364.
From the outset, computer simulations helped us to choose among the competing theories proposed to explain the experiment. Thus the first molecular dynamics simulation of water by Rahman and Stillinger6 clearly displayed that water is correctly described by a distorted hydrogen-bonding network or "continuum" model and not by its rival, the mixture/interstitial model. One may also hope that computer simulations will assist us in determining which model among the three described above can explain the nature of anomalous fluctuations in density. Indeed, simulations that address the question of the nature of density fluctuations at low temperature are found in the literature. The results are controversial at least. To support Stanley and Teixeira's theory Geiger and Stanley performed a MD simulation on ST2 water,' in which they observed the existence of low-density patches of water molecules.8 Rapaport performed a MD simulation of MCY water9 at -22 "C and did not observe any correlation between local density and the number of hydrogen bonds associated with each water molecule.'0 Speedy and Mezei in their work aimed to support Speedy's theory, reporting an increase in five-member rings concentration with a temperature decrease.]' Therefore they concluded that water anomalies may be related to self-replicating propensity of pentagons. They also have correlated the increase in pentagon concentration observed in "computer" water with the increase in the compressibility in "real" water. But what if the compressibility of the "computer" water (which of course may depend on the model used in the simulation) is not increasing when the temperature is decreased? To reach definite conclusions on the correlation between the polygon concentration increase and the compressibility increase, one has to calculate the compressibility of the "computer" water. We know of only one study where the direct calculation of compressibility and its temperature dependence was performed. This was done by Jorgensen and Madura in their Monte Carlo simulation performed on TIP4P potential.'* The reported values for KT in ref 12 are (67 f 13) X IO-'' Pa-' at T = 25 OC (the experimental value is 45.2 X IO-'' Pa-]), (24 3) X lo-" Pa-' at T = 0 OC (50.9 X lo-" Pa-' in (6) Rahman, A,; Stillinger, F. H. J . Cfiem. Phys. 1971, 55, 3336. (7) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 60, 1545. (8) Geiger, A.; Stanley, H. E. Pfiys. Rev. Left. 1982, 49, 1749. (9) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Pfiys. 1976,64, 1351.
(IO) Rapaport, D. C. Mol. Phys. 1983, 48, 23. ( 1 1 ) Speedy, R. J.; Mezei, M. J . Pfiys. Chem. 1985, 89, 171. (It)Jorgensen, W.L.; Madura, J. D. Mol. f h y s . 1985, 56, 1381
0022-3654/90/2094-8359$02.50/00 1990 American Chemical Society
