1c 2c $Id$ 3c 4c Modified from the HONDO dipole integral code 5c 6 subroutine hnd_dipder(xyzi, expi, coefi,i_nprim,i_ngen, Li, 7 * xyzj, expj, coefj, j_nprim, j_ngen, Lj, 8 * center, buf, lbuf, nint, scr, lscr) 9c 10c This routine calculates the dipole derivative integrals. This is a 11c wrapper to the work routine hnd_ddipint. 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 lscr ! [input] size of scratch array 22 integer lbuf ! [input] size of any integral 23 integer nint ! [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 center(3) ! [input] center for expansion 31 double precision scr(lscr) ! [scratch] scratch buffer 32 double precision buf(lbuf) ! [output] dipole derivative integrals 33c 34 double precision zero 35 data zero /0.0d+00/ 36c 37# include "hnd_pointers.fh" 38c 39c Zero out the integral 40c 41 call dcopy(nint*18,zero,0,buf,1) 42c 43c Use scratch for temps in integral calculation 44c Scratch needs are 6*(Li+2)*(Lj+2)+12*(Li+1)*(Lj+1) 45c 46c The pointers are defined in hdn_pointers.fh 47c 48 call hnd_ddipint(xyzi, expi, coefi,i_nprim,i_ngen, Li, 49 * xyzj, expj, coefj, j_nprim, j_ngen, Lj, 50 * center, buf, buf((nint*9)+1),nint, 51 * scr(tvp(1)),scr(tvp(2)),scr(tvp(3)),scr(tvp(4)),scr(tvp(5)), 52 * scr(tvp(6)),scr(tvp(7)),scr(tvp(8)),scr(tvp(9)), 53 * scr(tvp(10)),scr(tvp(11)),scr(tvp(12)),scr(tvp(13)), 54 * scr(tvp(14)),scr(tvp(15)),scr(tvp(16)),scr(tvp(17)), 55 * scr(tvp(18))) 56c 57 return 58 end 59c 60 subroutine hnd_ddipint(xyzi, expi, coefi,i_nprim,i_ngen, Li, 61 * xyzj, expj, coefj, j_nprim, j_ngen, Lj, center, didij, djdij, 62 * nint, xin, yin,zin, xxin, yyin, zzin, dxsdi, dysdi, dzsdi, 63 * dxxdi, dyydi, dzzdi, dxsdj, dysdj, dzsdj, dxxdj, dyydj, dzzdj) 64c 65c This is the main work routine for the dipole derivative integrals. 66c It is assumed the didij and djdij have been zeroed before entering. 67c 68 implicit none 69#include "hnd_tol.fh" 70#include "stdio.fh" 71#include "errquit.fh" 72 integer i_nprim ! [input] num. prims on function i 73 integer i_ngen ! [input] num general conts on func. i 74 integer Li ! [input] angular momentum of func. i 75 integer j_nprim ! [input] num. prims on function j 76 integer j_ngen ! [input] num general conts on func. j 77 integer Lj ! [input] angular momentum of func. j 78 integer nint ! [input] number of base integrals 79 double precision xyzi(3) ! [input] position of center i 80 double precision expi(i_nprim) ! [input] exponents on i 81 double precision coefi(i_nprim,i_ngen) ! [input] i coeffs 82 double precision xyzj(3) ! [input] position of center j 83 double precision expj(j_nprim) ! [input] exponents on j 84 double precision coefj(j_nprim,j_ngen) ! [input] j coeffs 85 double precision center(3) ! [input] center for expansion 86 double precision didij(nint,9) ! [output] dip.der. integrals wrt i 87 double precision djdij(nint,9) ! [ouptut] dip.der. integrals wrt j 88c 89 integer ni, nj 90 integer ijx, ijy, ijz, ix, iy, iz, jx, jy, jz 91 integer nder, lit, ljt, litder, ljtder 92 integer maxi, maxj 93 integer ig, jg, ij, i, j 94 double precision zero, one 95 double precision rln10, tol 96 double precision rr, ai, arri, axi, ayi, azi 97 double precision aj, aa, aa1, dum, fac 98 double precision csi, cpi, cdi, cfi, cgi 99 double precision csj, cpj, cdj, cfj, cgj 100 double precision ax, ay, az, dum1, dum2 101 double precision dij 102 double precision xin, yin,zin, xxin, yyin, zzin 103 double precision dxsdi, dysdi, dzsdi, dxxdi, dyydi, dzzdi 104 double precision dxsdj, dysdj, dzsdj, dxxdj, dyydj, dzzdj 105 double precision dumxdx, dumydx, dumzdx 106 double precision dumxdy, dumydy, dumzdy 107 double precision dumxdz, dumydz, dumzdz 108 double precision xint, yint, zint, xintx, yinty, zintz 109 double precision t, xc, yc, zc, x0, y0, z0 110 double precision xi, yi, zi, xj, yj, zj 111 double precision t1 112c 113 common/hnd_xyzdip/xint,yint,zint,xintx,yinty,zintz,t,xc,yc,zc, 114 1 x0,y0,z0,xi,yi,zi,xj,yj,zj,ni,nj 115 dimension xin(Li+2,Lj+2), yin(Li+2,Lj+2), zin(Li+2,Lj+2) 116 dimension xxin(Li+2,Lj+2), yyin(Li+2,Lj+2), zzin(Li+2,Lj+2) 117 dimension dxsdi(Li+1,Lj+1),dysdi(Li+1,Lj+1),dzsdi(Li+1,Lj+1) 118 dimension dxsdj(Li+1,Lj+1),dysdj(Li+1,Lj+1),dzsdj(Li+1,Lj+1) 119 dimension dxxdi(Li+1,Lj+1),dyydi(Li+1,Lj+1),dzzdi(Li+1,Lj+1) 120 dimension dxxdj(Li+1,Lj+1),dyydj(Li+1,Lj+1),dzzdj(Li+1,Lj+1) 121 dimension Nxyz(3) 122 data zero /0.0d+00/ 123 data one /1.0e+00/ 124 data rln10 /2.30258d+00/ 125c 126 tol=rln10*itol 127 nder = 1 128c 129 xc=center(1) 130 yc=center(2) 131 zc=center(3) 132c 133c ----- ishell ----- 134c 135 xi=xyzi(1) 136 yi=xyzi(2) 137 zi=xyzi(3) 138 lit = Li + 1 139 maxi=lit*(lit+1)/2 140c 141 litder=lit+nder 142c 143c ----- jshell ----- 144c 145 xj=xyzj(1) 146 yj=xyzj(2) 147 zj=xyzj(3) 148 ljt = Lj + 1 149 maxj=ljt*(ljt+1)/2 150c 151 ljtder=ljt+nder 152c 153 rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 154 nroots=(lit+ljt-2)/2+1 155 if(nroots.gt.maxrys) then 156 write(luout,9997) maxrys,lit,ljt,nroots 157 call errquit('hnd_dipder: angular momentum too high',555, 158 & BASIS_ERR) 159 endif 160c 161c ----- i primitive ----- 162c 163 do 7000 ig=1,i_nprim 164 ai=expi(ig) 165 arri=ai*rr 166 axi=ai*xi 167 ayi=ai*yi 168 azi=ai*zi 169 csi=coefi(ig,i_ngen) 170c 171c ----- j primitive ----- 172c 173 do 6000 jg=1,j_nprim 174 aj=expj(jg) 175 aa=ai+aj 176 aa1=one/aa 177 dum=aj*arri*aa1 178 if(dum.gt.tol) go to 6000 ! the integral is zero 179 fac= exp(-dum) 180 csj=coefj(jg,j_ngen) 181 ax=(axi+aj*xj)*aa1 182 ay=(ayi+aj*yj)*aa1 183 az=(azi+aj*zj)*aa1 184c 185c ----- density factor ----- 186c 187 dij=fac*csi*csj 188c 189c ----- dipole moment integrals ----- 190c 191 t = sqrt(aa) 192 t1=one/t 193 x0=ax 194 y0=ay 195 z0=az 196 do 370 j=1,ljtder 197 nj=j 198 do 370 i=1,litder 199 ni=i 200 call hnd_dipxyz 201 xin(i,j)=xint*t1 202 yin(i,j)=yint*t1 203 zin(i,j)=zint*t1 204 xxin(i,j)=xintx*t1 205 yyin(i,j)=yinty*t1 206 zzin(i,j)=zintz*t1 207 370 continue 208c 209 call hnd_deriaj(dxsdi,dysdi,dzsdi,dxsdj,dysdj,dzsdj, 210 * xin,yin,zin, 211 * lit,ljt,ai,aj) 212 call hnd_deriaj(dxxdi,dyydi,dzzdi,dxxdj,dyydj,dzzdj, 213 * xxin,yyin,zzin, 214 * lit,ljt,ai,aj) 215 ij=0 216 do 390 i=1,maxi 217 call getNxyz(Li,i,Nxyz) 218 ix = Nxyz(1) + 1 219 iy = Nxyz(2) + 1 220 iz = Nxyz(3) + 1 221 do 390 j=1,maxj 222 call getNxyz(Lj,j,Nxyz) 223 jx = Nxyz(1) + 1 224 jy = Nxyz(2) + 1 225 jz = Nxyz(3) + 1 226 ij=ij+1 227c 228c First do derivative wrt the first atom 229c 230 dumxdx = dxxdi(ix,jx)* yin(iy,jy)* zin(iz,jz) 231 dumydx = dxsdi(ix,jx)* yyin(iy,jy)* zin(iz,jz) 232 dumzdx = dxsdi(ix,jx)* yin(iy,jy)* zzin(iz,jz) 233 dumxdy = xxin(ix,jx)*dysdi(iy,jy)* zin(iz,jz) 234 dumydy = xin(ix,jx)*dyydi(iy,jy)* zin(iz,jz) 235 dumzdy = xin(ix,jx)*dysdi(iy,jy)* zzin(iz,jz) 236 dumxdz = xxin(ix,jx)* yin(iy,jy)*dzsdi(iz,jz) 237 dumydz = xin(ix,jx)* yyin(iy,jy)*dzsdi(iz,jz) 238 dumzdz = xin(ix,jx)* yin(iy,jy)*dzzdi(iz,jz) 239c 240 didij(ij,1) = didij(ij,1) + dij*dumxdx 241 didij(ij,2) = didij(ij,2) + dij*dumydx 242 didij(ij,3) = didij(ij,3) + dij*dumzdx 243 didij(ij,4) = didij(ij,4) + dij*dumxdy 244 didij(ij,5) = didij(ij,5) + dij*dumydy 245 didij(ij,6) = didij(ij,6) + dij*dumzdy 246 didij(ij,7) = didij(ij,7) + dij*dumxdz 247 didij(ij,8) = didij(ij,8) + dij*dumydz 248 didij(ij,9) = didij(ij,9) + dij*dumzdz 249c 250c Next do derivative wrt the second atom 251c 252 dumxdx = dxxdj(ix,jx)* yin(iy,jy)* zin(iz,jz) 253 dumydx = dxsdj(ix,jx)* yyin(iy,jy)* zin(iz,jz) 254 dumzdx = dxsdj(ix,jx)* yin(iy,jy)* zzin(iz,jz) 255 dumxdy = xxin(ix,jx)*dysdj(iy,jy)* zin(iz,jz) 256 dumydy = xin(ix,jx)*dyydj(iy,jy)* zin(iz,jz) 257 dumzdy = xin(ix,jx)*dysdj(iy,jy)* zzin(iz,jz) 258 dumxdz = xxin(ix,jx)* yin(iy,jy)*dzsdj(iz,jz) 259 dumydz = xin(ix,jx)* yyin(iy,jy)*dzsdj(iz,jz) 260 dumzdz = xin(ix,jx)* yin(iy,jy)*dzzdj(iz,jz) 261c 262 djdij(ij,1) = djdij(ij,1) + dij*dumxdx 263 djdij(ij,2) = djdij(ij,2) + dij*dumydx 264 djdij(ij,3) = djdij(ij,3) + dij*dumzdx 265 djdij(ij,4) = djdij(ij,4) + dij*dumxdy 266 djdij(ij,5) = djdij(ij,5) + dij*dumydy 267 djdij(ij,6) = djdij(ij,6) + dij*dumzdy 268 djdij(ij,7) = djdij(ij,7) + dij*dumxdz 269 djdij(ij,8) = djdij(ij,8) + dij*dumydz 270 djdij(ij,9) = djdij(ij,9) + dij*dumzdz 271 390 continue 272c 273 6000 continue 274 7000 continue 275 return 276c 277 9997 format(' in -hnd_stvint- the rys quadrature is not implemented', 278 1 ' beyond -nroots- = ',i2,/, 279 2 ' lit,ljt,nroots= ',3i3) 280 end 281