1 subroutine dim_elfder(xyzi,expi,coefi,i_nprim,i_ngen, Li, 2 1 xyzj,expj,coefj,j_nprim,j_ngen,Lj,nder,nint,elfder,scr,lscr, 3 2 xyzpt, npt) 4c 5c $Id: hnd_elfder.F 21348 2011-10-31 22:51:25Z d3p852 $ 6c 7 implicit double precision (a-h,o-z) 8#include "hnd_pointers.fh" 9 dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 10 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 11 dimension xyzpt(3,npt) 12 dimension scr(lscr) 13 dimension elfder(*) 14c 15c 16c ----- Wrapper routine that sets the sizes of scratch blocks ----- 17c 18 idim = max(nder*3,1) 19c 20 call dim_elfder1(xyzi,expi,coefi,i_nprim,i_ngen, Li, xyzj,expj, 21 1 coefj, j_nprim, j_ngen, Lj, nder, nint, elfder, xyzpt, npt, idim, 22 2 scr(elpt(1)),scr(elpt(2)),scr(elpt(3)),scr(elpt(4)),scr(elpt(5)), 23 3 scr(elpt(6)),scr(elpt(7)),scr(elpt(8)),scr(elpt(9))) 24c 25 return 26 end 27c 28 subroutine dim_elfder1(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj, 29 1 expj,coefj, j_nprim, j_ngen, Lj, nder, nint, elfder, xyzpt, npt, 30 2 idim,xv,yv,zv,dxv,dyv,dzv,ddxv,ddyv,ddzv) 31c 32 implicit double precision (a-h,o-z) 33#include "nwc_const.fh" 34#include "hnd_rys.fh" 35#include "stdio.fh" 36#include "hnd_tol.fh" 37#include "dimqm.fh" 38 common/hnd_xyzder/xint,yint,zint,tx,x0,y0,z0,xi,yi,zi,xj,yj,zj, 39 1 ni,nj,cx,cy,cz 40 dimension w2(maxrys),w4(maxrys) 41 dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 42 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 43 dimension xyzpt(3,npt),elfder(nint,idim,npt) 44c 45 dimension xv(Li+1,Lj+1,*) 46 dimension yv(Li+1,Lj+1,*) 47 dimension zv(Li+1,Lj+1,*) 48 dimension dxv(Li+1,Lj+1,*) 49 dimension dyv(Li+1,Lj+1,*) 50 dimension dzv(Li+1,Lj+1,*) 51 dimension ddxv(Li+1,Lj+1,*) 52 dimension ddyv(Li+1,Lj+1,*) 53 dimension ddzv(Li+1,Lj+1,*) 54c 55c Routine calculates the following integrals for a given shell: 56c 57c nder =-2 : Electronic wave function 58c nder =-1 : Electronic density 59c nder = 0 : Electrostatic potential 60c nder = 1 : Electric field 61c nder = 2 : Electric field gradient 62c 63c Maximum scratch size needs to be Max(nder*3,1))*(Li+1)*(Lj+1)*((Li+Lj)/2+1) 64c Data format elfder will be (nints,ipts,Max(nder*3,1)) 65c 66 data rln10 /2.30258d+00/ 67 data pi212 /1.1283791670955d+00/ ! 2/sqrt(pi) 68c 69 dtol=rln10*itol 70c 71c Zero integral array 72c 73 call dcopy(nint*max(nder*3,1)*npt,0.0d0,0,elfder,1) 74 xintmax = 0.0d0 75 yintmax = 0.0d0 76 zintmax = 0.0d0 77c 78c ----- ishell ----- 79c 80 xi=xyzi(1) 81 yi=xyzi(2) 82 zi=xyzi(3) 83 lit = Li + 1 84 maxi = lit*(lit+1)/2 85c 86c ----- jshell ----- 87c 88 xj=xyzj(1) 89 yj=xyzj(2) 90 zj=xyzj(3) 91 ljt = Lj + 1 92 maxj = ljt*(ljt+1)/2 93c write(luout,*) "i shell:", xi, yi, zi 94c write(luout,*) "j shell:", xj, yj, zj 95c 96 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 97 nroots=(lit+ljt+nder-2)/2+1 98 99 if(nroots.gt.maxrys) then 100 write(luout,9997) maxrys,lit,ljt,nroots 101 call errquit('hnd_elfgrd: need higher Rys root',nroots,INT_ERR) 102 endif 103c 104c ----- i primitive ----- 105c 106 do 7000 ig=1,i_nprim 107 ai=expi(ig) 108 arri=ai*rr 109 axi=ai*xi 110 ayi=ai*yi 111 azi=ai*zi 112 csi=coefi(ig,i_ngen) 113c write(luout,*) "Li:", Li 114c write(luout,*) "ai:", ai 115c 116c ----- j primitive ----- 117c 118 do 6000 jg=1,j_nprim 119 aj=expj(jg) 120c write(luout,*) "Lj:", Lj 121c write(luout,*) "aj:", aj 122 aa=ai+aj 123c write(luout,*) "aa:", aa 124 aa1=1.0d0/aa 125 dum=aj*arri*aa1 126 if(dum.gt.dtol) go to 6000 127 fac= exp(-dum) 128 csj=coefj(jg,j_ngen) 129 ax=(axi+aj*xj)*aa1 130 ay=(ayi+aj*yj)*aa1 131 az=(azi+aj*zj)*aa1 132c 133c ----- density factor ----- 134c 135 dum1 = csi * fac 136 dij = dum1 * csj 137c 138c Electronic density integrals and electronic wave 139c function are done differently 140c then the electric field, electric field gradient, and 141c electrostatic potential integrals 142c Switch between the two classes here: 143c 144 if (nder.ge.0) goto 399 145c 146c ----- electronic wave function / density integral ----- 147c 148c For wave function we simply set j to 0.0d0 and reset the 149c density factor appropriately 150c 151 if (nder.eq.-2) then 152 dij = csi 153 ax = xi 154 ay = yi 155 az = zi 156 aa = ai 157 ljt = 1 158 maxj = 1 159 endif 160c 161 do 300 ipt=1,npt 162 x0=xyzpt(1,ipt) 163 y0=xyzpt(2,ipt) 164 z0=xyzpt(3,ipt) 165 dum = aa*((x0-ax)**2+(y0-ay)**2+(z0-az)**2) 166 if(dum.gt.dtol) go to 300 167 fac = exp(-dum) 168c 169c ----- density values ----- 170c 171 do 270 j=1,ljt 172 nj=j 173 do 270 i=1,lit 174 ni=i 175 call hnd_denxyz 176 xv(i,j,1)=xint 177 yv(i,j,1)=yint 178 zv(i,j,1)=zint 179 270 continue 180c 181c ----- combining the pieces together ----- 182c 183 ij=0 184 do 290 j=1,maxj 185 call getNxyz(Lj,j,Nxyz) 186 jx = Nxyz(1) + 1 187 jy = Nxyz(2) + 1 188 jz = Nxyz(3) + 1 189 do 290 i=1,maxi 190 call getNxyz(Li,i,Nxyz) 191 ix = Nxyz(1) + 1 192 iy = Nxyz(2) + 1 193 iz = Nxyz(3) + 1 194 ij=ij+1 195 elfder(ij,1,ipt)=elfder(ij,1,ipt)+fac*dij*xv(ix,jx,1)* 196 1 yv(iy,jy,1)*zv(iz,jz,1) 197 290 continue 198 300 continue 199c 200 goto 6000 201c 202c ----- electric field (gradient) term ----- 203c 204 399 dij = dij * pi212 * aa1 205 aax=aa*ax 206 aay=aa*ay 207 aaz=aa*az 208c write(luout,*) "dij:", dij 209c aa = aa/erf(aa) 210 do 500 ipt=1,npt 211 znuc=1.0d0 212 cx=xyzpt(1,ipt) 213 cy=xyzpt(2,ipt) 214 cz=xyzpt(3,ipt) 215 cu = 1.0d0/(2.7306542620197547)**2 216c xci = exp(-ai*(xi-cx)**2/(ai + 1)) 217c yci = exp(-ai*(yi-cy)**2/(ai + 1)) 218c zci = exp(-ai*(zi-cz)**2/(ai + 1)) 219c xcj = exp(-aj*(xj-cx)**2/(aj + 1)) 220c ycj = exp(-aj*(yj-cy)**2/(aj + 1)) 221c zcj = exp(-aj*(zj-cz)**2/(aj + 1)) 222c xci = exp(-(aa*cu/(aa+cu))*(ax-cx)**2) 223c yci = exp(-(aa*cu/(aa+cu))*(ay-cy)**2) 224c zci = exp(-(aa*cu/(aa+cu))*(az-cz)**2) 225c temp = xci*yci*zci 226c damp = erfc(temp) 227c write(luout,*) xci, yci, zci 228c write(luout,*) xci*yci*zci 229c if(temp .ge. 1.0d-2) then 230c write(luout,*) "overlap:", temp 231c write(luout,*) "scale: " , erfc(temp) 232c end if 233 234c write(luout,*) ipt, 'ci' 235c write(luout,*) xci, yci, zci 236c write(luout,*) ipt, 'cj' 237c write(luout,*) xcj, ycj, zcj 238c 239c Testing DIM atom as Gaussian 240c rr = (ax-cx)**2+(ay-cy)**2+(az-cz)**2 241c aac = aa+ac 242c dum = aa*rr*ac 243c dijc = dij * exp(-dum/aac) 244 245c temp = sqrt((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 246c if(temp .le. 10.0) write(luout,*) temp 247 yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 248c if(yy .le. 20.0) write(luout,*) "yy:", yy 249c if(yy .le. 20.0) write(luout,*) "aa:", aa 250c if(yy .le. 20.0) write(luout,*) "temp:", temp 251c if(yy .le. 20.0) write(luout,*) "nroots:", nroots 252 call hnd_droot 253 do 420 iroot=1,nroots 254 uu=u9(iroot)*aa 255c write(luout,*) "uu:", uu 256c uu = erfc(uu) * uu 257c uu = 1 - exp(-uu) 258c ddfac = 1.0d0 259c if(uu .ge. 0.1d0) then 260c write(luout,*) "uu:", uu 261c ddfac = 0.0d0 262c end if 263 u4=uu*uu 264 ww=w9(iroot)*znuc 265 w2(iroot)=ww*uu*2.0d0 266 w4(iroot)=ww*u4*4.0d0 267 w9(iroot)=ww 268 tt=1.0d0/(aa+uu) 269 tx= sqrt(tt) 270 x0=(aax+uu*cx)*tt 271 y0=(aay+uu*cy)*tt 272 z0=(aaz+uu*cz)*tt 273c xci = exp(-(aa*uu*tt)*(ax-cx)**2) 274c yci = exp(-(aa*uu*tt)*(ay-cy)**2) 275c zci = exp(-(aa*uu*tt)*(az-cz)**2) 276c temp = erfc(xci*yci*zci) 277c if(temp .le. 0.5) then 278c write(luout,*) "overlap:", xci*yci*zci 279c end if 280c write(luout,*) x0, y0, z0 281c temp = sqrt((x0-xi)**2+(y0-yi)**2+(z0-zi)**2) 282c if(temp .le. 10.0) write(luout,*) "Ri0:", temp 283c temp = sqrt((x0-xj)**2+(y0-yj)**2+(z0-zj)**2) 284c if(temp .le. 10.0) write(luout,*) "Rj0:", temp 285c write(luout,*) temp 286 287 do 410 i=1,lit 288 ni=i 289 do 410 j=1,ljt 290 nj=j 291 goto (402,401) nder+1 292 call dim_dervxyz(2) 293 ddxv(i,j,iroot)=xint 294 ddyv(i,j,iroot)=yint 295 ddzv(i,j,iroot)=zint 296 401 call dim_dervxyz(1) 297 dxv(i,j,iroot)=xint 298 dyv(i,j,iroot)=yint 299 dzv(i,j,iroot)=zint 300 402 call dim_sxyz 301 xv(i,j,iroot)=xint 302 yv(i,j,iroot)=yint 303 zv(i,j,iroot)=zint 304c write(luout,*) xint, yint, zint 305 410 continue 306 420 continue 307c 308c ----- combining the pieces together ----- 309c 310 ij=0 311 do 440 j=1,maxj 312 call getNxyz(Lj,j,Nxyz) 313 jx = Nxyz(1) + 1 314 jy = Nxyz(2) + 1 315 jz = Nxyz(3) + 1 316 do 450 i=1,maxi 317 call getNxyz(Li,i,Nxyz) 318 ix = Nxyz(1) + 1 319 iy = Nxyz(2) + 1 320 iz = Nxyz(3) + 1 321 ij=ij+1 322 if (nder.eq.2) then 323 dumxx=0.0d0 324 dumyy=0.0d0 325 dumzz=0.0d0 326 dumxy=0.0d0 327 dumxz=0.0d0 328 dumyz=0.0d0 329 do 430 iroot=1,nroots 330 dum=xv(ix,jx,iroot)*yv(iy,jy,iroot)* 331 1 zv(iz,jz,iroot)*w2(iroot) 332 dumxx=dumxx - dum + ddxv(ix,jx,iroot)* 333 1 yv(iy,jy,iroot)*zv(iz,jz,iroot)*w4(iroot) 334 dumyy=dumyy - dum + ddyv(iy,jy,iroot)* 335 1 xv(ix,jx,iroot)*zv(iz,jz,iroot)*w4(iroot) 336 dumzz=dumzz - dum + ddzv(iz,jz,iroot)* 337 1 xv(ix,jx,iroot)*yv(iy,jy,iroot)*w4(iroot) 338 dumxy=dumxy + dxv(ix,jx,iroot)*dyv(iy,jy,iroot)* 339 1 zv(iz,jz,iroot)*w4(iroot) 340 dumxz=dumxz + dxv(ix,jx,iroot)*dzv(iz,jz,iroot)* 341 1 yv(iy,jy,iroot)*w4(iroot) 342 dumyz=dumyz + dyv(iy,jy,iroot)*dzv(iz,jz,iroot)* 343 1 xv(ix,jx,iroot)*w4(iroot) 344 430 continue 345 346 elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumxx*dij 347 elfder(ij,2,ipt) = elfder(ij,2,ipt) + dumyy*dij 348 elfder(ij,3,ipt) = elfder(ij,3,ipt) + dumzz*dij 349 elfder(ij,4,ipt) = elfder(ij,4,ipt) + dumxy*dij 350 elfder(ij,5,ipt) = elfder(ij,5,ipt) + dumxz*dij 351 elfder(ij,6,ipt) = elfder(ij,6,ipt) + dumyz*dij 352 elseif (nder.eq.1) then 353 dumx=0.0d0 354 dumy=0.0d0 355 dumz=0.0d0 356 do 431 iroot=1,nroots 357 dumx=dumx + dxv(ix,jx,iroot)*yv(iy,jy,iroot)* 358 1 zv(iz,jz,iroot)*w2(iroot) 359 dumy=dumy + xv(ix,jx,iroot)*dyv(iy,jy,iroot)* 360 1 zv(iz,jz,iroot)*w2(iroot) 361 dumz=dumz + xv(ix,jx,iroot)*yv(iy,jy,iroot)* 362 1 dzv(iz,jz,iroot)*w2(iroot) 363 431 continue 364 elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumx*dij 365 elfder(ij,2,ipt) = elfder(ij,2,ipt) + dumy*dij 366 elfder(ij,3,ipt) = elfder(ij,3,ipt) + dumz*dij 367 else 368 dumx=0.0d0 369 do 432 iroot=1,nroots 370 dumx=dumx + xv(ix,jx,iroot)*yv(iy,jy,iroot)* 371 1 zv(iz,jz,iroot)*w9(iroot) 372 432 continue 373 elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumx*dij 374 endif 375 450 continue 376 440 continue 377c 378 500 continue 379c 380 6000 continue 381 7000 continue 382c 383 return 384 9997 format(' in -elfgrd- , the rys quadrature is not implemented', 385 1 ' beyond -nroots- = ',i3,/, 386 2 ' lit,ljt,nroots = ',3i3) 387 end 388