1      Subroutine xc_pw91ldag(rs, i, e ,d1e, d2e)
2      implicit none
3      double precision rs ! [in]
4      integer i ! [in]
5c     i=1 eps_c(rs,0)
6c     i=2 eps_c(rs,1)
7c     i=3 alpha_c(rs)
8      double precision e ! [out]
9      double precision d1e ! [out]
10      double precision d2e ! [out]
11c
12c
13      double precision h1, h2, d1h1,d2h1,d2h2,d1x,d2x,
14     D     d1h2 ,x
15c
16c     Functional Parameters
17c
18      double precision A(3), alp(3), b(4,3)
19      save A, alp, b
20      integer j,n, initial
21      save initial
22      data A   / 0.0310907d0, 0.01554535d0, 0.0168869d0 /
23      data alp / 0.21370d0, 0.20548d0, 0.11125d0 /
24      data b   / 7.5957d0, 3.5876d0, 1.6382d0, 0.49294d0,
25     &          14.1189d0, 6.1977d0, 3.3662d0, 0.62517d0,
26     &          10.357d0, 3.6231d0, 0.88026d0, 0.49671d0 /
27      data initial /1/
28c
29c     Define miscellaneous parameters.
30c
31      if (initial.eq.1)then
32         initial = 0
33c        For convenience store -2A as A and multiply betas by 2A
34         do j = 1, 3
35            A(j) = -2d0*A(j)
36            do n = 1, 4
37               b(n,j) = -A(j)*b(n,j)
38            enddo
39         enddo
40c        Finally, change the sign on A for spin stiffness since
41c        the negative of that is fitted in the PW'91 paper.  We can't
42c        just take the negative of A at the start since A also contributes
43c        to the argument of the ln function.
44         A(3) = -A(3)
45      endif
46      x = sqrt(rs)
47      d1x = 0.5d0/x
48      d2x = -0.5d0*d1x/rs
49c
50      h2 = x*(b(1,i) + x*(b(2,i) + x*(b(3,i) + x*b(4,i))))
51      d1h2 = b(1,i)
52     &     + x*(2d0*b(2,i) + x*(3d0*b(3,i) + 4d0*x*b(4,i)))
53
54      d2h2 = 2d0*b(2,i) + x*(6d0*b(3,i) + 12d0*x*b(4,i))
55
56c
57      h1 = DLOG(1d0+1d0/h2)
58      d1h1 = -d1h2/(h2*(h2+1d0))
59
60      d2h1 = d1h1*d1h1*(2d0*h2+1d0) - d2h2/(h2*(h2+1d0))
61
62c
63      e = A(i)*(1d0+alp(i)*rs)*h1
64
65      d1e = A(i)*(2d0*alp(i)*x*h1+(1d0+alp(i)*rs)*d1h1)
66
67      d2e = A(i)*(2d0*alp(i)*h1+4d0*alp(i)*x*d1h1
68     &     +(1d0+alp(i)*rs)*d2h1)
69
70c
71c           Transform derivatives wrt x to derivatives wrt rs
72c
73
74c           Do 2nd derivative first so the x first derivative in d1e
75c           is not lost
76      d2e = d2e*d1x*d1x + d1e*d2x
77
78      d1e = d1e*d1x
79
80      return
81      end
82c $Id$
83