1c 2c $Id$ 3c 4c Modified from HONDO 5c 6 subroutine hnd_vstat(xyzi,expi,coefi,i_nprim,i_ngen, Li, 7 1 xyzj,expj,coefj, j_nprim, j_ngen, Lj, 8 1 xyz,zan,nat,dens,s,ti,vi,lstv,doS,doT,doV,scr,lscr) 9c 10c This is a wrapper routine, setting up scratch blocks used in actual 11c integral routine 12c 13 implicit none 14c 15 integer i_nprim ! [input] num. prims on function i 16 integer i_ngen ! [input] num general conts on func. i 17 integer Li ! [input] angular momentum of func. i 18 integer j_nprim ! [input] num. prims on function j 19 integer j_ngen ! [input] num general conts on func. j 20 integer Lj ! [input] angular momentum of func. j 21 integer nat ! [input] number of atoms 22 integer lscr ! [input] size of scratch array 23 integer lstv ! [input] size of any integral buffer 24 double precision xyzi(3) ! [input] position of center i 25 double precision expi(i_nprim) ! [input] exponents on i 26 double precision coefi(i_nprim,i_ngen) ! [input] i coeffs 27 double precision xyzj(3) ! [input] position of center j 28 double precision expj(j_nprim) ! [input] exponents on j 29 double precision coefj(j_nprim,j_ngen) ! [input] j coeffs 30 double precision xyz(3,nat) ! [input] all atom positions 31 double precision zan(nat) ! [input] charges on all atoms 32 double precision dens(*) ! [input] (transition)density 33 double precision scr(lscr) ! [scratch] scratch buffers 34 double precision s(lstv) ! [output] overlap integrals 35 double precision ti(lstv) ! [output] kinetic energy integrals 36 double precision vi(lstv) ! [output] potential integrals 37 logical doS ! [input] compute overlap (True/False) 38 logical doT ! [input] compute kinetic (True/False) 39 logical doV ! [input] compute potential (True/False) 40c 41 integer ss,st,tt,tv,vv, dd 42c 43c Use scratch for temps in integral calculation 44c Scratch needs are dependent on integral types: 45c doS : 3*(Li+1)*(Lj+3) 46c doT : 3*(Li+1)*(Lj+3)+3*(Li+1)*(Lj+1) 47c doV : 3*(Li+1)*(Lj+1)*((Li+Lj)/2+1) 48c 49 ss = 0 50 tt = 0 51 vv = 0 52 st = 1 53 tv = 1 54 dd = (Li+1)*(Li+2)*(Lj+1)*(Lj+2)/4 55 if (doS.or.doT) then 56 ss = dd + (Li+1) * (Lj+3) 57 st = st + 3 * ss 58 if (doV) tv = st 59 if (doT) then 60 tt = (Li+1) * (Lj+1) 61 if (doV) tv = tv + 3 * tt 62 endif 63 endif 64 if (doV) vv = (Li+1) * (Lj+1) * ((Li+Lj)/2+1) 65c 66 call hnd_vstat1(xyzi,expi,coefi,i_nprim,i_ngen, Li, 67 1 xyzj,expj,coefj, j_nprim, j_ngen, Lj, 68 1 xyz,zan,nat,dens,s,ti,vi,lstv,doS,doT,doV,scr(1), 69 3 scr(1+dd),scr(1+ss),scr(1+2*ss),scr(st),scr(st+tt), 70 4 scr(st+2*tt),scr(tv),scr(tv+vv),scr(tv+2*vv)) 71c 72 return 73 end 74c 75 subroutine hnd_vstat1(xyzi,expi,coefi,i_nprim,i_ngen, Li, 76 1 xyzj,expj,coefj, j_nprim, j_ngen, Lj, 77 1 xyz,zan,nat,dens,s,ti,vi,lstv,doS,doT,doV,dij, 78 1 xs,ys,zs,xt,yt,zt,xv,yv,zv) 79c 80 implicit none 81#include "stdio.fh" 82#include "hnd_rys.fh" 83#include "hnd_tol.fh" 84#include "errquit.fh" 85 integer i_nprim ! [input] num. prims on function i 86 integer i_ngen ! [input] num general conts on func. i 87 integer Li ! [input] angular momentum of func. i 88 integer j_nprim ! [input] num. prims on function j 89 integer j_ngen ! [input] num general conts on func. j 90 integer Lj ! [input] angular momentum of func. j 91 integer nat ! [input] number of atoms 92 integer lscr ! [input] size of scratch array 93 integer lstv ! [input] size of any integral buffer 94 double precision xyzi(3) ! [input] position of center i 95 double precision expi(i_nprim) ! [input] exponents on i 96 double precision coefi(i_nprim,i_ngen) ! [input] i coeffs 97 double precision xyzj(3) ! [input] position of center j 98 double precision expj(j_nprim) ! [input] exponents on j 99 double precision coefj(j_nprim,j_ngen) ! [input] j coeffs 100 double precision xyz(3,nat) ! [input] all atom positions 101 double precision zan(nat) ! [input] charges on all atoms 102 double precision dens(*) ! [input] (transition)density 103 double precision s(lstv) ! [output] overlap integrals 104 double precision ti(lstv) ! [output] kinetic energy integrals 105 double precision vi(lstv) ! [output] potential integrals 106 logical doS ! [input] compute overlap (True/False) 107 logical doT ! [input] compute kinetic (True/False) 108 logical doV ! [input] compute potential (True/False) 109 double precision tol, aa, aa1, rr, ai, arri, axi, ayi, azi 110 double precision csi, cpi, cdi, cfi, cgi, aj, fac 111 double precision csj, cpj, cdj, cfj, cgj, ax, ay, az 112 double precision xs, ys, zs, xt, yt, zt, dum, dum1, dum2 113 double precision xint, yint, zint, t, x0, y0, z0, xi, yi, zi 114 double precision xj, yj, zj, ccx, ccy, ccz 115 double precision rln10, zero, one, pi212 116 double precision aax, aay, aaz, znuc, cx, cy, cz, uu, ww, tt 117 double precision xv, yv, zv, dij 118 integer lit, ljt, maxi, maxj, ljtmod, jg 119 integer ijx, ijy, ijz, ij, i, j 120 integer ni, nj, ig, ic, iroot 121 integer ix, iy, iz, jx, jy, jz, Nxyz(3) 122 logical some 123 common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj, 124 1 ni,nj,ccx,ccy,ccz 125 dimension xs(Li+1,Lj+3) ,ys(Li+1,Lj+3) ,zs(Li+1,Lj+3) 126 dimension xt(Li+1,Lj+1) ,yt(Li+1,Lj+1) ,zt(Li+1,Lj+1) 127 dimension xv(Li+1,Lj+1,*),yv(Li+1,Lj+1,*),zv(Li+1,Lj+1,*) 128 dimension dij(*) 129 data rln10 /2.30258d+00/ 130 data zero /0.0d+00/ 131 data one /1.0d+00/ 132 data pi212 /1.1283791670955d+00/ 133c 134c write(*,*)"basis function information" 135c write(*,*)"xyz=", xyz, xyzi,xyzj 136c write(*,*)"expi,coefi,i_nprim,i_ngen, Li", 137c & expi, coefi,i_nprim,i_ngen, Li 138c write(*,*)"expj,coefj,j_nprim,j_ngen, Lj", 139c & expj, coefj,j_nprim,j_ngen, Lj 140 141c write(*,*)"vints in hnd_vstat, before" 142c write(*,"(4f12.8)")((xyz(ig,jg),ig=1,3),vi(jg),jg=1,nat) 143c 144 tol=rln10*itol 145c 146c ----- calculate -s-, -t-, and -v- integrals ----- 147c 148 some = .false. 149 if(some) write(luout,9999) 150c 151c ----- ishell ----- 152c 153 xi=xyzi(1) 154 yi=xyzi(2) 155 zi=xyzi(3) 156 lit = Li + 1 157 maxi=lit*(lit+1)/2 158c 159c ----- jshell ----- 160c 161 xj=xyzj(1) 162 yj=xyzj(2) 163 zj=xyzj(3) 164 ljt = Lj + 1 165 maxj=ljt*(ljt+1)/2 166 ljtmod=ljt+2 167c 168 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 169 nroots=(lit+ljt-2)/2+1 170 if(nroots.gt.maxrys) then 171 write(luout,9997) maxrys,lit,ljt,nroots 172 call errquit('hnd_stvint: angular momentum too high',555, 173 & INT_ERR) 174 endif 175c 176 if (doS) call dcopy(lstv,zero,0,s,1) 177 if (doT) call dcopy(lstv,zero,0,ti,1) 178c 179c if (doV) call dcopy(lstv,zero,0,vi,1) 180c 181c do not fill vi with zero's since the contributions are summed over 182c all basis pairs 183c 184c ----- i primitive ----- 185c 186 do 7000 ig=1,i_nprim 187 ai=expi(ig) 188 arri=ai*rr 189 axi=ai*xi 190 ayi=ai*yi 191 azi=ai*zi 192 csi=coefi(ig,i_ngen) 193c 194c ----- j primitive ----- 195c 196 do 6000 jg=1,j_nprim 197 aj=expj(jg) 198 aa=ai+aj 199 aa1=one/aa 200 dum=aj*arri*aa1 201 if(dum.gt.tol) go to 6000 ! the integral is zero 202 fac= exp(-dum) 203 csj=coefj(jg,j_ngen) 204 ax=(axi+aj*xj)*aa1 205 ay=(ayi+aj*yj)*aa1 206 az=(azi+aj*zj)*aa1 207c 208c ----- density factor ----- 209c 210 ij=0 211 dum2=fac*csi*csj 212 do 3600 i=1,maxi 213 do 360 j=1,maxj 214 ij=ij+1 215 dij(ij)=dum2*dens(ij) 216 360 continue 217 3600 continue 218c 219c ----- overlap and kinetic energy ----- 220c 221 if (doS.or.doT) then 222 t = sqrt(aa1) 223 x0=ax 224 y0=ay 225 z0=az 226 do 3700 j=1,ljtmod 227 nj=j 228 do 370 i=1,lit 229 ni=i 230 call hnd_sxyz 231 xs(i,j)=xint*t 232 ys(i,j)=yint*t 233 zs(i,j)=zint*t 234 370 continue 235 3700 continue 236c 237 if (doT) call hnd_txyz(xt,yt,zt,xs,ys,zs,lit,ljt,aj) 238c 239 ij=0 240 do 390 i=1,maxi 241 call getNxyz(Li,i,Nxyz) 242 ix = Nxyz(1) + 1 243 iy = Nxyz(2) + 1 244 iz = Nxyz(3) + 1 245 do 380 j=1,maxj 246 call getNxyz(Lj,j,Nxyz) 247 jx = Nxyz(1) + 1 248 jy = Nxyz(2) + 1 249 jz = Nxyz(3) + 1 250 ij=ij+1 251 dum =xs(ix,jx)*ys(iy,jy)*zs(iz,jz) 252 if (doS) s(ij)= s(ij)+ dum*dij(ij) 253 if (doT) then 254 dum =xt(ix,jx)*ys(iy,jy)*zs(iz,jz) 255 1 +xs(ix,jx)*yt(iy,jy)*zs(iz,jz) 256 2 +xs(ix,jx)*ys(iy,jy)*zt(iz,jz) 257 ti(ij)=ti(ij)+(dum*dij(ij)) 258 endif 259 380 continue 260 390 continue 261 endif ! of overlap and kinetic energy integral 262c 263c ----- nuclear attraction ----- 264c 265 if (doV) then 266 aax=aa*ax 267 aay=aa*ay 268 aaz=aa*az 269 do 500 ic=1,nat 270 znuc=-zan(ic) 271 cx=xyz(1,ic) 272 cy=xyz(2,ic) 273 cz=xyz(3,ic) 274 yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 275 call hnd_droot 276 do 420 iroot=1,nroots 277 uu=u9(iroot)*aa 278 ww=w9(iroot)*znuc 279 tt=one/(aa+uu) 280 t = sqrt(tt) 281 x0=(aax+uu*cx)*tt 282 y0=(aay+uu*cy)*tt 283 z0=(aaz+uu*cz)*tt 284 do 4100 j=1,ljt 285 nj=j 286 do 410 i=1,lit 287 ni=i 288 call hnd_sxyz 289 xv(i,j,iroot)=xint 290 yv(i,j,iroot)=yint 291 zv(i,j,iroot)=zint*ww 292 410 continue 293 4100 continue 294 420 continue 295c 296 ij=0 297 do 450 i=1,maxi 298 call getNxyz(Li,i,Nxyz) 299 ix = Nxyz(1) + 1 300 iy = Nxyz(2) + 1 301 iz = Nxyz(3) + 1 302 do 440 j=1,maxj 303 call getNxyz(Lj,j,Nxyz) 304 jx = Nxyz(1) + 1 305 jy = Nxyz(2) + 1 306 jz = Nxyz(3) + 1 307 dum=zero 308 do 430 iroot=1,nroots 309c write(*,*)"ic,ix,iy,iz,jx,jy,jz,xv,yy,zv", 310c & ic,ix,iy,iz,jx,jy,jz, 311c & xv(ix,jx,iroot),yv(iy,jy,iroot), 312c & zv(iz,jz,iroot) 313 dum=dum+xv(ix,jx,iroot) 314 & *yv(iy,jy,iroot)*zv(iz,jz,iroot) 315 430 continue 316 dum=dum*(aa1*pi212) 317 ij=ij+1 318c write(*,"('den(ij)=', f10.4,' ai=',f10.4, 319c & ' aj=',f10.4 320c & ' ci=',f10.4,' aj=',f10.4,' soo=', 321c & f10.4)") 322c & dens(ij),ai,aj,dum1/fac,dum2/dum1,dij(ij) 323c write(*,"('i=',i5,' j=',i5,' ij=',i5, 324c & ' dum=',f12.8,' dij=',f12.8, 325c & ' dij*dum=',f12.8)") 326c & i,j,ij,dum,dij(ij),dum*dij(ij) 327 vi(ic)=vi(ic)+dum*dij(ij) 328 440 continue 329 450 continue 330c 331 500 continue 332 endif ! of nuclear part 333c 334 6000 continue 335 7000 continue 336c 337 if(some) write(luout,9998) 338 return 339 9999 format(/,10x,20(1h-),/,10x,'1 electron integrals', 340 2 /,10x,20(1h-)) 341 9998 format(' ...... end of one-electron integrals ......') 342 9997 format(' in -hnd_stvint- the rys quadrature is not implemented', 343 1 ' beyond -nroots- = ',i2,/, 344 2 ' lit,ljt,nroots= ',3i3) 345 end 346 347 348 349 350 351 352