- Email: [email protected]
Convergence of the bertaut series and the calculation of the electrostatic energy of an extended crystalline lattice
VoIume 62, number 3
CHEhlICAL PHYSICS LETTERS
CONVERGENCE
OF THE RERTAUT
15 April 1979
SERIES AND THE CALCULATION
OF THE ELECTROSTATIC ENERGY OF AN EXTENDED CRYSTALLINE
LATTICE
H.D.B. JENKINS and K.F. PRATT ScItooi
Of Cknistiy and Moleculnr Sciences. UIziiersiry of Warwick. Co&enrr_r.WestMidiands, UAF
Received 10 Jnmury 1979; in tin21 form 16 hiarch1979
Computations are reported wbicb enable the prediction of the necessary truncation of the Bertaut series for urious cbrge density functions consistent with various degrees of conwrgence. The results are sufficiently e.rtensive to be appliable to nlmost all httices-
I_ Introduction Studies designed to obtain accurate lattice energies and related thermochemical data for various ionic salts having complex ions and subsequent extension of this pork into the mineraIogica1 field here led to a situation where the studies of Jenkins [ 1,2] on the truncation of the Bertaut Madelung series have become insufficiently extensive to predict truncation limits. The present work seeks to update and extend these studies so as to apply to almost alI crystai lattices_ The series developed by Bertaut i3] adopts a summation technique which operates entirely over the reciprocal, rather than the real, lattice. The computations necessary to evaluate this summation are extensive for all but the simplest lattice- It is therefore desirable to be able to predict, with reasonable accuracy, the truncation point which can be used to terminate this series, for particular degrees of convergence. Templeton [4] began investigations by consideration of the truncation of the series when using 3 uniform charge density function to represent the point charges (see section 2). His work however covered a range of lattices which included only simple ionic salts- Further work by Jones and Templeton [S] included consideration of other charge density functions, but stili on!y for this limited range of lattices- The later work by Jenkins [I .2 J covered a range of lattices greater than that of Templeton, but still insufficient to include many lattices which are commonly encountered- This work examined truncaticn properties specifically for a “parzlbolic” charge density function_ The development of a completely general computer program for Iattice energy calculation (LATEN [6]) has necessitated 3 general approach to this subject_ Different charge density functions are considered for a range of possible series convergence and results cttlculated to cover virtualiy all possible ionic lattices- Convergence is investigated in terms of the etectrostatic energy rather than the Rladelung constant, as has been the case previously IlA51.
The electrostatic
energy of an ionic lattice (Uelec) is given by the Bertaut approach
[3] as
N u e,ec = W+2Z) 416
,g
nic$
- UWnVZ)
F
I FI~129;J~z2 - U,, •i- i.J,, ,
(1)
Volume 62, number 3
CHEMICAL
PHYSICS LlZl’-lTRS
15 April 1979
where K is the conversion Factor from e2A-1
to lcJ mol-’ (taken as 1389.30 in this work), g1 is a constant dependent on the charge density function chosen to represent the spherically symmetrical chrrrge distribution (see below), R is the cut-off radius of the charge density function, 2 is the number of molecules per unit cell. N is the number of types of atom in the lattice, ni is the number of the jth atom type in the unit cell, of atomic charge qi. r/is the unit cell volume, h is a reciprocal lattice vector, Fjr is the coulombic structure factor, orb the Fourier transform of the charge density function, U,, the energy of overlap of the atomic charge distributions (chosen to be zero for all but the gaussian function) and Use, the self-energy of the compIax ions in the iattice (simpIy the electrostatic energy terms internal to the ions). A basic law of electrostatics [7] allows a point charge to be replaced by a spherically symmetrical charge density cloud at its centre for the purposes of representing an attraction from a point outside the charge cloud. Various charge density functions (uIr) have been employed in relation to the above equation. All but one can be symbolised
by ~~=li~(R-r)~,
R>r;
(2) =0
>
R
where Xc,is a notrnalisation constant. Jones and Templeton ES] state that functions where x = 3 or greater do not improve convergence over those of lower order, and so the forms considered here are x = 0 (I = U) a zmifum distribution, x = I (I = L) it linc~r distribution, and x = 2 (t = P) .x parabolic distribution- A gatcssiatz distrtbution. corresponding to aG, = k& exp(---k&irr2)
,
all r ,
(3
is also employed_ In (2) R is chosen as half the shortest interatomic distance in the lattice in order to avoid overlap between charge clouds. For the gaussian case overlap cannot be avoided, zmd the evaluation of the energy of this overlap has been discussed by 3ones and Templeton [5] _Values of g1 and kt, and the function elf, have been listed by Jones and Templeton in their paper ISIThe error introduced into the electrostatic energy (AU,,,,) b,v t runcntion of the infinite summation in (1) at a reciprocal lrrttice vector of value nz, is given by
Using the approximations introduced by Templeton [4] (valid where the number of terms in a small iz interval is Iage). which consist of taking the average vAe for the coutombic structure factor and replacing the summation by an integral, we have
Simplifying snd introducing
the Isttice psrameter
(P)
A’
we obtain AU&c
= KpQ&iR
,
(7)
where Q,=Pfs&da.
Cr=2iihR
and
fl= ‘2mz~R_
@I 417
Volume
62. number 3
CHEMICAL
PHYSICS
LElTERS
15 April 1979
The evaluation of the integraIs in the various forms of (8) each require, for the power charge density functions, severa! stages of integration by parts. The specific resuhs are, for a “umfonn" chargedensity function:
in agreement with that quoted by Templeton Si(x) =
f
x F
141, where
dri _
0
The result for the ‘%zed7 Q,=144
I
4(1
function is ! 0 sin f3(i - cos p) + (II+
- cos p)’ 7p7
-
+(cos2~P5cosp) ___63003
cos p)( 1 - cos P) + sin S(5 i- cos P)
2206
(5t4cosp)sinP 1260/i”
21Op4
1osps (4cos3,p+5cosP) 126Op -
(11) +&[&
- Si(S)]+
&[$r
- Si(2/3)]]
,
and for the “‘panzbolic” is sin o(3 f cos B) f (9 + 5 cos P) ----f5 sin P zp8 14p7 84p6
Qp = 3600
f--- cos p I O0803
sin p
cos p
20 L6fl’
--i-&&-Si(P)] 20166
cos P
sin P
84f15
336@
1
,
which is equivaIent to the form of the result of Jenkins [ 1,2]_ Thegaussion charge density function gives an expression Qc = rr3/2 erfc(p/&)
,
erfc(x)
exp(-&)
(13)
where = (~/KJ/~)
J x
For a particular Iattice parameter, convergence we wish to predict a in the summation), i.e. of p_ This actually most convenient). For n ~lJ,Iec < 0.5 X IO+
drr _
(14)
“cut-off” radius of the charge density function and number of decimal places of value of m (the maximum reciprocal !attice vector which it is necessary to include is most easily achieved by reference to graphs of @ against P ($/la against P/R is decimal places of convergence:
kJ mol-l
_
(15)
The cases for n = -2 (i-50 kJ mol-l) to n = 3 (kO.0005kJ mol;‘) are considered here (except for the function where JZ = 3 is not considered due to the very large nz v&es that this value of Jr necessitates)_
“tc~zifom”
3s Results The curves for the “zuzifom”
charge density function are given in fig. 1, those for the “‘linear” case in fig. 2,
for the “parabolic" in fig. 3,and For the “gmsskm” in fig_4. For a particular ionic lattice a value of PJR is readily computed, 418
using this a vahre of /J/a may be interpolated
for the desired degree of convergence_ The maximum
re-
Volume 62, numbct 3
0
20
CHEIMICAL PHYSICS LlXlXLS
+0
60
80 P/R&)
tw
20
--
QO
60
80 P/R (A’, -
40
60
80 P
Fig. f. Corwrgence using “‘ttniform” chrse density function
0
20
-
fern---2co2.
OL-
a
15 April 1979
RCA)
100 -
Fig. 2. Conuersen& using “Zirwar”charge density function for ?Z= -2 to 3.
J 100
Fig. 3. Convergence using “‘prubolic” charge density function for n = -2 to 3.
Fig_ 4_ Ccmer~ence for n = -2 to 3.
I.&I& “Qgmatin” charge density function
419
Volume 62, number 3
15 April 1979
CHEMICAL PHYSICS LEL-TERS
ciprocaI Iattice vector which it is necessary to include in the summation
(m) is then given by:
m =/3/21-R _
(16)
The integrals in (IO) and (14) were evahrated using the NAG [S] routines SISADF and SISAi)F respectively. Figs_ I4 show that, except for very simple lattices, having lattice parameters less than four, or for very low degrees of convergence, the “parabolrir”and “CgatfisiaKcharge density functions are generally most efficient, i.e. imply smaIIer m (@) values. The choice between these two is very dependent on lattice parameter and degree of convergence, however generally for compiex Iattices and “accurate” results the gaussian function is the most effi-
cient_
4. Discussion
In the computer program LATEN [6], m _ order to facilitate its use by chemists unfamiliar with the precise methods of lattice energy cakulation, a routine has been built into the program to enable the truncation point to be estimated during the in, rather than input as data. For the most commonly used, “‘parabolic” charge density function. obtained
such a prediction of the truncation point involves soiution of eq. (12) to obtain a value offi (Qp being from (7))_ The soIution of such an equation is obviousiy best avoided, where possible, for the production
of an efficient computer program. The prediction in LATEN is made therefore by reference to the curves in Iigs. I-4, each of which have been approximated by series of straight Iines- These sets of straight lines are chosen so as to most cIoseIy appro?cimate the curves, without passing beneath them at any point_ Using this approach the program estimzttes the truncation point in a mere fraction of the time that would be required to solve equation of the form of (12), with only a small loss in accuracy for the establishment of the truncation point (negligible in most cases)_ In the calculation of the electrostatic energy of a lattice it is often desirable to produce a resuIt which is a function of a charge distribution in a cornpIes ion (see for example ref. [9])_ The lattice parameter(P) should then be computed using the mztximum wIue of Sions ~7; for which the resuhs are to be relied on (usually the classica oxidation numbers of the component atoms) (i-e_ the maximum P)_ The ability to be able to predict accurately the possible truncation point of the Bertaut summation leads to quite appreciable savings in computer time, as well as greater reliability of results_ The computer time required for the electrostatic energy calculation for comp1e.u lsttices and minerals. is often large and so such optimisation is highly desirable_
AcknovvIedgement The computations were performed on the Burroughs 6700 computer at the University of Warwick Computer Centre_ KFP wishes to ncknowI+dge the provision of a Science Research Council Studentship_ References [I I H-D-B. Jenkins. Chem. Phys. Letters 9 (1971) 473. ItI H-D-B. Jenhins.Chem. Phys- Letters 20 (1973) 155. I31 E-F- Bertzrut,J- Phys. Radium 13 (1953) 499_ [jl [Sl 161 [7j [al 191
420
D-if-Templeton, J- Chem- Phls_ 23 (19.55) 167-P. R.E_ Jones and D.H. Templeton, J. Chem_ Phks. 25 (1956) 106X H.D.B. Jenkins and K-F. Pratt. Computer PII~S.Commun.. in preptimrion. B-J- Blay and B- Bleany. Electricity and mqnetism. 2nd Ed- (Ckuendon Press, Oxford, 1965). NAG FORTRAN Routines, .hlk_6 (FLM 6). Numeric=ilAlgorithms Croup. Oxford. UK (1977). H.D.B. Jenkins and K-F_ Pratt. J. Chem_ Sot. Faraday II 74 (1978) 968.
CHEhlICAL PHYSICS LETTERS
CONVERGENCE
OF THE RERTAUT
15 April 1979
SERIES AND THE CALCULATION
OF THE ELECTROSTATIC ENERGY OF AN EXTENDED CRYSTALLINE
LATTICE
H.D.B. JENKINS and K.F. PRATT ScItooi
Of Cknistiy and Moleculnr Sciences. UIziiersiry of Warwick. Co&enrr_r.WestMidiands, UAF
Received 10 Jnmury 1979; in tin21 form 16 hiarch1979
Computations are reported wbicb enable the prediction of the necessary truncation of the Bertaut series for urious cbrge density functions consistent with various degrees of conwrgence. The results are sufficiently e.rtensive to be appliable to nlmost all httices-
I_ Introduction Studies designed to obtain accurate lattice energies and related thermochemical data for various ionic salts having complex ions and subsequent extension of this pork into the mineraIogica1 field here led to a situation where the studies of Jenkins [ 1,2] on the truncation of the Bertaut Madelung series have become insufficiently extensive to predict truncation limits. The present work seeks to update and extend these studies so as to apply to almost alI crystai lattices_ The series developed by Bertaut i3] adopts a summation technique which operates entirely over the reciprocal, rather than the real, lattice. The computations necessary to evaluate this summation are extensive for all but the simplest lattice- It is therefore desirable to be able to predict, with reasonable accuracy, the truncation point which can be used to terminate this series, for particular degrees of convergence. Templeton [4] began investigations by consideration of the truncation of the series when using 3 uniform charge density function to represent the point charges (see section 2). His work however covered a range of lattices which included only simple ionic salts- Further work by Jones and Templeton [S] included consideration of other charge density functions, but stili on!y for this limited range of lattices- The later work by Jenkins [I .2 J covered a range of lattices greater than that of Templeton, but still insufficient to include many lattices which are commonly encountered- This work examined truncaticn properties specifically for a “parzlbolic” charge density function_ The development of a completely general computer program for Iattice energy calculation (LATEN [6]) has necessitated 3 general approach to this subject_ Different charge density functions are considered for a range of possible series convergence and results cttlculated to cover virtualiy all possible ionic lattices- Convergence is investigated in terms of the etectrostatic energy rather than the Rladelung constant, as has been the case previously IlA51.
The electrostatic
energy of an ionic lattice (Uelec) is given by the Bertaut approach
[3] as
N u e,ec = W+2Z) 416
,g
nic$
- UWnVZ)
F
I FI~129;J~z2 - U,, •i- i.J,, ,
(1)
Volume 62, number 3
CHEMICAL
PHYSICS LlZl’-lTRS
15 April 1979
where K is the conversion Factor from e2A-1
to lcJ mol-’ (taken as 1389.30 in this work), g1 is a constant dependent on the charge density function chosen to represent the spherically symmetrical chrrrge distribution (see below), R is the cut-off radius of the charge density function, 2 is the number of molecules per unit cell. N is the number of types of atom in the lattice, ni is the number of the jth atom type in the unit cell, of atomic charge qi. r/is the unit cell volume, h is a reciprocal lattice vector, Fjr is the coulombic structure factor, orb the Fourier transform of the charge density function, U,, the energy of overlap of the atomic charge distributions (chosen to be zero for all but the gaussian function) and Use, the self-energy of the compIax ions in the iattice (simpIy the electrostatic energy terms internal to the ions). A basic law of electrostatics [7] allows a point charge to be replaced by a spherically symmetrical charge density cloud at its centre for the purposes of representing an attraction from a point outside the charge cloud. Various charge density functions (uIr) have been employed in relation to the above equation. All but one can be symbolised
by ~~=li~(R-r)~,
R>r;
(2) =0
>
R
where Xc,is a notrnalisation constant. Jones and Templeton ES] state that functions where x = 3 or greater do not improve convergence over those of lower order, and so the forms considered here are x = 0 (I = U) a zmifum distribution, x = I (I = L) it linc~r distribution, and x = 2 (t = P) .x parabolic distribution- A gatcssiatz distrtbution. corresponding to aG, = k& exp(---k&irr2)
,
all r ,
(3
is also employed_ In (2) R is chosen as half the shortest interatomic distance in the lattice in order to avoid overlap between charge clouds. For the gaussian case overlap cannot be avoided, zmd the evaluation of the energy of this overlap has been discussed by 3ones and Templeton [5] _Values of g1 and kt, and the function elf, have been listed by Jones and Templeton in their paper ISIThe error introduced into the electrostatic energy (AU,,,,) b,v t runcntion of the infinite summation in (1) at a reciprocal lrrttice vector of value nz, is given by
Using the approximations introduced by Templeton [4] (valid where the number of terms in a small iz interval is Iage). which consist of taking the average vAe for the coutombic structure factor and replacing the summation by an integral, we have
Simplifying snd introducing
the Isttice psrameter
(P)
A’
we obtain AU&c
= KpQ&iR
,
(7)
where Q,=Pfs&da.
Cr=2iihR
and
fl= ‘2mz~R_
@I 417
Volume
62. number 3
CHEMICAL
PHYSICS
LElTERS
15 April 1979
The evaluation of the integraIs in the various forms of (8) each require, for the power charge density functions, severa! stages of integration by parts. The specific resuhs are, for a “umfonn" chargedensity function:
in agreement with that quoted by Templeton Si(x) =
f
x F
141, where
dri _
0
The result for the ‘%zed7 Q,=144
I
4(1
function is ! 0 sin f3(i - cos p) + (II+
- cos p)’ 7p7
-
+(cos2~P5cosp) ___63003
cos p)( 1 - cos P) + sin S(5 i- cos P)
2206
(5t4cosp)sinP 1260/i”
21Op4
1osps (4cos3,p+5cosP) 126Op -
(11) +&[&
- Si(S)]+
&[$r
- Si(2/3)]]
,
and for the “‘panzbolic” is sin o(3 f cos B) f (9 + 5 cos P) ----f5 sin P zp8 14p7 84p6
Qp = 3600
f--- cos p I O0803
sin p
cos p
20 L6fl’
--i-&&-Si(P)] 20166
cos P
sin P
84f15
336@
1
,
which is equivaIent to the form of the result of Jenkins [ 1,2]_ Thegaussion charge density function gives an expression Qc = rr3/2 erfc(p/&)
,
erfc(x)
exp(-&)
(13)
where = (~/KJ/~)
J x
For a particular Iattice parameter, convergence we wish to predict a in the summation), i.e. of p_ This actually most convenient). For n ~lJ,Iec < 0.5 X IO+
drr _
(14)
“cut-off” radius of the charge density function and number of decimal places of value of m (the maximum reciprocal !attice vector which it is necessary to include is most easily achieved by reference to graphs of @ against P ($/la against P/R is decimal places of convergence:
kJ mol-l
_
(15)
The cases for n = -2 (i-50 kJ mol-l) to n = 3 (kO.0005kJ mol;‘) are considered here (except for the function where JZ = 3 is not considered due to the very large nz v&es that this value of Jr necessitates)_
“tc~zifom”
3s Results The curves for the “zuzifom”
charge density function are given in fig. 1, those for the “‘linear” case in fig. 2,
for the “parabolic" in fig. 3,and For the “gmsskm” in fig_4. For a particular ionic lattice a value of PJR is readily computed, 418
using this a vahre of /J/a may be interpolated
for the desired degree of convergence_ The maximum
re-
Volume 62, numbct 3
0
20
CHEIMICAL PHYSICS LlXlXLS
+0
60
80 P/R&)
tw
20
--
QO
60
80 P/R (A’, -
40
60
80 P
Fig. f. Corwrgence using “‘ttniform” chrse density function
0
20
-
fern---2co2.
OL-
a
15 April 1979
RCA)
100 -
Fig. 2. Conuersen& using “Zirwar”charge density function for ?Z= -2 to 3.
J 100
Fig. 3. Convergence using “‘prubolic” charge density function for n = -2 to 3.
Fig_ 4_ Ccmer~ence for n = -2 to 3.
I.&I& “Qgmatin” charge density function
419
Volume 62, number 3
15 April 1979
CHEMICAL PHYSICS LEL-TERS
ciprocaI Iattice vector which it is necessary to include in the summation
(m) is then given by:
m =/3/21-R _
(16)
The integrals in (IO) and (14) were evahrated using the NAG [S] routines SISADF and SISAi)F respectively. Figs_ I4 show that, except for very simple lattices, having lattice parameters less than four, or for very low degrees of convergence, the “parabolrir”and “CgatfisiaKcharge density functions are generally most efficient, i.e. imply smaIIer m (@) values. The choice between these two is very dependent on lattice parameter and degree of convergence, however generally for compiex Iattices and “accurate” results the gaussian function is the most effi-
cient_
4. Discussion
In the computer program LATEN [6], m _ order to facilitate its use by chemists unfamiliar with the precise methods of lattice energy cakulation, a routine has been built into the program to enable the truncation point to be estimated during the in, rather than input as data. For the most commonly used, “‘parabolic” charge density function. obtained
such a prediction of the truncation point involves soiution of eq. (12) to obtain a value offi (Qp being from (7))_ The soIution of such an equation is obviousiy best avoided, where possible, for the production
of an efficient computer program. The prediction in LATEN is made therefore by reference to the curves in Iigs. I-4, each of which have been approximated by series of straight Iines- These sets of straight lines are chosen so as to most cIoseIy appro?cimate the curves, without passing beneath them at any point_ Using this approach the program estimzttes the truncation point in a mere fraction of the time that would be required to solve equation of the form of (12), with only a small loss in accuracy for the establishment of the truncation point (negligible in most cases)_ In the calculation of the electrostatic energy of a lattice it is often desirable to produce a resuIt which is a function of a charge distribution in a cornpIes ion (see for example ref. [9])_ The lattice parameter(P) should then be computed using the mztximum wIue of Sions ~7; for which the resuhs are to be relied on (usually the classica oxidation numbers of the component atoms) (i-e_ the maximum P)_ The ability to be able to predict accurately the possible truncation point of the Bertaut summation leads to quite appreciable savings in computer time, as well as greater reliability of results_ The computer time required for the electrostatic energy calculation for comp1e.u lsttices and minerals. is often large and so such optimisation is highly desirable_
AcknovvIedgement The computations were performed on the Burroughs 6700 computer at the University of Warwick Computer Centre_ KFP wishes to ncknowI+dge the provision of a Science Research Council Studentship_ References [I I H-D-B. Jenkins. Chem. Phys. Letters 9 (1971) 473. ItI H-D-B. Jenhins.Chem. Phys- Letters 20 (1973) 155. I31 E-F- Bertzrut,J- Phys. Radium 13 (1953) 499_ [jl [Sl 161 [7j [al 191
420
D-if-Templeton, J- Chem- Phls_ 23 (19.55) 167-P. R.E_ Jones and D.H. Templeton, J. Chem_ Phks. 25 (1956) 106X H.D.B. Jenkins and K-F. Pratt. Computer PII~S.Commun.. in preptimrion. B-J- Blay and B- Bleany. Electricity and mqnetism. 2nd Ed- (Ckuendon Press, Oxford, 1965). NAG FORTRAN Routines, .hlk_6 (FLM 6). Numeric=ilAlgorithms Croup. Oxford. UK (1977). H.D.B. Jenkins and K-F_ Pratt. J. Chem_ Sot. Faraday II 74 (1978) 968.