1      subroutine teigd(pcap,qcap,u,t,dt,dos)
2C$Id$
3      implicit none
4      logical klnemn
5      real *8  pcap(*),qcap(*),u(*),t(*),dos(*),dt(*)
6c.......................................................................
7c
8c     two-electron integral routine for s,p,d, and f functions.
9c.......................................................................
10#include "cguess.fh"
11c
12      integer i, k, l, im, m, n
13      integer j, nstep1, nstep2, kin, kmx, kl, mmx, mmxp, mn
14      integer nmx, ntest, min
15      logical open
16c
17      real *8 twopow(14), pifac(14)
18      real *8 pi, on6, on15, on70, on35, pot, potn, factkl
19      real *8 prfac1, prfac2, prfac3, prfac4
20      real *8 zp, zq, zpq, zr, zpr, zqr, zs, zps, zrs, zqs, zpqrs
21      real *8 zpqrs2, zprzqs, zpszqr
22      real *8 xterm, factmn, term
23      real *8 xfac1, xfac2, xfac3, xfac11, xfac21, xfsum
24      real *8 x0, y0, pj, qj, y1, y2, y3, x2, x4, y01, y02
25      real *8 y21, y22, y4, y5, y6
26c
27c.......................................................................
28c     angular factors for exchange integrals. obtained as sum of squares
29c     of slater coefficients c(kappa;l1,m1;l2,m2) divided by
30c     2*(2*l1+1)*(2*l2+1) for any kappa, l1, and l2.
31c.......................................................................
32      real *8 ss0, sp1, pp0, pp2, sd2, pd1, pd3, dd0, dd2
33      real *8 sf3, pf2, pf4, df1, df3, df5, ff0, ff2, ff4, ff6
34      real *8 a, b
35      real *8 f0pol, f2pol, df1pol, df3pol
36      data ss0,sp1,pp0,pp2,sd2,pd1,pd3,dd0,dd2,
37     x     sf3,pf2,pf4,df1,df3,df5,ff0,ff2,ff4,ff6
38     x/ .50000000000d+00, .16666666667d+00, .16666666667d+00,
39     x  .66666666667d-01, .10000000000d+00, .66666666667d-01,
40     x  .42857142857d-01, .10000000000d+00, .28571428571d-01,
41     x  .71428571429d-01, .42857142857d-01,
42     x  .31746031746d-01, .42857142857d-01, .19047619048d-01,
43     x  .21645021645d-01, .71428571429d-01, .19047619048d-01,
44     x  .12987012987d-01, .16650016650d-01/
45      f0pol(a,b) = 3*(16*a**6+104*a**5*b+286*a**4*b**2+429*(a*b)
46     +             **3+286*a**2*b**4+104*a*b**5+16*b**6)
47      f2pol(a,b) = 8*a**4 + 52*a**3*b + 143*(a*b)**2 + 52*a*b**3 +
48     +             8*b**4
49      df1pol(a,b) = 8*a**4 + 44*a**3*b + 99*(a*b)**2 + 44*a*b**3 +
50     +              8*b**4
51      df3pol(a,b) = 2*a**2 + 11*a*b + 2*b**2
52c.......................................................................
53c     two-electron integral routine for lcgo atom scf.
54c     restricted to principal quantum numbers 1,2 and 3 for respectively
55c     s,p, and d orbitals.
56c.......................................................................
57      pi = 3.14159265d0
58      on6 = 1.0d0/6.0d0
59      on15 = 1.0d0/15.d0
60      on70 = 1.0d0/70.d0
61      on35 = 2*on70
62      pifac(1) = sqrt(3.14159265d0)
63      twopow(1) = 1
64      do i = 2 , 14
65         twopow(i) = twopow(i-1)*0.5d0
66         pifac(i) = pifac(i-1)*(2*i-1)
67      enddo
68c.......................................................................
69c
70c     this part sets up the coefficients lambda,p,q and mu,r,s.
71c.......................................................................
72      j = 0
73      nstep1 = 0
74      kmx = 0
75      factkl = 1
76      kl = 0
77      pot = 0.0d0
78      potn = 0.0d0
79      cin = 0.0d0
80      do i = 1 , nsym
81         prfac1 = twopow(i+1)*factkl
82         kin = kmx + 1
83         kmx = kin + nbas(i) - 1
84         do k = kin , kmx
85            zp = zeta(k)
86            do l = kin , k
87               kl = kl + 1
88               pcap(kl) = 0.0d0
89               qcap(kl) = 0.0d0
90               zq = zeta(l)
91               zpq = zp + zq
92               prfac2 = prfac1*u(kl)
93               xfac1 = prfac2*zpq**i*pifac(i)
94               nstep2 = 0
95               mmx = 0
96               factmn = 1.0d0
97               mn = 0
98               do im = 1 , i
99                  open = (nosh(i).ne.0 .and. nosh(im).ne.0)
100                  prfac3 = prfac2*pifac(im)*twopow(im)*factmn
101                  xfac2 = xfac1*factmn*twopow(im+1)
102                  min = mmx + 1
103                  mmx = min - 1 + nbas(im)
104                  mmxp = mmx
105                  if (im.eq.i) mmxp = k
106                  do m = min , mmxp
107                     zr = zeta(m)
108                     zpr = zp + zr
109                     zqr = zq + zr
110                     nmx = l
111                     if (m.lt.k) nmx = m
112                     do n = min , nmx
113                        mn = mn + 1
114                        klnemn = (kl.ne.mn)
115                        j = j + 1
116c.......................................................................
117c
118c     j is the number label of the matrix elements to be calculated
119c     i=lambda+1,k=p,l=q,im=mu+1,m=r,n=s
120c.......................................................................
121                        zs = zeta(n)
122                        zrs = zr + zs
123                        zqs = zq + zs
124                        zps = zp + zs
125                        zpqrs = zpq + zrs
126                        zpqrs2 = 2*zpqrs**2
127                        zprzqs = zpr*zqs
128                        zpszqr = zps*zqr
129                        xterm = (1.0d0/sqrt(zpqrs))**(2*(i+im)-3)
130                        prfac4 = prfac3*u(mn)*xterm
131                        xfac3 = xfac2*u(mn)*xterm
132                        xfac11 = xfac3*(zrs/zprzqs)**im
133                        xfac21 = xfac3*(zrs/zpszqr)**im
134                        xfsum = xfac11 + xfac21
135                        ntest = i*(i-1)/2 + im
136                        go to (30,40,50,60,70,80,90,100,110,120,140) ,
137     +                         ntest
138c.......................................................................
139c
140c     i=1,im=1,(ss)-loop. x0=j0(ss),y0=k0(ss)
141c.......................................................................
142 30                     x0 = prfac4
143                        y0 = xfsum
144                        pj = x0 - y0*ss0
145                        qj = -ajmn(1)*y0
146                        go to 130
147c.......................................................................
148c
149c     i=2,im=1,(sp)-loop. x0=j0(sp),y1=k1(sp)
150c.......................................................................
151 40                     x0 = prfac4*(3*zpq+2*zrs)
152                        y1 = xfsum
153                        pj = x0 - y1*sp1
154                        qj = -ajmn(2)*y1
155                        go to 130
156c.......................................................................
157c
158c     i=2,im=2,(pp)-loop. x0=j0(pp),y0=k0(pp),y2=k2(pp)
159c.......................................................................
160 50                     x0 = prfac4*(zpqrs2+zpq*zrs)
161                        y0 = xfsum*zpqrs2
162                        xfsum = xfac11*zpr*zqs + xfac21*zps*zqr
163                        y0 = y0 + xfsum
164                        y2 = xfsum*5
165                        pj = x0 - y0*pp0 - y2*pp2
166                        qj = -(ajmn(3)*y0+ajmn(4)*y2)
167                        go to 130
168c.......................................................................
169c
170c     i=3,im=1,(sd)-loop. x0=j0(sd),y2=k2(sd)
171c.......................................................................
172 60                     x0 = prfac4*(15*zpq**2+20*zpq*zrs+8*zrs**2)
173                        y2 = xfsum
174                        pj = x0 - y2*sd2
175                        qj = -ajmn(5)*y2
176                        go to 130
177c.......................................................................
178c
179c     i=3,im=2,(pd)-loop. x0=j0(pd),y1=k1(pd),y3=k3(pd),x2=j2(pd)
180c.......................................................................
181 70                     x0 = prfac4*(10*zpq**3+35*zpq**2*zrs+
182     +                       28*zpq*zrs**2+8*zrs**3)
183                        y1 = xfsum*zpqrs2
184                        xfsum = xfac11*zpr*zqs + xfac21*zps*zqr
185                        y1 = y1 + 3*xfsum
186                        y3 = 7*xfsum
187                        pj = x0 - y1*pd1 - y3*pd3
188                        qj = 0.d0
189                        if (open) then
190                           x2 = prfac4*5*zpq*zrs*(7*zpq+2*zrs)
191                           qj = ajmn(21)*x2 - (ajmn(6)*y1+ajmn(7)*y3)
192                        end if
193                        go to 130
194c.......................................................................
195c
196c     i=3,im=3,(dd)-loop. x0=j0(dd),y0=k0(dd),y2=k2(dd),y4=k4(dd)
197c.......................................................................
198 80                     zprzqs = zpr*zqs
199                        zpszqr = zps*zqr
200                        x0 = prfac4*((zpqrs2+zpq*zrs)
201     +                       *zpqrs2*2+7*(zpq*zrs)**2)
202                        y01 = xfac11*((zpqrs2+zprzqs)
203     +                        *zpqrs2*2+7*zprzqs**2)
204                        y02 = xfac21*((zpqrs2+zpszqr)
205     +                        *zpqrs2*2+7*zpszqr**2)
206                        y0 = y01 + y02
207                        xfac11 = xfac11*7*zprzqs
208                        xfac21 = xfac21*7*zpszqr
209                        y21 = xfac11*(zpqrs2+5*zprzqs)
210                        y22 = xfac21*(zpqrs2+5*zpszqr)
211                        y2 = y21 + y22
212                        y4 = (xfac11*zprzqs+xfac21*zpszqr)*9
213                        pj = x0 - y0*dd0 - (y2+y4)*dd2
214                        qj = -ajmn(8)*y0 - ajmn(9)*y2 - ajmn(10)*y4
215                        go to 130
216c.......................................................................
217c
218c     i=4,im=1,(sf)-loop. x0=j0(sf),y3=k3(sf)
219c.......................................................................
220 90                     x0 = prfac4*(105*zpq**3+210*zpq**2*zrs+
221     +                       168*zpq*zrs**2+48*zrs**3)
222                        y3 = xfsum
223                        pj = x0 - y3*sf3
224                        qj = -ajmn(11)*y3
225                        go to 130
226c.......................................................................
227c
228c     i=4,im=2,(pf)-loop. x0=j0(pf),y2=k2(pf),y4=k4(pf),x2=j2(pf)
229c.......................................................................
230 100                    x0 = prfac4*(70*zpq**4+315*zpq**3*zrs+
231     +                       378*(zpq*zrs)**2+216*zpq*zrs**3+48*zrs**4)
232c.......................................................................
233c                       y2 = xfsum*zpqrs2
234c.......................................................................
235                        y2 = xfac11*(2*zpr**2+9*zpr*zqs+2*zqs**2)
236     +                       + xfac21*(2*zps**2+9*zps*zqr+2*zqr**2)
237                        xfsum = xfac11*zpr*zqs + xfac21*zps*zqr
238c.......................................................................
239c                       y2 = y2 + 5*xfsum
240c.......................................................................
241                        y4 = 9*xfsum
242                        pj = x0 - y2*pf2 - y4*pf4
243                        qj = 0.0d0
244                        if (open) then
245                           x2 = prfac4*5*zpq*zrs*(63*zpq**2+26*zpq*zrs+
246     +                          8*zrs**2)
247                           qj = ajmn(22)*x2 - ajmn(12)*y2 - ajmn(13)*y4
248                        end if
249                        go to 130
250c.......................................................................
251c
252c     i=4,im=3,(df)-loop. x0=j0(df),y1=k1(df),y3=k3(df),y5=k5(df)
253c                         x2=j2(df),x4=j4(df).
254c.......................................................................
255 110                    x0 = prfac4*(56*zpq**5+308*zpq**4*zrs+
256     +                       693*zpq**3*zrs**2+594*zpq**2*zrs**3+
257     +                       264*zpq*zrs**4+48*zrs**5)
258                        y1 = xfac11*df1pol(zpr,zqs)
259     +                       + xfac21*df1pol(zps,zqr)
260                        xfac11 = 9*zpr*zqs*xfac11
261                        xfac21 = 9*zps*zqr*xfac21
262                        y3 = xfac11*df3pol(zpr,zqs)
263     +                       + xfac21*df3pol(zps,zqr)
264                        y5 = 11*(xfac11*zpr*zqs+xfac21*zps*zqr)
265                        pj = x0 - y1*df1 - y3*df3 - y5*df5
266                        qj = 0.0d0
267                        if (open) then
268                           prfac4 = prfac4*7*zpq*zrs
269                           x2 = prfac4*(18*zpq**3+99*zpq**2*zrs+
270     +                          44*zpq*zrs**2+8*zrs**3)
271                           x4 = 9*prfac4*(11*zpq+2*zrs)*zpq*zrs
272                           qj = x2*ajmn(23) + x4*ajmn(24) - y1*ajmn(14)
273     +                          - y3*ajmn(15) - y5*ajmn(16)
274                        end if
275                        go to 130
276c.......................................................................
277c
278c     i=4,im=4,(ff)-loop. x0=j0(ff),y0=k0(ff),y2=k2(ff),y4=k4(ff)
279c              y6=k6(ff)
280c.......................................................................
281 120                    x0 = prfac4*f0pol(zpq,zrs)
282                        y0 = xfac11*f0pol(zpr,zqs)
283     +                       + xfac21*f0pol(zps,zqr)
284                        xfac11 = xfac11*zpr*zqs*9
285                        xfac21 = xfac21*zps*zqr*9
286                        y2 = xfac11*f2pol(zpr,zqs)
287     +                       + xfac21*f2pol(zps,zqr)
288                        xfac11 = xfac11*zpr*zqs*11
289                        xfac21 = xfac21*zps*zqr*11
290                        y4 = xfac11*(2*zpr**2+13*zpr*zqs+2*zqs**2)
291     +                       + xfac21*(2*zps**2+13*zps*zqr+2*zqr**2)
292                        y6 = 13*(xfac11*zpr*zqs+xfac21*zps*zqr)
293                        pj = x0 - y0*ff0 - y2*ff2 - y4*ff4 - y6*ff6
294                        qj = -y0*ajmn(17) - y2*ajmn(18) - y4*ajmn(19)
295     +                       - y6*ajmn(20)
296 130                    pcap(kl) = pcap(kl) + pj*dt(mn)
297                        if (klnemn) pcap(mn) = pcap(mn) + pj*dt(kl)
298                        term = dt(kl)*pj*dt(mn)
299                        if (open) then
300                           qcap(kl) = qcap(kl) + qj*dos(mn)
301                           if (klnemn) qcap(mn) = qcap(mn) + qj*dos(kl)
302                           term = term - dos(kl)*qj*dos(mn)
303                        end if
304                        if (.not.klnemn) term = term*0.5d0
305                        pot = pot + term
306 140                 continue
307                     enddo
308                  enddo
309                  factmn = factmn/im
310               enddo
311               potn = potn + u(kl)*dt(kl)
312               cin = cin + t(kl)*dt(kl)
313            enddo
314         enddo
315         factkl = factkl/i
316      enddo
317      pot = pot - zn*potn
318      energ = cin + pot
319      vir = pot/cin
320      return
321      end
322