1c
2c $Id$
3c
4
5*     ******************************************************
6*     *                                                    *
7*     *             nwpw_WGaussian                         *
8*     *                                                    *
9*     ******************************************************
10*
11*     Calculates the two electron two center Gaussian integral
12*
13*                                            //
14*    WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
15*                                            || ------------------------------------  dr dr'
16*                                            //                |r-r'|
17*
18*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
19*
20*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
21*
22*     The normalization constant C_l is defined such at
23*            /
24*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
25*            /
26*
27      real*8 function nwpw_WGaussian(la,ma,sa,lb,mb,sb,Rab)
28      implicit none
29      integer la,ma,lb,mb
30      real*8  sa,sb,Rab(3)
31
32*     *** local variables ***
33      integer l,m,fac,fac2
34      real*8 c,x,y,pi,tmp,mtmp,alpha
35      real*8 cos_theta,phi,R
36
37*     **** external functions ****
38      integer  nwpw_doublefactorial
39      external nwpw_doublefactorial
40      real*8   nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
41      external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
42
43      call nwpw_timing_start(49)
44      pi = 4.0d0*datan(1.0d0)
45      x = dble(nwpw_doublefactorial(2*la+1))
46      y = dble(nwpw_doublefactorial(2*lb+1))
47      alpha = dsqrt(0.25d0*(sa*sa + sb*sb))
48      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
49      cos_theta = Rab(3)/R
50
51      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
52         phi = 0.0d0
53      else
54         phi = datan2(Rab(2),Rab(1))
55      end if
56
57      if (mod(2*la+lb,2).eq.1) then
58         c = -32.0d0*pi/(x*y)
59      else
60         c = 32.0d0*pi/(x*y)
61      end if
62
63      if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
64         fac = -1
65      else
66         fac = 1
67      end if
68
69      tmp = 0.0d0
70      do l = abs(la-lb), (la+lb), 2
71         mtmp = 0.0d0
72         do m=-l,l
73            mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb)
74     >                   *Tesseral_lm(l,m,cos_theta,phi)
75         end do
76         tmp = tmp + fac * mtmp * nwpw_GaussBessel(la+lb,l,alpha,R)
77         fac = -fac
78      end do
79      call nwpw_timing_end(49)
80
81      nwpw_WGaussian = c * tmp
82      return
83      end
84
85
86*     ******************************************************
87*     *                                                    *
88*     *             nwpw_dWGaussian                        *
89*     *                                                    *
90*     ******************************************************
91*
92*     Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab
93*
94*                                            //
95*   dWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
96*                                            || ------------------------------------  dr dr'
97*                                            //                |r-r'|
98*
99*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
100*
101*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
102*
103*     The normalization constant C_l is defined such at
104*            /
105*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
106*            /
107*
108      subroutine nwpw_dWGaussian(la,ma,sa,lb,mb,sb,Rab,W,dW)
109      implicit none
110      integer la,ma,lb,mb
111      real*8  sa,sb,Rab(3)
112      real*8  W,dW(3)
113
114*     *** local variables ***
115      integer l,m,fac
116      real*8 c,x,y,pi,tmp,mtmp,mtmpx,mtmpy,mtmpz,alpha
117      real*8 cos_theta,phi,R,gg1,gg2,gg3,Tx,Ty,Tz
118
119*     **** external functions ****
120      integer  nwpw_doublefactorial
121      external nwpw_doublefactorial
122      real*8   nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm
123      external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm
124
125      call nwpw_timing_start(55)
126      pi = 4.0d0*datan(1.0d0)
127      x = dble(nwpw_doublefactorial(2*la+1))
128      y = dble(nwpw_doublefactorial(2*lb+1))
129      alpha = dsqrt(0.25d0*(sa*sa + sb*sb))
130
131      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
132      cos_theta = Rab(3)/R
133      phi       = datan2(Rab(2),Rab(1))
134
135      if (mod(2*la+lb,2).eq.1) then
136         c = -32.0d0*pi/(x*y)
137      else
138         c = 32.0d0*pi/(x*y)
139      end if
140
141      if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
142         fac = -1
143      else
144         fac = 1
145      end if
146
147      W = 0.0d0
148      dW(1) = 0.0d0
149      dW(2) = 0.0d0
150      dW(3) = 0.0d0
151      do l = abs(la-lb), (la+lb), 2
152         mtmp  = 0.0d0
153         mtmpx = 0.0d0
154         mtmpy = 0.0d0
155         mtmpz = 0.0d0
156         do m=-l,l
157            gg1 = nwpw_gaunt(.false.,l,m,la,ma,lb,mb)
158            call dTesseral_lm(l,m,cos_theta,phi,Tx,Ty,Tz)
159            mtmp = mtmp   + gg1*Tesseral_lm(l,m,cos_theta,phi)
160            mtmpx = mtmpx + gg1*Tx
161            mtmpy = mtmpy + gg1*Ty
162            mtmpz = mtmpz + gg1*Tz
163         end do
164         gg2 = nwpw_GaussBessel(la+lb,l,alpha,R)
165         gg3 = nwpw_dGaussBessel(la+lb,l,alpha,R)
166         W     = W     + fac*(mtmp*gg2)
167         dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3)
168         dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3)
169         dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3)
170         fac = -fac
171      end do
172      W = W*c
173      dW(1) = dW(1)*c
174      dW(2) = dW(2)*c
175      dW(3) = dW(3)*c
176      call nwpw_timing_end(55)
177
178      return
179      end
180
181*     ******************************************************
182*     *                                                    *
183*     *             nwpw_UGaussian                         *
184*     *                                                    *
185*     ******************************************************
186*
187*     Calculates the two electron one center Gaussian integral
188*
189*                                            //
190*    WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r) * g(lb,mb,sb;r')
191*                                            || ------------------------------------  dr dr'
192*                                            //                |r-r'|
193*
194*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
195*
196*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
197*
198*     The normalization constant C_l is defined such at
199*            /
200*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
201*            /
202*
203*
204*   Note - this routine is equivalent to find_self_energy_coeff in the paw code
205*
206      real*8 function nwpw_UGaussian(la,ma,sa,lb,mb,sb)
207      implicit none
208      integer la,ma,lb,mb
209      real*8  sa,sb
210
211*     *** local variables ***
212      real*8 x,y,twopi,s,U
213
214*     **** external functions ****
215      integer  nwpw_doublefactorial
216      external nwpw_doublefactorial
217
218      call nwpw_timing_start(48)
219      U = 0.0d0
220      if ((la.eq.lb).and.(ma.eq.mb)) then
221         twopi = 8.0d0*datan(1.0d0)
222         x = dble((2*la+1)*nwpw_doublefactorial(2*la+1))
223         y = (dsqrt(0.5d0*(sa*sa+sb*sb)))**(2*la+1)
224         U = 4.0d0*dsqrt(twopi)/(x*y)
225      end if
226      call nwpw_timing_end(48)
227
228      nwpw_UGaussian = U
229      return
230      end
231
232
233*************************** real combo versions ***********************************
234
235*     ******************************************************
236*     *                                                    *
237*     *             nwpw_WGaussian3                        *
238*     *                                                    *
239*     ******************************************************
240*
241*     Calculates the 4 terms sum
242*
243*        WGaussian(la,ma,sa,lb,mb,sb,Rab)
244*      - WGaussian(la,ma,sa,lb,mb,sm,Rab)
245*      - WGaussian(la,ma,sm,lb,mb,sb,Rab)
246*      + WGaussian(la,ma,sm,lb,mb,sm,Rab)
247*
248*    of the two electron two center Gaussian integral
249*
250*                                            //
251*    WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
252*                                            || ------------------------------------  dr dr'
253*                                            //                |r-r'|
254*
255*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
256*
257*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
258*
259*     The normalization constant C_l is defined such at
260*            /
261*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
262*            /
263*
264      real*8 function nwpw_WGaussian3(la,ma,sa,lb,mb,sb,sm,Rab)
265      implicit none
266      integer la,ma,lb,mb
267      real*8  sa,sb,sm,Rab(3)
268
269*     *** local variables ***
270      integer l,m,fac,fac2
271      real*8 c,x,y,pi,mtmp
272      real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm
273      real*8 tmp
274      real*8 cos_theta,phi,R,gg2
275
276*     **** external functions ****
277      integer  nwpw_doublefactorial
278      external nwpw_doublefactorial
279      real*8   nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
280      external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
281
282      call nwpw_timing_start(49)
283      pi = 4.0d0*datan(1.0d0)
284      x = dble(nwpw_doublefactorial(2*la+1))
285      y = dble(nwpw_doublefactorial(2*lb+1))
286      alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb))
287      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
288      alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb))
289      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
290
291      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
292      cos_theta = Rab(3)/R
293
294      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
295         phi = 0.0d0
296      else
297         phi = datan2(Rab(2),Rab(1))
298      end if
299
300      if (mod(2*la+lb,2).eq.1) then
301         c = -32.0d0*pi/(x*y)
302      else
303         c = 32.0d0*pi/(x*y)
304      end if
305
306      if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
307         fac = -1
308      else
309         fac = 1
310      end if
311
312      tmp = 0.0d0
313      do l = abs(la-lb), (la+lb), 2
314         mtmp = 0.0d0
315         do m=-l,l
316            mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb)
317     >                   *Tesseral_lm(l,m,cos_theta,phi)
318         end do
319         gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R)
320     >       - nwpw_GaussBessel(la+lb,l,alpha_am,R)
321     >       - nwpw_GaussBessel(la+lb,l,alpha_mb,R)
322     >       + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
323
324         tmp = tmp + fac*mtmp*gg2
325         fac = -fac
326      end do
327      call nwpw_timing_end(49)
328
329      nwpw_WGaussian3 = c*tmp
330      return
331      end
332
333*     ******************************************************
334*     *                                                    *
335*     *             nwpw_WGaussian2                        *
336*     *                                                    *
337*     ******************************************************
338*
339*     Calculates the 4 term sum
340*
341*        WGaussian(la,ma,sa,lb,mb,sa,Rab)
342*      - WGaussian(la,ma,sa,lb,mb,sm,Rab)
343*      - WGaussian(la,ma,sm,lb,mb,sa,Rab)
344*      + WGaussian(la,ma,sm,lb,mb,sm,Rab)
345*
346*    of the two electron two center Gaussian integral
347*
348*                                            //
349*    WGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
350*                                            || ------------------------------------  dr dr'
351*                                            //                |r-r'|
352*
353*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
354*
355*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
356*
357*     The normalization constant C_l is defined such at
358*            /
359*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
360*            /
361*
362      real*8 function nwpw_WGaussian2(la,ma,sa,lb,mb,sm,Rab)
363      implicit none
364      integer la,ma,lb,mb
365      real*8  sa,sb,sm,Rab(3)
366
367*     *** local variables ***
368      integer l,m,fac,fac2
369      real*8 c,x,y,pi,mtmp
370      real*8 alpha_aa,alpha_am,alpha_mm
371      real*8 tmp
372      real*8 cos_theta,phi,R,gg2
373
374*     **** external functions ****
375      integer  nwpw_doublefactorial
376      external nwpw_doublefactorial
377      real*8   nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
378      external nwpw_gaunt,nwpw_GaussBessel,Tesseral_lm
379
380      call nwpw_timing_start(49)
381      pi = 4.0d0*datan(1.0d0)
382      x = dble(nwpw_doublefactorial(2*la+1))
383      y = dble(nwpw_doublefactorial(2*lb+1))
384      alpha_aa = dsqrt(0.25d0*(sa*sa + sa*sa))
385      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
386      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
387
388      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
389      cos_theta = Rab(3)/R
390
391      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
392         phi = 0.0d0
393      else
394         phi = datan2(Rab(2),Rab(1))
395      end if
396
397      if (mod(2*la+lb,2).eq.1) then
398         c = -32.0d0*pi/(x*y)
399      else
400         c = 32.0d0*pi/(x*y)
401      end if
402
403      if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
404         fac = -1
405      else
406         fac = 1
407      end if
408
409      tmp = 0.0d0
410      do l = abs(la-lb), (la+lb), 2
411         mtmp = 0.0d0
412         do m=-l,l
413            mtmp = mtmp + nwpw_gaunt(.false.,l,m,la,ma,lb,mb)
414     >                   *Tesseral_lm(l,m,cos_theta,phi)
415         end do
416
417         gg2 = nwpw_GaussBessel(la+lb,l,alpha_aa,R)
418     >       + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
419     >       - 2.0d0*nwpw_GaussBessel(la+lb,l,alpha_am,R)
420
421         tmp = tmp + fac*mtmp*gg2
422
423         fac = -fac
424      end do
425      call nwpw_timing_end(49)
426
427      nwpw_WGaussian2 = c*tmp
428      return
429      end
430
431*     ******************************************************
432*     *                                                    *
433*     *             nwpw_dWGaussian3                       *
434*     *                                                    *
435*     ******************************************************
436*
437*     Calculates the 4 terms sum
438*
439*        WGaussian(la,ma,sa,lb,mb,sb,Rab)
440*      - WGaussian(la,ma,sa,lb,mb,sm,Rab)
441*      - WGaussian(la,ma,sm,lb,mb,sb,Rab)
442*      + WGaussian(la,ma,sm,lb,mb,sm,Rab)
443*
444*     of the two electron two center Gaussian integral and it's derivative wrt to Rab
445*
446*                                            //
447*   dWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
448*                                            || ------------------------------------  dr dr'
449*                                            //                |r-r'|
450*
451*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Tlm(rhat)
452*
453*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
454*
455*     The normalization constant C_l is defined such at
456*            /
457*            | g(l,m,s;r) * |r|**l *Tlm(rhat) dr = 1
458*            /
459*
460      subroutine nwpw_dWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab,W,dW)
461      implicit none
462      integer la,ma,lb,mb
463      real*8  sa,sb,sm,Rab(3)
464      real*8  W,dW(3)
465
466*     *** local variables ***
467      integer l,m,fac
468      real*8 c,x,y,pi,tmp,mtmp,mtmpx,mtmpy,mtmpz,alpha
469      real*8 cos_theta,phi,R,gg1,gg2,gg3,Tx,Ty,Tz
470      real*8 W_ab,W_am,W_mb,W_mm
471      real*8 dW_ab(3),dW_am(3),dW_mb(3),dW_mm(3)
472      real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm
473
474*     **** external functions ****
475      integer  nwpw_doublefactorial
476      external nwpw_doublefactorial
477      real*8   nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm
478      external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel,Tesseral_lm
479
480      call nwpw_timing_start(55)
481      pi = 4.0d0*datan(1.0d0)
482      x = dble(nwpw_doublefactorial(2*la+1))
483      y = dble(nwpw_doublefactorial(2*lb+1))
484      alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb))
485      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
486      alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb))
487      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
488
489      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
490      cos_theta = Rab(3)/R
491      phi       = datan2(Rab(2),Rab(1))
492
493      if (mod(2*la+lb,2).eq.1) then
494         c = -32.0d0*pi/(x*y)
495      else
496         c = 32.0d0*pi/(x*y)
497      end if
498
499      if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
500         fac = -1
501      else
502         fac = 1
503      end if
504
505      W = 0.0d0
506      dW(1) = 0.0d0
507      dW(2) = 0.0d0
508      dW(3) = 0.0d0
509
510      do l = abs(la-lb), (la+lb), 2
511         mtmp  = 0.0d0
512         mtmpx = 0.0d0
513         mtmpy = 0.0d0
514         mtmpz = 0.0d0
515         do m=-l,l
516            gg1 = nwpw_gaunt(.false.,l,m,la,ma,lb,mb)
517            call dTesseral_lm(l,m,cos_theta,phi,Tx,Ty,Tz)
518            mtmp = mtmp   + gg1*Tesseral_lm(l,m,cos_theta,phi)
519            mtmpx = mtmpx + gg1*Tx
520            mtmpy = mtmpy + gg1*Ty
521            mtmpz = mtmpz + gg1*Tz
522         end do
523
524         gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R)
525     >       - nwpw_GaussBessel(la+lb,l,alpha_am,R)
526     >       - nwpw_GaussBessel(la+lb,l,alpha_mb,R)
527     >       + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
528
529         gg3 = nwpw_dGaussBessel(la+lb,l,alpha_ab,R)
530     >       - nwpw_dGaussBessel(la+lb,l,alpha_am,R)
531     >       - nwpw_dGaussBessel(la+lb,l,alpha_mb,R)
532     >       + nwpw_dGaussBessel(la+lb,l,alpha_mm,R)
533
534         W     = W     + fac*(mtmp*gg2)
535         dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3)
536         dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3)
537         dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3)
538
539         fac = -fac
540      end do
541      W = W*c
542      dW(1) = dW(1)*c
543      dW(2) = dW(2)*c
544      dW(3) = dW(3)*c
545      call nwpw_timing_end(55)
546
547      return
548      end
549
550
551*************************** real combo versions ***********************************
552
553
554
555*************************** complex versions ***********************************
556
557
558*     ******************************************************
559*     *                                                    *
560*     *             nwpw_CWGaussian                        *
561*     *                                                    *
562*     ******************************************************
563*
564*     Calculates the two electron two center Gaussian integral
565*
566*                                             //
567*    CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
568*                                             || ------------------------------------  dr dr'
569*                                             //                |r-r'|
570*
571*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat)
572*
573*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
574*
575*     The normalization constant C_l is defined such at
576*            /
577*            | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1
578*            /
579*
580      complex*16 function nwpw_CWGaussian(la,ma,sa,lb,mb,sb,Rab)
581      implicit none
582      integer la,ma,lb,mb
583      real*8  sa,sb,Rab(3)
584
585*     *** local variables ***
586      integer l,m,fac,fac2
587      real*8 c,x,y,pi,alpha
588      real*8 cos_theta,phi,R
589      complex*16 mtmp,tmp
590
591*     **** external functions ****
592      integer  nwpw_doublefactorial
593      external nwpw_doublefactorial
594      real*8   nwpw_gaunt,nwpw_GaussBessel
595      external nwpw_gaunt,nwpw_GaussBessel
596      complex*16 Yspherical_lm
597      external   Yspherical_lm
598      real*8   nwpw_gaunt_sub,gen_gaunt_coeff_sub
599      external nwpw_gaunt_sub,gen_gaunt_coeff_sub
600
601
602      call nwpw_timing_start(49)
603      pi = 4.0d0*datan(1.0d0)
604      x = dble(nwpw_doublefactorial(2*la+1))
605      y = dble(nwpw_doublefactorial(2*lb+1))
606      alpha = dsqrt(0.25d0*(sa*sa + sb*sb))
607      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
608      cos_theta = Rab(3)/R
609
610      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
611         phi = 0.0d0
612      else
613         phi = datan2(Rab(2),Rab(1))
614      end if
615
616      if (mod(2*la+lb,2).eq.1) then
617         c = -32.0d0*pi/(x*y)
618      else
619         c = 32.0d0*pi/(x*y)
620      end if
621
622      !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1)
623      !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
624      !if (mod(lb+mb,2).eq.1) then
625
626
627      !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
628      if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
629         fac = -1
630      else
631         fac = 1
632      end if
633
634      m = ma + mb
635      tmp = 0.0d0
636      do l = abs(la-lb), (la+lb), 2
637         if (abs(m).le.l) then
638c            write(*,*) "l,m=",l,m,la,ma,lb,mb,
639c     >                 nwpw_gaunt(.true.,l,m,la,ma,lb,-mb),
640c     >                 nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb),
641c     >                 gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb)
642            mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb)
643     >            *YSpherical_lm(l,m,cos_theta,phi)
644            tmp = tmp + fac * mtmp * nwpw_GaussBessel(la+lb,l,alpha,R)
645         end if
646         fac = -fac
647      end do
648      call nwpw_timing_end(49)
649
650      nwpw_CWGaussian = c * tmp
651      return
652      end
653
654
655*     ******************************************************
656*     *                                                    *
657*     *             nwpw_dCWGaussian                       *
658*     *                                                    *
659*     ******************************************************
660*
661*     Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab
662*
663*                                             //
664*   dCWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
665*                                             || ------------------------------------  dr dr'
666*                                             //                |r-r'|
667*
668*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat)
669*
670*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
671*
672*     The normalization constant C_l is defined such at
673*            /
674*            | g(l,m,s;r) * |r|**l * congj(Ylm(rhat)) dr = 1
675*            /
676*
677      subroutine nwpw_dCWGaussian(la,ma,sa,lb,mb,sb,Rab,W,dW)
678      implicit none
679      integer la,ma,lb,mb
680      real*8  sa,sb,Rab(3)
681      complex*16  W,dW(3)
682
683*     *** local variables ***
684      integer l,m,fac
685      real*8 c,x,y,pi,alpha
686      real*8 cos_theta,phi,R,gg1,gg2,gg3
687      complex*16 mtmp,mtmpx,mtmpy,mtmpz,Tx,Ty,Tz
688
689*     **** external functions ****
690      integer  nwpw_doublefactorial
691      external nwpw_doublefactorial
692      real*8   nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel
693      external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel
694      complex*16 Yspherical_lm
695      external   Yspherical_lm
696
697      call nwpw_timing_start(55)
698      pi = 4.0d0*datan(1.0d0)
699      x = dble(nwpw_doublefactorial(2*la+1))
700      y = dble(nwpw_doublefactorial(2*lb+1))
701      alpha = dsqrt(0.25d0*(sa*sa + sb*sb))
702
703      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
704      cos_theta = Rab(3)/R
705      phi       = datan2(Rab(2),Rab(1))
706
707      if (mod(2*la+lb,2).eq.1) then
708         c = -32.0d0*pi/(x*y)
709      else
710         c = 32.0d0*pi/(x*y)
711      end if
712
713      !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
714      !if (mod(lb+mb,2).eq.1) then
715
716      if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
717         fac = -1
718      else
719         fac = 1
720      end if
721
722      W = dcmplx(0.0d0,0.0d0)
723      dW(1) = dcmplx(0.0d0,0.0d0)
724      dW(2) = dcmplx(0.0d0,0.0d0)
725      dW(3) = dcmplx(0.0d0,0.0d0)
726      m = ma + mb
727      do l = abs(la-lb), (la+lb), 2
728         if (abs(m).le.l) then
729            gg1 = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb)
730            call dYspherical_lm(l,m,cos_theta,phi,Tx,Ty,Tz)
731            mtmp  = gg1*Yspherical_lm(l,m,cos_theta,phi)
732            mtmpx = gg1*Tx
733            mtmpy = gg1*Ty
734            mtmpz = gg1*Tz
735            gg2 = nwpw_GaussBessel(la+lb,l,alpha,R)
736            gg3 = nwpw_dGaussBessel(la+lb,l,alpha,R)
737            W     = W     + fac*(mtmp*gg2)
738            dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3)
739            dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3)
740            dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3)
741         end if
742         fac = -fac
743      end do
744      W = W*c
745      dW(1) = dW(1)*c
746      dW(2) = dW(2)*c
747      dW(3) = dW(3)*c
748      call nwpw_timing_end(55)
749
750      return
751      end
752
753
754*************************** complex versions ***********************************
755
756*************************** complex combo versions ***********************************
757
758*     ******************************************************
759*     *                                                    *
760*     *             nwpw_CWGaussian3                       *
761*     *                                                    *
762*     ******************************************************
763*
764*     Calculates the two electron two center Gaussian integral
765*
766*                                             //
767*    CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
768*                                             || ------------------------------------  dr dr'
769*                                             //                |r-r'|
770*
771*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat)
772*
773*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
774*
775*     The normalization constant C_l is defined such at
776*            /
777*            | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1
778*            /
779*
780      complex*16 function nwpw_CWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab)
781      implicit none
782      integer la,ma,lb,mb
783      real*8  sa,sb,sm,Rab(3)
784
785*     *** local variables ***
786      integer l,m,fac,fac2
787      real*8 c,x,y,pi
788      real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm
789      real*8 cos_theta,phi,R,gg2
790      complex*16 mtmp,tmp
791
792*     **** external functions ****
793      integer  nwpw_doublefactorial
794      external nwpw_doublefactorial
795      real*8   nwpw_gaunt,nwpw_GaussBessel
796      external nwpw_gaunt,nwpw_GaussBessel
797      complex*16 Yspherical_lm
798      external   Yspherical_lm
799      real*8   nwpw_gaunt_sub,gen_gaunt_coeff_sub
800      external nwpw_gaunt_sub,gen_gaunt_coeff_sub
801
802
803      call nwpw_timing_start(49)
804      pi = 4.0d0*datan(1.0d0)
805      x = dble(nwpw_doublefactorial(2*la+1))
806      y = dble(nwpw_doublefactorial(2*lb+1))
807      alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb))
808      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
809      alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb))
810      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
811      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
812      cos_theta = Rab(3)/R
813
814      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
815         phi = 0.0d0
816      else
817         phi = datan2(Rab(2),Rab(1))
818      end if
819
820      if (mod(2*la+lb,2).eq.1) then
821         c = -32.0d0*pi/(x*y)
822      else
823         c = 32.0d0*pi/(x*y)
824      end if
825
826      !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1)
827      !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
828      !if (mod(lb+mb,2).eq.1) then
829
830
831      !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
832      if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
833         fac = -1
834      else
835         fac = 1
836      end if
837
838      m = ma + mb
839      tmp = dcmplx(0.0d0,0.0d0)
840      do l = abs(la-lb), (la+lb), 2
841         if (abs(m).le.l) then
842c            write(*,*) "l,m=",l,m,la,ma,lb,mb,
843c     >                 nwpw_gaunt(.true.,l,m,la,ma,lb,-mb),
844c     >                 nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb),
845c     >                 gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb)
846            mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb)
847     >            *YSpherical_lm(l,m,cos_theta,phi)
848
849            gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R)
850     >          - nwpw_GaussBessel(la+lb,l,alpha_am,R)
851     >          - nwpw_GaussBessel(la+lb,l,alpha_mb,R)
852     >          + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
853
854            tmp = tmp + fac*mtmp*gg2
855
856         end if
857         fac = -fac
858      end do
859      call nwpw_timing_end(49)
860
861      nwpw_CWGaussian3 = c*tmp
862      return
863      end
864
865
866*     ******************************************************
867*     *                                                    *
868*     *             nwpw_CWGaussian2                       *
869*     *                                                    *
870*     ******************************************************
871*
872*     Calculates the two electron two center Gaussian integral
873*
874*                                             //
875*    CWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
876*                                             || ------------------------------------  dr dr'
877*                                             //                |r-r'|
878*
879*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat)
880*
881*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
882*
883*     The normalization constant C_l is defined such at
884*            /
885*            | g(l,m,s;r) * |r|**l *conjg(Ylm(rhat)) dr = 1
886*            /
887*
888      complex*16 function nwpw_CWGaussian2(la,ma,sa,lb,mb,sm,Rab)
889      implicit none
890      integer la,ma,lb,mb
891      real*8  sa,sm,Rab(3)
892
893*     *** local variables ***
894      integer l,m,fac,fac2
895      real*8 c,x,y,pi
896      real*8 alpha_aa,alpha_am,alpha_mm
897      real*8 cos_theta,phi,R,gg2
898      complex*16 mtmp,tmp
899
900*     **** external functions ****
901      integer  nwpw_doublefactorial
902      external nwpw_doublefactorial
903      real*8   nwpw_gaunt,nwpw_GaussBessel
904      external nwpw_gaunt,nwpw_GaussBessel
905      complex*16 Yspherical_lm
906      external   Yspherical_lm
907      real*8   nwpw_gaunt_sub,gen_gaunt_coeff_sub
908      external nwpw_gaunt_sub,gen_gaunt_coeff_sub
909
910
911      call nwpw_timing_start(49)
912      pi = 4.0d0*datan(1.0d0)
913      x = dble(nwpw_doublefactorial(2*la+1))
914      y = dble(nwpw_doublefactorial(2*lb+1))
915      alpha_aa = dsqrt(0.25d0*(sa*sa + sa*sa))
916      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
917      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
918      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
919      cos_theta = Rab(3)/R
920
921      if ((dabs(Rab(2)).lt.1.0d-9).and.(dabs(Rab(1)).lt.1.0d-9)) then
922         phi = 0.0d0
923      else
924         phi = datan2(Rab(2),Rab(1))
925      end if
926
927      if (mod(2*la+lb,2).eq.1) then
928         c = -32.0d0*pi/(x*y)
929      else
930         c = 32.0d0*pi/(x*y)
931      end if
932
933      !phase_factor = (-1)**(m1+l1)/sigma**(l1+l2+1)
934      !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
935      !if (mod(lb+mb,2).eq.1) then
936
937
938      !if (mod(lb+mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
939      if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
940         fac = -1
941      else
942         fac = 1
943      end if
944
945      m = ma + mb
946      tmp = dcmplx(0.0d0,0.0d0)
947      do l = abs(la-lb), (la+lb), 2
948         if (abs(m).le.l) then
949c            write(*,*) "l,m=",l,m,la,ma,lb,mb,
950c     >                 nwpw_gaunt(.true.,l,m,la,ma,lb,-mb),
951c     >                 nwpw_gaunt_sub(.true.,l,m,la,ma,lb,-mb),
952c     >                 gen_gaunt_coeff_sub(l,m,la,ma,lb,-mb)
953            mtmp = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb)
954     >            *YSpherical_lm(l,m,cos_theta,phi)
955
956            gg2 = nwpw_GaussBessel(la+lb,l,alpha_aa,R)
957     >          + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
958     >          - 2.0d0*nwpw_GaussBessel(la+lb,l,alpha_am,R)
959
960            tmp = tmp + fac*mtmp*gg2
961
962         end if
963         fac = -fac
964      end do
965      call nwpw_timing_end(49)
966
967      nwpw_CWGaussian2 = c*tmp
968      return
969      end
970
971
972
973*     ******************************************************
974*     *                                                    *
975*     *             nwpw_dCWGaussian3                      *
976*     *                                                    *
977*     ******************************************************
978*
979*     Calculates the two electron two center Gaussian integral and it's derivative wrt to Rab
980*
981*                                             //
982*   dCWGaussian(la,ma,sa,lb,Ra,mb,sb,Rb)   =  || g(la,ma,sa;r-Ra) * g(lb,mb,sb;r'-Rb)
983*                                             || ------------------------------------  dr dr'
984*                                             //                |r-r'|
985*
986*     where g(l,m,s; r) = C_l * |r|**l * exp(-(r/s)**2) * Ylm(rhat)
987*
988*          and C_l = 2**(l+2) / (sqrt(pi) (2*l+1)!! s**(2*l+3) )
989*
990*     The normalization constant C_l is defined such at
991*            /
992*            | g(l,m,s;r) * |r|**l * congj(Ylm(rhat)) dr = 1
993*            /
994*
995      subroutine nwpw_dCWGaussian3(la,ma,sa,lb,mb,sb,sm,Rab,W,dW)
996      implicit none
997      integer la,ma,lb,mb
998      real*8  sa,sb,sm,Rab(3)
999      complex*16  W,dW(3)
1000
1001*     *** local variables ***
1002      integer l,m,fac
1003      real*8 c,x,y,pi
1004      real*8 alpha_ab,alpha_am,alpha_mb,alpha_mm
1005      real*8 cos_theta,phi,R,gg1,gg2,gg3
1006      complex*16 mtmp,mtmpx,mtmpy,mtmpz,Tx,Ty,Tz
1007
1008*     **** external functions ****
1009      integer  nwpw_doublefactorial
1010      external nwpw_doublefactorial
1011      real*8   nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel
1012      external nwpw_gaunt,nwpw_GaussBessel,nwpw_dGaussBessel
1013      complex*16 Yspherical_lm
1014      external   Yspherical_lm
1015
1016      call nwpw_timing_start(55)
1017      pi = 4.0d0*datan(1.0d0)
1018      x = dble(nwpw_doublefactorial(2*la+1))
1019      y = dble(nwpw_doublefactorial(2*lb+1))
1020      alpha_ab = dsqrt(0.25d0*(sa*sa + sb*sb))
1021      alpha_am = dsqrt(0.25d0*(sa*sa + sm*sm))
1022      alpha_mb = dsqrt(0.25d0*(sm*sm + sb*sb))
1023      alpha_mm = dsqrt(0.25d0*(sm*sm + sm*sm))
1024
1025      R = dsqrt(Rab(1)*Rab(1) + Rab(2)*Rab(2) + Rab(3)*Rab(3))
1026      cos_theta = Rab(3)/R
1027      phi       = datan2(Rab(2),Rab(1))
1028
1029      if (mod(2*la+lb,2).eq.1) then
1030         c = -32.0d0*pi/(x*y)
1031      else
1032         c = 32.0d0*pi/(x*y)
1033      end if
1034
1035      !if (mod((abs(la-lb)+la+lb)/2,2).eq.1) then
1036      !if (mod(lb+mb,2).eq.1) then
1037
1038      if (mod(mb+(abs(la-lb)+la+lb)/2,2).eq.1) then
1039         fac = -1
1040      else
1041         fac = 1
1042      end if
1043
1044      W = dcmplx(0.0d0,0.0d0)
1045      dW(1) = dcmplx(0.0d0,0.0d0)
1046      dW(2) = dcmplx(0.0d0,0.0d0)
1047      dW(3) = dcmplx(0.0d0,0.0d0)
1048
1049      m = ma + mb
1050      do l = abs(la-lb), (la+lb), 2
1051         if (abs(m).le.l) then
1052            gg1 = nwpw_gaunt(.true.,l,m,la,ma,lb,-mb)
1053            call dYspherical_lm(l,m,cos_theta,phi,Tx,Ty,Tz)
1054            mtmp  = gg1*Yspherical_lm(l,m,cos_theta,phi)
1055            mtmpx = gg1*Tx
1056            mtmpy = gg1*Ty
1057            mtmpz = gg1*Tz
1058
1059            gg2 = nwpw_GaussBessel(la+lb,l,alpha_ab,R)
1060     >          - nwpw_GaussBessel(la+lb,l,alpha_am,R)
1061     >          - nwpw_GaussBessel(la+lb,l,alpha_mb,R)
1062     >          + nwpw_GaussBessel(la+lb,l,alpha_mm,R)
1063
1064            gg3 = nwpw_dGaussBessel(la+lb,l,alpha_ab,R)
1065     >          - nwpw_dGaussBessel(la+lb,l,alpha_am,R)
1066     >          - nwpw_dGaussBessel(la+lb,l,alpha_mb,R)
1067     >          + nwpw_dGaussBessel(la+lb,l,alpha_mm,R)
1068
1069            W     = W     + fac*(mtmp*gg2)
1070            dW(1) = dW(1) + fac*(mtmpx*gg2 + mtmp*(Rab(1)/R)*gg3)
1071            dW(2) = dW(2) + fac*(mtmpy*gg2 + mtmp*(Rab(2)/R)*gg3)
1072            dW(3) = dW(3) + fac*(mtmpz*gg2 + mtmp*(Rab(3)/R)*gg3)
1073         end if
1074         fac = -fac
1075      end do
1076      W = W*c
1077      dW(1) = dW(1)*c
1078      dW(2) = dW(2)*c
1079      dW(3) = dW(3)*c
1080      call nwpw_timing_end(55)
1081
1082      return
1083      end
1084
1085
1086*************************** complex combo versions ***********************************
1087