1 subroutine cart_dens_trans_prod_sph( 2 $ l1, x1, y1, z1, 3 $ l2, x2, y2, z2, 4 $ a, b, c, work, dinv, ld, dens, ndens) 5* 6* $Id$ 7* 8 implicit none 9 integer l1, l2 ! The angular momenta 10 integer ld ! Dimension parameter for dinv 11 double precision x1, y1, z1, x2, y2, z2, a, b, c 12 double precision dinv(((ld+1)*(ld+2)*(ld+3))/6,-ld:ld,0:ld) 13 double precision dens(((l1+1)*(l1+2))/2,((l2+1)*(l2+2))/2) 14 double precision ndens(*), work(*) 15c 16c Combine the function of cart_dens_translate/product/to_sph. 17c 18c Given in dens a density multiplying cartesian shells of given 19c angular momentum and centers, 20c 21c d(x,y,z) = sum(i,j,k) sum(p,q,r) dens(i,j,k; p,q,r) * 22c . (x-z1)^i(y-y1)^j(z-z1)^k * 23c . (x-z2)^p(y-y2)^q(z-z2)^r * 24c 25c (where i+j+k = l1 and p+q+r = l2) return the density in the form 26c 27c d(x,y,z) = sum(n=0..l1+l2) sum(l=n,0,-2) sum(m=-l,l) 28c . ndens(n,l,m) Xlm(x-a,y-b,z-c) r^(n-l) 29c 30c (where r is the distance of (x,y,z) from (a,b,c) 31c 32c The output order in ndens is as described in cart_dens_to_sph 33c 34c Both work and ndens need to be dimensioned 35c (((lmax+1)*(lmax+2)*(lmax+3))/6)**2 where lmax=max(l1,l2). 36c The output data in ndens is of length 37c (((l12+1)*(l12+2)*(l12+3))/6) where l12 = l1 + l2. 38c 39c Dinv are the cartesian to spherical transformation coeffs 40c generated by calling xlm_coeff_inv, and ld>=l1+l2. 41c 42c The general case ... code special cases for optimization later. 43c 44c This routine ONLY translates the polynomial or angular components 45c ... any rescaling due to the gaussian factors must be done 46c separately. 47c 48c The input density in dens is preserved. 49c 50 call cart_dens_translate(l1, x1, y1, z1, l2, x2, y2, z2, 51 $ a, b, c, dens, ndens, work) 52 call cart_dens_product(l1,l2,ndens,work) 53 call cart_dens_to_sph(l1+l2,work,ndens,dinv,ld) 54c 55 end 56 subroutine cart_dens_translate( 57 $ l1, x1, y1, z1, 58 $ l2, x2, y2, z2, 59 $ a, b, c, dens, ndens, work) 60 implicit none 61 integer l1, l2 62 double precision x1, y1, z1, x2, y2, z2, a, b, c 63 double precision dens(((l1+1)*(l1+2))/2,((l2+1)*(l2+2))/2) 64 double precision ndens(*), work(*) 65c 66c Given a block of a density matrix for two cartesian shells 67c of angular momentum l1 and l2, respectively, translate 68c the center of expansion for the polynomials from the separate 69c centers of the shells to the common center (a, b, c). 70c 71c Work should be the same size as ndens which is 72c ndens((l1+1)*(l1+2)*(l1+3))/6,((l2+1)*(l2+2)*(l2+3))/6) 73c 74c The input density d(ijk,pqr) is defined just for 75c i+j+k=l1 and p+q+r=l2 in the standard order for a cartesian 76c shell. The output density of necessity will include all 77c functions with i+j+k<=l1 and p+q+r<=l2. See cart_poly_translate 78c for the ordering and addressing of these functions. 79c 80c The input density is preserved. 81c 82 integer num1, num2, i1, i2, loff1, loff2 83c 84 integer itri, ind, loff, i, j, l 85c 86 itri(i,j) = ((i*(i-1))/2 + j) 87 ind(i,j,l) = itri(l-i+1,l-i-j+1) ! Index of x^i*y^j*z^(l-i-j) in 88c . ! cartesian shell of rank l 89 loff(l) = ((l*(l+1)*(l+2))/6) 90c 91c 92c Translate shell 1. First transpose into work(2,1), then 93c translate 1. 94c 95 num1 = ((l1+1)*(l1+2))/2 96 num2 = ((l2+1)*(l2+2))/2 97c 98 if (l1 .gt. 0) then 99 call dfill(loff(l1+1)*num2,0.0d0,work,1) 100 loff1= loff(l1)*num2 101 do i1 = 1, num1 102 do i2 = 1, num2 103 work(i2 + loff1) = dens(i1,i2) 104 enddo 105 loff1 = loff1 + num2 106 enddo 107 call cart_poly_translate(l1,num2,work,ndens,a-x1,b-y1,c-z1) 108c 109c Transpose back to work(1,2) and translate shell 2 110c 111 num1 = loff(l1+1) 112 num2 = ((l2+1)*(l2+2))/2 113c 114 call dfill(loff(l1+1)*loff(l2+1),0.0d0,work,1) 115 loff2= loff(l2)*num1 116 do i2 = 1, num2 117 loff1 = 0 118 do i1 = 1, num1 119 work(i1 + loff2) = ndens(i2 + loff1) 120 loff1 = loff1 + num2 121 enddo 122 loff2 = loff2 + num1 123 enddo 124 else 125 call dfill(loff(l1+1)*loff(l2+1),0.0d0,work,1) 126 loff2= loff(l2)*num1 127 do i2 = 1, num2 128 do i1 = 1, num1 129 work(i1 + loff2) = dens(i1, i2) 130 enddo 131 loff2 = loff2 + num1 132 enddo 133 endif 134c 135 if (l2 .gt. 0) then 136 call cart_poly_translate(l2,num1,work,ndens,a-x2,b-y2,c-z2) 137 else 138 call dcopy(loff(l1+1)*loff(l2+1),work,1,ndens,1) 139 endif 140c 141 end 142 subroutine cart_poly_translate(lmax,npoly,coeff,ncoeff,a,b,c) 143 implicit none 144 integer lmax, npoly 145 double precision coeff(npoly,*) 146 double precision ncoeff(npoly,*) 147 double precision a, b, c 148c 149c Given a set of polynomials (m=1,..,npoly) defined as 150c 151c Pm = sum(i,j,k) coeff(m,i,j,k) x^i y^j z^k 152c 153c recompute the coefficients for the polynomial being 154c expanded about the center (a, b, c) 155c 156c Use the standard binomial expansion 157c 158c (x-a+a)^i = (x'+a)^i = a^i + i*a^(i-1)*x' + ... + x'^i 159c 160c coeff of a^p*x^(i-p) is i!/((i-p)!p!) 161c 162c ncoeff returns the result, the input coeffs are DESTROYED 163c 164c You should have at all points Pm(x,y,z) = Pmnew(x-a,y-b,z-c) 165c 166c Included in coeff are all terms x^i y^j z^k with 0 <= i+j+k <= lmax 167c with all terms for each value of l stored contiguously. Within 168c a given l value (i+j+k=l) the offset of x^i y^j z^k is given 169c by ind(i,j,l). The number of preceeding coefficients for all 170c polynomials up to, but not including, l is loff(l) = l*(l+1)*(l+2)/6. 171c 172c So the location of (i,j,k) is loff(l)+ind(i,j,l) where l=i+j+k. 173c 174c i j k 175c 0 0 0 176c 1 0 0, 0 1 0, 0 0 1, 177c 2 0 0, 1 1 0, 1 0 1, 0 2 0, 0 1 1, 0 0 2 178c ... 179c 180 integer i, j, k, l, m, p, nijk, ijk, ijkp 181 double precision fac 182 integer itri, ind, loff 183c 184 itri(i,j) = ((i*(i-1))/2 + j) 185 ind(i,j,l) = itri(l-i+1,l-i-j+1) ! Index of x^i*y^j*z^(l-i-j) in 186c . ! cartesian shell of rank l 187 loff(l) = ((l*(l+1)*(l+2))/6) 188c 189 nijk = loff(lmax+1) 190c 191c First translate x 192c 193 call dfill(npoly*nijk, 0.0d0, ncoeff, 1) 194 do l = 0, lmax 195 do i = l, 0, -1 196 do j = l-i,0,-1 197 k = l-i-j 198 ijk = loff(l) + ind(i,j,l) 199 fac = 1.0d0 200 do p = 0, i 201 if (abs(fac) .gt. 0.0d0) then 202 ijkp = loff(l-p) + ind(i-p,j,l-p) 203 do m = 1, npoly 204 ncoeff(m,ijkp) = ncoeff(m,ijkp) + 205 $ fac*coeff(m,ijk) 206 enddo 207 fac = fac * a * dble(i-p) / dble(p+1) 208 endif 209 enddo 210 enddo 211 enddo 212 enddo 213c 214c Then y 215c 216 call dfill(npoly*nijk, 0.0d0, coeff, 1) 217c 218 do l = 0, lmax 219 do i = l, 0, -1 220 do j = l-i,0,-1 221 k = l-i-j 222 ijk = loff(l) + ind(i,j,l) 223 fac = 1.0d0 224 do p = 0, j 225 if (abs(fac) .gt. 0.0d0) then 226 ijkp = loff(l-p) + ind(i,j-p,l-p) 227 do m = 1, npoly 228 coeff(m,ijkp) = coeff(m,ijkp) + 229 $ fac*ncoeff(m,ijk) 230 enddo 231 fac = fac * b * dble(j-p) / dble(p+1) 232 endif 233 enddo 234 enddo 235 enddo 236 enddo 237c 238c Then z 239c 240 call dfill(npoly*nijk, 0.0d0, ncoeff, 1) 241 do l = 0, lmax 242 do i = l, 0, -1 243 do j = l-i,0,-1 244 k = l-i-j 245 ijk = loff(l) + ind(i,j,l) 246 fac = 1.0d0 247 do p = 0, k 248 if (abs(fac) .gt. 0.0d0) then 249 ijkp = loff(l-p) + ind(i,j,l-p) 250 do m = 1, npoly 251 ncoeff(m,ijkp) = ncoeff(m,ijkp) + 252 $ fac*coeff(m,ijk) 253 enddo 254 fac = fac * c * dble(k-p) / dble(p+1) 255 endif 256 enddo 257 enddo 258 enddo 259 enddo 260c 261 end 262 subroutine cart_dens_product(l1,l2,dens,ndens) 263 implicit none 264 integer l1, l2 265 double precision dens(*), ndens(*) 266c 267c Given in dens(ijk,pqr) the coefficients multiplying 268c x^(i+p)y^(j+q)z^(k+r) for i+j+k <= l1 and p+q+r <= l2 269c in the ordering generated by cart_{dens,poly}_translate. 270c Resum everything so that we return in ndens(ijk) the equivalent 271c coefficients for x^i*y^j*z^k for (i+j+k) <= (l1+l2) again 272c in the order described in cart_poly_translate 273c 274c The arrays should be dimensioned 275c 276c dens(loff(l1+1),loff(l2+1)) 277c ndens(loff(l1+l2+1)) 278c 279 integer i, j, l, l1p, l2p, i1, i2, j1, j2, k1, k2, ipt, ipt12 280 integer i12, j12, l12, loff12 281 integer itri, ind, loff 282c 283 itri(i,j) = ((i*(i-1))/2 + j) 284 ind(i,j,l) = itri(l-i+1,l-i-j+1) ! Index of x^i*y^j*z^(l-i-j) in 285c . ! cartesian shell of rank l 286 loff(l) = (((l)*(l+1)*(l+2))/6) 287c 288 call dfill(loff(l1+l2+1),0.0d0,ndens,1) 289c 290 ipt = 0 291 do l2p = 0, l2 292 do i2 = l2p,0,-1 293 do j2 = l2p-i2,0,-1 294 k2 = l2p-i2-j2 295 do l1p = 0, l1 296 l12 = l1p + l2p 297 loff12 = loff(l12) 298 do i1 = l1p,0,-1 299 i12 = i1 + i2 300 do j1 = l1p-i1,0,-1 301 k1 = l1p-i1-j1 302 ipt = ipt + 1 303c 304 j12 = j1 + j2 305 ipt12 = loff12 + ind(i12,j12,l12) 306 ndens(ipt12) = ndens(ipt12) + dens(ipt) 307c 308 enddo 309 enddo 310 enddo 311 enddo 312 enddo 313 enddo 314c 315 end 316 subroutine cart_dens_to_sph(lmax,dens,ndens,dinv,ld) 317 implicit none 318 integer lmax, ld 319 double precision dens(*),ndens(*) 320 double precision dinv(((ld+1)*(ld+2)*(ld+3))/6,-ld:ld,0:ld) 321#include "errquit.fh" 322c 323c Given in dens a set of density matrix coefficients 324c for a cartesian basis of the form resulting from 325c cart_dens_product, transform these into the real solid 326c spherical harmonic basis returning them in ndens. 327c 328c dinv should contain the transformation coeffs resulting 329c from call xlm_coeff_inv(ld, d, dinv) (ld >= lmax) 330c 331c The input order of coeffs for x^i y^j z^k in dens is 332c 333c do l = 0, lmax 334c . do i = l,0,-1 335c . do j = l-i,0,-1 336c . k = l-i-j 337c . d(i,j,k) 338c 339c the output order of coeffs for r^(n-l)Xlm in ndens is 340c 341c do n = 0, lmax 342c . do l = n,0,-2 343c . do m = -l,l 344c . d(n,l,m) 345c 346c Both dens and ndens are the same size 347c 348c (((lmax+1)*(lmax+2)*(lmax+3))/6) 349c 350 integer ind, n, nbase, loff, l, m, ijk 351 double precision sum 352 loff(l) = ((l*(l+1)*(l+2))/6) 353c 354 if (ld .lt. lmax) call errquit( 355 $ 'carttrans: ld lt lmax ',ld+lmax*10000, CALC_ERR) 356 ind = 0 357 do n = 0, lmax 358 nbase = loff(n) 359 do l = n,0,-2 360 do m = -l, l 361 sum = 0.0d0 362 do ijk = 1,(n+1)*(n+2)/2 363 sum = sum + dens(ijk+nbase)*dinv(ijk+nbase,m,l) 364 enddo 365 ind = ind + 1 366 ndens(ind) = sum 367 enddo 368 enddo 369 if (ind .ne. loff(n+1)) stop 9997 370 enddo 371c 372 end 373 subroutine gaussian_product( 374 $ zeta1, x1, y1, z1, 375 $ zeta2, x2, y2, z2, 376 $ zeta, a, b, c, prefac) 377 implicit none 378 double precision 379 $ zeta1, x1, y1, z1, 380 $ zeta2, x2, y2, z2, 381 $ zeta, a, b, c, prefac 382c 383c The Gaussian product theorem 384c 385c exp(-zeta1*(r-r1)^2)*exp(-zeta2*(r-r2)^2) = 386c . exp(-q*(r1-r2)^2)*exp(-zeta*(r-P)^2) 387c 388c where 389c 390c q = (zeta1*zeta2)/(zeta1+zeta2) 391c zeta = (zeta1+zeta2) 392c P = (zeta1*r1+ zeta2*r2)/(zeta1+zeta2) 393c 394c Return zeta, P=(a,b,c), and prefac=exp(-q*(r1-r2)^2) 395c 396 double precision rzeta, q, dx, dy, dz, r12sq 397c 398 zeta = (zeta1+zeta2) 399 rzeta = 1.0d0/zeta 400 q = (zeta1*zeta2)*rzeta 401 a = (zeta1*x1 + zeta2*x2)*rzeta 402 b = (zeta1*y1 + zeta2*y2)*rzeta 403 c = (zeta1*z1 + zeta2*z2)*rzeta 404 dx = (x1-x2) 405 dy = (y1-y2) 406 dz = (z1-z2) 407 r12sq = dx*dx + dy*dy + dz*dz 408 prefac= exp(-q*r12sq) 409c 410 end 411