1      Subroutine xc_murakn(r,w,nr,range,lreset)
2
3C$Id$
4      implicit none
5
6      double precision r(*) ! grid pts coord [output]
7      double precision w(*) ! grid pts weights [output]
8c
9      double precision rm
10      parameter (rm=3d0)
11c
12      integer nr,i
13      double precision fmn,qi,ri,wi,alpha
14      double precision range ! max extent of basis f
15c
16      logical lreset
17cedo      double precision alphad(36)
18cedo      data alphad/ 5d0, 5d0,
19cedoc Li-N
20cedo     *       7d0,  7d0,  5d0,  5d0,  5d0,  5d0,  5d0,  5d0,
21cedoC Na-Ar
22cedo     *       7d0,  7d0,  5d0,  5d0,  5d0,  5d0,  5d0,  5d0,
23cedoC K-Kr
24cedo     *       7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0,
25cedoC Cu-Kr
26cedo     *       7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0,  7d0/
27c
28c Definition of Euler-Maclaurin numerical quadrature points and weights
29c for radial integrals.
30c Transformation from 0<r<infty to 0<x<1 according to
31c ME Mura and PJ Knowles, J Chem Phys 104, 9848 (1996)
32c
33c***************************************************************************
34      if(lreset) then
35         alpha=-range/log(1d0-((1d0+nr)/(2d0+nr))**3)
36      else
37c
38c     huub recipe
39c
40         alpha=3.3d0*range
41c         alpha=5d0
42c         if(znumber.lt.37) alpha=alphad(znumber)
43      endif
44      fmn = rm/(1d0+nr)
45
46      do 10 i = 1,nr
47
48       qi = dble(i)/(nr+1d0)
49       ri = -alpha*log(1.D0 - qi**rm)
50       wi = fmn*alpha*(ri*ri)/(1.D0 - qi**rm)*qi**(rm-1.d0)
51
52       r(i) = ri
53       w(i) = wi
54
55   10 continue
56
57      end
58