1 implicit none 2* 3* $Id$ 4* 5 integer lmax, npoly, numl, lmax2 6 double precision util_random 7 external util_random 8 parameter (lmax = 3, npoly = 7) 9 parameter (lmax2 = lmax+lmax) 10 parameter (numl = (lmax+1)*(lmax+2)*(lmax+3)/6) 11 integer i, j, k, l, m, ijk, ind, nlm, n 12 integer i1, i2, j1, j2, k1, k2, l1, l2, l1p, l2p 13 double precision x1, x2, y1, y2, z1, z2 14 double precision test1, test2, test3, r, factor 15 double precision coeff(npoly,numl) 16 double precision scoeff(npoly,numl) 17 double precision ncoeff(npoly,numl) 18 double precision dens(numl*numl),ndens(numl*numl),work(numl*numl) 19 double precision junk, a, b, c, x, y, z, test, testn, err 20 double precision d(((lmax2+1)*(lmax2+2))/2, -lmax2:lmax2, 0:lmax2) 21 double precision 22 $ dinv(((lmax2+1)*(lmax2+2)*(lmax2+3))/6, -lmax2:lmax2,0:lmax2) 23 double precision q(-lmax2:lmax2,0:lmax2) 24c 25 junk = util_random(55512121) 26 call xlm_init 27 call xlm_coeff_inv(lmax2,d,dinv) 28c 29 do l = 0, lmax 30 do m = 1, npoly 31 ijk = (l*(l+1)*(l+2))/6 32 do i = l,0,-1 33 do j = l-i,0,-1 34 k = l-i-j 35 ijk = ijk + 1 36 scoeff(m,ijk) = util_random(0) - 0.5d0 37 coeff(m,ijk) = scoeff(m,ijk) 38 enddo 39 enddo 40 enddo 41 enddo 42c 43 a = util_random(0)-0.5d0 44 b = util_random(0)-0.5d0 45 c = util_random(0)-0.5d0 46 x = util_random(0)-0.5d0 47 y = util_random(0)-0.5d0 48 z = util_random(0)-0.5d0 49 write(6,1) ' New center ', a, b, c 50 write(6,1) ' Target ', x, y, z 51 1 format(a,3f12.6) 52c 53 call cart_poly_translate(lmax,npoly,coeff,ncoeff,a,b,c) 54c 55 err = 0.0d0 56 do m = 1, npoly 57 test = 0.0d0 58 testn= 0.0d0 59 ijk = 0 60 do l = 0, lmax 61 do i = l,0,-1 62 do j = l-i,0,-1 63 k = l-i-j 64 ijk = ijk + 1 65 test = test + scoeff(m,ijk)*x**i*y**j*z**k 66 testn = testn + ncoeff(m,ijk)* 67 $ (x-a)**i*(y-b)**j*(z-c)**k 68 enddo 69 enddo 70 enddo 71 err = max(err,abs(test-testn)) 72** write(6,*) m, test, testn, test-testn 73 enddo 74 write(6,*) ' Err from translation of poly ', err 75c 76 do i = 1, numl**2 77 dens(i) = util_random(0) - 0.5d0 78 enddo 79 x1 = util_random(0)-0.5d0 80 y1 = util_random(0)-0.5d0 81 z1 = util_random(0)-0.5d0 82 x2 = util_random(0)-0.5d0 83 y2 = util_random(0)-0.5d0 84 z2 = util_random(0)-0.5d0 85c 86 write(6,1) ' New center ', a, b, c 87 write(6,1) ' Center1 ', x1, y1, z1 88 write(6,1) ' Center2 ', x2, y2, z2 89 err = 0.0d0 90 do l1 = 0, lmax 91 do l2 = 0, lmax 92c 93 test1 = 0.0d0 94 ind = 0 95 do i2 = l2,0,-1 96 do j2 = l2-i2,0,-1 97 k2 = l2-i2-j2 98 do i1 = l1,0,-1 99 do j1 = l1-i1,0,-1 100 k1 = l1-i1-j1 101 ind = ind + 1 102 test1 = test1 + dens(ind)* 103 $ (x-x1)**i1*(y-y1)**j1*(z-z1)**k1* 104 $ (x-x2)**i2*(y-y2)**j2*(z-z2)**k2 105 enddo 106 enddo 107 enddo 108 enddo 109c 110 call cart_dens_translate(l1,x1,y1,z1,l2,x2,y2,z2, 111 $ a,b,c,dens,ndens,work) 112c 113 test2 = 0.0d0 114 ind = 0 115 do i2 = l2,0,-1 116 do j2 = l2-i2,0,-1 117 k2 = l2-i2-j2 118 do i1 = l1,0,-1 119 do j1 = l1-i1,0,-1 120 k1 = l1-i1-j1 121 ind = ind + 1 122 test2 = test2 + dens(ind)* 123 $ (x-x1)**i1*(y-y1)**j1*(z-z1)**k1* 124 $ (x-x2)**i2*(y-y2)**j2*(z-z2)**k2 125 enddo 126 enddo 127 enddo 128 enddo 129c 130 test3 = 0.0d0 131 ind = 0 132 do l2p = 0, l2 133 do i2 = l2p,0,-1 134 do j2 = l2p-i2,0,-1 135 k2 = l2p-i2-j2 136 do l1p = 0, l1 137 do i1 = l1p,0,-1 138 do j1 = l1p-i1,0,-1 139 k1 = l1p-i1-j1 140 ind = ind + 1 141 test3 = test3 + ndens(ind)* 142 $ (x-a)**i1*(y-b)**j1*(z-c)**k1* 143 $ (x-a)**i2*(y-b)**j2*(z-c)**k2 144 enddo 145 enddo 146 enddo 147 enddo 148 enddo 149 enddo 150c 151 err = max(err,abs(test3-test1)) 152** write(6,*) l1, l2, test1, test3, test3-test1 153c 154 call cart_dens_product(l1,l2,ndens,work) 155 ijk = 0 156 test = 0.0d0 157 do l = 0, l1+l2 158 do i = l,0,-1 159 do j = l-i,0,-1 160 k = l-i-j 161 ijk = ijk + 1 162 test = test + work(ijk)* 163 $ (x-a)**i*(y-b)**j*(z-c)**k 164 enddo 165 enddo 166 enddo 167c 168 err = max(err,abs(test-test1)) 169** write(6,*) test1-test 170c 171 call cart_dens_to_sph(l1+l2,work,ndens,dinv,lmax2) 172 test = 0.0d0 173 r = sqrt((x-a)**2 + (y-b)**2 + (z-c)**2) 174 call xlm(l1+l2,x-a,y-b,z-c,q,lmax2) 175 nlm = 0 176 do n = 0, l1+l2 177 do l = n, 0, -2 178 factor = r**(n-l) 179 do m = -l, l 180 nlm = nlm + 1 181 test = test + factor*q(m,l)*ndens(nlm) 182 enddo 183 enddo 184 enddo 185 err = max(err,abs(test-test1)) 186** write(6,*) test1-test 187c 188 call cart_dens_trans_prod_sph( 189 $ l1, x1, y1, z1, 190 $ l2, x2, y2, z2, 191 $ a, b, c, work, dinv, lmax2, dens, ndens) 192 193 test = 0.0d0 194 r = sqrt((x-a)**2 + (y-b)**2 + (z-c)**2) 195 call xlm(l1+l2,x-a,y-b,z-c,q,lmax2) 196 nlm = 0 197 do n = 0, l1+l2 198 do l = n, 0, -2 199 factor = r**(n-l) 200 do m = -l, l 201 nlm = nlm + 1 202 test = test + factor*q(m,l)*ndens(nlm) 203 enddo 204 enddo 205 enddo 206 err = max(err,abs(test-test1)) 207** write(6,*) test1-test 208 enddo 209 enddo 210 write(6,*) ' Error from translation of density ', err 211c 212 end 213