1 subroutine hnd_giaoh11(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj, 2 & expj,coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,scr,lscr, 3 & para,dia) 4c 5c $Id$ 6c 7c ----- Wrapper routine that sets the sizes of scratch blocks ----- 8c 9 implicit double precision (a-h,o-z) 10#include "hnd_pointers.fh" 11 dimension scr(lscr) 12 logical para,dia 13 dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 14 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 15 dimension e11(nint,3,3,*) 16 dimension catms(3,nat) 17c 18c Use scratch for temps in integral calculation 19c Scratch needs are 20c 11*3*(Li+2)*(Lj+2)*((Li+Lj+3)/2+1) 21c 22c The pointers are defined in hdn_pointers.fh 23c 24 call hnd_giaoh111(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,expj, 25 1 coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,para,dia, 26 2 scr(gh01(1)),scr(gh01(2)),scr(gh01(3)),scr(gh01(4)), 27 3 scr(gh01(5)),scr(gh01(6)),scr(gh01(7)),scr(gh01(8)), 28 4 scr(gh01(9)),scr(gh01(10)),scr(gh01(11))) 29c 30 return 31 end 32c 33 subroutine hnd_giaoh111(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj, 34 1 expj,coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,para,dia, 35 4 v0,dv0,v0d,v1,dv1,v11,dv11,v1d,v12,dv12,v11d) 36c 37 implicit double precision (a-h,o-z) 38#include "nwc_const.fh" 39#include "hnd_rys.fh" 40#include "hnd_tol.fh" 41#include "stdio.fh" 42 common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj, 43 1 ni,nj,cx,cy,cz 44 dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 45 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 46 dimension e11(nint,3,3,*) 47 dimension tm(21),catms(3,nat) 48 dimension v0(3,Li+2,Lj+2,*) ! blocks for vx0,vy0,vz0 49 dimension dv0(3,Li+2,Lj+2,*) ! blocks for dvx0,dvy0,dvz0 50 dimension v0d(3,Li+2,Lj+2,*) ! blocks for vx0d,vy0d,vz0d 51 dimension v1(3,Li+2,Lj+2,*) ! blocks for vx1,vy1,vz1 52 dimension dv1(3,Li+2,Lj+2,*) ! blocks for dvx1,dvy1,dvz1 53 dimension v11(3,Li+2,Lj+2,*) ! blocks for vx1x,vy1y,vz1z 54 dimension dv11(3,Li+2,Lj+2,*) ! blocks for dvx1x,dvy1y,dvz1z 55 dimension v1d(3,Li+2,Lj+2,*) ! blocks for vx1d,vy1d,vz1d 56 dimension v12(3,Li+2,Lj+2,*) ! blocks for xvx1,yvy1,zvz1 57 dimension dv12(3,Li+2,Lj+2,*) ! blocks for xdvx1,ydvy1,zdvz1 58 dimension v11d(3,Li+2,Lj+2,*) ! blocks for xvx1d,yvy1d,zvz1d 59 logical para,dia 60 data rln10 /2.30258d+00/ 61 data pi212 /1.1283791670955d+00/ 62c 63 tol=rln10*itol 64c 65c Zero integral array 66c 67 call dcopy(nint*9*nat,0.0d0,0,e11,1) 68c 69c ----- ishell ----- 70c 71 xi=xyzi(1) 72 yi=xyzi(2) 73 zi=xyzi(3) 74 lit = Li + 1 75 maxi = lit*(lit+1)/2 76 litmod=lit+1 77c 78c ----- jshell ----- 79c 80 xj=xyzj(1) 81 yj=xyzj(2) 82 zj=xyzj(3) 83 ljt = Lj + 1 84 maxj = ljt*(ljt+1)/2 85 ljtmod=ljt+1 86c 87 qijx=yi*zj-zi*yj 88 qijy=zi*xj-xi*zj 89 qijz=xi*yj-yi*xj 90 tijx=xi - xj 91 tijy=yi - yj 92 tijz=zi - zj 93c 94 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 95c 96 nroots=(lit+ljt+3-2)/2+1 97 if(nroots.gt.maxrys) then 98 write(luout,9999) maxrys,lit,ljt,nroots 99 call errquit('hnd_giaoh11: need higher Rys rt',nroots,INT_ERR) 100 endif 101c 102c ----- i primitive ----- 103c 104 do ig=1, i_nprim 105 ai=expi(ig) 106 arri=ai*rr 107 axi=ai*xi 108 ayi=ai*yi 109 azi=ai*zi 110 csi=coefi(ig,i_ngen) 111c 112c ----- j primitive ----- 113c 114 do jg=1,j_nprim 115 aj=expj(jg) 116 aa=ai+aj 117 aa1=1.0d0/aa 118 dum=aj*arri*aa1 119 if(dum.gt.tol) goto 1000 120 fac= exp(-dum) 121 csj=coefj(jg,j_ngen) 122 ax=(axi+aj*xj)*aa1 123 ay=(ayi+aj*yj)*aa1 124 az=(azi+aj*zj)*aa1 125c 126c ----- density factor ----- 127c 128 cij=csi*csj*fac*pi212*aa1 129c 130c ----- -h11- dia + paramagnetic terms ----- 131c 132 aax=aa*ax 133 aay=aa*ay 134 aaz=aa*az 135 do kat=1,nat 136 cx=catms(1,kat) 137 cy=catms(2,kat) 138 cz=catms(3,kat) 139 xx=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 140 yy=xx 141 call hnd_droot 142 do ir=1,nroots 143 uu=u9(ir)*aa 144 ww=w9(ir) 145 vv=ww 146 ww=ww*(uu+uu) 147 tt=1.0d0/(aa+uu) 148 t = sqrt(tt) 149 x0=(aax+uu*cx)*tt 150 y0=(aay+uu*cy)*tt 151 z0=(aaz+uu*cz)*tt 152 do j=1,ljtmod 153 nj=j 154 do i=1,litmod 155 ni=i 156 call hnd_sxyz 157c 158c ----- for 1/r ----- 159c 160 if(kat.eq.1) then 161 v0(1,i,j,ir)=xint ! vx0 162 v0(2,i,j,ir)=yint ! vy0 163 v0(3,i,j,ir)=zint*vv ! vz0 164 endif 165c 166c ----- for x/r**3 ----- 167c 168 v1(1,i,j,ir)=xint ! vx1 169 v1(2,i,j,ir)=yint ! vy1 170 v1(3,i,j,ir)=zint*ww ! vz1 171 call hnd_dervxyz(1) 172 dv1(1,i,j,ir)=xint ! dvx1 173 dv1(2,i,j,ir)=yint ! dvy1 174 dv1(3,i,j,ir)=zint*ww ! dvz1 175 enddo ! i-loop 176 enddo ! j-loop 177 enddo ! ir-loop 178c 179 do ir=1,nroots 180 do j=1,ljt 181 do i=1,lit 182 v11(1,i,j,ir)= v1(1,i,j+1,ir) ! vx1x 183 v11(2,i,j,ir)= v1(2,i,j+1,ir) ! vy1y 184 v11(3,i,j,ir)= v1(3,i,j+1,ir) ! vz1z 185 dv11(1,i,j,ir)=dv1(1,i,j+1,ir) ! dvx1x 186 dv11(2,i,j,ir)=dv1(2,i,j+1,ir) ! dvy1y 187 dv11(3,i,j,ir)=dv1(3,i,j+1,ir) ! dvz1z 188 enddo ! i-loop 189 enddo ! j-loop 190 do j=1,ljtmod 191 do i=1,lit 192 v12(1,i,j,ir)= v1(1,i+1,j,ir) ! xvx1 193 v12(2,i,j,ir)= v1(2,i+1,j,ir) ! yvy1 194 v12(3,i,j,ir)= v1(3,i+1,j,ir) ! zvz1 195 dv12(1,i,j,ir)=dv1(1,i+1,j,ir) ! xdvx1 196 dv12(2,i,j,ir)=dv1(2,i+1,j,ir) ! ydvy1 197 dv12(3,i,j,ir)=dv1(3,i+1,j,ir) ! zdvz1 198 enddo ! i-loop 199 enddo ! j-loop 200 enddo ! ir-loop 201c 202 do ir=1,nroots 203c 204c ----- for 1/r derivatives wrt. centers -i- and -j- ----- 205c 206 if (kat.eq.1) then 207c 208c ----- derivatives with respect to xj ... ----- 209c 210 do i=1,lit 211 v0d(1,i,1,ir)=-(-(aj+aj)*v0(1,i,2,ir)) ! vx0d and vx0 212 v0d(2,i,1,ir)=-(-(aj+aj)*v0(2,i,2,ir)) ! vy0d and vy0 213 v0d(3,i,1,ir)=-(-(aj+aj)*v0(3,i,2,ir)) ! vz0d and vz0 214 do j=2,ljt 215 v0d(1,i,j,ir)=-(dble(j-1)*v0(1,i,j-1,ir)- ! vx0d and vx0 216 & (aj+aj)*v0(1,i,j+1,ir)) 217 v0d(2,i,j,ir)=-(dble(j-1)*v0(2,i,j-1,ir)- ! vy0d and vy0 218 & (aj+aj)*v0(2,i,j+1,ir)) 219 v0d(3,i,j,ir)=-(dble(j-1)*v0(3,i,j-1,ir)- ! vz0d and vz0 220 & (aj+aj)*v0(3,i,j+1,ir)) 221 enddo ! j-loop 222 enddo ! i-loop 223c 224c ----- derivatives with respect to xi ... ----- 225c 226 do j=1,ljt 227 dv0(1,1,j,ir)=-(-(ai+ai)*v0(1,2,j,ir)) ! dvx0 and vx0 228 dv0(2,1,j,ir)=-(-(ai+ai)*v0(2,2,j,ir)) ! dvy0 and vy0 229 dv0(3,1,j,ir)=-(-(ai+ai)*v0(3,2,j,ir)) ! dvz0 and vz0 230 do i=2,lit 231 dv0(1,i,j,ir)=-(dble(i-1)*v0(1,i-1,j,ir)- ! dvx0 and vx0 232 & (ai+ai)*v0(1,i+1,j,ir)) 233 dv0(2,i,j,ir)=-(dble(i-1)*v0(2,i-1,j,ir)- ! dvy0 and vy0 234 & (ai+ai)*v0(2,i+1,j,ir)) 235 dv0(3,i,j,ir)=-(dble(i-1)*v0(3,i-1,j,ir)- ! vz0d and vz0 236 & (ai+ai)*v0(3,i+1,j,ir)) 237 enddo ! i-loop 238 enddo ! j-loop 239 endif 240c 241c ----- d/dx ... operators ----- 242c 243 do i=1,lit 244 v1d(1,i,1,ir)= (-(aj+aj)*v1(1,i,2,ir)) ! vx1d and vx1 245 v1d(2,i,1,ir)= (-(aj+aj)*v1(2,i,2,ir)) ! vy1d and vy1 246 v1d(3,i,1,ir)= (-(aj+aj)*v1(3,i,2,ir)) ! vz1d and vz1 247 do j=2,ljt 248 v1d(1,i,j,ir)= (dble(j-1)*v1(1,i,j-1,ir)- ! vx1d and vx1 249 & (aj+aj)*v1(1,i,j+1,ir)) 250 v1d(2,i,j,ir)= (dble(j-1)*v1(2,i,j-1,ir)- ! vy1d and vy1 251 & (aj+aj)*v1(2,i,j+1,ir)) 252 v1d(3,i,j,ir)= (dble(j-1)*v1(3,i,j-1,ir)- ! vz1d and vz1 253 & (aj+aj)*v1(3,i,j+1,ir)) 254 enddo ! j-loop 255 enddo ! i-loop 256c 257 do i=1,lit 258 v11d(1,i,1,ir)= (-(aj+aj)*v12(1,i,2,ir)) ! xvx1d and xvx1 259 v11d(2,i,1,ir)= (-(aj+aj)*v12(2,i,2,ir)) ! yvy1d and yvy1 260 v11d(3,i,1,ir)= (-(aj+aj)*v12(3,i,2,ir)) ! zvz1d and zvz1 261 do j=2,ljt 262 v11d(1,i,j,ir)= (dble(j-1)*v12(1,i,j-1,ir)- ! xvx1d and xvx1 263 & (aj+aj)*v12(1,i,j+1,ir)) 264 v11d(2,i,j,ir)= (dble(j-1)*v12(2,i,j-1,ir)- ! yvy1d and yvy1 265 & (aj+aj)*v12(2,i,j+1,ir)) 266 v11d(3,i,j,ir)= (dble(j-1)*v12(3,i,j-1,ir)- ! zvz1d and zvz1 267 & (aj+aj)*v12(3,i,j+1,ir)) 268 enddo ! j-loop 269 enddo ! i-loop 270c 271 enddo ! ir-loop 272c 273 ij=0 274 do j=1,maxj 275 call getNxyz(Lj,j,Nxyz) 276 jx = Nxyz(1) + 1 277 jy = Nxyz(2) + 1 278 jz = Nxyz(3) + 1 279 do i=1,maxi 280 call getNxyz(Li,i,Nxyz) 281 ix = Nxyz(1) + 1 282 iy = Nxyz(2) + 1 283 iz = Nxyz(3) + 1 284 transx=0.0d0 285 transy=0.0d0 286 transz=0.0d0 287 call dcopy(21,0.0d0,0,tm,1) 288 do ir=1,nroots 289c 290 if(kat.eq.1) then 291c 292c ----- translation invariance for nuclear derivatives of 1/r ----- 293c 294 transx=transx 295 1 +dv0(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0(3,iz,jz,ir) 296 2 +dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1(3,iz,jz,ir) 297 3 +v0d(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0(3,iz,jz,ir) 298 transy=transy 299 1 +v0(1,ix,jx,ir)*dv0(2,iy,jy,ir)*v0(3,iz,jz,ir) 300 2 +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir) 301 3 +v0(1,ix,jx,ir)*v0d(2,iy,jy,ir)*v0(3,iz,jz,ir) 302 transz=transz 303 1 +v0(1,ix,jx,ir)*v0(2,iy,jy,ir)*dv0(3,iz,jz,ir) 304 2 +v1(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir) 305 3 +v0(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0d(3,iz,jz,ir) 306 endif 307c 308c ----- for h(1,1) diamagnetic ----- 309c 310c xx = xx + dvx1x * vy1 * vz1 311c xy = xy + vx1x * dvy1 * vz1 312c xz = xz + vx1x * vy1 * dvz1 313c yx = yx + dvx1 * vy1y * vz1 314c yy = yy + vx1 * dvy1y * vz1 315c yz = yz + vx1 * vy1y * vz1 316c zx = zx + dvx1 * vy1 * dvz1 317c zy = zy + vx1 * dvy1 * vz1z 318c zz = zz + vx1 * vy1 * dvz1z 319c 320 if (dia) then 321 tm(1)=tm(1)+dv11(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1(3,iz,jz,ir) 322 tm(2)=tm(2)+v11(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir) 323 tm(3)=tm(3)+v11(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir) 324 tm(4)=tm(4)+dv1(1,ix,jx,ir)*v11(2,iy,jy,ir)*v1(3,iz,jz,ir) 325 tm(5)=tm(5)+v1(1,ix,jx,ir)*dv11(2,iy,jy,ir)*v1(3,iz,jz,ir) 326 tm(6)=tm(6)+v1(1,ix,jx,ir)*v11(2,iy,jy,ir)*dv1(3,iz,jz,ir) 327 tm(7)=tm(7)+dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v11(3,iz,jz,ir) 328 tm(8)=tm(8)+v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v11(3,iz,jz,ir) 329 tm(9)=tm(9)+v1(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv11(3,iz,jz,ir) 330 endif 331c 332c ----- for h(1,1) paramagnetic ----- 333c 334c t10 = t10 + vx1 * dvy1 * vz1d - vx1 * vy1d * dvz1 335c t11 = t11 + vx1d * vy1 * vz1d - dvx1 * vy1 * vz1d 336c t12 = t12 + dvx1 * vy1d * vz1 - vx1d * dvy1 * vz1 337c 338 if (para) then 339 tm(10)=tm(10) 340 1 +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1d(3,iz,jz,ir) 341 2 -v1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv1(3,iz,jz,ir) 342 tm(11)=tm(11) 343 1 +v1d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir) 344 2 -dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1d(3,iz,jz,ir) 345 tm(12)=tm(12) 346 1 +dv1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v1(3,iz,jz,ir) 347 2 -v1d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir) 348c 349c t13 = t10 + xvx1 * dvy1 * vz1d - xvx1 * vy1d * dvz1 350c t14 = t11 + xvx1d * vy1 * dvz1 - xdvx1 * vy1 * vz1d 351c t15 = t12 + xdvx1 * vy1d * vz1 - xvx1d * dvy1 * vz1 352c t16 = t12 + vx1 * ydvy1 * vz1d - vx1 * yvy1d * dvz1 353c t17 = t12 + vx1d * yvy1 * dvz1 - dvx1 * yvy1 * vz1d 354c t18 = t12 + dvx1 * yvy1d * vz1 - vx1d * ydvy1 * vz1 355c t19 = t12 + vx1 * dvy1 * zvz1d - vx1 * vy1d * zdvz1 356c t20 = t12 + vx1d * vy1 * zdvz1 - dvx1 * vy1 * zvz1d 357c t21 = t12 + dvx1 * vy1d * zvz1 - vx1d * dvy1 * zvz1 358c 359 tm(13)=tm(13) 360 1 +v12(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1d(3,iz,jz,ir) 361 2 -v12(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv1(3,iz,jz,ir) 362 tm(14)=tm(14) 363 1 +v11d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir) 364 2 -dv12(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1d(3,iz,jz,ir) 365 tm(15)=tm(15) 366 1 +dv12(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v1(3,iz,jz,ir) 367 2 -v11d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir) 368 tm(16)=tm(16) 369 1 +v1(1,ix,jx,ir)*dv12(2,iy,jy,ir)*v1d(3,iz,jz,ir) 370 2 -v1(1,ix,jx,ir)*v11d(2,iy,jy,ir)*dv1(3,iz,jz,ir) 371 tm(17)=tm(17) 372 1 +v1d(1,ix,jx,ir)*v12(2,iy,jy,ir)*dv1(3,iz,jz,ir) 373 2 -dv1(1,ix,jx,ir)*v12(2,iy,jy,ir)*v1d(3,iz,jz,ir) 374 tm(18)=tm(18) 375 1 +dv1(1,ix,jx,ir)*v11d(2,iy,jy,ir)*v1(3,iz,jz,ir) 376 2 -v1d(1,ix,jx,ir)*dv12(2,iy,jy,ir)*v1(3,iz,jz,ir) 377 tm(19)=tm(19) 378 1 +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v11d(3,iz,jz,ir) 379 2 -v1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv12(3,iz,jz,ir) 380 tm(20)=tm(20) 381 1 +v1d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv12(3,iz,jz,ir) 382 2 -dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v11d(3,iz,jz,ir) 383 tm(21)=tm(21) 384 1 +dv1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v12(3,iz,jz,ir) 385 2 -v1d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v12(3,iz,jz,ir) 386 endif 387c 388 enddo 389 ij=ij+1 390c 391c ----- h(1,1) diamagnetic ----- 392c 393 if (dia) then 394 e11(ij,1,1,kat)=e11(ij,1,1,kat)+(tm(5)+tm(9))*cij 395 e11(ij,1,2,kat)=e11(ij,1,2,kat)-tm(2)*cij 396 e11(ij,1,3,kat)=e11(ij,1,3,kat)-tm(3)*cij 397 e11(ij,2,1,kat)=e11(ij,2,1,kat)-tm(4)*cij 398 e11(ij,2,2,kat)=e11(ij,2,2,kat)+(tm(1)+tm(9))*cij 399 e11(ij,2,3,kat)=e11(ij,2,3,kat)-tm(6)*cij 400 e11(ij,3,1,kat)=e11(ij,3,1,kat)-tm(7)*cij 401 e11(ij,3,2,kat)=e11(ij,3,2,kat)-tm(8)*cij 402 e11(ij,3,3,kat)=e11(ij,3,3,kat)+(tm(1)+tm(5))*cij 403 endif 404c 405c ----- h(1,1) paramagnetic ----- 406c 407c ----- a bit unclear ... yx ... instead of ... xy ... ----- 408c 409 if (para) then 410 tm(1)=(qijx*tm(10)+tijy*tm(19)-tijz*tm(16))*cij 411 tm(2)=(qijx*tm(11)+tijy*tm(20)-tijz*tm(17))*cij 412 tm(3)=(qijx*tm(12)+tijy*tm(21)-tijz*tm(18))*cij 413 tm(4)=(qijy*tm(10)+tijz*tm(13)-tijx*tm(19))*cij 414 tm(5)=(qijy*tm(11)+tijz*tm(14)-tijx*tm(20))*cij 415 tm(6)=(qijy*tm(12)+tijz*tm(15)-tijx*tm(21))*cij 416 tm(7)=(qijz*tm(10)+tijx*tm(16)-tijy*tm(13))*cij 417 tm(8)=(qijz*tm(11)+tijx*tm(17)-tijy*tm(14))*cij 418 tm(9)=(qijz*tm(12)+tijx*tm(18)-tijy*tm(15))*cij 419c 420 e11(ij,1,1,kat)=e11(ij,1,1,kat)+tm(1) 421 e11(ij,1,2,kat)=e11(ij,1,2,kat)+tm(4) 422 e11(ij,1,3,kat)=e11(ij,1,3,kat)+tm(7) 423 e11(ij,2,1,kat)=e11(ij,2,1,kat)+tm(2) 424 e11(ij,2,2,kat)=e11(ij,2,2,kat)+tm(5) 425 e11(ij,2,3,kat)=e11(ij,2,3,kat)+tm(8) 426 e11(ij,3,1,kat)=e11(ij,3,1,kat)+tm(3) 427 e11(ij,3,2,kat)=e11(ij,3,2,kat)+tm(6) 428 e11(ij,3,3,kat)=e11(ij,3,3,kat)+tm(9) 429 endif 430c 431c if(kat.eq.1) then 432c if(abs(transx*cij).ge.1.0d-08.or.abs(transy*cij) 433c & .ge.1.0d-08.or.abs(transz*cij).ge.1.0d-08) then 434c write(luout,9993) 435c call errquit('hnd_giaoh11: no transl inv',0,INT_ERR) 436c endif 437c endif 438c 439 enddo ! j-loop final summation 440 enddo ! i-loop final summation 441c 442 enddo ! kat-loop 443c 444 1000 continue 445 enddo ! jprim loop 446 enddo ! iprim loop 447c 448 return 449 9999 format(' in -giaoh11- , the rys quadrature is not implemented', 450 1 ' beyond -nroots- = ',i3,/,' lit,ljt,nroots= ',3i3) 451 9993 format(' something wrong with translational', 452 1 ' invariance in -giaoh11- . stop. ') 453 end 454