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