1c 2c $Id$ 3c 4c Modified from HONDO 5c 6 subroutine hnd_hlf_ij2(xyzi,expi,coefi,i_nprim,i_ngen, Li, 7 1 xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij, 8 2 zan, xyz, nder, nint, scr,lscr) 9c 10c This is a wrapper routine, setting up scratch blocks used in actual 11c integral routine 12c 13 implicit none 14 integer i_nprim ! [input] num. prims on function i 15 integer i_ngen ! [input] num general conts on func. i 16 integer Li ! [input] angular momentum of func. i 17 integer j_nprim ! [input] num. prims on function j 18 integer j_ngen ! [input] num general conts on func. j 19 integer Lj ! [input] angular momentum of func. j 20 integer nder ! [input] 1=1rst der; 2=2nd der 21 integer nint ! [input] number of base integrals 22 integer kat ! [input] lexical number of an atom 23 integer lscr ! [input] length of scratch space 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,*) ! [input] all atom positions 31 double precision dvij(nint,3) ! [output] 1rst. der integrals 32 double precision ddvij(nint,3,3) ! [output] 2nd der integrals 33 double precision zan(*) ! [input] nuclear charges 34 double precision scr(lscr) ! [input] scratch buffer 35c 36# include "hnd_pointers.fh" 37c 38c Use scratch for temps in integral calculation 39c Scratch needs are dependent on nder: 40c nder=1 (6 * (Li+1)*(Lj+1)*((Li+Lj+2)/2+1) + 3 0.0d0 length blocks 41c nder=2 (9 * (Li+1)*(Lj+1)*((Li+Lj+2)/2+1) 42c 43c The pointers are defined in hnd_pointers.fh 44c 45 call hnd_hlf_ij21(xyzi,expi,coefi,i_nprim,i_ngen, Li, 46 2 xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij, 47 3 zan, xyz, nder, nint, 48 4 scr(tvp(22)),scr(tvp(23)),scr(tvp(24)), 49 5 scr(tvp(25)),scr(tvp(26)),scr(tvp(27)), 50 6 scr(tvp(46)),scr(tvp(47)),scr(tvp(48))) 51c 52 return 53 end 54c 55 subroutine hnd_hlf_ij21(xyzi,expi,coefi,i_nprim,i_ngen, Li, 56 1 xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij, 57 2 zan, xyz, nder, nint,xv,yv,zv,dxv,dyv,dzv,ddxv,ddyv,ddzv) 58c 59c This is the routine that actually computes the 1rst and 2nd 60c derivatives for the Helman-Feynman term. It is assumed that the buffers 61c for the integrals have been 0.0d0ed before entering this routine. 62c 63 implicit none 64#include "stdio.fh" 65#include "hnd_rys.fh" 66#include "hnd_tol.fh" 67#include "errquit.fh" 68 integer i_nprim ! [input] num. prims on function i 69 integer i_ngen ! [input] num general conts on func. i 70 integer Li ! [input] angular momentum of func. i 71 integer j_nprim ! [input] num. prims on function j 72 integer j_ngen ! [input] num general conts on func. j 73 integer Lj ! [input] angular momentum of func. j 74 integer nder ! [input] 1=1rst der; 2=2nd der 75 integer nint ! [input] number of base integrals 76 integer kat ! [input] lexical number of an atom 77 double precision xyzi(3) ! [input] position of center i 78 double precision expi(i_nprim) ! [input] exponents on i 79 double precision coefi(i_nprim,i_ngen) ! [input] i coeffs 80 double precision xyzj(3) ! [input] position of center j 81 double precision expj(j_nprim) ! [input] exponents on j 82 double precision coefj(j_nprim,j_ngen) ! [input] j coeffs 83 double precision xyz(3,*) ! [input] all atom positions 84 double precision dvij(nint,3) ! [output] 1rst. der integrals 85 double precision ddvij(nint,3,3) ! [output] 2nd der integrals 86 double precision zan(*) ! [input] nuclear charges 87c 88 character*8 errmsg 89 common/hnd_xyzder/xint,yint,zint,tx,x0,y0,z0,xi,yi,zi,xj,yj,zj, 90 1 ni,nj,cx,cy,cz 91 double precision xint, yint, zint, tx, x0, y0, z0, xi, yi, zi 92 double precision xj, yj, zj, cx, cy, cz 93 double precision ijx, ijy, ijz 94 double precision rln10, tol, rr, ai, aj, arri 95 double precision axi, ayi, azi, csi, cpi, cdi, cfi, cgi 96 double precision aa, aa1, dum, fac, csj, cpj, cdj, cfj, cgj 97 double precision ax, ay, az, dum1, dum2, pij 98 double precision dumx, dumy, dumz, dumxx, dumyy, dumzz 99 double precision dumxy, dumxz, dumyz 100 double precision pi212, aax, aay, aaz, znuc 101 double precision uu, u2, u4, ww, w2, w4, tt, xv, yv, zv 102 double precision dxv, dyv, dzv, ddxv, ddyv, ddzv 103 integer ni, nj, iroot, Nxyz 104 integer lit,maxi, ljt, maxj 105 integer ig, jg, ij, i, j, ix, iy, iz, jx, jy, jz 106 dimension xv(Li+1,Lj+1,*), yv(Li+1,Lj+1,*), zv(Li+1,Lj+1,*) 107 dimension dxv(Li+1,Lj+1,*), dyv(Li+1,Lj+1,*), dzv(Li+1,Lj+1,*) 108 dimension ddxv(Li+1,Lj+1,*),ddyv(Li+1,Lj+1,*),ddzv(Li+1,Lj+1,*) 109 dimension Nxyz(3) 110 dimension w2(maxrys),w4(maxrys) 111 dimension errmsg(3) 112 data errmsg /'program ','stop in ','-hlfspd-'/ 113 data rln10 /2.30258d+00/ 114 data pi212 /1.1283791670955d+00/ 115c 116 tol =rln10*itol 117c 118c ----- calculate -helfey- term ----- 119c 120c ----- ishell ----- 121c 122 xi=xyzi(1) 123 yi=xyzi(2) 124 zi=xyzi(3) 125 lit = Li + 1 126 maxi=lit*(lit+1)/2 127c 128c ----- jshell ----- 129c 130 xj=xyzj(1) 131 yj=xyzj(2) 132 zj=xyzj(3) 133 ljt = Lj + 1 134 maxj=ljt*(ljt+1)/2 135c 136 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 137 nroots=(lit+ljt+nder-2)/2+1 138 if(nroots.gt.maxrys) then 139 write(luout,9997) maxrys,lit,ljt,nroots 140 call hnd_hnderr(3,errmsg) 141 endif 142c 143c ----- i primitive ----- 144c 145 do 7000 ig=1,i_nprim 146 ai=expi(ig) 147 arri=ai*rr 148 axi=ai*xi 149 ayi=ai*yi 150 azi=ai*zi 151 csi=coefi(ig,i_ngen) 152c 153c ----- j primitive ----- 154c 155 do 6000 jg=1,j_nprim 156 aj=expj(jg) 157 aa=ai+aj 158 aa1=1.0d0/aa 159 dum=aj*arri*aa1 160 if(dum.gt.tol) go to 6000 161 fac= exp(-dum) 162 csj=coefj(jg,j_ngen) 163 ax=(axi+aj*xj)*aa1 164 ay=(ayi+aj*yj)*aa1 165 az=(azi+aj*zj)*aa1 166c 167c ----- density factor ----- 168c 169 pij=fac*csi*csj*pi212*aa1 170c 171c ----- hellmann-feynman term ----- 172c 173 aax=aa*ax 174 aay=aa*ay 175 aaz=aa*az 176c 177c ----- kat ----- 178c 179 znuc=-zan(kat) 180 cx=xyz(1,kat) 181 cy=xyz(2,kat) 182 cz=xyz(3,kat) 183 yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2) 184 call hnd_droot 185 do 420 iroot=1,nroots 186 uu=u9(iroot)*aa 187 u2=uu 188 u4=uu*uu 189 ww=w9(iroot)*znuc 190 w2(iroot)=ww*u2*2.0d0 191 w4(iroot)=ww*u4*4.0d0 192 tt=1.0d0/(aa+uu) 193 tx= sqrt(tt) 194 x0=(aax+uu*cx)*tt 195 y0=(aay+uu*cy)*tt 196 z0=(aaz+uu*cz)*tt 197 do 410 j=1,ljt 198 nj=j 199 do 410 i=1,lit 200 ni=i 201 202 call hnd_sxyz 203 204 xv(i,j,iroot)=xint 205 yv(i,j,iroot)=yint 206 zv(i,j,iroot)=zint 207 208 call hnd_dervxyz(1) 209 210 dxv(i,j,iroot)=xint 211 dyv(i,j,iroot)=yint 212 dzv(i,j,iroot)=zint 213 214 if (nder .gt. 1) then 215 call hnd_dervxyz(2) 216 217 ddxv(i,j,iroot)=xint 218 ddyv(i,j,iroot)=yint 219 ddzv(i,j,iroot)=zint 220 endif 221 410 continue 222 420 continue 223c 224 ij=0 225 do 450 i=1,maxi 226 call getNxyz(Li,i,Nxyz) 227 ix = Nxyz(1) + 1 228 iy = Nxyz(2) + 1 229 iz = Nxyz(3) + 1 230 do 440 j=1,maxj 231 call getNxyz(Lj,j,Nxyz) 232 jx = Nxyz(1) + 1 233 jy = Nxyz(2) + 1 234 jz = Nxyz(3) + 1 235 ij=ij+1 236 if (nder.eq.1) then 237 dumx=0.0d0 238 dumy=0.0d0 239 dumz=0.0d0 240 do 430 iroot=1,nroots 241 dumx = dumx + dxv(ix,jx,iroot)* yv(iy,jy,iroot)* 242 1 zv(iz,jz,iroot)*w2(iroot) 243 dumy = dumy + xv(ix,jx,iroot)* dyv(iy,jy,iroot)* 244 1 zv(iz,jz,iroot)*w2(iroot) 245 dumz = dumz + xv(ix,jx,iroot)* yv(iy,jy,iroot)* 246 1 dzv(iz,jz,iroot)*w2(iroot) 247 430 continue 248 dumx =dumx *pij 249 dumy =dumy *pij 250 dumz =dumz *pij 251 dvij(ij,1)=dvij(ij,1)+dumx 252 dvij(ij,2)=dvij(ij,2)+dumy 253 dvij(ij,3)=dvij(ij,3)+dumz 254 elseif (nder.eq.2) then 255 dumxx=0.0d0 256 dumyy=0.0d0 257 dumzz=0.0d0 258 dumxy=0.0d0 259 dumxz=0.0d0 260 dumyz=0.0d0 261 do 431 iroot=1,nroots 262 dum=xv(ix,jx,iroot)*yv(iy,jy,iroot)*zv(iz,jz,iroot)*w2(iroot) 263 dumxx=dumxx-dum+ ddxv(ix,jx,iroot)* yv(iy,jy,iroot)* 264 1 zv(iz,jz,iroot)*w4(iroot) 265 dumyy=dumyy-dum+ xv(ix,jx,iroot)*ddyv(iy,jy,iroot)* 266 1 zv(iz,jz,iroot)*w4(iroot) 267 dumzz=dumzz-dum+ xv(ix,jx,iroot)* yv(iy,jy,iroot)* 268 1 ddzv(iz,jz,iroot)*w4(iroot) 269 dumxy=dumxy+ dxv(ix,jx,iroot)* dyv(iy,jy,iroot)* 270 1 zv(iz,jz,iroot)*w4(iroot) 271 dumxz=dumxz+ dxv(ix,jx,iroot)* yv(iy,jy,iroot)* 272 1 dzv(iz,jz,iroot)*w4(iroot) 273 dumyz=dumyz+ xv(ix,jx,iroot)* dyv(iy,jy,iroot)* 274 1 dzv(iz,jz,iroot)*w4(iroot) 275 431 continue 276 dumxx=dumxx*pij 277 dumyy=dumyy*pij 278 dumzz=dumzz*pij 279 dumxy=dumxy*pij 280 dumxz=dumxz*pij 281 dumyz=dumyz*pij 282 ddvij(ij,1,1)=ddvij(ij,1,1)+dumxx 283 ddvij(ij,1,2)=ddvij(ij,1,2)+dumxy 284 ddvij(ij,1,3)=ddvij(ij,1,3)+dumxz 285 ddvij(ij,2,1)=ddvij(ij,2,1)+dumxy 286 ddvij(ij,2,2)=ddvij(ij,2,2)+dumyy 287 ddvij(ij,2,3)=ddvij(ij,2,3)+dumyz 288 ddvij(ij,3,1)=ddvij(ij,3,1)+dumxz 289 ddvij(ij,3,2)=ddvij(ij,3,2)+dumyz 290 ddvij(ij,3,3)=ddvij(ij,3,3)+dumzz 291 endif 292c 293 440 continue 294 450 continue 295c 296 6000 continue 297 7000 continue 298c 299 return 300 9997 format(' in -hlf- , the rys quadrature is not implemented', 301 1 ' beyond -nroots- = ',i3,/, 302 2 ' lit,ljt,nroots = ',3i3) 303 end 304