1 subroutine hnd_giaotv10(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj, 2 & expj,coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint,tv10,scr,lscr) 3c 4c $Id$ 5c 6c ----- Wrapper routine that sets the sizes of scratch blocks ----- 7c 8 implicit double precision (a-h,o-z) 9#include "hnd_pointers.fh" 10 dimension scr(lscr) 11 dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 12 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 13 dimension tv10(nint,3),coord(3,nat),zan(nat) 14c 15c Use scratch for temps in integral calculation 16c Scratch needs are 17c 6*(Li+2)*(Lj+3) + 6*(Li+2)*(Lj+1)*((Li+Lj+1)/2+1) 18c 19c The pointers are defined in hdn_pointers.fh 20c 21 call hnd_giaotv101(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,expj, 22 1 coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint,tv10,scr(gh01(1)), 23 2 scr(gh01(12)),scr(gh01(2)),scr(gh01(3)),scr(gh01(4)), 24 3 scr(gh01(5)),scr(gh01(6)),scr(gh01(7)),scr(gh01(8)),scr(gh01(9)), 25 4 scr(gh01(10)),scr(gh01(11))) 26c 27 return 28 end 29c 30 subroutine hnd_giaotv101(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj, 31 1 expj,coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint, 32 2 tv10,xs,ys,zs,xxs,yys,zzs,xv,yv,zv,xxv,yyv,zzv) 33 implicit double precision (a-h,o-z) 34#include "nwc_const.fh" 35#include "hnd_rys.fh" 36#include "hnd_tol.fh" 37#include "stdio.fh" 38 integer NoKinetic ! FA-11-08-10 39 common /skipKinetic/NoKinetic ! FA-11-08-10 40 41 common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj, 42 1 ni,nj,ccx,ccy,ccz 43 dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim) 44 dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen) 45 dimension tv10(nint,3),coord(3,nat),zan(nat) 46 dimension xs(Li+2,Lj+3),ys(Li+2,Lj+3),zs(Li+2,Lj+3) 47 dimension xv(Li+2,Lj+1,*),yv(Li+2,Lj+1,*),zv(Li+2,Lj+1,*) 48 dimension xxs(Li+2,Lj+3),yys(Li+2,Lj+3),zzs(Li+2,Lj+3) 49 dimension xxv(Li+2,Lj+1,*),yyv(Li+2,Lj+1,*),zzv(Li+2,Lj+1,*) 50 data rln10 /2.30258d+00/ 51 data pi212 /1.1283791670955d+00/ 52c 53 tol=rln10*itol 54c 55c Zero integral array 56c 57 call dcopy(nint*3,0.0d0,0,tv10,1) 58c 59c ----- ishell ----- 60c 61 xi=xyzi(1) 62 yi=xyzi(2) 63 zi=xyzi(3) 64 lit = Li + 1 65 maxi = lit*(lit+1)/2 66 litmod=lit+1 67c 68c ----- jshell ----- 69c 70 xj=xyzj(1) 71 yj=xyzj(2) 72 zj=xyzj(3) 73 ljt = Lj + 1 74 maxj = ljt*(ljt+1)/2 75 ljtmod=ljt+2 76c 77 qijx=yi*zj-zi*yj 78 qijy=zi*xj-xi*zj 79 qijz=xi*yj-yi*xj 80 tijx=xi-xj 81 tijy=yi-yj 82 tijz=zi-zj 83c 84 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 85c 86 nroots=(1+lit+ljt-2)/2+1 87 if(nroots.gt.maxrys) then 88 write(luout,9997) maxrys,lit,ljt,nroots 89 call errquit('hnd_giaotv10: need higher Rys rt',nroots,INT_ERR) 90 endif 91c 92c ----- i primitive ----- 93c 94 do ig=1, i_nprim 95 ai=expi(ig) 96 arri=ai*rr 97 axi=ai*xi 98 ayi=ai*yi 99 azi=ai*zi 100 csi=coefi(ig,i_ngen) 101c 102c ----- j primitive ----- 103c 104 do jg=1,j_nprim 105 aj=expj(jg) 106 aa=ai+aj 107 aa1=1.0d0/aa 108 dum=aj*arri*aa1 109 if(dum.gt.tol) goto 1000 110 fac= exp(-dum) 111 csj=coefj(jg,j_ngen) 112 ax=(axi+aj*xj)*aa1 113 ay=(ayi+aj*yj)*aa1 114 az=(azi+aj*zj)*aa1 115c 116c ----- density factor ----- 117c 118 cij=csi*csj*fac 119c 120c ----- kinetic energy ----- 121c 122 if (NoKinetic .eq. 1) goto 20 ! SKIP KINETIC ENERGY-FA-11-08-10 123 t = sqrt(aa1) 124 x0=ax 125 y0=ay 126 z0=az 127 do j=1,ljtmod 128 nj=j 129 do i=1,litmod 130 ni=i 131 call hnd_sxyz 132 xs(i,j)=xint*t 133 ys(i,j)=yint*t 134 zs(i,j)=zint*t 135 enddo 136 do i=1,lit 137 xxs(i,j)=xs(i+1,j) 138 yys(i,j)=ys(i+1,j) 139 zzs(i,j)=zs(i+1,j) 140 enddo 141 enddo 142c 143 do i=1,lit 144 xv(i,1,1)=(xs(i,1)-xs(i,3)*(aj+aj))*aj 145 yv(i,1,1)=(ys(i,1)-ys(i,3)*(aj+aj))*aj 146 zv(i,1,1)=(zs(i,1)-zs(i,3)*(aj+aj))*aj 147 xxv(i,1,1)=(xxs(i,1)-xxs(i,3)*(aj+aj))*aj 148 yyv(i,1,1)=(yys(i,1)-yys(i,3)*(aj+aj))*aj 149 zzv(i,1,1)=(zzs(i,1)-zzs(i,3)*(aj+aj))*aj 150 if (ljt.gt.1) then 151 xv(i,2,1)=(xs(i,2)*dble(2+2-1)-xs(i,4)*(aj+aj))*aj 152 yv(i,2,1)=(ys(i,2)*dble(2+2-1)-ys(i,4)*(aj+aj))*aj 153 zv(i,2,1)=(zs(i,2)*dble(2+2-1)-zs(i,4)*(aj+aj))*aj 154 xxv(i,2,1)=(xxs(i,2)*dble(2+2-1)-xxs(i,4)*(aj+aj))*aj 155 yyv(i,2,1)=(yys(i,2)*dble(2+2-1)-yys(i,4)*(aj+aj))*aj 156 zzv(i,2,1)=(zzs(i,2)*dble(2+2-1)-zzs(i,4)*(aj+aj))*aj 157 endif 158 do j=3,ljt 159 xv(i,j,1)=(xs(i,j)*dble(j+j-1)-xs(i,j+2)*(aj+aj))*aj 160 1 -xs(i,j-2)*dble(((j-1)*(j-2))/2) 161 yv(i,j,1)=(ys(i,j)*dble(j+j-1)-ys(i,j+2)*(aj+aj))*aj 162 1 -ys(i,j-2)*dble(((j-1)*(j-2))/2) 163 zv(i,j,1)=(zs(i,j)*dble(j+j-1)-zs(i,j+2)*(aj+aj))*aj 164 1 -zs(i,j-2)*dble(((j-1)*(j-2))/2) 165 xxv(i,j,1)=(xxs(i,j)*dble(j+j-1)-xxs(i,j+2)*(aj+aj))*aj 166 1 -xxs(i,j-2)*dble(((j-1)*(j-2))/2) 167 yyv(i,j,1)=(yys(i,j)*dble(j+j-1)-yys(i,j+2)*(aj+aj))*aj 168 1 -yys(i,j-2)*dble(((j-1)*(j-2))/2) 169 zzv(i,j,1)=(zzs(i,j)*dble(j+j-1)-zzs(i,j+2)*(aj+aj))*aj 170 1 -zzs(i,j-2)*dble(((j-1)*(j-2))/2) 171 enddo 172 enddo 173c 174 ij=0 175 do j=1,maxj 176 call getNxyz(Lj,j,Nxyz) 177 jx = Nxyz(1) + 1 178 jy = Nxyz(2) + 1 179 jz = Nxyz(3) + 1 180 do i=1,maxi 181 call getNxyz(Li,i,Nxyz) 182 ix = Nxyz(1) + 1 183 iy = Nxyz(2) + 1 184 iz = Nxyz(3) + 1 185 dum = xv(ix,jx,1)* ys(iy,jy) * zs(iz,jz) 186 1 + xs(ix,jx) * yv(iy,jy,1)* zs(iz,jz) 187 2 + xs(ix,jx) * ys(iy,jy) * zv(iz,jz,1) 188 dumx=xxv(ix,jx,1)* ys(iy,jy) * zs(iz,jz) 189 1 +xxs(ix,jx) * yv(iy,jy,1)* zs(iz,jz) 190 2 +xxs(ix,jx) * ys(iy,jy) * zv(iz,jz,1) 191 dumy= xv(ix,jx,1)*yys(iy,jy) * zs(iz,jz) 192 1 + xs(ix,jx) *yyv(iy,jy,1)* zs(iz,jz) 193 2 + xs(ix,jx) *yys(iy,jy) * zv(iz,jz,1) 194 dumz= xv(ix,jx,1)* ys(iy,jy) *zzs(iz,jz) 195 1 + xs(ix,jx) * yv(iy,jy,1)*zzs(iz,jz) 196 2 + xs(ix,jx) * ys(iy,jy) *zzv(iz,jz,1) 197 ij=ij+1 198 tv10(ij,1)=tv10(ij,1)+cij*(qijx*dum+tijy*dumz-tijz*dumy) 199 tv10(ij,2)=tv10(ij,2)+cij*(qijy*dum+tijz*dumx-tijx*dumz) 200 tv10(ij,3)=tv10(ij,3)+cij*(qijz*dum+tijx*dumy-tijy*dumx) 201 enddo 202 enddo 203 20 continue ! SKIP KINETIC ENERGY-FA-11-08-10 204c 205c ----- nuclear attraction ----- 206c 207 aax=aa*ax 208 aay=aa*ay 209 aaz=aa*az 210 do ic=1,nat 211 znuc=-zan(ic) 212 cx=coord(1,ic) 213 cy=coord(2,ic) 214 cz=coord(3,ic) 215 yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 216 call hnd_droot 217 do iroot=1,nroots 218 uu=u9(iroot)*aa 219 ww=w9(iroot)*znuc 220 tt=1.0d0/(aa+uu) 221 t = sqrt(tt) 222 x0=(aax+uu*cx)*tt 223 y0=(aay+uu*cy)*tt 224 z0=(aaz+uu*cz)*tt 225 do j=1,ljt 226 nj=j 227 do i=1,litmod 228 ni=i 229 call hnd_sxyz 230 xv(i,j,iroot)=xint 231 yv(i,j,iroot)=yint 232 zv(i,j,iroot)=zint*ww 233 enddo 234 do i=1,lit 235 xxv(i,j,iroot)=xv(i+1,j,iroot) 236 yyv(i,j,iroot)=yv(i+1,j,iroot) 237 zzv(i,j,iroot)=zv(i+1,j,iroot) 238 enddo 239 enddo 240 enddo 241c 242 ij=0 243 do j=1,maxj 244 call getNxyz(Lj,j,Nxyz) 245 jx = Nxyz(1) + 1 246 jy = Nxyz(2) + 1 247 jz = Nxyz(3) + 1 248 do i=1,maxi 249 call getNxyz(Li,i,Nxyz) 250 ix = Nxyz(1) + 1 251 iy = Nxyz(2) + 1 252 iz = Nxyz(3) + 1 253 dum =0.0d0 254 dumx=0.0d0 255 dumy=0.0d0 256 dumz=0.0d0 257 do iroot=1,nroots 258 dum =dum + xv(ix,jx,iroot)* yv(iy,jy,iroot)* zv(iz,jz,iroot) 259 dumx=dumx+xxv(ix,jx,iroot)* yv(iy,jy,iroot)* zv(iz,jz,iroot) 260 dumy=dumy+ xv(ix,jx,iroot)*yyv(iy,jy,iroot)* zv(iz,jz,iroot) 261 dumz=dumz+ xv(ix,jx,iroot)* yv(iy,jy,iroot)*zzv(iz,jz,iroot) 262 enddo 263 dum =dum *(aa1*pi212) 264 dumx=dumx*(aa1*pi212) 265 dumy=dumy*(aa1*pi212) 266 dumz=dumz*(aa1*pi212) 267 ij=ij+1 268 tv10(ij,1)=tv10(ij,1)+cij*(qijx*dum+tijy*dumz-tijz*dumy) 269 tv10(ij,2)=tv10(ij,2)+cij*(qijy*dum+tijz*dumx-tijx*dumz) 270 tv10(ij,3)=tv10(ij,3)+cij*(qijz*dum+tijx*dumy-tijy*dumx) 271 enddo 272 enddo 273c 274c ----- end loop over atoms ----- 275c 276 enddo 277c 278c ----- end loop over primitives ----- 279c 280 1000 continue 281 enddo 282 enddo 283c 284 return 285 9997 format(' in -gitv10- the rys quadrature is not implemented', 286 1 ' beyond -nroots- = ',i2,/,' lit,ljt,nroots= ',3i3) 287 end 288