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