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