Acceleration of Convergence for Lattlce Sums

7320

J. Phys. Chem. 1989, 93, 7320-7327

Acceleration of Convergence for Lattlce Sums Naoki Karasawa and William A. Goddard III* Arthur Amos Noyes Laboratory of Chemical Physics,? California Institute of Technology, Pasadena, California 91 I25 (Received: January 6, 1989: In Final Form: May 10, 1989)

The lattice sums of nonbond interactions (electrostatic and dispersion) for computer simulations of periodic systems typically coverge very slowly. Here, we examine the accelerated convergence of these sums based on Ewald procedures and derive equations for the energies, forces, stresses, and curvatures. A method is proposed and tested for selecting the convergence acceleration parameter 1) on the basis of minimum calculation time. As an illustration and check, the properties of argon and NaCl crystals are calculated by using these equations and compared with values obtained analytically.

I. Introduction To carry out molecular mechanics and molecular dynamics calculations, it is necessary to sum various nonbonded interactions over all pairs of atoms. Thus, the electrostatic energy for a collection of point charges, {qi},is

where Qij = (Cunit/e)qAj and qi is the charge of the atom i. Here the prime indicates that i = j terms are excluded, e is the dielectric constant (which we take as l.O), and Cunit = 332.0647 puts the final energy in kcal/mol if distances are in angstroms (Cunit = 14.400 if the energy is in electronvolts). Similarly, the dispersion part of the vdw interaction has the form

where Bij is negative. These sums are notoriously slow to converge. This is illustrated in Table I for NaCl and polyethylene (PE) crystals. To obtain an electrostatic energy good to 0.01 kcal/mol by this procedure would require calculating all terms larger than 0.001 kcal/mol, which for unit charges would require a cutoff distance of 332 000 atoms/cm3 this would require A! With a typical density of 10ls atoms! Clearly such large cutoffs are untenable. The convergence is speeded by grouping together all atoms of a cell and summing complete unit cells as illustrated in Table 11. However, the convergence is still far too slow. The general solution to this problem originated with Ewald in 1921 using convergence functions for 1 / R interacti0ns.l Nijboer and de Wet@ generalized this approach to include all cases where the interactions are proportional to negative powers of distance for a single atom in a cell. Williams3 extended the formulas to allow multiple atoms in a cell. Consider a general lattice sum of the type 1 s,=-c 2 Iri - rj - RLlm

(4) converges much faster than (3). The second term converges slowly, but by taking the Fourier transform, the resulting sum (over the reciprocal lattice vectors) converges much faster. Following Nijboer and de Wette2 and Williams,3 we choose 4,,, as

11. Coulomb Sums For the Coulomb case ( 5 ) becomes &(r) = erfc ( r / g ) = 1 - erf ( r / g )

(6) where erf is the error function. In this case, the Coulomb interaction with infinite range is replaced by an interaction of range g and summed in the real space, while the long range corrections are summed in the reciprocal space. The parameter determines how much of the real space sum is converted to the reciprocal space sum (large g leads to a larger real space sum). A . Total Energy. The energy sum is

(7)

where it is assumed that each cell is neutral i

and the prime indicates that the term at the origin is excluded. Here h is the reciprocal lattice vector, is the volume of the unit cell, a = Iri - rj - RLl/g, b = 1/2hg,and h = lhl. The first term in (7) is just the first term of (4). The last term in (7) arises from the exclusion of i = j terms when L = 0, since in the reciprocal space sum, these terms are included. The second term of (7) arises from the expansion of the second term of (4) in terms of Fourier transforms of the point charge 6 function distribution 1 cQi,4(r - ri + rj + RL) = -CS(h)S(-h)Pr (sa)

L,ij

where ri and r are the basis vectors of atoms i and j in the cell and RL is the {attice translation vector. The sums over i and j each go over all atoms inside the cell except that i # j when L = 0. The total electrostatic energy is given by S1, while the dispersion term is given by S6. Multiplying every term in (3) by the convergence function 4, and then by (1 - +,), we obtain

s,

=

1

-I3 2

L.ij

Aij4m(Iri

- rj - RLI)

Irj - ri - RLlm

1

+ ?LiJ X

a b

iJ,L

Aij

(3)

Aij[l - 4 A r i - rj - RA)] Iri - r, - RLlm (4)

If &(r) is a rapidly decreasing function, then the first term of 'Contribution No. 7902.

0022-3654/89/2093-7320$01.50/0

(8)

cqj = 0

where

is referred to as the structure factor. The quantity S(h)S(-h) in (7) has the form S(h)S(-h) = c Q i j cos [h.(ri - rj)] (10) ij

which Cowley et aL4 have shown can be rewritten as Cunit S(h)S(-h) = -([xqi cos (h-ri)]* + [Eqisin (h.ri)12] (1 1) t

i

i

(1) Tosi, M. P . Solid Stare Phys. 1964, 16, 107.

(2) Nijboer, B. R. A.; de Wette, F. W. Physica 1957, 23, 309. (3) (a) Williams, D. E. Acta Crysrallogr., Sect. A 1971, 27, 452. (b) Williams, D. E. In Ctysral Cohesion and Conformational Energies; Metzger, R. M., Ed.;Springer-Verlag: Berlin-Heidelberg-New York, 1981; pp 3-40.

0 1989 American Chemical Society

Acceleration of Convergence for Lattice Sums

The Journal of Physical Chemistry, Vol. 93, No. 21, 1989 7321

TABLE I: Convergence of Nonbond Interactions Using Atom-Based Cutoffs" NaCld

8,

A

Rinner? 8 9 10 14 15 19 20 24 25

10 10 10 15 15 20 20 25 25 EwaldC

polyethylene

terms

EQ

Edisp

608 608 608 2152 21 52 5296 5296 10320 10320

-98 1.65 -1556.20 -1492.67 1446.83 1281.61 -1743.32 1207.33 1210.98 109.11 -824.59

-20.330 -20.41 3 -20.445 -20.726 -20.735 -20.794 -20.800 -20.821 -20.822 -20.846

terms 3072 3072 3072 10394 10394 246 15 246 15 48015 48015

EQ

Edisp

-5 1.870 -2.343 -111.998 -40.272 -52.972 -79.905 -63.080 -107.974 -54.830 -65.685

-797.687 -797.798 -791.863 -798.184 -798.198 -798.275 -798.280 -798.306 -798.308 -798.338

'Energies are in kcal/mol. bEach term of (1) or (2) is multiplied by the cubic cutoff function S(Rij), where S(R) = 1.0 if R < Rinoor,S(R) = 0.0 if R > R,,, and S(R) = [R20uter- R2]2[R20ul, 2R2 - 3R2i,,,,er]/[R20u,r - RZinnJ3, otherwise. c v = 2.5 for NaCl and q = 2.0 for polyethylene. mol-'. qNa = e, qcl = -e, A = 5.63 A. dParameters used are BNaNa = -24.180, BNaCI= -161.20, Bclcl = -1669.58 kcal

+

Ab

TABLE II: Convergence of Nonbond Interactions Using Cell-Based Cutoffs' NaCP

'4

'CUI.

EQ -824.5563 -824.5563 -824.5563 -824.5768 -824.5836 -824.5852 -824.5861

termsb 3500 (2 2 2) 3500 (2 2 2) 3500 (2 2 2) 9604 (3 3 3) 20412 (4 4 4) 37268 (5 5 5) 512 224od 296 + 2240"

6 8 10

15 20 25 EwaldC

polyethylene

+

Edap

-20.761 5 -20.7615 -20.7615 -20.8 161 -20.8323 -20.8387 -20.8463

termsb 6930 (1 3 2) 14850 (2 4 2) 20790 (2 4 3) 54054 (3 6 4) 86394 (3 8 5) 162162 (4 10 6) 1748 1584d 1044 6204d

+

+

EO

E,+*"

-65.6032 -65.6402 -65.6521 -65.6655 -65.6669 -65.6725 -65.6836

-798.0302 -798.1933 -798.23 12 -798.3005 -798.3 176 -798.3274 -798.3383

'Energies are in kcal/mol. bThe numbers in parentheses show the ( a b c ) index for the last shell considered. = 2.5 for NaCl and v = 2.0 for polyethylene. 89 = 0.001 kcal/mol. = IO4 kcal/mol. dThe number of terms in real space and reciprocal space sums, respectively (eq 10 is used). #Parameters are given in Table I.

Since the summation in (1 1) runs only over single atoms rather than pairs, the computations for the reciprocal space sum are significantly reduced. B. Force and Stress. The force a t atom p is given by

--as1

reciprocal lattice and real lattice vectors, and H i s the transpose of H . The stress IIaB is given by the relationS

where E is the total (potential) energy of the system. For the electrostatic stress, we have

a r.

4* ' ,Fh[EQpi

sin (h-(rp - ri))]h-2e-b'(12)

i

where up = Irp - ri - RLI/q. By use of (1 l), the second term can be rewritten as Cunit 4~ ' -q,Eh(sin t f l h

(h.r,)[Cqi cos (h-r,)]- cos (her,)

1 = -(&-~HHH~-I

X

i

[Eqisin (h.ri)])h-2e-b'(13) i

To calculate the stress, consider the matrix H that contains the real space unit cell vectors in Cartesian coordinates, i.e., H = [a,b,c],where a,b,c are the unit cell vectors.5 Then we have

fl = det(H) ri - rj - RL = H(si - sj - L)

(14)

h = 2rk1n

(16)

(15)

and

where si is the fractional coordinate vector of the atom i, and n and L are vectors whose elements are the integers specifying the (4) Cowley, E. R.;Jacucci, G.; Klein, M. L.; McDonald, I. R.Phys. Rev. B 1976, 14, 1758. In this paper, A,(T) in eq 3 must be multiplied by a factor (1/4.

(5) Nose, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055.

The stresses can also be obtained by differentiating the energy with respect to strains. The strain tensor is defined by6

2

- 1)

(19)

where H o contains the original cell vectors and H contains the deformed cell vectors. t is symmetric so we define six independent strain components ei such that e, = e l l , e2 = t22, e3 = t33re4 = 2q3, es = 2e3,,and e6 = 2tI2. Then we have

where II, = 1111, 112 = 1122,113 = II,,, 114 = II23, IIs = II,,, and n6 = II12. To derive this equation, we use the relations

= SiaHiy

(ei = 0 )

(21a)

(6) Parrinello, M.; Rahman, A. J . Appl. Phys. 1981, 52, 7182. (7) Kittel, C. Introduction to Solid State Physics, 5th ed.; John Wiley & Sons: New York, 1976. (8) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: London, 1954.

Karasawa and Goddard

7322 The Journal of Physical Chemistry, Vol. 93, No. 21, 1989

-m, - - CHkL'(HOiyH0jk+ HOjyHoik) aeij

k

= Gj,EliT

+ GiaHj7

in the reciprocal space sum of (24) to reduce the computing time. B. Force and Stress. The force on atom p is

(ei = 0 )

(21b)

aS6

= -d'P 1 = - CCBpi(rp- ri - RL) X

f6,

derived from (19). C. Curvatures. The second derivatives with respect to XX and YY strains have the form

q8

+

-

(9

L i

( 6 ~ ; ~ 6apd

- rj - RL)14+ 2F(ri - rj -

sin [h-(rp - ri)]]h3[ n1I2erfc (b)

1

4HhI2 + G (22a) J

nn,,

- C.CBij(6a-8 + 6ad

-

(27)

+ 3 c 4 + a-2)e-a2(ri- rj - RL),

X

L ij

(ri - rj - RL)o

1

(5

=

2q8

+ H(h12+ h,Z) + G

+

j

The stress term is 1

(ii)

2n '

,912 + 3ap4 + ap-2)e-u2+ xh[CB X 12R h "

n3~2

+CCBijcos [h.(ri - rj)] X 24Q b i j 1

(22b)

where F = [erfc ( a ) + 2~-1/2ae-0~]/~3

G = v2e-@/4b2 H = -q4(l

Here we have used

+ b2)e-b2/8b4 where ri - rj - RL is given in (15) and

111. Dispersion Sums For the dispersion sum with m = 6 in (3), the form for & is $6

[

= 1

+ (r/q)2 + i(r/v)4]e-(r/q)2

(23)

C. Curvatures. The second derivatives with respect to XX and YY strains have the form (9

A. Energy. This leads to an energy of

a2s6 =-- 1 &$e, erfc (b) l)e-@] b

1 + 6Qq3 CBij- -12$ xBii ij

I

1 + (G -

- rj -RL)14

2q8 L , i j RL)12

+ 2F(ri - rj -

n3~2 +CCB,, cos [h.(ri - rj)] 240 b

X

j j

(24)

i

where b = 1/2hq, a = Iri - rj - RLl/q, and h = IbI. The first term in (24) arises from the first term in (4) and can be written as

(ii) d2S6 =-- 1

%de,

2$

- rj -RL)12(ri- rj - RL)Z2

L,ij

n3/2

'

CCBijcos [h-(ri - rj)] 24Q h

X

j j

where a, = lRLI/7. The second term arises from the second term in (4) using (9) and (10). The third term is from h = 0 in the second term, while the fourth term is from exclusion of i = j terms when L = 0. Often the dispersion terms are assumed to satisfy the geometric combination rules -Bij = (BiiBjj)l12

(26)

and expressions similar to (24) are based on this form.3 However, in the above derivation, no such assumptions are made. When (26) is satisfied, we can use XBij cos [h-(ri - rj)] = i j

- [ x ( [ B i i l ) 1 /cos 2 ( h ~ ~-) [Z(lBii1)1/2 ] ~ sin (h.ri)I2 (26') i

i

i$h12h22H'+ H(h12+ h2')

+G

where F=(6d

+ 6ad + 3a4 + u-2)e-u2 erfc (b)

G = q3

+

(&

e(n1l2erfc (b)

H= 9

dF F'= da'

-

- i)e-b2}

e-bz) b

H I = -d H

db2

1

+



we obtain

TABLE 111: Predicted Values of tiw (Using (34)) crystal NaCl (8 atoms/cell) NaCl (2 atoms/cell) Ar (4 atomslcell) ‘ Ar (1 atom/cell) polyethylene (12 atoms/cell) orthorhombic poly(oxymethy1ene) (16 atoms/cell)

IRL,minl,

A

Ihminl,

5.63 3.98 5.31 3.755 2.54 3.56

A-’

1.116 1.933 1.183 2.050 0.8656 0.8212

?opt*

A

3.18 2.03 3.00 1.91 2.42 2.94

IV. Selection of tl A . Choice of q from Cell Parameters. The parameter q should be chosen to optimize the convergence of both the real lattice and reciprocal lattice sums. If 9 is too small, the reciprocal lattice sum does not converge quickly since b = 1/2h7remains small when h becomes large. Similarly, if 7 is too large, the real space sum does not converge quickly since a = Iri - rj - RLI/? remains small for large IRLl. Since both sums are multiplied by terms like e-$ or erfc (x) (where x = a for real space and x = b for reciprocal space), a simple estimate of the optimum parameter 9 for rapid convergence of both sums can be obtained as follows. These terms decrease rapidly for large x , so that both sums converge similarly when a = b. This leads to

which determines 71 in terms of L and h. The simplest choice of is to use the minimum length of RL and of h so that we have the same orders of magnitude of a and b when these have the minimum values (of course, IRLl # 0 and h # 0). Then

Cunit

xCQ

v2 erfc

(5)

(36)

(b) Reciprocal Space Sum. The error in the reciprocal space sum due to a cutoff at H,,, is given by

2a

I

e-b2

E Q . ~=~OI CS(h)S(-h)--B(h h h2 - H,,J

(37)

Replacing the sum by an integral and replacing S(h)S(-h) by , obtain Cunit p ( q 2 ) / c we

Using (36) and (38), we can evaluate the cutoff distances R,,, and H,,, to obtain a given accuracy 6Q (for given q ) . Because of the neutrality of the cell, there will be a great deal of cancellation in (35) and (37) that is ignored in (36) and (38). Consequently, eq (36) and (38) generally overestimate the error by a factor of 10 or so. 2. Dispersion Sums. ( a ) Real Space Sum. The total error in the real space dispersion sum is given by

7

?opt2

lRL,minl =-

(33)

(39)

Mmin

Assuming a C b,c and h, qopt2

h,,hb (33) leads to aQ = Q =.rrla X bJ ?rb sin y

Approximating this by an integral, we have

C

-

(34)

This determines ‘lop,in terms of the crystal structure. The ‘lop, thus obtained may not be the optimum g for minimizing computation time, since we assumed that the computation time per each term is the same for both real and reciprocal space sums. Table 111 shows the calculated topt for various systems. B. Accuracy Specified Cutofls. With q specified there is still an infinite number of terms in the sums over the real space and reciprocal space lattices, and we use an accuracy criteria to specify limits on these sums. To this end we specify a tolerance 6 and carry out the sums until the neglected terms have a total contribution smaller than this tolerance. 1 . Electrostatic Sums. ( a ) Real Space Sum. Using a cutoff distance R,,, introduces an error in the total energy for the real space sum of

where p ( B ) = CiJBijI is the sum of coefficients of dispersion terms. We approximate the above equation as

( b ) Reciprocal Space Sum. The total error in the reciprocal space sum is given by

where 0(x) is the step function [0(x) = 1 when x > 0 and 0(x) = 0 when x C 01. To estimate this error, we replace discrete sum (35) by a continuous integral. Defining the average interaction as ( q 2 ) = q?/N, we obtain

EL,rccip = - CCIBijlcos [h.(r, - rj)]h3 24Q b ij

Since

Approximating this sum as an integral leads to

xi

a3/2

’

7324 The Journal of Physical Chemistry, Vol. 93, No. 21, 1989


By using all2 erfc (b) I (l/b)e-b2, we have

n

Using (40) and (42), we can calculate the cutoff distances R,,, and H,,, that lead to a total error 6v for given 7. 3. Repulsive Terms. The repulsive terms in simulations generally have the form ' EI2= -1 xA..R.,-12 2 i j 11 '1

(43)

1 ' E,, = - ZAi,e-'URf/ 2 ij

(44)

or

We can estimate the R,,, for these cases in a similar way. The results are (45)

where M ( A ) = &Aij and ( A ) is determined by the relation M(A)e-(V&t = xAije-hi/Rcut

Q

and similarly for Lb, L,. The sum over the reciprocal lattice has limits of Ha = (a/Zn)H,, and similarly for Hband H,. Within this set of cells, we eliminate any terms for which R > &,,. This leads to a considerable reduction in effort, as indicated in Table IV. C. Time Minimized 7. Since the choice of 7 affects the R,,, and H,,, and thereby the number of terms in real and reciprocal space sums, we can choose 7 so as to minimize calculation time while retaining prespecified computational accuracy. This is illustrated in Table V for the case of poly(oxymethy1ene). In this table the electrostatic energy using (1 1) (a2) as well as those using (10) ( a l ) are shown. On the basis of several test calculations, we estimate that trcal/trccipi= 4 in ( a l ) and (b) and treal/trccip= 10 in (a2), where trealis the time per term for the real space calculation, while tAp is the time per term for the reciprocal space calculation. Of course, these values will depend on the actual program and computer used. The reason that the evaluation of the reciprocal space term is much faster than the real space term is that the reciprocal space lattice sum involves a factor that only depends on h, which is calculated only once, while in the real space sum, all factors depend on atomic distances and must be calculated for each pair. Excluding the interactions between the same atoms in different cells (since there are easy to evaluate) and using (lo), the number of terms in the real space and reciprocal space sums can be estimated as follows (assuming constant atomic density) (47) (9) Calculations optimizing the cell parameters and internal coordinates of crystals were carried out using the PolyGraf polymer simulation program from Molecular Simulations, Inc. (BioDesign Division), Pasadena, CA 91 101.

(48)

where N is the number of atoms in the unit cell. When (1 1) is used, N(N - 1)/2 in (48) should be replaced by 2N. Since there is a cancellation of terms in the electrostatic sum due to the charge neutrality in the unit cell, we use 6, = lobdispto obtain similar accuracy in both sums. By using the estimated number of terms and the values of tr,,/trccip, we select 7 so that N,,,, + Nr,/m is minimized where m is the ratio treal/treci. In Table V we show the results for poly(oxymethy1ene). In this system, there are 24 centers of which 16 are charged. In ( a l ) and (b), the 7 leading to the fastest calculation is 7 =2.5 A, while in (a2), the 7 leading to the fastest calculation is r ) i= 2.0 A. The r) predicted using (33) is 2.94 A, in reasonable agreement with ( a l ) and (b). In (a2), it is shown that the computing times are shorter than those in ( a l ) for any 7 due to the smaller number of terms in the reciprocal space sum. Also, the computing time for each term in the reciprocal space sum is shorter in (a2) than in (al). Therefore, the time-minimized 7 in this case is smaller than for other cases. The accuracy here is 10" kcal/mol for electrostatic and lo4 kcal/mol for dispersion. Note that the error in the energy is independent of 7 in these calculations, showing the effectiveness of accuracy-specified cutoffs.

V. Application to Ar Crystal As a test case for the Ewald dispersion sum formlas, we will calculate the structure and properties of argon face-centered cubic (fcc) crystal analytically and compare with those obtained by using the lattice sums. A. Model. We will describe the Ar-Ar interactions with the Lennard-Jones 12-6 form

ij

4. Implementations of the Limits. The limit on the real space sum, L, is chosen so that all cells containing any atoms within R,,, of any atom in the unit cell are i n c l ~ d e d . ~This leads to Rcu,bcsin CY +1 La =

N ( N - 1) 2

U(R) = DO[

(

$)12

- 2( $)6]

(49)

where Ro = 3.8666 A, and Do= 0.2351 kcal/mol are adjusted to reproduce the experimental lattice spacing and heat of vaporization. For an fcc crystal structure, the total energy E,,, is7

where, R is the nearest-neighbor distance, N = 4 is the number of atoms in the conventional fcc cell, p12= 12.13 1 88, and p6 = 14.453 92. B. Equilibrium Lattice Constants. Requiring aE,,,/aR = 0 leads to the equilibrium nearest-neighbor distance

The lattice constant at equilibrium becomes A , = 2'I2R, = 1.37353Ro = 5.3109

A

(51b) the experimental valuelo at 0 K. C. Cohesive Energy. Substituting (51) into (50) leads to 1 P62 2 Pl2 so that the cohesive energy (per atom) is UtJR,)

= - -NDo-

U,, = 8.6102D0 (53) This leads to U,, = 2.0244 kcal/mol (of atoms), the experimental heat of vaporization at 0 K (after correcting for zero point energy).' ( I O ) Donohue, J. Structures of the Elemenzs R. E. Krieger Publishing Co.: Malabar, FL, 1982. For argon, A = 5.3109 A at 0 K. (1 1) Hultgren, R.; et al. Selected Values of the Thermodynamic Properties ofthe Elements; American Society for Metals: Metals Park, OH, 1973. For argon, AHv = 1.848 kcal/mol at 0 K. Estimating the zero point energy as 0.1764 kcal/mol (from our calculations) leads to Umh= 2.0244 kcal/mol.



TABLE I V Number of Terms in the Sum for Given Ewald ( v ) and Accuracy (6) Parameters

7,

8,

aQ, kcal/mol

R,,,, 8,

10-2 10-3 10-4 10-2 10-3 10-4 10-2

4.645 5.152 5.6 17 8.127 8.938 9.686 11.721 12.828 13.852

1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5

I 0-3 10-4

H,,,

(a) Electrostatic Sum NaCl real terms recip termsu

4.568 4.983 5.366 2.683 2.934 3.173 1.891 2.073 2.243

72 104 104 296 512 608 1088 1328 1616

7168 10192 12880 1568 2240 2576 504 728 896

polyethylene H,,, A-1 real terms

R,,,, 8, 4.226 4.773 5.269 7.460 8.332 9.126 10.814 11.999 13.084

4.090 4.549 4.964 2.394 2.670 2.922 1.678 1.879 2.061

recip terms'

212 344 468 1264 1748 2342 3820 5322 6824

7128 9372 11088 924 1584 2904 528 660 660

(b) DisDersion Sum 7,

8,

1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5

hap,

kcal/(mol 8,) R,,,, 8, 10-3 10-4 I 0-5

lo-' 10-4 10-5 10-3 10-4 10-5

H,,,,

4.819 5.304 5.752 7.453 8.310 9.096 9.872 11.122 12.262

NaCl real terms

5.089 5.485 5.856 2.695 2.966 3.217 1.740 1.954 2.143

72 104 104 296 296 512 608 896 1232

recip termsd

R,,, 8,

Hat, 8,-l

IO864 I3552 I7304 1568 2240 2576 504 728 728

4.867 5.348 5.793 7.539 8.387 9.168 9.997 11.235 12.366

5.001 5.410 5.787 2.639 2.915 3.167 1.690 1.910 2.105

polyethylene real terms 372 480 628 1302 1816 2390 3066 4326 5854

recip terms' 1 I880 16104 20328 1584 2376 3432 528 660 660

'The number of atom pairs in the sum. D. Bulk Modulus. The bulk modulus, a measure of the stiffness of the crystal, is defined by (54)

+

Here the sum is over all integer values of 11, 12, l3 where Il l2 + l3 must be even. The constants s1,8 = 0.031 84001, s2,8 = 0.015673 32, s ~=,0.320449, ~ and s2,5= 0.140899 were obtained by numerical computations. The result is

Writing the volume of the crystal as V = NR3/Z1l2, we have

CI1= 148.762(

(55)

CI2= C4,= 85.123(

where b12 = (l/8)p12NsD&~'2and b6 = (l/2)p@D&06. The equilibrium volume is V, = ( 1/21/2)NRo'(P12/p6)11z, leading to The Poisson ratio

106.328($)

-

(56)

E. Elastic Constants. The elastic constants have the form8 I

CI,= -

$) $)

u

for this cubic system el

=

= 0.3640

'I2

CI1 +

Cl2

is independent of Ro and Do. If Do is expressed in kcal/mol and Ro is expressed in A, then a conversion factor of 6.947 80 is needed to obtain pressures in GPa. Thus the above expressions lead to

CI1= 4.2031 GPa CI2= 2.4051 GPa

and and (For a two-body potential such as (49), the shear constant C, is equal to Clz.) This leads to

1 B = -(CI1 + 2CI2) = 3.0044 GPa 3 These can be compared to the experimental values at 0 K ofI2

CI1= 4.0 GPa CI2= 2.0 GPa

C4,= 2.0 GPa where lattice sum constants s ~ and , s~ 1 4

~ are , given ~

by

F. Results. In Table VI we compare the results of optimizing the unit cell parameters of Ar with various ways of carrying out the lattice sums (all of which used periodic boundary conditions). In these calculations we used conjugate gradient minimization to optimize both cell parameters and atomic coordinates inside (12) Meixner, H.; Leideren, P.; Berberich, P.; Liischer, E. Phys. Left. A 1972, 40, 257.

7326 The Journal of Physical Chemistry, Vol. 93, No. 21, 1989 TABLE V Dependence of Timing on Choice of r) for Orthorhombic Poly(oxy methylene)' ( a l ) Electrostatic Sum. 6, = 0.001 kcal/mol. Equation IO is used. Use of 6Q = kcal/mol leads to E = -51.19144 kcal/mol. 7, A

1.0 1.5 2.0 2.25 2.5 3.0 3.5 4.0 4.5 5.0

E, kcal/mol -51.1916 -51.1914 -51.1916 -51.1915 -51.1915 -51.1914 -51.1915 -51.1914 -51.1914 -51.1914

N I,, 116 443 1154 1656 2306 4232 6928 10830 15795 22136

N,,

trcalltrsip

tmt

tEO~c, s

92160 27360 10320 6720 4800 2880 1440 960 480 480

8.3 6.9 4.4 5.9 5.1 4.6 4.3 3.9 4.4 4.0

23156 7283 3734 3336 3506 4952 7288 11070 15915 22256

7.38 2.25 1.08 1.11 1.13 1.47 2.02 2.67 4.35 5.55

(a2) Electrostatic Sum. 6Q = 0.001 kcal/mol. Equation 11 is used. 7,

A

0.75 1.0 1.25 1.5 1.75 2.0 2.25 2.5 3.0

E , kcal/mol -51.1919 -51.1916 -51.1915 -51.1914 -51.1914 -51.1916 -51.1915 -51.1915 -51.1914

N,I 36 116 234 443 741 1154 1656 2306 4232

N,,

tr,~/trsip

Lt

tCP~c, s

61312 24576 12608 7296 4160 2752 1792 1280 768

35.3 13.6 12.7 12.1 8.5 7.3 9.4 7.8 7.8

6168 2573 1494 1172 1 I57 1429 1835 2434 4308

3.10 1.32 0.74 0.57 0.49 0.50 0.75 0.86 1.32

(b) Dispersion Sum. 6djrp = IO4 kcal/mol. Use of &p = IOw5 kcal/mol leads to E = -8085.19405 kcal/mol. 7,A

E , kcal/mol

1.0 1.5 2.0 2.25 2.5 3.0 3.5 4.0 4.5 5.0

-8085.1939 -8085.1939 -8085.1940 -8085.1940 -8085.1939 -8085.1940 -8085.1940 -8085.1940 -8085.1940 -8085.1940

N I,, 440 1479 3209 4261 5656 9136 13741 19501 26581 34738

N,,

tdtrecip

tmt

416760 104880 35880 23736 15456 8832 4416 2208 1104 552

6.6 4.0 2.6 3.8 3.2 3.1 2.8 2.5 2.5 2.7

104630 27699 12179 10195 9520 11344 14845 20053 26857 34876

tcalC,s 33.93 8.39 3.17 1.78 2.29 2.41 2.83 3.15 4.04 6.03

' E is the energy, N,,]and Nra are the numbers of terms in the real space and reciprocal space sums,tr,l/trsip is the actual ratio of calculation time per term for real and reciprocal space, tmt is the estimated relative time, and talc is the actual calculation time. The estimated N,,,/m, where m = 4 in ( a l ) relative time is calculated from and (b) and m = IO in (a2).

+

the celL9 In each case the final root mean square (RMS) force per atomic degree of freedom is less than 0.001 kcal/mol, and the RMS stress for the six cell parameters is less than 0.0001 kcal/mol. Using traditional distance cutoffs with R,, = 9 8, leads to an error in the lattice constant of 0.006 8, or 0.1%, an error in the cohesive energy of 0.68 kcal/mol (of cells) or 8.4%, and an error in the bulk modulus of 0.09 GPa or 3.0%. In the Ewald calculations, we evaluated the repulsive terms by using direct sums and used the error bounding procedures of section 1V.B with several energy criteria, where 6, = 106,. We also carried out calculations as a function of strain by introducing a finite strain to the system, optimizing the atomic coordinates, and calculating the stresses9 The highly repulsive nature of the atomic interactions for short R (due to the Pauli principle that requires orthogonalization of the overlapping atomic orbitals) leads to a quite nonlinear stressstrain relation. Using the stress versus strain results for small strains (&0.005, &0.01), we obtained numerical estimates of CIIand CI2to compare with the analytic values. The results agreed to two decimal places in each case. VI. Application to NaCI Crystal As a test case for the Ewald Coulomb sum formula, we will calculate the structure and properties of NaCl crystal analytically and compare with those obtained by using the lattice sums. A . Model. For the atom-atom interactions, we use the form


uij= Xe-R/p - " R

f-

q2

(nearest neighbors)

(62)

(otherwise)

PijR

The total energy E,,, (relative to free ions) is given by7

Here R is the nearest neighbor distance, N = 4 is the number of pairs of NaCl atoms in the unit cell, and z = 6 is the number of nearest neighbors. The Madelung constant is a = 1.747565, while unit charges on each ion lead to q2 = 332.0647 A kcal/mol (allowing R to be in 8, and E in kcal/mol). B. Parameters. We will choose the parameters X and p in (62) to obtain the experimental lattice constant13 A = 5.578 8,

and cohesive energyi3 E,,,(R,) = -740.0 kcal/mol

at 0 K. The equilibrium nearest neighbor distance Re = (1 / 2 ) 4 is obtained by requiring that dE,,,/dR = 0, leading to

Substituting Re into (63) leads to

Using these results leads to p

= 0.309 223 8,

X = 3.84485@c/p = 31 765.8 kcal/mol

C. Bulk Modulus and Elastic Constants. The volume of the crystal Vis given by V = 2NR3; hence

(

-(

E,,,( V) = N(zX exp[ & ) 1 ' 3 / p ]

- $)'I3

aq")

(66)

Using eq 54, we have B = ( VdZEtot z ) ,

=

$(f i ) -

= 25.9858 GPa

(67)

where the equilibrium volume V, is 2NRd and the conversion factor 1 kcal/(mol A3) = 6.947 80 GPa was used. This can be compared with the experimental value ofL326.60 GPa. Since each atom in the NaCl crystal is at a point of inversion symmetry, no internal strain is induced when external stress is applied. Hence, the two independent elastic constants are given by

where e, and e2 are strains (XX and YY components) and the first term in CI1arises from the repulsive part of the potential. The lattice sums for the second partial derivatives of the elec(13) (a) Sangster, M. J. L.; Schrder, U.; Atwood, R. M.J. Phys. C 1978, Zl. 1523. (b) Sangster, M. J. L.; Atwood, R. M. J . Phys. C 1978, ZZ, 1541.



TABLE VI: Equilibrium Properties of Argon Crystal (Using Four Atoms per Cell)

A, A U, kcal/mol B, GPa C,,, GPa C12, GPa C4,, GPa d

exper' 5.3109 -8.0975 2.67 4.0 2.0 2.0 0.32

0.00001b

exact 5.3109 -8.0976 3.0044 4.2031 2.405 1 2.405 1 0.364

5.3109 -8.0970 3.0044 4.2031 2.4051 2.4051 0.3640 0.171 120 480

time! s pairs R spaced pairs K spaced

Ewald 0.001b 5.3108 -8.098 1 3.0054 4.2048 2.4057 2.4057 0.3639 0.143 72 336

0.01b 5.3090 -8.1171 3.0164 4.2216 2.4138 2.4138 0.3638 0.084 24 156

9Ac 5.3172 -7.4202 3.0940 4.4123 2.4348 2.4348 0.3556 0.045 120

direct summation 15 Ac 5.3176 -7.9756 2.9482 4.1309 2.3569 2.3569 0.3633 0.124 480

30 Ac 5.3111 -8.0830 3.0050 4.2052 2.4049 2.4049 0.3638 0.866 4320

'See ref 10 (A), 11 (U), and 12 (C,,, C,z, and C4 ). bAccuracy parameters in kcal/mol. 'Cutoff distances (Itc,,). Using a cubic spline function to decrease the potential from full value at R,,,-1 to zero at R,,,-0.5 A. dOnly attractive terms are considered.

1

trostatic energy with respect to the strain components converge very slowly, requiring Ewald sums. By using (22), we have

-a2s -I -

I l:

(ri - rj - RL)I2

+ - CS(h)S(-h)

"See ref 13. bAccuracy parameters in kcal/mol.

X

2* ' (ri - rj - RL)12(ri- rj - RL)2Z+ - zS(h)S(-h) o

2-

(1

b

+ bz) h2

In this model, C,, = CI2. The calculations of these sums give Cll = 51.4366 GPa C12= C4, = 13.2598 GPa 1 B = -(Cll 3

TABLE VII: Equilibrium Properties of NaCl Crystal (Using Four Molecules per Cell) Ewald exuef exact O.lb 0.01b 0.001b 5.5116 5.5780 A, A 5.578 5.5780 5.5776 U, kcal/mol -740.0 -740.00 -739.9938 -740.0184 -739.9985 B, GPa' 26.60 25.9858 25.97 25.963 25.9854 C,,, GPa 57.33 51.46 51.433 51.4366 CI2,GPa 11.23 13.23 13.228 13.2598 C,, GPa 13.31 13.23 13.228 13.2598 U 0.164 0.205 0.205 0.2050

+ 2CI2) = 25.9854 GPa

while the Poisson ratio (61) is c = 0.2050 D. Result. In Table VII, we show the results of optimizing the unit cell parameters of NaCl with three different accuracy parameters.

VII. Summary and Discussion Using Ewald-inspired approximations for accelerated convergence of lattice sums, we developed and tested equations for the energies, forces, stresses, and second derivatives for both electrostatic (1/R) and dispersion (1/R6) lattice sums. In addition, we developed an approach for estimating the convergence parameter 7 in order to minimize the computation time while retaining a fixed level of accuracy. With these accuracy specification procedures, the costs of carrying out accurate lattice sums is less than that for normal direct sums (at the same level of accuracy) despite the more complicated formulas for accelerated convergence. We suggest that these procedures may also prove equally valuable for biological systems. A molecule such as hemoglobin with -6000 atoms would lead to 18000000 pairwise interactions, which are truncated to -500000 by using energy cutoffs of 9 A. We suspect that such cutoffs lead to errors in the Coulomb interactions similar to those shown for NaC1. By considering the hemoglobin to be in a unit cell sufficiently large that interactions between cells is small, one could use the Ewald procedures to generate a given level of accuracy. This approach would be especially valuable for including explicit solvent (water) in the calculation. Water has large charges (0.4 e on each H and -0.8 e on 0) so that use of a finite solvent shell leads to very large surface effects. With a periodic cell containing H b and H 2 0 , one could eliminate surface effects while obtaining accurate energies. Acknowledgment. This work was partially supported by a grant from Imperial Chemical Industries, Cleveland, England, and by a grant from the Air Force Office of Scientific Research (No. AFOSR-88-005 1). We thank Molecular Simulation, Inc. (Biodesign, Inc.), Pasadena, CA for use of the PolyGraf polymer simulation program. Registry No. NaC1, 7647-14-5.

Acceleration of Convergence for Lattlce Sums

Recommend Documents