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