1c
2C$Id$
3c
4      subroutine xc_kop(tol_rho,whichf,x,
5     &     kalpha, dkadxa)
6      implicit none
7c
8      character*4 whichf
9      double precision tol_rho,x
10      double precision kalpha,dkadxa
11      if(whichf.eq.'be88') then
12         call xc_kbecke88(tol_rho,x,
13     &     kalpha, dkadxa)
14      endif
15      if(whichf.eq.'pb96') then
16         call xc_kpbe96(tol_rho,x,
17     &     kalpha, dkadxa)
18      endif
19      return
20      end
21      subroutine xc_kbecke88(tol_rho,x,
22     &     kalpha, dkadxa)
23      implicit none
24c
25      double precision tol_rho,x
26      double precision kalpha,dkadxa
27c
28      double precision BETA, C
29      Parameter (BETA = 0.0042D0)
30      double precision g,gdenom,dgdenom,dg
31      double precision arcsinh, darcsinh
32      arcsinh(x)=log(x+dsqrt(1d0+x*x))
33      darcsinh(x)=1d0/dsqrt(1d0+x*x)
34c
35c
36c     Uniform electron gas constant
37c
38      C =  3d0*(0.75d0/acos(-1d0))**(1d0/3d0)
39
40      if (x.gt.0d0)then
41         gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
42         g = 2d0*BETA*x*x / gdenom
43         dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
44         dg = g*(2d0/x-dgdenom/gdenom)
45
46         kalpha= C + g
47         dkadxa = dg
48
49      else
50         kalpha= C
51         dkadxa = 0d0
52      endif
53      return
54      end
55      subroutine xc_kpbe96(tol_rho,x,
56     &     kalpha, dkadxa)
57      implicit none
58c
59      double precision tol_rho,x
60      double precision kalpha,dkadxa
61c
62      double precision pi,um, uk, umk
63      parameter(um=0.21951d0, uk=0.804d0, umk=um/uk)
64      double precision C
65      double precision forty8,deno
66c
67c
68c     Uniform electron gas constant
69c
70      pi = acos(-1.d0)
71      C =  3d0*(0.75d0/pi)**(1d0/3d0)
72
73      if (x.gt.0d0)then
74         forty8=1d0/((48d0*pi*pi)**(2d0/3d0))
75         deno=1d0/(1d0+um*x*x*forty8/uk)
76         kalpha= C * (1d0 + uk - uk *deno)
77         dkadxa = C * (2d0*um*x*deno*deno*
78     *        forty8)
79
80      else
81         kalpha= C
82         dkadxa = 0d0
83      endif
84      return
85      end
86