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