1      subroutine xc_tarad(znumber,r,w,nr,origin,which,rlast)
2
3C$Id$
4      implicit none
5
6      double precision r(*) ! grid pts coord [output]
7      double precision w(*) ! grid pts weights [output]
8      double precision origin ! origin of the grid [input]
9      integer znumber !   atomic number [input]
10c
11      integer nr,i
12      double precision xi,wi
13      double precision ri,zetad(36),zeta,rlast
14      character*8 which
15      double precision alpha,rsc(400),wsc(400)
16      parameter (alpha=0.6d0)
17      data zetad/ 0.8d0, 0.9d0,
18c Li-N
19     *     1.8d0,1.4d0,1.3d0,1.1d0,0.9d0,0.9d0,0.9d0,0.9d0,
20C Na-Ar
21     *     1.4d0,1.3d0,1.3d0,1.2d0,1.1d0,1.0d0,1.0d0,1.0d0,
22C K-Kr
23     *     1.5d0,1.4d0,1.3d0,1.2d0,1.2d0,1.2d0,1.2d0,1.2d0,1.2d0,1.1d0,
24C Cu-Kr
25     *     1.1d0,1.1d0,1.1d0,1.0d0,0.9d0,0.9d0,0.9d0,0.9d0/
26c
27c
28c Definition of Gauss-Chebyshev numerical quadrature points and weights
29c for radial integrals.
30c Transformation from 0<r<infty to -1<x<1 according to
31c O Treutler and R Alrhichs, J Chem Phys 102, 346 (1995)
32c
33c***************************************************************************
34      zeta=1.d0
35      if(znumber.lt.37) zeta = zetad(znumber)
36
37      if(which.eq.'chebyshe') then
38      do  i = 1,nr
39        call grid_gausscheb(wi,xi,i,nr)
40        ri = zeta/log(2.D0)*(1.d0+xi)**alpha*log(2d0/(1d0-xi))
41        r(i) = ri+origin
42        w(i)=wi*r(i)*r(i)*ri*(alpha/(1d0+xi)+
43     +       1d0/(1d0-xi)/log(2d0/(1d0-xi)))
44c
45c       becke transf
46c        ri=(1d0+xi)/(1d0-xi)
47c        w(i) = wi*ri*ri/(1d0-xi)**2*2d0
48      enddo
49      elseif(which.eq.'chebyshr') then
50c
51c     get zeta from rmax
52c
53         call grid_gausscheb(wi,xi,nr,nr)
54c M4
55c         zeta=rlast*log(2d0)/
56c     /        ((1.d0+xi)**alpha*log(2d0/(1d0-xi)))
57c        ri = zeta/log(2.D0)*(1.d0+xi)**alpha*log(2d0/(1d0-xi))
58c M3
59         zeta=(rlast-origin)*log(2d0)/
60     /        (log(2d0/(1d0-xi)))
61        ri = zeta/log(2.D0)*log(2d0/(1d0-xi))
62        r(nr) = ri+origin
63        w(nr)=wi*r(nr)*r(nr)*ri*(alpha/(1d0+xi)+
64     +       1d0/(1d0-xi)/log(2d0/(1d0-xi)))
65
66      do  i = 1,nr-1
67         call grid_gausscheb(wi,xi,i,nr)
68        ri = zeta/log(2.D0)*(1.d0+xi)**alpha*log(2d0/(1d0-xi))
69        r(i) = ri+origin
70        w(i)=wi*r(i)*r(i)*ri*(alpha/(1d0+xi)+
71     +       1d0/(1d0-xi)/log(2d0/(1d0-xi)))
72c
73c       becke transf
74c        ri=(1d0+xi)/(1d0-xi)
75c        w(i) = wi*ri*ri/(1d0-xi)**2*2d0
76      enddo
77      elseif(which.eq.'legendre') then
78        call grid_gaussleg(nr, rsc,
79     &       wsc)
80        do i=1,nr
81          xi=rsc(nr-i+1)
82          ri = zeta/log(2.D0)*(1.d0+xi)**alpha*log(2d0/(1d0-xi))
83          r(i) = ri+origin
84          r(i) = ri+origin
85          wi = wsc(nr-i+1)
86          w(i)=wi*r(i)*r(i)*ri*(alpha/(1d0+xi)+
87     +         1d0/(1d0-xi)/log(2d0/(1d0-xi)))
88        enddo
89      endif
90      return
91      end
92      subroutine xc_interv(rmin,rmax,r,w,nr,which)
93      implicit none
94#include "errquit.fh"
95c
96c     compute Gauss-Chebyshev quadrature for the
97c     interval rmin<r<rmax
98c
99      double precision rmin ! interval minimum [input]
100      double precision rmax ! interval maximum [input]
101      double precision r(*) ! grid pts coord [output]
102      double precision w(*) ! grid pts weights [output]
103      integer nr            ! [input]
104      character*8 which
105c
106      double precision slope,intercept
107      double precision wi,xi,ri,wsc(300),rsc(300)
108      integer i
109c
110      if(nr.gt.300) call errquit(
111     (    ' xc_interv: too many radial pts ',nr, UNKNOWN_ERR)
112      slope=(rmax-rmin)*.5d0
113      intercept=(rmin+rmax)*.5d0
114      if(which.eq.'chebyshe') then
115        do i=1,nr
116          call grid_gausscheb(wi,xi,i,nr)
117          ri = slope*xi + intercept
118          r(i) = ri
119          w(i) = wi*ri*ri*slope
120        enddo
121      elseif(which.eq.'legendre') then
122        call grid_gaussleg(nr, rsc,
123     &       wsc)
124        do i=1,nr
125          ri = slope*rsc(nr-i+1) + intercept
126          r(i) = ri
127          w(i) = wsc(nr-i+1)*ri*ri*slope
128        enddo
129          else
130        call errquit('unrecognized quadr. ',999, INPUT_ERR)
131      endif
132      return
133      end
134
135      subroutine grid_gausscheb(wi,xi,i_in,n)
136      implicit none
137      double precision wi,xi ! [output]
138      integer i_in,n ! [input]
139c
140c     Gauss Chebyshev quadrature scheme
141c     JM Perez-Yorda and E San-Fabian, CPC 77, 46 (1993)
142c
143      double precision pi
144      integer i
145c
146      pi=acos(-1d0)
147      i=n-i_in+1
148      wi=sin(pi*i/(1d0+n))**4/(3d0*(1d0+n))*16d0
149      xi=(1.d0+n-i-i)/(1d0+n)+
150     +     2d0/pi*(1d0+2d0/3d0*sin(pi*i/(1d0+n))**2)*
151     *     cos(pi*i/(1d0+n))*sin(pi*i/(1d0+n))
152      return
153      end
154