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