- Email: [email protected]
The convergence of the direct lattice sums in the random phase approximation method applied to periodic infinite systems
Volume 2 IO, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
The convergence of the direct lattice sums in the random phase approximation method applied to periodic infinite systems Benoit Champagne a,b*’ , David H. Mosley a and Jean-Marie Andre a a Laboratoire de Chrmie ThPorique Apphquie, Facultk Universitaires Notre-Dame de la Paix, 6 I, rue de Eruxelles~ 5000 Namur, Belgium b Quantum Theory Projects Williamson Hall, 362, University of Florida, Gainesville, FL 32611.2085,
USA
Received 8 February 1993; in final form 14 May 1993
The random phase approximation procedure applied to stereoregular polymers requires the evaluation of two types of twoelectron integrals between crystal orbitals. Proofs are given in this work that for both types, the direct lattice sums converge. Infinite hydrogen chains are used as model examples.
1. Introduction Polymer quantum chemistry is rapidly expanding and is proving an efficient tool to describe, understand, and characterize the electronic structure of polymers that are, for instance, candidates in nonlinear optics [ I]. In practice, in order to reduce the computational effort, the one-dimensional periodicity of the direct lattice is taken into account and leads to an electronic characterization of these systems in the form of a band structure and its corresponding crystalline orbitals. Unlike finite molecular systems, the evaluation of the electrostatic interactions in infinite systems is a difficult task because of their long-range behaviour. At the restricted Hartree-Fock level, numerous studies have dealt with the evaluation of the Coulomb and exchange terms and with the convergence of their direct lattice summations [ 2- 111. Recently we have proposed methods [ 12-141 to compute the longitudinal polarizabilities per unit cell of infinite periodic systems. Coupled Hartree-Fock values are obtained via the random phase approximation (RPA) method adapted to these polymeric systems. The method requires the computation of numerous twoelectron integrals between crystalline orbitals. As a consequence of this, an obvious question arises: what behaviour do the direct lattice sums involved in these computations exhibit? In addressing this question, this work is organized as follows: - As the polarizability results are influenced by the computed band structure, the first part summarizes the key points concerning the lattice sums that occur in the calculation of the band structure of stereoregular polymers. This summary focuses attention onto the truncation of the lattice sum in the evaluation of the exchange terms. Numerical results on model molecular hydrogen chain, demonstrating the effect of improper truncation on the band structure, are used to define suitable extents of the lattice sums in order to obtain converged crystal orbital energies. - The second part consists of a brief description of the polymeric RPA method and highlights the direct lattice summations. -The last part demonstrates the converging behaviour of the different lattice sums involved in the evaluation ’ Research Assistant of the National Fund for Scientific Research (Belgium).
232
0009-2614/93/$
06.00 0 1993 Elsevier Science Publishers B.V. All rights reserved,
Volume 2 IO, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July I993
of the two-electron integrals between crystalline orbitals. Numerical results on the molecular hydrogen chain support these theoretical establishments,
2. The band structure calculation In restricted Hartree-Fock theory, the many-electron wavefunctions of the closed-shell one-dimensional periodic systems are approximated by Slater determinants which are constructed from doubly occupied crystalline orbitals. These crystalline orbitals are single-particle states defined in the linear combination of atomic orbitals (LCAO ) formalism by
where n and k are respectively the band index and the wave vector or quasimomentum associated with the particle. Thus, two indices label a crystalline orbital. 2NS 1 is the (odd) number of unit cells considered (N+co) in the band structure calculation which corresponds also to the number of k states in one band or to the periodicity [ (2N+ 1)a] of the crystalline orbitals imposed by the Born-Karman cyclic boundary conditions. I/ dm is then the normalization factor, a is the unit cell length. The C,,(k) terms are the k dependent LCAO coefficients. &(r-R,-jue,), abbreviated to Xi(r), is the pth atomic orbital centered in thejth unit cell. In practice, these xJp(r) functions are chosen to be contractions of Gaussian functions. e, is the unit vector in the periodicity direction. The polymeric LCAO coefficients C,,(k) and their associated energies t,(k) are, respectively, the eigenvectors and eigenvalues of the matrix equation F(k)C(k)=S(k)C(k)~(k),
(2)
where F(k) and S(k) are the k-dependent Fock and overlap matrices. The standard theory of band structure calculations is described in several papers [ 15-2 I] and is reviewed in two recent monographs [ 1,221. However, due to the effect of the band structure onto the RPA calculation of the polarizability per unit cell and due to the similarity of the lattice sums in both the SCF and RPA procedures, we find it necessary to briefly describe the procedure of evaluation of the k-dependent overlap and Fock matrices here. Their elements are defined by Fourier transforms, e ikjaFOJ P4 ) j=-N
$,(k)=
2 1=--N
e ikjasOj P4 9
where the cell-dependent matrix elements between the pth atomic orbital in the reference cell (j=O) and the qth atomic orbital in the jth cell read as Fj$= (x; ]hSCFIX$ ,
(5)
sg=(x;Ixg
(6)
3
hsc’ is the monoelectronic Fock operator which describes the motion of a single electron in the field of the fixed nuclei and the average Coulomb and Pauli fields due to the presence of the other electrons. The cell-
dependent Fock matrix elements may be partitioned into its different components, namely, the kinetic (T), the electron-nuclei attraction ( I’), the Coulomb (C) and the exchange (X) terms:
233
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
(11) where Q is the number of atoms per subunit of which the nuclei of charge Z, have positions A + hue,. The twoelectron integrals between atomic orbitals are defined by
(12) The density matrix elements 0% are the Fourier transform of the summation over the no doubly occupied bands of the products of the LCAO coefficients, (13) In theory, the lattice sums range from --co to +co, though this is impractical computationally. Two types of lattice sums exist: the Fourier transform to compute S(k) and F(k) and the internal lattice sums entering in the evaluation of FFi. In our study, the calculation of the band structure is performed by the PLH program [ 23,241. An important feature of this program is the division of the cell-index summations occurring in eqs. (3 )-( 1I ) into three regions as shown in fig. 1. - The Sz{ and F$, terms are evaluated in an interval ranging from -N to N. This short-range region (see fig. 1) is defined such that the overlap Sz{ and kinetic T$$ terms are negligible at the extremities. Using Gaussian functions as atomic orbitals ensures an exponential decrease of these integrals as the interatomic distance increases. In consequence, they rapidly tend to zero. Moreover, as every two-electron integral
that corresponds to the electrostatic interaction between the charge distributions $xJ, and #““, evolves proportionally to Sj; and S?:, the cell index summations over n in C$, and over h and n in X2 converge and are limited to the short-range region. - In the computation of F?i, the Coulombic interactions ( VFi+C$) are grouped to be exactly computed because individually, the sums over h are divergent in character and behave like the logarithmically divergent 2N+l mm
-M 4
-WN
+M
2MtI
>
Fig. 1. Divisions of the cell-index summations.
234
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
harmonic series. In the short- and intermediate-range regions [ -M, Mj, all the integrals are computed exactly whereas, for the regions going from -A4 to -co and from +M to co, a truncated multipole expansion technique, including all dipole-dipole and monopole-quadrupole interactions, is used in order to take into account the long-range Coulombic interactions [ 2,3]. A4 is defined such that the two charge distributions (xix’, and &,x,““) do not overlap, which allows the use of the multipole expansion technique. - All the Q may be evaluated rapidly with a high precision due to the convergence of the h and n lattice sums. However, they behave like l/ ]j] as j increases. Thus, the density matrix elements determine the convcrgencc of the j sum. For non-metallic systems, since D,,J k) is continuous and infinitely differentiable evcrywhere in the first Brillouin zone, its Fourier coefficients Drq decrease asymptotically with increasing j faster than any negative power of lj] [ 5,7,25]. (In the metallic case, due to the discontinuity in the D,,(k) function, the 0% Fourier coefficients decay like 1 / ]jl . As a consequence, the j cell index summation is very slow and needs as many terms as millions to converge. However, this situation leads to an infinite polarizability and is thus beyond the scope of this work.) Thus, thej summation which appears in the evaluation ofthe k-dependent Fock matrix elements converges everywhere in the first Brillouin zone. However, very many terms are often needed to stabilize the exchange term. This is particularly the case of systems with a small band gap for which the decrease of Dgn+lmh with increasing j is slow. Practically, this lattice sum is truncated to include just the short-range region and thus creates errors in the evaluation of the Fock matrix elements. In some cases, the errors in F(k) can lead to convergence problems in the self-consistent procedure which is used to solve the Hartree-Fock equation [ 26 1. The problem of evaluating these sums correctly to infinity is considered in many studies [4-l 11. As the polarizability calculations are based on the band structure results, we firstly emphasize the effect of this lattice sum truncation on the computed band structure. For this purpose, we use the molecular hydrogen chain model. The intramolecular distances are fixed at 2.0 au and the intermolecular distance is fixed either at 3.0 au or at 4.0 au in order to visualize the effect of the bond length alternation. For these systems, no convergence problems appear in the SCF procedure. In tables 1 and 2, the energy values of the top and the bottom of the valence and conduction bands are given according to the number of unit cells considered in the short-range region obtained using the PLH program [23,24]. We have used a split valence basis set improved with p orbitals. It is the (3)-21G(*)* atomic basis set [ 27,281. The parentheses indicate the gratuitousness of these symbols since there is no second-row atom here. As would be expected, an increase in the bond length alternation accelerates the convergence of the energies of the top and bottom of the valence and conduction bands with respect to 2N+ 1 [ 4-61. The results on Table 1 Convergence of the electronic band structure of a model molecular hydrogen chain (intermolecular distance 3.0 au) with the number of interacting unit cells (2N-t 1) in the short-range region. 2M+ I is fixed at 45. The energy values in atomic units concern the bottom and top ofthe valence and conduction bnndsobtained with the (3)-21G(*)* basis set (1.0 au ofenergy=4.3598 lo-” J~27.211 eV) 2N+ 1
tvd (bottom)
Gal (top)
cm,,,, (bottom
3 5 7 9 11 13 1.5 17 19 21
-0.60013 -0.6009 I -0.60055 - 0.60062 -0.60060 -0.60060 -0.60060 - 0.60060 -0.60060 - 0.60060
- 0.42372 -0.42619 -0.42582 - 0.42579 -0.42578 -0.42578 - 0.42578 - 0.42578 - 0.42578 - 0.42578
0.01186 0.02958 0.03 128 0.03 I63 0.03 I72 0.03 I74 0.03175 0.03 I75 0.03 I75 0.03 I75
)
tcond ttop
)
0.38022 0.53529 0.52485 0.52716 0.52658 0.52674 0.52670 0.52671 0.5267 1 0.5267 I
235
Volume 2 10, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
Table 2 Convergence of the electronic band structure ofa model molecular hydrogen chain (intermolecular distance 4.0 au) with the number of interacting unit cells (2N+ 1) in the short-range region. 2Mt 1 is fixed at 45. The energy values in atomic units concern the bottom and top of the valence and conduction bands obtained with the (3)-21G( *)* basis set 2Ni- 1
eval(bottom )
&I
3 5 7 9 11 13 15 17 19 21
-0.55029 -0.55039 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038
-0.47820 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826
(top)
&md
(bottom)
0.06584 0.06970 -0.06992 -0.06994 -0.06995 -0.06995 -0.06995 -0.06995 -0.06995 -0.06995
hmd
(top)
0.34659 0.35579 0.35528 0.35533 0.35533 0.35533 0.35533 0.35533 0.35533 0.35533
the hydrogen chain model justify our initial choice of 2N+ 1= 2 1 in the forthcoming calculations of the polarizability per unit cell.
3. The RPA method applied to polymers Recently, we have shown [ 13) that the evaluation of the longitudinal polarizabilities per unit cell of infinite periodic systems at the coupled Hartree-Fock level needs the use of the RPA method due to the unbounded character of the dipole moment operator. This method, based on propagator techniques [ 29-321, is equally eflicient at computing the static polarizability and the frequency-dependent polarizabilities and may be adapted to include electron correlation at the desired level. The frequency-independent longitudinal polarizability of a closed-shell infinite periodic system is given by (14)
where the asterisk (*) indicates the complex conjugated values. The different matrix elements are defined by
(15)
(A(k,k’))ai.~j=[~a(k)-~;(k)I~ab~~~~~~+2(~~.(k)~;(k)I~j(k’)~~(k’))-(#~,(k)4~(kf)I~j(k’)~;(k)) (16) 2
(B(k,k’))ai.6j=(9i(k)9~(k’)I~ji(~)~=,(k))-2(~i(k)~=(k)IQ),(k’)~b(k’)) *
(17)
The subscript labels a, b, ...(i. i, ...) correspond to unoccupied (occupied) bands. Chemist’s or Mulliken notation has been chosen to define the two-electron integrals between crystalline orbitals, (~=(k)4,(k)l~,(k’)46(k’))=
J‘/ dr,drz @X(k,r&W,
rl) &G(Y
r&%(k’, rZ) .
(18)
The first step consists of solving the following linear equations system: Ax+EW=n, B*X+A*Y
236
= Sit.
(19)
CHEMICAL PHYSICS
Volume2IO,number1,2,3
23 July 1993
LETTERS
The matrices A, B, Q are of infinite dimension since there are an infinite number (2N+ 1, N+w) of k values and thus an infinite number of k states in each band. Moreover, there is an equation for each triplet formed by a k value and two band indices. This triplet represents a particle-hole excitation that is vertical in order to preserve the momentum. As is the case in many polymeric techniques, the infinite sum over k is transformed into an integration in the first Brillouin zone, nla
(2N2~1’“$‘~ j (A(k,k’))rr;,bjXbj(k’)+(B(k,k’))ai,bjYbl(k’)dk’=R,i(k) 7
-n/n
x/a (2N2i1)a
TUT
j (B*(k, k’)),,,Xb(k’)+ -n/a
(A*(k, k’))oi,bjYtj(k’)
dk’=!&(k)
*
(20)
At this state, trapezoidal quadrature is used to obtain X(k) and Y (k). The second step consists of the following scalar product that provides the longitudinal polarizability of the polymer
(21) Again, the infinite sum over k is replaced by integration in the first Brillouin zone, This straightforwardly provides the longitudinal polarizability per unit cell, ~“~R~i(k)~~i(k)+~~,,,(k)Y~i(k)
dk >
1 a
(22)
Computing the polarizabilities per unit cell consists in a post-SCF calculation which makes use of the band structure obtained by the PLH program [ 23,241.
4. Convergence behaviour of the direct lattice sums involved in the crystal orbital two-electron integral
evaluation As the two-electron integrals correspond to repulsion between two delocalized one-electron densities represented by products of crystalline orbitals, we do not expect these integrals to be divergent. However, it is important to consider a formal but necessary analysis of the behaviour and nature of the lattice sums occurring in the evaluation of the elements of A and B. By using the LCAO development of the crystalline orbitals ( 1), the four different two-electron integrals needed to evaluate the matrix elements of A (k, k’) and B (k, k’) read
1 (@a,(k)@i(k)I&(k’)@b(k’) )= 2N+
1
pg, E, ,z, Jz, C~(k)C,,(k)CX(k’)C,b(k’) i
(23)
eikha
I h=-N pz, (Qi(k)@a,(k)
I@(k’)h6(k’)
I=
f,
ri,
sg,
C~i(k)C,,(k)C~(k’)C,b(k’)
(24)
& f
eikhn
h=-N
231
Volume 210, number
1,2,3
CHEMICAL
PHYSICS LETTERS
23 July I993
(25)
(26)
The summations over h, n, m are formally infinite, but their extent Iimited in practice to the intervals [ -A’, N] and [ -P, P] will become clear in the discussion which follows. Moreover, the need to distinguish between two cell-index domains will become apparent in the following part. The multiplicative factor 1/ (2NS 1) coming from the normalization factor is exactly compensated by the 2N+ 1 factor which multiplies the integrals in eq. (20) possessing a diverging character when N tends to infinity. Hence, the discussion about the convergence of the direct lattice sums concerns the contents of the large brackets of integrals (23)-(26). We distinguish between two types of integral according to the matching of the quasimomenta k and k’ in the description of the repulsive one-electron distributions. Indeed, by permuting the crystal orbitals & and #,, integral (24) can be obtained from integral (23) and integral (26) from integral (25). It is well-known that the integrals e
:I:
‘7)
and
describe the Coulombic repulsion between the electron clouds defined by the products of two atomic functions (respectively x:x$ with x?x~‘” and ,$x?+~ with x:x$). It comes straightforwardly that their values depend upon the overlap between these functions (S$Spm”‘+n=SighS:: and Sz:+mS$). As the overlap integrals decrease as exp( -v2), where r is the distance between the atomic centers, the sum over m converges rapidly in the integrals (25 ) and (26 ). The remaining part of this two-electron integral to evaluate consists of summations over the set of the atomic orbitals and of a 2D Fourier transform which also converge rapidly due to the behaviour of the overlap terms as the distance increases. Thus the summations in eqs. (25) and (26) can be restricted to a short-range region, going from -N to N, in the same way as for the evaluation of the k-dependent overlap matrix in section 2. For the second type of integral ( (23) and (24) ), it seems more difficult because its evaluation involves a sum over the cell index m which has a diverging character. Indeed, its elements evolve like l/ 1m I. In order to demonstrate the converging behaviour of these sums, the integral (23) is written in a new form, Pz, $, $, $i, C~i(k)C,,(k)C~(~)c,b(k’)
(2N+ 1)(Ci(k)@a,(k) I @;,,_(k’)@bCk’) )= i m=-P
e ikba
I
5
f h=
-N
n---N
eik%a(~
“,I:
mbn)
I .
(27)
The sum over m can be decomposed into two regions: the first going from -P to P and the second from P to cc and from -P to - 00. Then, we evaluate the repulsion integral between atomic orbitals for ( m I 1 P. Whatever the angular momentum of Gaussian-type orbitals (s, p, d, ...). the integral may be written in the following general form [ 33 ] : 238
SOhsOn m ‘s
Ima+A(p,q,r,.5h,n)l
+
C(P,4, I maSA(p,
23 July 1993
CHEMICAL PHYSICS LETTERS
Volume 2 10, number1,2,3
m, 9, r, $7h, n1 + Ima+A(p,q,r,.hkn)IZ
f-9s, A,n)
(28)
q, r, s, h, n) 13 + ... .
A(p, q, r, s, h, n), B(p, q, T, s, h, n) and C(p, q, r, s, h, n) are constants and depend upon the centers and the forms of the contracted Gaussian functions. For example, B(p, q, r, s, h, n), C(p, q, r, s, h, n), ... are zero if all the functions are of s-type [ 331. The second, third, ... terms present a converging character as m increases. Indeed, the summations of the type I%:= _-ml/ I m Ix converge when X> 2. Our discussion is then restricted to the first term which may create a problem due to the m-cell summation. As a consequence of the small value of A(p, q, r, s, h, n) with respect to mu, the denominator of the first term of eq. (28) may be expanded in a MacLaurin series, 1
1
+ A (P, 4,5 3, k fl)
Ima+A(a,q,r,s,h,n)l =-----lmal -
lma12
T
A@, 4, r,s, h, n)* +A(~,q,r,f,kfl)~_
Ima13
-
lH4
+ .,. .
(29)
Again, the summations over m for the second, third, ... terms present a converging behaviour. Hence, our discussion is finally focused on the first term of eq. (29). As it only depends upon m, for large m, the bracket of eq. (27) reads
hpz, f, f, f,
c,‘,(k)C,,(k)C~(k’)C,,(k’)lf eikhu5 eik’naS%?f . h=-N
tl-N
(30)
After reordering the different sums and terms as (31) there appears in the equation the definition of the k-dependent overlap matrix elements,
(32) and finally, by using the orthonormalization condition between the crystalline orbitals we find that the term is equal to zero. As for the [ - M, M] domain in the Coulombic terms, [ -P, P] domain is expanded such that the one-electron distributions do not overlap, and the MacLaurin series of formula (29) may be applied. Thus, the infinite sum may be restricted to the intermediate-range region [ -P, P] and these two-electron integrals over crystalline orbitals will converge. However, it is necessary to make a distinction between the sum over m and the sums over h and n entering in the evaluation of the 2D Fourier transforms. It is necessary to distinguish between the effects of truncation of the cell-index summations in the evaluation of the two-electron integrals between crystalline orbitals from the truncation of these sums in the SCF procedure which influence the crystalline orbitals and their energies and thus the elements of matrices A and B. Firstly, we compute the band structure with 2N+ 1 and 2M+ 1 respectively fixed at 2 1 and 45 using the (3)2 1G(*)* basis set [ 27,281, as it was concluded at the end of the last section that such cell-index domains ensure the band structure to have converged. Then, we evaluate the matrix elements of A and B by considering different cell-index summation domains [ -IV, N] and [ -P, P]. The RPA longitudinal polarizabilities per unit cell obtained by solving eqs. (20) and (22) are listed in tables 3 and 4 for the molecular hydrogen chains of two different intermolecular distances defined previously according to the number of interacting unit cells in the [ -N, N] and [ -P, P] domains. For a fixed value of 2N+ 1, the polarizability per unit cell converges towards the asymptotic value as P increases showing the good behaviour of the m summation in eqs. (23)-(26) that we have demonstrated. It is, 239
Volume 2 10, number l,2,3
CHEMICAL PHYSICS LETTERS
Table 3 Effect of the size of the intermediate-range region (2Pt 1) on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au (1.Oau of polarizability=1,6488xlO-*‘C*m*J-‘=0.14818A3)
Table 4 Effect of the size of the short-range region (2N+ I ) on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au. P is fixed at 2N+ 1 2N-k 1
2N+ I
2PfI
23 July 1993
%(RPA)I(2N+
1)
dRPAlI(2NtI) &,,,=3.0au
di,,,=4.0au
3 3 3 3
5 1 11 25
27.098 27.064 27.064 27.064
17.220 17.209 17.209 17.209
7 7 1 7 7 7
9 II 13 15 17 25
28.26 1 28.361 28.420 28.411 28.411 28.411
17.443 17.415 17.494 17.492 17.492 17.492
13 13 13 13 13 13 13
15 17 19 21 23 25 27
28.464 28.489 28.506 28.5 19 28.528 28.536 28.532
17.505 17.513 17.519 17.523 17.526 17.528 17.521
13
29
28.532
17.521
3 5 7 9 II 13 15 17 19 21 23 25 27
&,,,= 3.0 au
&,..= 4.0 au
27.064 28.182 28.411 28.482 28.514 28.532 28.543 28.550 28.555 28.559 28.562 28.563 28.563
17.209 17.439 17.492 17.512 17.522 17.527 17.531 17.533 17.534 17.535 17.536 17.536 11.536
however, important to mention that proper convergence needs correctly balanced summations, i.e. Pa 2N+ 1. Ensuring the correct balance of the summations, table 4 shows the convergence of the polarizability values obtained with increasingly large [ -N, N] domains. In practice, depending upon the level of accuracy desired, the extent of the summation domains [ -N, N] and [ -P, P] could be taken to be [ - 6, 61 and [ - 13, 131 or [ - 3, 31 and [ - 7, 71 in order to obtain an acceptable precision. In the particular case of hydrogen chain model, you obtain respectively 99.9O/band 99.5% of the asymptotic value. As would be expected, the more regular the chain or the smaller the bond length alternation, the slower the convergence of the longitudinal polarizability per unit cell is. Finally, in order to highlight the influence of the band structure calculation on the polarizability, we have computed the polarizability per unit cell by fixing N and P respectively at 6 and 12 in the RPA procedure but we vary the [ -N, N] short-range region in the band structure calculation. These results are presented in table 5. The use of 2N+ 1 = 7 in the SCF procedure already provides polarizability values stabilized to the fifth digit whereas in section 2, we have shown that the saturation of the crystal orbital energies up to the same accuracy needs a larger short-range region. Hence, polarizability calculations do not appear to be as sensitive as the band structure to the extent of the short-range domain used in the SCF procedure. Our results demonstrate the good convergence properties of the two-electron integrals between crystalline orbitals which appear in the RPA method applied to periodic infinite systems. In addition, our results show the need to distinguish between two regions: the cell-index summation domain [ -N, N] concerning the 2D Fourier transform and the [-P, P] domain. In order to get proper convergence, correctly balanced summations have to be performed i.e. Pa 2NS 1, Moreover, the need to find adequate methods to evaluate correctly the
240
CHEMICAL PHYSICS LETTERS
Volume 210, number 1,2,3
23 July 1993
Table 5 Effect of the size of the short-range region (2N+ 1) in the SCF procedure on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au. The RPA calculations are performed by fixing Nand Pat 6 and 12 respectively 2N+ 1
%(RPA)I(2N+
1)
&,,,=3.0 au
di.,.,=4.0 au
30.107 28.551 28.536 28.536 28.536 28.536
17.569 17.528 17.528 17.528 17.528 17.528
exchange terms does not seem to be as important in the calculation of the polarizability as in the calculation of the band structure.
Acknowledgement The authors would like to thank Professor J. Delhalle and Dr. J.G. Fripiat for their interest brought to this work and their comments in the field of polymer quantum chemistry. BC thanks the Belgian National Fund for Scientific Research (FNRS) for his Research Assistant position. DHM thanks the Services de la Programmation de la Politique Scientifique (SPPS) for the grant received in the framework of the ELSAM (Electronic Large Scale Computational System for Advanced Materials) project, a part of the Belgian National Program of Impulsion in Information Technology. All calculations reported here have been performed on the NamurScientific Computing Facility (Namur-SCF), a result of a cooperation between the Belgian National Fund for Scientific Research (FNRS), IBM-Belgium, and the Facultes Universitaires Notre-Dame de la Paix (FUNDP).
References [ I] J.M. And@ J. Delhalle and J.L. Br&las, Quantum chemistry aided design of organic polymers for molecular electronics (World Scientific, Singapore, 1991) [2] L. Piela and J. Delhalle, Intern. J. Quantum Chem. 13 (1978) 605. [ 31 J. Delhalle, L. Piela, J.L. Bredas and J.M. Andr6, Phys. Rev. B 22 ( 1980) 6254. [ 41 L. Piela, J.M. And&, J.G. Fripiat and J. Delhalle, Chem. Phys. Letters 77 ( 1981) 143. [ 51I. Delhalle, Intern. J. Quantum Chem. 26 ( 1984) 7 17. [6] J.G. Fripiat, J. Delhalle, J.M. Andre and J.L. Calais, Intern. J. Quantum Chem. Quantum Chem. Symp. S24 ( 1990) 593. [ 71 J. Delhalle and J.L. Calais, J. Chem. Phys. 85 ( 1986) 5286. [ 81 J. Delhalle and J.L. Calais, Intern. J. Quantum Chem. Quantum Chem. Symp. S21 (1987) 115. 191J. Delhalle, M.H. Delvaux, J.G. Fripiat, J.M. Andre and J.L. Calais, J. Chem. Phys. 88 ( 1988) 3 141. [IO] L.Z. Stolarczyk, M. Jeziorskaand H.J. Monkhont, Phys. Rev. B 37 (1988) 10646. [ II] M. Jeziorska, L.Z. Stolarczyk, J. Paldus and H.J. Monkhorst, Phys. Rev. B 41 ( 1990) 12473. [ 121 B. Champagne and J.M. Andre, Intern. J. Quantum Chem. 42 (1992) 1009. [ 131 B. Champagne, J.G. Fripiat and J.M. Andrk, J. Chem. Phys. 96 (1992) 8330. [ 141 B. Champagne, D.H. Mosley, J.G. Fripiat and J.M. And&, Intern. J. Quantum Chem. 46 ( 1993) I. [ 151J.M. AndrC, L. Gouverneur and G. Leroy, Intern. J. Quantum Chem. I (1967) 427. [ 161J.M. And&, L. Gouvemeur and G. Leroy, Intern. J. Quantum Chem. 1 ( 1967) 45 I. [ 171J.M. Andre and G. Leroy, Theoret. Chim. Acta 9 ( 1967) 123. [ 181 G. Del Re, J. Ladik and G. Biczo, Phys. Rev. 155 (1967) 997. [ 191J. Ladik, Advan. Quantum Chem. 7 ( 1973) 397.
241
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
[20 ] J.M. Andre, Advan. Quantum Chem. I2 ( 1980) 65. [21]M.Kertesz,Advan.QuantumChem. 15 (1982) 161. [22] J. Ladik, Quantum theory of polymers as solids (Plenum Press, New York, 1988). [23] J.M. Andre, J.L. Brtdas, J. Delhalle, D.J. Vanderveken, D.P. Vercauteren and J.G. Fripiat, in: MOTECC-9 I, modern techniques in computational chemistry, ed. E. Clementi (ESCOM, Leiden, 1991) p. 793. [24] J.G. Fripiat, D.J. Vanderveken, D.P. Vercauteren and J.M. Andre, in: MOTECC-91, modem techniques in computational chemistry, PLH-91 from MOTECC-91, Input/Output Documentation, ed. E. Clementi, part 3 (ESCOM, Leiden, 1991 ) p. 195. [25] N.K. Bary, Treatise on trigonometric series (Pergamon Press, Oxford, 1964). [26] S. Suhai, P.S. Bagus and J. Ladik, Chem. Phys. 68 (1982) 467. (271 PC. Hariharan and J.A. Pople, Theoret. Chim. Acta 28 (1973) 213. [28] W.J. Pietro, M.M. Francl, W.J. Hehre, D.J. DeFrees, J.A. Pople and J.S. Binkley, J. Am. Chem. Sot. 104 (1982) 5039. [29] J. Linderbergand Y. ohm, Propagators in quantum chemistry (Academic Press, New York, 1973). [30] J. Oddershede, Advan. Quantum Chem. 11 (1978) 257. [31] P. Jorgensen and J. Simons, Second quantization-based methods in quantum chemistry (Academic Press, New York, 1981). [32] J. Oddershede, P. Jorgensen and D.L. Yeager, Computer Phys. Rept. 2 ( 1984) 33. [33] E. Clementi and D.R. Davies, J. Comput. Phys. 1 ( 1966) 86.
242
CHEMICAL PHYSICS LETTERS
23 July 1993
The convergence of the direct lattice sums in the random phase approximation method applied to periodic infinite systems Benoit Champagne a,b*’ , David H. Mosley a and Jean-Marie Andre a a Laboratoire de Chrmie ThPorique Apphquie, Facultk Universitaires Notre-Dame de la Paix, 6 I, rue de Eruxelles~ 5000 Namur, Belgium b Quantum Theory Projects Williamson Hall, 362, University of Florida, Gainesville, FL 32611.2085,
USA
Received 8 February 1993; in final form 14 May 1993
The random phase approximation procedure applied to stereoregular polymers requires the evaluation of two types of twoelectron integrals between crystal orbitals. Proofs are given in this work that for both types, the direct lattice sums converge. Infinite hydrogen chains are used as model examples.
1. Introduction Polymer quantum chemistry is rapidly expanding and is proving an efficient tool to describe, understand, and characterize the electronic structure of polymers that are, for instance, candidates in nonlinear optics [ I]. In practice, in order to reduce the computational effort, the one-dimensional periodicity of the direct lattice is taken into account and leads to an electronic characterization of these systems in the form of a band structure and its corresponding crystalline orbitals. Unlike finite molecular systems, the evaluation of the electrostatic interactions in infinite systems is a difficult task because of their long-range behaviour. At the restricted Hartree-Fock level, numerous studies have dealt with the evaluation of the Coulomb and exchange terms and with the convergence of their direct lattice summations [ 2- 111. Recently we have proposed methods [ 12-141 to compute the longitudinal polarizabilities per unit cell of infinite periodic systems. Coupled Hartree-Fock values are obtained via the random phase approximation (RPA) method adapted to these polymeric systems. The method requires the computation of numerous twoelectron integrals between crystalline orbitals. As a consequence of this, an obvious question arises: what behaviour do the direct lattice sums involved in these computations exhibit? In addressing this question, this work is organized as follows: - As the polarizability results are influenced by the computed band structure, the first part summarizes the key points concerning the lattice sums that occur in the calculation of the band structure of stereoregular polymers. This summary focuses attention onto the truncation of the lattice sum in the evaluation of the exchange terms. Numerical results on model molecular hydrogen chain, demonstrating the effect of improper truncation on the band structure, are used to define suitable extents of the lattice sums in order to obtain converged crystal orbital energies. - The second part consists of a brief description of the polymeric RPA method and highlights the direct lattice summations. -The last part demonstrates the converging behaviour of the different lattice sums involved in the evaluation ’ Research Assistant of the National Fund for Scientific Research (Belgium).
232
0009-2614/93/$
06.00 0 1993 Elsevier Science Publishers B.V. All rights reserved,
Volume 2 IO, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July I993
of the two-electron integrals between crystalline orbitals. Numerical results on the molecular hydrogen chain support these theoretical establishments,
2. The band structure calculation In restricted Hartree-Fock theory, the many-electron wavefunctions of the closed-shell one-dimensional periodic systems are approximated by Slater determinants which are constructed from doubly occupied crystalline orbitals. These crystalline orbitals are single-particle states defined in the linear combination of atomic orbitals (LCAO ) formalism by
where n and k are respectively the band index and the wave vector or quasimomentum associated with the particle. Thus, two indices label a crystalline orbital. 2NS 1 is the (odd) number of unit cells considered (N+co) in the band structure calculation which corresponds also to the number of k states in one band or to the periodicity [ (2N+ 1)a] of the crystalline orbitals imposed by the Born-Karman cyclic boundary conditions. I/ dm is then the normalization factor, a is the unit cell length. The C,,(k) terms are the k dependent LCAO coefficients. &(r-R,-jue,), abbreviated to Xi(r), is the pth atomic orbital centered in thejth unit cell. In practice, these xJp(r) functions are chosen to be contractions of Gaussian functions. e, is the unit vector in the periodicity direction. The polymeric LCAO coefficients C,,(k) and their associated energies t,(k) are, respectively, the eigenvectors and eigenvalues of the matrix equation F(k)C(k)=S(k)C(k)~(k),
(2)
where F(k) and S(k) are the k-dependent Fock and overlap matrices. The standard theory of band structure calculations is described in several papers [ 15-2 I] and is reviewed in two recent monographs [ 1,221. However, due to the effect of the band structure onto the RPA calculation of the polarizability per unit cell and due to the similarity of the lattice sums in both the SCF and RPA procedures, we find it necessary to briefly describe the procedure of evaluation of the k-dependent overlap and Fock matrices here. Their elements are defined by Fourier transforms, e ikjaFOJ P4 ) j=-N
$,(k)=
2 1=--N
e ikjasOj P4 9
where the cell-dependent matrix elements between the pth atomic orbital in the reference cell (j=O) and the qth atomic orbital in the jth cell read as Fj$= (x; ]hSCFIX$ ,
(5)
sg=(x;Ixg
(6)
3
hsc’ is the monoelectronic Fock operator which describes the motion of a single electron in the field of the fixed nuclei and the average Coulomb and Pauli fields due to the presence of the other electrons. The cell-
dependent Fock matrix elements may be partitioned into its different components, namely, the kinetic (T), the electron-nuclei attraction ( I’), the Coulomb (C) and the exchange (X) terms:
233
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
(11) where Q is the number of atoms per subunit of which the nuclei of charge Z, have positions A + hue,. The twoelectron integrals between atomic orbitals are defined by
(12) The density matrix elements 0% are the Fourier transform of the summation over the no doubly occupied bands of the products of the LCAO coefficients, (13) In theory, the lattice sums range from --co to +co, though this is impractical computationally. Two types of lattice sums exist: the Fourier transform to compute S(k) and F(k) and the internal lattice sums entering in the evaluation of FFi. In our study, the calculation of the band structure is performed by the PLH program [ 23,241. An important feature of this program is the division of the cell-index summations occurring in eqs. (3 )-( 1I ) into three regions as shown in fig. 1. - The Sz{ and F$, terms are evaluated in an interval ranging from -N to N. This short-range region (see fig. 1) is defined such that the overlap Sz{ and kinetic T$$ terms are negligible at the extremities. Using Gaussian functions as atomic orbitals ensures an exponential decrease of these integrals as the interatomic distance increases. In consequence, they rapidly tend to zero. Moreover, as every two-electron integral
that corresponds to the electrostatic interaction between the charge distributions $xJ, and #““, evolves proportionally to Sj; and S?:, the cell index summations over n in C$, and over h and n in X2 converge and are limited to the short-range region. - In the computation of F?i, the Coulombic interactions ( VFi+C$) are grouped to be exactly computed because individually, the sums over h are divergent in character and behave like the logarithmically divergent 2N+l mm
-M 4
-WN
+M
2MtI
>
Fig. 1. Divisions of the cell-index summations.
234
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
harmonic series. In the short- and intermediate-range regions [ -M, Mj, all the integrals are computed exactly whereas, for the regions going from -A4 to -co and from +M to co, a truncated multipole expansion technique, including all dipole-dipole and monopole-quadrupole interactions, is used in order to take into account the long-range Coulombic interactions [ 2,3]. A4 is defined such that the two charge distributions (xix’, and &,x,““) do not overlap, which allows the use of the multipole expansion technique. - All the Q may be evaluated rapidly with a high precision due to the convergence of the h and n lattice sums. However, they behave like l/ ]j] as j increases. Thus, the density matrix elements determine the convcrgencc of the j sum. For non-metallic systems, since D,,J k) is continuous and infinitely differentiable evcrywhere in the first Brillouin zone, its Fourier coefficients Drq decrease asymptotically with increasing j faster than any negative power of lj] [ 5,7,25]. (In the metallic case, due to the discontinuity in the D,,(k) function, the 0% Fourier coefficients decay like 1 / ]jl . As a consequence, the j cell index summation is very slow and needs as many terms as millions to converge. However, this situation leads to an infinite polarizability and is thus beyond the scope of this work.) Thus, thej summation which appears in the evaluation ofthe k-dependent Fock matrix elements converges everywhere in the first Brillouin zone. However, very many terms are often needed to stabilize the exchange term. This is particularly the case of systems with a small band gap for which the decrease of Dgn+lmh with increasing j is slow. Practically, this lattice sum is truncated to include just the short-range region and thus creates errors in the evaluation of the Fock matrix elements. In some cases, the errors in F(k) can lead to convergence problems in the self-consistent procedure which is used to solve the Hartree-Fock equation [ 26 1. The problem of evaluating these sums correctly to infinity is considered in many studies [4-l 11. As the polarizability calculations are based on the band structure results, we firstly emphasize the effect of this lattice sum truncation on the computed band structure. For this purpose, we use the molecular hydrogen chain model. The intramolecular distances are fixed at 2.0 au and the intermolecular distance is fixed either at 3.0 au or at 4.0 au in order to visualize the effect of the bond length alternation. For these systems, no convergence problems appear in the SCF procedure. In tables 1 and 2, the energy values of the top and the bottom of the valence and conduction bands are given according to the number of unit cells considered in the short-range region obtained using the PLH program [23,24]. We have used a split valence basis set improved with p orbitals. It is the (3)-21G(*)* atomic basis set [ 27,281. The parentheses indicate the gratuitousness of these symbols since there is no second-row atom here. As would be expected, an increase in the bond length alternation accelerates the convergence of the energies of the top and bottom of the valence and conduction bands with respect to 2N+ 1 [ 4-61. The results on Table 1 Convergence of the electronic band structure of a model molecular hydrogen chain (intermolecular distance 3.0 au) with the number of interacting unit cells (2N-t 1) in the short-range region. 2M+ I is fixed at 45. The energy values in atomic units concern the bottom and top ofthe valence and conduction bnndsobtained with the (3)-21G(*)* basis set (1.0 au ofenergy=4.3598 lo-” J~27.211 eV) 2N+ 1
tvd (bottom)
Gal (top)
cm,,,, (bottom
3 5 7 9 11 13 1.5 17 19 21
-0.60013 -0.6009 I -0.60055 - 0.60062 -0.60060 -0.60060 -0.60060 - 0.60060 -0.60060 - 0.60060
- 0.42372 -0.42619 -0.42582 - 0.42579 -0.42578 -0.42578 - 0.42578 - 0.42578 - 0.42578 - 0.42578
0.01186 0.02958 0.03 128 0.03 I63 0.03 I72 0.03 I74 0.03175 0.03 I75 0.03 I75 0.03 I75
)
tcond ttop
)
0.38022 0.53529 0.52485 0.52716 0.52658 0.52674 0.52670 0.52671 0.5267 1 0.5267 I
235
Volume 2 10, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
Table 2 Convergence of the electronic band structure ofa model molecular hydrogen chain (intermolecular distance 4.0 au) with the number of interacting unit cells (2N+ 1) in the short-range region. 2Mt 1 is fixed at 45. The energy values in atomic units concern the bottom and top of the valence and conduction bands obtained with the (3)-21G( *)* basis set 2Ni- 1
eval(bottom )
&I
3 5 7 9 11 13 15 17 19 21
-0.55029 -0.55039 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038 -0.55038
-0.47820 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826 -0.47826
(top)
&md
(bottom)
0.06584 0.06970 -0.06992 -0.06994 -0.06995 -0.06995 -0.06995 -0.06995 -0.06995 -0.06995
hmd
(top)
0.34659 0.35579 0.35528 0.35533 0.35533 0.35533 0.35533 0.35533 0.35533 0.35533
the hydrogen chain model justify our initial choice of 2N+ 1= 2 1 in the forthcoming calculations of the polarizability per unit cell.
3. The RPA method applied to polymers Recently, we have shown [ 13) that the evaluation of the longitudinal polarizabilities per unit cell of infinite periodic systems at the coupled Hartree-Fock level needs the use of the RPA method due to the unbounded character of the dipole moment operator. This method, based on propagator techniques [ 29-321, is equally eflicient at computing the static polarizability and the frequency-dependent polarizabilities and may be adapted to include electron correlation at the desired level. The frequency-independent longitudinal polarizability of a closed-shell infinite periodic system is given by (14)
where the asterisk (*) indicates the complex conjugated values. The different matrix elements are defined by
(15)
(A(k,k’))ai.~j=[~a(k)-~;(k)I~ab~~~~~~+2(~~.(k)~;(k)I~j(k’)~~(k’))-(#~,(k)4~(kf)I~j(k’)~;(k)) (16) 2
(B(k,k’))ai.6j=(9i(k)9~(k’)I~ji(~)~=,(k))-2(~i(k)~=(k)IQ),(k’)~b(k’)) *
(17)
The subscript labels a, b, ...(i. i, ...) correspond to unoccupied (occupied) bands. Chemist’s or Mulliken notation has been chosen to define the two-electron integrals between crystalline orbitals, (~=(k)4,(k)l~,(k’)46(k’))=
J‘/ dr,drz @X(k,r&W,
rl) &G(Y
r&%(k’, rZ) .
(18)
The first step consists of solving the following linear equations system: Ax+EW=n, B*X+A*Y
236
= Sit.
(19)
CHEMICAL PHYSICS
Volume2IO,number1,2,3
23 July 1993
LETTERS
The matrices A, B, Q are of infinite dimension since there are an infinite number (2N+ 1, N+w) of k values and thus an infinite number of k states in each band. Moreover, there is an equation for each triplet formed by a k value and two band indices. This triplet represents a particle-hole excitation that is vertical in order to preserve the momentum. As is the case in many polymeric techniques, the infinite sum over k is transformed into an integration in the first Brillouin zone, nla
(2N2~1’“$‘~ j (A(k,k’))rr;,bjXbj(k’)+(B(k,k’))ai,bjYbl(k’)dk’=R,i(k) 7
-n/n
x/a (2N2i1)a
TUT
j (B*(k, k’)),,,Xb(k’)+ -n/a
(A*(k, k’))oi,bjYtj(k’)
dk’=!&(k)
*
(20)
At this state, trapezoidal quadrature is used to obtain X(k) and Y (k). The second step consists of the following scalar product that provides the longitudinal polarizability of the polymer
(21) Again, the infinite sum over k is replaced by integration in the first Brillouin zone, This straightforwardly provides the longitudinal polarizability per unit cell, ~“~R~i(k)~~i(k)+~~,,,(k)Y~i(k)
dk >
1 a
(22)
Computing the polarizabilities per unit cell consists in a post-SCF calculation which makes use of the band structure obtained by the PLH program [ 23,241.
4. Convergence behaviour of the direct lattice sums involved in the crystal orbital two-electron integral
evaluation As the two-electron integrals correspond to repulsion between two delocalized one-electron densities represented by products of crystalline orbitals, we do not expect these integrals to be divergent. However, it is important to consider a formal but necessary analysis of the behaviour and nature of the lattice sums occurring in the evaluation of the elements of A and B. By using the LCAO development of the crystalline orbitals ( 1), the four different two-electron integrals needed to evaluate the matrix elements of A (k, k’) and B (k, k’) read
1 (@a,(k)@i(k)I&(k’)@b(k’) )= 2N+
1
pg, E, ,z, Jz, C~(k)C,,(k)CX(k’)C,b(k’) i
(23)
eikha
I h=-N pz, (Qi(k)@a,(k)
I@(k’)h6(k’)
I=
f,
ri,
sg,
C~i(k)C,,(k)C~(k’)C,b(k’)
(24)
& f
eikhn
h=-N
231
Volume 210, number
1,2,3
CHEMICAL
PHYSICS LETTERS
23 July I993
(25)
(26)
The summations over h, n, m are formally infinite, but their extent Iimited in practice to the intervals [ -A’, N] and [ -P, P] will become clear in the discussion which follows. Moreover, the need to distinguish between two cell-index domains will become apparent in the following part. The multiplicative factor 1/ (2NS 1) coming from the normalization factor is exactly compensated by the 2N+ 1 factor which multiplies the integrals in eq. (20) possessing a diverging character when N tends to infinity. Hence, the discussion about the convergence of the direct lattice sums concerns the contents of the large brackets of integrals (23)-(26). We distinguish between two types of integral according to the matching of the quasimomenta k and k’ in the description of the repulsive one-electron distributions. Indeed, by permuting the crystal orbitals & and #,, integral (24) can be obtained from integral (23) and integral (26) from integral (25). It is well-known that the integrals e
:I:
‘7)
and
describe the Coulombic repulsion between the electron clouds defined by the products of two atomic functions (respectively x:x$ with x?x~‘” and ,$x?+~ with x:x$). It comes straightforwardly that their values depend upon the overlap between these functions (S$Spm”‘+n=SighS:: and Sz:+mS$). As the overlap integrals decrease as exp( -v2), where r is the distance between the atomic centers, the sum over m converges rapidly in the integrals (25 ) and (26 ). The remaining part of this two-electron integral to evaluate consists of summations over the set of the atomic orbitals and of a 2D Fourier transform which also converge rapidly due to the behaviour of the overlap terms as the distance increases. Thus the summations in eqs. (25) and (26) can be restricted to a short-range region, going from -N to N, in the same way as for the evaluation of the k-dependent overlap matrix in section 2. For the second type of integral ( (23) and (24) ), it seems more difficult because its evaluation involves a sum over the cell index m which has a diverging character. Indeed, its elements evolve like l/ 1m I. In order to demonstrate the converging behaviour of these sums, the integral (23) is written in a new form, Pz, $, $, $i, C~i(k)C,,(k)C~(~)c,b(k’)
(2N+ 1)(Ci(k)@a,(k) I @;,,_(k’)@bCk’) )= i m=-P
e ikba
I
5
f h=
-N
n---N
eik%a(~
“,I:
mbn)
I .
(27)
The sum over m can be decomposed into two regions: the first going from -P to P and the second from P to cc and from -P to - 00. Then, we evaluate the repulsion integral between atomic orbitals for ( m I 1 P. Whatever the angular momentum of Gaussian-type orbitals (s, p, d, ...). the integral may be written in the following general form [ 33 ] : 238
SOhsOn m ‘s
Ima+A(p,q,r,.5h,n)l
+
C(P,4, I maSA(p,
23 July 1993
CHEMICAL PHYSICS LETTERS
Volume 2 10, number1,2,3
m, 9, r, $7h, n1 + Ima+A(p,q,r,.hkn)IZ
f-9s, A,n)
(28)
q, r, s, h, n) 13 + ... .
A(p, q, r, s, h, n), B(p, q, T, s, h, n) and C(p, q, r, s, h, n) are constants and depend upon the centers and the forms of the contracted Gaussian functions. For example, B(p, q, r, s, h, n), C(p, q, r, s, h, n), ... are zero if all the functions are of s-type [ 331. The second, third, ... terms present a converging character as m increases. Indeed, the summations of the type I%:= _-ml/ I m Ix converge when X> 2. Our discussion is then restricted to the first term which may create a problem due to the m-cell summation. As a consequence of the small value of A(p, q, r, s, h, n) with respect to mu, the denominator of the first term of eq. (28) may be expanded in a MacLaurin series, 1
1
+ A (P, 4,5 3, k fl)
Ima+A(a,q,r,s,h,n)l =-----lmal -
lma12
T
A@, 4, r,s, h, n)* +A(~,q,r,f,kfl)~_
Ima13
-
lH4
+ .,. .
(29)
Again, the summations over m for the second, third, ... terms present a converging behaviour. Hence, our discussion is finally focused on the first term of eq. (29). As it only depends upon m, for large m, the bracket of eq. (27) reads
hpz, f, f, f,
c,‘,(k)C,,(k)C~(k’)C,,(k’)lf eikhu5 eik’naS%?f . h=-N
tl-N
(30)
After reordering the different sums and terms as (31) there appears in the equation the definition of the k-dependent overlap matrix elements,
(32) and finally, by using the orthonormalization condition between the crystalline orbitals we find that the term is equal to zero. As for the [ - M, M] domain in the Coulombic terms, [ -P, P] domain is expanded such that the one-electron distributions do not overlap, and the MacLaurin series of formula (29) may be applied. Thus, the infinite sum may be restricted to the intermediate-range region [ -P, P] and these two-electron integrals over crystalline orbitals will converge. However, it is necessary to make a distinction between the sum over m and the sums over h and n entering in the evaluation of the 2D Fourier transforms. It is necessary to distinguish between the effects of truncation of the cell-index summations in the evaluation of the two-electron integrals between crystalline orbitals from the truncation of these sums in the SCF procedure which influence the crystalline orbitals and their energies and thus the elements of matrices A and B. Firstly, we compute the band structure with 2N+ 1 and 2M+ 1 respectively fixed at 2 1 and 45 using the (3)2 1G(*)* basis set [ 27,281, as it was concluded at the end of the last section that such cell-index domains ensure the band structure to have converged. Then, we evaluate the matrix elements of A and B by considering different cell-index summation domains [ -IV, N] and [ -P, P]. The RPA longitudinal polarizabilities per unit cell obtained by solving eqs. (20) and (22) are listed in tables 3 and 4 for the molecular hydrogen chains of two different intermolecular distances defined previously according to the number of interacting unit cells in the [ -N, N] and [ -P, P] domains. For a fixed value of 2N+ 1, the polarizability per unit cell converges towards the asymptotic value as P increases showing the good behaviour of the m summation in eqs. (23)-(26) that we have demonstrated. It is, 239
Volume 2 10, number l,2,3
CHEMICAL PHYSICS LETTERS
Table 3 Effect of the size of the intermediate-range region (2Pt 1) on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au (1.Oau of polarizability=1,6488xlO-*‘C*m*J-‘=0.14818A3)
Table 4 Effect of the size of the short-range region (2N+ I ) on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au. P is fixed at 2N+ 1 2N-k 1
2N+ I
2PfI
23 July 1993
%(RPA)I(2N+
1)
dRPAlI(2NtI) &,,,=3.0au
di,,,=4.0au
3 3 3 3
5 1 11 25
27.098 27.064 27.064 27.064
17.220 17.209 17.209 17.209
7 7 1 7 7 7
9 II 13 15 17 25
28.26 1 28.361 28.420 28.411 28.411 28.411
17.443 17.415 17.494 17.492 17.492 17.492
13 13 13 13 13 13 13
15 17 19 21 23 25 27
28.464 28.489 28.506 28.5 19 28.528 28.536 28.532
17.505 17.513 17.519 17.523 17.526 17.528 17.521
13
29
28.532
17.521
3 5 7 9 II 13 15 17 19 21 23 25 27
&,,,= 3.0 au
&,..= 4.0 au
27.064 28.182 28.411 28.482 28.514 28.532 28.543 28.550 28.555 28.559 28.562 28.563 28.563
17.209 17.439 17.492 17.512 17.522 17.527 17.531 17.533 17.534 17.535 17.536 17.536 11.536
however, important to mention that proper convergence needs correctly balanced summations, i.e. Pa 2N+ 1. Ensuring the correct balance of the summations, table 4 shows the convergence of the polarizability values obtained with increasingly large [ -N, N] domains. In practice, depending upon the level of accuracy desired, the extent of the summation domains [ -N, N] and [ -P, P] could be taken to be [ - 6, 61 and [ - 13, 131 or [ - 3, 31 and [ - 7, 71 in order to obtain an acceptable precision. In the particular case of hydrogen chain model, you obtain respectively 99.9O/band 99.5% of the asymptotic value. As would be expected, the more regular the chain or the smaller the bond length alternation, the slower the convergence of the longitudinal polarizability per unit cell is. Finally, in order to highlight the influence of the band structure calculation on the polarizability, we have computed the polarizability per unit cell by fixing N and P respectively at 6 and 12 in the RPA procedure but we vary the [ -N, N] short-range region in the band structure calculation. These results are presented in table 5. The use of 2N+ 1 = 7 in the SCF procedure already provides polarizability values stabilized to the fifth digit whereas in section 2, we have shown that the saturation of the crystal orbital energies up to the same accuracy needs a larger short-range region. Hence, polarizability calculations do not appear to be as sensitive as the band structure to the extent of the short-range domain used in the SCF procedure. Our results demonstrate the good convergence properties of the two-electron integrals between crystalline orbitals which appear in the RPA method applied to periodic infinite systems. In addition, our results show the need to distinguish between two regions: the cell-index summation domain [ -N, N] concerning the 2D Fourier transform and the [-P, P] domain. In order to get proper convergence, correctly balanced summations have to be performed i.e. Pa 2NS 1, Moreover, the need to find adequate methods to evaluate correctly the
240
CHEMICAL PHYSICS LETTERS
Volume 210, number 1,2,3
23 July 1993
Table 5 Effect of the size of the short-range region (2N+ 1) in the SCF procedure on the longitudinal polarizability per unit cell of the molecular hydrogen chain models. All values are given in au. The RPA calculations are performed by fixing Nand Pat 6 and 12 respectively 2N+ 1
%(RPA)I(2N+
1)
&,,,=3.0 au
di.,.,=4.0 au
30.107 28.551 28.536 28.536 28.536 28.536
17.569 17.528 17.528 17.528 17.528 17.528
exchange terms does not seem to be as important in the calculation of the polarizability as in the calculation of the band structure.
Acknowledgement The authors would like to thank Professor J. Delhalle and Dr. J.G. Fripiat for their interest brought to this work and their comments in the field of polymer quantum chemistry. BC thanks the Belgian National Fund for Scientific Research (FNRS) for his Research Assistant position. DHM thanks the Services de la Programmation de la Politique Scientifique (SPPS) for the grant received in the framework of the ELSAM (Electronic Large Scale Computational System for Advanced Materials) project, a part of the Belgian National Program of Impulsion in Information Technology. All calculations reported here have been performed on the NamurScientific Computing Facility (Namur-SCF), a result of a cooperation between the Belgian National Fund for Scientific Research (FNRS), IBM-Belgium, and the Facultes Universitaires Notre-Dame de la Paix (FUNDP).
References [ I] J.M. And@ J. Delhalle and J.L. Br&las, Quantum chemistry aided design of organic polymers for molecular electronics (World Scientific, Singapore, 1991) [2] L. Piela and J. Delhalle, Intern. J. Quantum Chem. 13 (1978) 605. [ 31 J. Delhalle, L. Piela, J.L. Bredas and J.M. Andr6, Phys. Rev. B 22 ( 1980) 6254. [ 41 L. Piela, J.M. And&, J.G. Fripiat and J. Delhalle, Chem. Phys. Letters 77 ( 1981) 143. [ 51I. Delhalle, Intern. J. Quantum Chem. 26 ( 1984) 7 17. [6] J.G. Fripiat, J. Delhalle, J.M. Andre and J.L. Calais, Intern. J. Quantum Chem. Quantum Chem. Symp. S24 ( 1990) 593. [ 71 J. Delhalle and J.L. Calais, J. Chem. Phys. 85 ( 1986) 5286. [ 81 J. Delhalle and J.L. Calais, Intern. J. Quantum Chem. Quantum Chem. Symp. S21 (1987) 115. 191J. Delhalle, M.H. Delvaux, J.G. Fripiat, J.M. Andre and J.L. Calais, J. Chem. Phys. 88 ( 1988) 3 141. [IO] L.Z. Stolarczyk, M. Jeziorskaand H.J. Monkhont, Phys. Rev. B 37 (1988) 10646. [ II] M. Jeziorska, L.Z. Stolarczyk, J. Paldus and H.J. Monkhorst, Phys. Rev. B 41 ( 1990) 12473. [ 121 B. Champagne and J.M. Andre, Intern. J. Quantum Chem. 42 (1992) 1009. [ 131 B. Champagne, J.G. Fripiat and J.M. Andrk, J. Chem. Phys. 96 (1992) 8330. [ 141 B. Champagne, D.H. Mosley, J.G. Fripiat and J.M. And&, Intern. J. Quantum Chem. 46 ( 1993) I. [ 151J.M. AndrC, L. Gouverneur and G. Leroy, Intern. J. Quantum Chem. I (1967) 427. [ 161J.M. And&, L. Gouvemeur and G. Leroy, Intern. J. Quantum Chem. 1 ( 1967) 45 I. [ 171J.M. Andre and G. Leroy, Theoret. Chim. Acta 9 ( 1967) 123. [ 181 G. Del Re, J. Ladik and G. Biczo, Phys. Rev. 155 (1967) 997. [ 191J. Ladik, Advan. Quantum Chem. 7 ( 1973) 397.
241
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
[20 ] J.M. Andre, Advan. Quantum Chem. I2 ( 1980) 65. [21]M.Kertesz,Advan.QuantumChem. 15 (1982) 161. [22] J. Ladik, Quantum theory of polymers as solids (Plenum Press, New York, 1988). [23] J.M. Andre, J.L. Brtdas, J. Delhalle, D.J. Vanderveken, D.P. Vercauteren and J.G. Fripiat, in: MOTECC-9 I, modern techniques in computational chemistry, ed. E. Clementi (ESCOM, Leiden, 1991) p. 793. [24] J.G. Fripiat, D.J. Vanderveken, D.P. Vercauteren and J.M. Andre, in: MOTECC-91, modem techniques in computational chemistry, PLH-91 from MOTECC-91, Input/Output Documentation, ed. E. Clementi, part 3 (ESCOM, Leiden, 1991 ) p. 195. [25] N.K. Bary, Treatise on trigonometric series (Pergamon Press, Oxford, 1964). [26] S. Suhai, P.S. Bagus and J. Ladik, Chem. Phys. 68 (1982) 467. (271 PC. Hariharan and J.A. Pople, Theoret. Chim. Acta 28 (1973) 213. [28] W.J. Pietro, M.M. Francl, W.J. Hehre, D.J. DeFrees, J.A. Pople and J.S. Binkley, J. Am. Chem. Sot. 104 (1982) 5039. [29] J. Linderbergand Y. ohm, Propagators in quantum chemistry (Academic Press, New York, 1973). [30] J. Oddershede, Advan. Quantum Chem. 11 (1978) 257. [31] P. Jorgensen and J. Simons, Second quantization-based methods in quantum chemistry (Academic Press, New York, 1981). [32] J. Oddershede, P. Jorgensen and D.L. Yeager, Computer Phys. Rept. 2 ( 1984) 33. [33] E. Clementi and D.R. Davies, J. Comput. Phys. 1 ( 1966) 86.
242