1
2*    ************************************
3*    *					*
4*    *	    gen_BEEF_BW_unrestricted	*
5*    *					*
6*    ************************************
7*
8*    This function returns the BEEF exchange-correlation
9*  energy density, xce, and its derivatives with respect
10*  to nup, ndn, |grad nup|, |grad ndn|, and |grad n|.
11*
12*   Entry - n2ft3d     : number of grid points
13*           dn_in(*,2) : spin densites nup and ndn
14*           agr_in(*,3): |grad nup|, |grad ndn|, and |grad n|
15*           x_parameter: scale parameter for exchange
16*           c_parameter: scale parameter for correlation
17*
18*   Exit - xce(*)  : PBE96 energy density
19*        - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn
20*        - fdn(*,3): d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|
21*                    d(n*xce)/d|grad n|
22
23      subroutine gen_BEEF_BW_unrestricted(n2ft3d,
24     >                           dn_in,agr_in,
25     >                           x_parameter,c_parameter,alphac,
26     >                           xce,fn,fdn)
27      implicit none
28
29      integer n2ft3d
30      real*8 dn_in(n2ft3d,2)
31      real*8 agr_in(n2ft3d,3)
32      real*8 x_parameter,c_parameter,alphac
33      real*8 xce(n2ft3d)
34      real*8 fn(n2ft3d,2)
35      real*8 fdn(n2ft3d,3)
36
37*     **** Density cutoff parameter ****
38      real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
39      parameter (DNS_CUT = 1.0d-20)
40      parameter (ETA=1.0d-20)
41      parameter (ETA2=1.0d-14)
42      parameter (alpha_zeta=(1.0d0-ETA2))
43      parameter (alpha_zeta2=(1.0d0-ETA2))
44
45c     ***** PBE96 GGA exchange constants ******
46      real*8 MU,KAPPA
47      parameter (MU    = 0.2195149727645171d0)
48      parameter (KAPPA = 0.8040000000000000d0)
49
50c     ****** PBE96 GGA correlation constants ******
51      real*8 GAMMA,BETA,BOG
52      parameter (GAMMA	= 0.031090690869655d0)
53      parameter (BETA	= 0.066724550603149d0)
54      !parameter (BETA	= 0.066725d0)
55      parameter (BOG    = BETA/GAMMA)
56
57
58c     ****** Perdew-Wang92 LDA correlation coefficients *******
59      real*8 GAM,iGAM,FZZ,iFZZ
60      parameter (GAM  	= 0.519842099789746329d0)
61      parameter (iGAM  	= 1.0d0/GAM)
62      parameter (FZZ    = (8.0d0/(9.0d0*GAM)) )
63      parameter (iFZZ    = 0.125d0*9.0d0*GAM)
64
65      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
66      parameter (A_1  = 0.0310907d0)
67      !parameter (A_1  = 0.031091d0)
68      parameter (A1_1 =	0.2137000d0)
69      parameter (B1_1 =	7.5957000d0)
70      parameter (B2_1 =	3.5876000d0)
71      parameter (B3_1 =	1.6382000d0)
72      parameter (B4_1 =	0.4929400d0)
73
74      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
75      parameter (A_2  =  0.01554535d0)
76      !parameter (A_2  =  0.015545d0)
77      parameter (A1_2 =	 0.20548000d0)
78      parameter (B1_2 =	14.11890000d0)
79      parameter (B2_2 =	 6.19770000d0)
80      parameter (B3_2 =	 3.36620000d0)
81      parameter (B4_2 =	 0.62517000d0)
82
83      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
84      parameter (A_3  =  0.0168869d0)
85      !parameter (A_3  =  0.016887d0)
86      parameter (A1_3 =	 0.1112500d0)
87      parameter (B1_3 =	10.3570000d0)
88      parameter (B2_3 =	 3.6231000d0)
89      parameter (B3_3 =	 0.8802600d0)
90      parameter (B4_3 =	 0.4967100d0)
91
92c     **** other constants ****
93      real*8 onethird,fourthird,fivethird,onesixthm
94      real*8 twothird,sevensixthm
95      real*8 onethirdm
96      parameter (onethird=1.0d0/3.0d0)
97      parameter (onethirdm=-1.0d0/3.0d0)
98      parameter (twothird=2.0d0/3.0d0)
99      parameter (fourthird=4.0d0/3.0d0)
100      parameter (fivethird=5.0d0/3.0d0)
101      parameter (onesixthm=-1.0d0/6.0d0)
102      parameter (sevensixthm=-7.0d0/6.0d0)
103
104*---- Beef expansion parameters given by Wellendorff et al-----------------*
105      real*8 malphac
106c      real*8 alphac,malphac
107c      parameter (alphac=0.6001664769d0)
108c      parameter (malphac=1.0d0-alphac)
109      real*8 am(0:29)
110      data am(0:29)/1.516501714d0,  4.413532099d-1,-9.182135241d-2,
111     >            -2.352754331d-2, 3.418828455d-2, 2.411870076d-3,
112     >            -1.416381352d-2, 6.975895581d-4, 9.859205137d-3,
113     >            -6.737855051d-3,-1.573330824d-3, 5.036146253d-3,
114     >            -2.569472453d-3,-9.874953976d-4, 2.033722895d-3,
115     >            -8.018718848d-4,-6.688078723d-4, 1.030936331d-3,
116     >            -3.673838660d-4,-4.213635394d-4,5.761607992d-4,
117     >            -8.346503735d-5,-4.458447585d-4,4.601290092d-4,
118     >            -5.231775398d-6,-4.239570471d-4,3.750190679d-4,
119     >             2.114938125d-5, -1.904911565d-4,7.384362421d-5/
120
121c     **** local variables ****
122      integer i,j
123      real*8 n,agr
124      real*8 nup,agrup
125      real*8 ndn,agrdn
126      real*8 kf,ks,s,n_onethird,pi,rs_scale
127      real*8 rs         ! Wigner radius
128      real*8 rss        ! rss  = sqrt(rs)
129      real*8 rs_n       ! rs_n = n*drs/dn
130      real*8 t,t2,t4,t6
131      real*8 t_nup      ! t_nup = n*dt/dnup
132      real*8 t_ndn      ! t_ndn = n*dt/dndn
133      real*8 t_agr      ! t_agr = n*dt/dagr
134      real*8 zet,twoksg
135      real*8 zet_nup    ! zet_nup = n*dzet/dnup
136      real*8 zet_ndn    ! zet_nup = n*dzet/dnup
137      real*8 zetp_1_3,zetm_1_3
138      real*8 zetpm_1_3,zetmm_1_3
139      real*8 phi,phi3,phi4
140      real*8 phi_zet
141      real*8 A,A2
142      real*8 A_phi,A_ec_lda
143      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8
144      real*8 PON,FZ,z4
145      real*8 tau
146      real*8 F
147      real*8 Fs                ! dF/ds
148      real*8 Hpbe
149      real*8 Hpbe_t            ! dHpbe/dt
150      real*8 Hpbe_phi          ! dHpbe/dphi
151      real*8 Hpbe_ec_lda       ! dHpbe/d(ec_lda)
152      real*8 Hpbe_nup,Hpbe_ndn ! n*dHpbe/dnup, n*dHpbe/dndn
153      real*8 Ipbe
154      real*8 Ipbe_t,Ipbe_A     ! dIpbe/dt, dIpbe/dA
155
156      real*8 exup,exdn,ex,ex_lda
157      real*8 ecu,ecp,eca,ec,ec_lda
158      real*8 ecu_rs,ecp_rs,eca_rs
159      real*8 ec_lda_rs,ec_lda_zet  ! d(ec_lda)/drs, d(ec_lda)/dzet
160      real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn
161      real*8 fnxup,fdnxup          ! d(n*ex)/dnup, d(n*ex)/dndn
162      real*8 fnxdn,fdnxdn          ! d(n*ex)/d|grad nup|, d(n*ex)/d|grad ndn|
163      real*8 fncup,fncdn           ! d(n*ec)/dnup, d(n*ec)/dndn
164      real*8 fdnx_const
165
166      real*8 p0,p1,p2,dp,s2,oneovers2,Ft,dp2,sgn
167
168      malphac=1.0d0-alphac
169
170      pi = 4.0d0*datan(1.0d0)
171      rs_scale = (0.75d0/pi)**onethird
172      fdnx_const = -3.0d0/(8.0d0*pi)
173
174
175!$OMP DO
176      do i=1,n2ft3d
177         nup     = dn_in(i,1)+ETA
178         agrup   = agr_in(i,1)
179
180         ndn     = dn_in(i,2)+ETA
181         agrdn   = agr_in(i,2)
182
183c        ****************************************************************
184c        ***** calculate polarized Exchange energies and potentials *****
185c        ****************************************************************
186
187c        ************
188c        **** up ****
189c        ************
190         n     = 2.0d0*nup
191         agr   = 2.0d0*agrup
192
193         n_onethird = (3.0d0*n/pi)**onethird
194         ex_lda     = -0.75d0*n_onethird
195
196         kf = (3.0d0*pi*pi*n)**onethird
197         s  = agr/(2.0d0*kf*n)
198
199         t  = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0
200         if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t)
201         F  = 0.0d0
202         Ft = 0.0d0
203
204         s2   = t*t - 1.0d0
205         if (dabs(s2).lt.1.0d-12) then
206            if (t.gt.0.0d0) then
207               do j=0,29
208                  F  = F  + am(j)
209                  Ft = Ft + am(j)*0.5d0*dble(j*(j+1))
210               end do
211            else
212               sgn = 1.0d0
213               do j=0,29
214                  F  = F  + sgn*am(j)
215                  Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1))
216                  sgn = -sgn
217               end do
218            end if
219         else
220            oneovers2 = 1.0d0/s2
221            p0 = 1.0d0
222            dp = 0
223            F  = F  + am(0)*p0
224            Ft = Ft + am(0)*dp
225            p1 = t
226            dp = 1.0d0
227            F  = F  + am(1)*t
228            Ft = Ft + am(1)*dp
229            do j=2,29
230               p2    = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0)
231               dp2   = dble(j)*oneovers2*(t*p2-p1)
232               F  = F + am(j)*p2
233               Ft = Ft + am(j)*dp2
234               p0 = p1
235               p1 = p2
236            end do
237         end if
238         Fs  = (16.0d0*s/(4.0d0+s*s)**2)*Ft
239
240         exup = ex_lda*F
241         fnxup = fourthird*(exup - ex_lda*Fs*s)
242         fdnxup = fdnx_const*Fs
243
244c        **************
245c        **** down ****
246c        **************
247         n     = 2.0d0*ndn
248         agr   = 2.0d0*agrdn
249
250         n_onethird = (3.0d0*n/pi)**onethird
251         ex_lda     = -0.75d0*n_onethird
252
253         kf = (3.0d0*pi*pi*n)**onethird
254         s  = agr/(2.0d0*kf*n)
255
256         t  = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0
257         if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t)
258         F  = 0.0d0
259         Ft = 0.0d0
260
261         s2   = t*t - 1.0d0
262         if (dabs(s2).lt.1.0d-12) then
263            if (t.gt.0.0d0) then
264               do j=0,29
265                  F  = F  + am(j)
266                  Ft = Ft + am(j)*0.5d0*dble(j*(j+1))
267               end do
268            else
269               sgn = 1.0d0
270               do j=0,29
271                  F  = F  + sgn*am(j)
272                  Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1))
273                  sgn = -sgn
274               end do
275            end if
276         else
277            oneovers2 = 1.0d0/s2
278            p0 = 1.0d0
279            dp = 0
280            F  = F  + am(0)*p0
281            Ft = Ft + am(0)*dp
282            p1 = t
283            dp = 1.0d0
284            F  = F  + am(1)*t
285            Ft = Ft + am(1)*dp
286            do j=2,29
287               p2    = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0)
288               dp2   = dble(j)*oneovers2*(t*p2-p1)
289               F  = F + am(j)*p2
290               Ft = Ft + am(j)*dp2
291               p0 = p1
292               p1 = p2
293            end do
294         end if
295         Fs  = (16.0d0*s/(4.0d0+s*s)**2)*Ft
296
297         exdn   = ex_lda*F
298         fnxdn  = fourthird*(exdn - ex_lda*Fs*s)
299         fdnxdn = fdnx_const*Fs
300
301         n = nup+ndn
302
303         ex = (exup*nup+ exdn*ndn)/ n
304
305c        *******************************************************************
306c        ***** calculate polarized correlation energies and potentials *****
307c        *******************************************************************
308         agr   = agr_in(i,3)
309
310         zet = (nup-ndn)/n
311c         if (zet.gt.0.0d0) zet = zet - ETA2
312c         if (zet.lt.0.0d0) zet = zet + ETA2
313c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n**2
314c        if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2
315c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n
316c        zet_nup =  2*ndn/n
317c        zet_ndn = -2*nup/n
318         zet_nup = -(zet - 1.0d0)
319         zet_ndn = -(zet + 1.0d0)
320         zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm
321         zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm
322         zetp_1_3  = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2
323         zetm_1_3  = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2
324
325
326         phi = 0.5d0*( zetp_1_3**2 + zetm_1_3**2)
327         phi_zet = alpha_zeta*( zetpm_1_3 - zetmm_1_3)/3.0d0
328         F =(  (1.0d0+zet*alpha_zeta)*zetp_1_3
329     >       + (1.0d0-zet*alpha_zeta)*zetm_1_3
330     >       - 2.0d0)*iGAM
331
332         FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM)
333
334
335
336*        **** calculate Wigner radius ****
337         rs    = rs_scale/(n**onethird)
338         rss   = dsqrt(rs)
339
340*        **** calculate n*drs/dn ****
341c        rs_n = onethirdm*rs/n
342         rs_n = onethirdm*rs
343
344
345
346c        **** calculate t ****
347         kf = (3.0d0*pi*pi*n)**onethird
348         ks = dsqrt(4.0d0*kf/pi)
349
350         twoksg = 2.0d0*ks*phi
351
352         t  = agr/(twoksg*n)
353
354*        *** calculate n*dt/dnup, n*dt/dndn, n*dt/d|grad n| ****
355         t_nup = sevensixthm*t - (phi_zet)*(zet_nup)*t/phi
356         t_ndn = sevensixthm*t - (phi_zet)*(zet_ndn)*t/phi
357         t_agr  = 1.0d0/(twoksg)
358
359
360
361
362c        **************************************************
363c        ***** compute LSDA correlation energy density ****
364c        **************************************************
365         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs)
366         call LSDT(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs)
367         call LSDT(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs)
368
369         z4 = zet**4
370
371         ec_lda = ecu*(1.0d0-F*z4)
372     >          + ecp*F*z4
373     >          - eca*F*(1.0d0-z4)/FZZ
374
375         ec_lda_rs = ecu_rs*(1.0d0-F*z4)
376     >             + ecp_rs*F*z4
377     >             - eca_rs*F*(1.0d0-z4)/FZZ
378
379         ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ)
380     >              - FZ*eca*iFZZ
381
382
383
384
385c        ********************************************
386c        **** calculate PBE96 correlation energy ****
387c        ********************************************
388         phi3 = phi**3
389         phi4 = phi3*phi
390         PON  = -ec_lda/(phi3*GAMMA)
391         tau  = DEXP(PON)
392
393         A = BOG/(tau-1.0d0+ETA)
394         A2 = A*A
395         t2 = t*t
396         t4 = t2*t2
397         t6 = t4*t2
398         Q4 = 1.0d0 + A*t2
399         Q5 = 1.0d0 + 2.0d0*A*t2
400         Q6 = 2.0d0 + A*t2
401         Q7 = 1.0d0+A*t2+A2*t4
402         Q8 = Q7*Q7
403
404         Ipbe = 1.0d0 + BOG*t2*Q4/Q7
405         Hpbe = GAMMA*phi3*DLOG(Ipbe)
406
407         Ipbe_t =  BOG*(2.0d0*t)*Q5/Q8
408         Ipbe_A = -BOG*(A*t6)   *Q6/Q8
409
410         A_ec_lda  = tau/(BETA*phi3)*A2
411         A_phi     = -3.0d0*ec_lda*tau/(BETA*phi4)*A2
412
413
414         Hpbe_ec_lda = (GAMMA*phi3/Ipbe)*Ipbe_A*A_ec_lda
415
416         Hpbe_phi    = 3.0d0*Hpbe/phi
417     >               + (GAMMA*phi3/Ipbe)*Ipbe_A*A_phi
418
419         Hpbe_t      = (GAMMA*phi3/Ipbe)*Ipbe_t
420
421         ec_lda_nup = ec_lda_zet
422     >              - zet * ec_lda_zet
423     >              + rs_n * ec_lda_rs
424         ec_lda_ndn = -ec_lda_zet
425     >              - zet  * ec_lda_zet
426     >              + rs_n * ec_lda_rs
427
428
429
430         Hpbe_nup  = ec_lda_nup   * Hpbe_ec_lda
431     >          + phi_zet*zet_nup * Hpbe_phi
432     >          + t_nup           * Hpbe_t
433
434         Hpbe_ndn  = ec_lda_ndn   * Hpbe_ec_lda
435     >          + phi_zet*zet_ndn * Hpbe_phi
436     >          + t_ndn           * Hpbe_t
437
438
439
440         ec = ec_lda + Hpbe
441
442         fncup  = ec + (ec_lda_nup + Hpbe_nup)
443         fncdn  = ec + (ec_lda_ndn + Hpbe_ndn)
444
445         xce(i)   = x_parameter*ex     + c_parameter*ec
446         fn(i,1)  = x_parameter*fnxup  + c_parameter*fncup
447         fn(i,2)  = x_parameter*fnxdn  + c_parameter*fncdn
448
449         fdn(i,1) = x_parameter*fdnxup
450         fdn(i,2) = x_parameter*fdnxdn
451         fdn(i,3) = c_parameter*t_agr*Hpbe_t
452
453      end do
454!$OMP END DO
455
456
457
458      return
459      end
460
461
462
463
464
465*    ************************************
466*    *                                  *
467*    *        gen_BEEF_BW_restricted    *
468*    *                                  *
469*    ************************************
470
471*   This routine calculates the non-Langreth terms of the BEEF-vdw exchange-correlation
472*   potential(xcp) and energy density(xce).
473*
474*
475*   Entry - n2ft3d     : number of grid points
476*           rho_in(*) :  density (nup+ndn)
477*           agr_in(*): |grad rho_in|
478*           x_parameter: scale parameter for exchange
479*           c_parameter: scale parameter for correlation
480*
481*     Exit  - xce(n2ft3d) : PBE96 exchange correlation energy density
482*             fn(n2ft3d)  : d(n*xce)/dn
483*             fdn(n2ft3d) : d(n*xce/d|grad n|
484*
485      subroutine gen_BEEF_BW_restricted(n2ft3d,rho_in,agr_in,
486     >                                x_parameter,c_parameter,alphac,
487     >                                xce,fn,fdn)
488*
489      implicit none
490      integer    n2ft3d
491      real*8     rho_in(n2ft3d)
492      real*8     agr_in(n2ft3d)
493      real*8     x_parameter,c_parameter,alphac
494      real*8     xce(n2ft3d)
495      real*8     fn(n2ft3d)
496      real*8     fdn(n2ft3d)
497
498*     **** Density cutoff parameter ****
499      real*8 DNS_CUT,ETA
500      parameter (DNS_CUT = 1.0d-20)
501      parameter (ETA     = 1.0d-20)
502
503c     ***** PBE96 GGA exchange constants ******
504      real*8 MU,KAPPA
505      parameter (MU    = 0.2195149727645171d0)
506      parameter (KAPPA = 0.8040000000000000d0)
507
508c     ****** PBE96 GGA correlation constants ******
509      real*8 GAMMA,BETA,BOG
510      parameter (GAMMA  = 0.031090690869655d0)
511      parameter (BETA   = 0.066724550603149d0)
512      parameter (BOG    = BETA/GAMMA)
513
514c     ****** Perdew-Wang92 LDA correlation coefficients *******
515      real*8 GAM,iGAM,FZZ,iFZZ
516      parameter (GAM    = 0.519842099789746329d0)
517      parameter (iGAM   = 1.0d0/GAM)
518      parameter (FZZ    = (8.0d0/(9.0d0*GAM)) )
519      parameter (iFZZ    = 0.125d0*9.0d0*GAM)
520
521      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
522      parameter (A_1  = 0.0310907d0)
523      !parameter (A_1  = 0.031091d0)
524      parameter (A1_1 = 0.2137000d0)
525      parameter (B1_1 = 7.5957000d0)
526      parameter (B2_1 = 3.5876000d0)
527      parameter (B3_1 = 1.6382000d0)
528      parameter (B4_1 = 0.4929400d0)
529
530      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
531      parameter (A_2  =  0.01554535d0)
532      !parameter (A_2  =  0.015545d0)
533      parameter (A1_2 =  0.20548000d0)
534      parameter (B1_2 = 14.11890000d0)
535      parameter (B2_2 =  6.19770000d0)
536      parameter (B3_2 =  3.36620000d0)
537      parameter (B4_2 =  0.62517000d0)
538
539      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
540      parameter (A_3  =  0.0168869d0)
541      !parameter (A_3  =  0.016887d0)
542      parameter (A1_3 =  0.1112500d0)
543      parameter (B1_3 = 10.3570000d0)
544      parameter (B2_3 =  3.6231000d0)
545      parameter (B3_3 =  0.8802600d0)
546      parameter (B4_3 =  0.4967100d0)
547
548c     **** other constants ****
549      real*8 onethird,fourthird,fivethird,onesixthm
550      real*8 twothird,sevensixthm,sevensixths
551      real*8 onethirdm
552      parameter (onethird=1.0d0/3.0d0)
553      parameter (onethirdm=-1.0d0/3.0d0)
554      parameter (twothird=2.0d0/3.0d0)
555      parameter (fourthird=4.0d0/3.0d0)
556      parameter (fivethird=5.0d0/3.0d0)
557      parameter (onesixthm=-1.0d0/6.0d0)
558      parameter (sevensixthm=-7.0d0/6.0d0)
559      parameter (sevensixths=7.0d0/6.0d0)
560
561
562*---- Beef expansion parameters given by Wellendorff et al-----------------*
563      real*8 malphac
564c      real*8 alphac,malphac
565c      parameter (alphac=0.6001664769d0)
566c      parameter (malphac=1.0d0-alphac)
567      real*8 am(0:29)
568      data am(0:29)/1.516501714d0,  4.413532099d-1,-9.182135241d-2,
569     >            -2.352754331d-2, 3.418828455d-2, 2.411870076d-3,
570     >            -1.416381352d-2, 6.975895581d-4, 9.859205137d-3,
571     >            -6.737855051d-3,-1.573330824d-3, 5.036146253d-3,
572     >            -2.569472453d-3,-9.874953976d-4, 2.033722895d-3,
573     >            -8.018718848d-4,-6.688078723d-4, 1.030936331d-3,
574     >            -3.673838660d-4,-4.213635394d-4,5.761607992d-4,
575     >            -8.346503735d-5,-4.458447585d-4,4.601290092d-4,
576     >            -5.231775398d-6,-4.239570471d-4,3.750190679d-4,
577     >             2.114938125d-5, -1.904911565d-4,7.384362421d-5/
578
579c     **** local variables ****
580      integer i,j
581      real*8 n,agr
582      real*8 kf,ks,s,n_onethird,pi,rs_scale
583      real*8 fdnx_const
584      real*8 rs,rss,t,t2,t4,t6
585      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q8,Q9,B
586      real*8 Ht
587      real*8 B_ec,Hrs,H_B
588      real*8 F,Fs
589
590      real*8 ex_lda,ec_lda
591      real*8 ec_lda_rs
592      real*8 ex,ec,H
593      real*8 fnx,fdnx,fnc,fdnc
594
595
596      real*8 p0,p1,p2,dp,s2,oneovers2,Ft,dp2,sgn
597
598
599      malphac=1.0d0-alphac
600      pi         = 4.0d0*datan(1.0d0)
601      rs_scale   = (0.75d0/pi)**onethird
602      fdnx_const = -3.0d0/(8.0d0*pi)
603
604!$OMP DO
605      do i=1,n2ft3d
606         n     = rho_in(i)+ETA
607         agr   = agr_in(i)
608
609c        ***** calculate unpolarized Exchange energies and potentials *****
610         n_onethird = (3.0d0*n/pi)**onethird
611         ex_lda     = -0.75d0*n_onethird
612
613         kf = (3.0d0*pi*pi*n)**onethird
614         s  = agr/(2.0d0*kf*n)
615
616         t  = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0
617         if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t)
618         F  = 0.0d0
619         Ft = 0.0d0
620         s2   = t*t - 1.0d0
621         if (dabs(s2).lt.1.0d-12) then
622            if (t.gt.0.0d0) then
623               do j=0,29
624                  F  = F  + am(j)
625                  Ft = Ft + am(j)*0.5d0*dble(j*(j+1))
626               end do
627            else
628               sgn = 1.0d0
629               do j=0,29
630                  F  = F  + sgn*am(j)
631                  Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1))
632                  sgn = -sgn
633               end do
634            end if
635         else
636            oneovers2 = 1.0d0/s2
637            p0 = 1.0d0
638            dp = 0
639            F  = F  + am(0)*p0
640            Ft = Ft + am(0)*dp
641            p1 = t
642            dp = 1.0d0
643            F  = F  + am(1)*t
644            Ft = Ft + am(1)*dp
645            do j=2,29
646               p2    = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0)
647               dp2   = dble(j)*oneovers2*(t*p2-p1)
648               F  = F + am(j)*p2
649               Ft = Ft + am(j)*dp2
650               p0 = p1
651               p1 = p2
652            end do
653         end if
654         Fs  = (16.0d0*s/(4.0d0+s*s)**2)*Ft
655
656         ex   = ex_lda*F
657         fnx  = fourthird*(ex - ex_lda*Fs*s)
658         fdnx = fdnx_const*Fs
659
660
661*        *********************************************************************
662c        ***** calculate unpolarized correlation energies and potentials *****
663*        *********************************************************************
664
665c        **** calculate rs and t ****
666         rs    = rs_scale/(n**onethird)
667         rss   = dsqrt(rs)
668
669         kf = (3.0d0*pi*pi*n)**onethird
670         ks = dsqrt(4.0d0*kf/pi)
671         t  = agr/(2.0d0*ks*n)
672
673
674c        **** unpolarized LDA correlation energy ****
675c        **** ec_p = correlation energy          ****
676c        ****   ec_p_rs = dec_p/drs              ****
677c        ****   uc_p    = dec_p/dn               ****
678         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs)
679c        **** PBE96 correlation energy  corrections ****
680         t2 = t*t
681         t4 = t2*t2
682         B = -ec_lda/GAMMA
683         B = BOG/(exp(B)-1.0d0+ETA)
684         Q4 = 1.0d0 + B*t2
685         Q5 = 1.0d0 + B*t2 + B*B*t4
686         H = GAMMA*dlog(1.0d0 + BOG*Q4*t2/Q5)
687
688
689c        **** PBE96 correlation fdn and fdnc derivatives ****
690         t6   = t4*t2
691
692         B_ec = (B/BETA)*(BOG+B)
693
694         Q8  = Q5*Q5+BOG*Q4*Q5*t2
695         Q9  = 1.0d0+2*B*t2
696         H_B  = -BETA*B*t6*(2.0d0+B*t2)/Q8
697         Hrs  = H_B*B_ec*ec_lda_rs
698
699         Ht  = 2.0d0*BETA*Q9/Q8*t
700
701         ec   = ec_lda + malphac*H
702         fnc  = ec  - (onethird*rs*ec_lda_rs)
703     >             - malphac*(onethird*rs*Hrs)
704     >             - malphac*(sevensixths*t*Ht)
705         fdnc = malphac*(0.5d0* Ht/ks)
706
707         xce(i) = x_parameter*ex   + c_parameter*ec
708         fn(i)  = x_parameter*fnx  + c_parameter*fnc
709         fdn(i) = x_parameter*fdnx + c_parameter*fdnc
710
711      end do
712!$OMP END DO
713
714      return
715      end
716
717