1*    ************************************
2*    *					*
3*    *	    gen_PBEsol_unrestricted	*
4*    *					*
5*    ************************************
6*
7*    This function returns the PBEsol exchange-correlation
8*  energy density, xce, and its derivatives with respect
9*  to nup, ndn, |grad nup|, |grad ndn|, and |grad n|.
10*
11*   Entry - n2ft3d     : number of grid points
12*           dn_in(*,2) : spin densites nup and ndn
13*           agr_in(*,3): |grad nup|, |grad ndn|, and |grad n|
14*           x_parameter: scale parameter for exchange
15*           c_parameter: scale parameter for correlation
16*
17*   Exit - xce(*)  : PBEsol energy density
18*        - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn
19*        - fdn(*,3): d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|
20*                    d(n*xce)/d|grad n|
21
22      subroutine gen_PBEsol_BW_unrestricted(n2ft3d,
23     >                           dn_in,agr_in,
24     >                           x_parameter,c_parameter,
25     >                           xce,fn,fdn)
26      implicit none
27
28      integer n2ft3d
29      real*8 dn_in(n2ft3d,2)
30      real*8 agr_in(n2ft3d,3)
31      real*8 x_parameter,c_parameter
32      real*8 xce(n2ft3d)
33      real*8 fn(n2ft3d,2)
34      real*8 fdn(n2ft3d,3)
35
36*     **** Density cutoff parameter ****
37      real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
38      parameter (DNS_CUT = 1.0d-20)
39      parameter (ETA=1.0d-20)
40      parameter (ETA2=1.0d-14)
41      parameter (alpha_zeta=(1.0d0-ETA2))
42      parameter (alpha_zeta2=(1.0d0-ETA2))
43
44c     ***** PBEsol GGA exchange constants ******
45      real*8 MU,KAPPA
46      parameter (MU    = (10.0d0/81.0d0))
47      parameter (KAPPA = 0.8040000000000000d0)
48
49c     ****** PBEsol GGA correlation constants ******
50      real*8 GAMMA,BETA,BOG
51      parameter (GAMMA	= 0.031090690869655d0)
52      !parameter (BETA	= 0.066724550603149d0)
53      parameter (BETA	= 0.046d0)
54      parameter (BOG    = BETA/GAMMA)
55
56
57c     ****** Perdew-Wang92 LDA correlation coefficients *******
58      real*8 GAM,iGAM,FZZ,iFZZ
59      parameter (GAM  	= 0.519842099789746329d0)
60      parameter (iGAM  	= 1.0d0/GAM)
61      parameter (FZZ    = (8.0d0/(9.0d0*GAM)) )
62      parameter (iFZZ    = 0.125d0*9.0d0*GAM)
63
64      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
65      parameter (A_1  = 0.0310907d0)
66      !parameter (A_1  = 0.031091d0)
67      parameter (A1_1 =	0.2137000d0)
68      parameter (B1_1 =	7.5957000d0)
69      parameter (B2_1 =	3.5876000d0)
70      parameter (B3_1 =	1.6382000d0)
71      parameter (B4_1 =	0.4929400d0)
72
73      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
74      parameter (A_2  =  0.01554535d0)
75      !parameter (A_2  =  0.015545d0)
76      parameter (A1_2 =	 0.20548000d0)
77      parameter (B1_2 =	14.11890000d0)
78      parameter (B2_2 =	 6.19770000d0)
79      parameter (B3_2 =	 3.36620000d0)
80      parameter (B4_2 =	 0.62517000d0)
81
82      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
83      parameter (A_3  =  0.0168869d0)
84      !parameter (A_3  =  0.016887d0)
85      parameter (A1_3 =	 0.1112500d0)
86      parameter (B1_3 =	10.3570000d0)
87      parameter (B2_3 =	 3.6231000d0)
88      parameter (B3_3 =	 0.8802600d0)
89      parameter (B4_3 =	 0.4967100d0)
90
91c     **** other constants ****
92      real*8 onethird,fourthird,fivethird,onesixthm
93      real*8 twothird,sevensixthm
94      real*8 onethirdm
95      parameter (onethird=1.0d0/3.0d0)
96      parameter (onethirdm=-1.0d0/3.0d0)
97      parameter (twothird=2.0d0/3.0d0)
98      parameter (fourthird=4.0d0/3.0d0)
99      parameter (fivethird=5.0d0/3.0d0)
100      parameter (onesixthm=-1.0d0/6.0d0)
101      parameter (sevensixthm=-7.0d0/6.0d0)
102
103c     **** local variables ****
104      integer i
105      real*8 n,agr
106      real*8 nup,agrup
107      real*8 ndn,agrdn
108      real*8 kf,ks,s,P0,n_onethird,pi,rs_scale
109      real*8 rs         ! Wigner radius
110      real*8 rss        ! rss  = sqrt(rs)
111      real*8 rs_n       ! rs_n = n*drs/dn
112      real*8 t,t2,t4,t6
113      real*8 t_nup      ! t_nup = n*dt/dnup
114      real*8 t_ndn      ! t_ndn = n*dt/dndn
115      real*8 t_agr      ! t_agr = n*dt/dagr
116      real*8 zet,twoksg
117      real*8 zet_nup    ! zet_nup = n*dzet/dnup
118      real*8 zet_ndn    ! zet_nup = n*dzet/dnup
119      real*8 zetp_1_3,zetm_1_3
120      real*8 zetpm_1_3,zetmm_1_3
121      real*8 phi,phi3,phi4
122      real*8 phi_zet
123      real*8 A,A2
124      real*8 A_phi,A_ec_lda
125      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8
126      real*8 PON,FZ,z4
127      real*8 tau
128      real*8 F
129      real*8 Fs                ! dF/ds
130      real*8 Hpbe
131      real*8 Hpbe_t            ! dHpbe/dt
132      real*8 Hpbe_phi          ! dHpbe/dphi
133      real*8 Hpbe_ec_lda       ! dHpbe/d(ec_lda)
134      real*8 Hpbe_nup,Hpbe_ndn ! n*dHpbe/dnup, n*dHpbe/dndn
135      real*8 Ipbe
136      real*8 Ipbe_t,Ipbe_A     ! dIpbe/dt, dIpbe/dA
137
138      real*8 exup,exdn,ex,ex_lda
139      real*8 ecu,ecp,eca,ec,ec_lda
140      real*8 ecu_rs,ecp_rs,eca_rs
141      real*8 ec_lda_rs,ec_lda_zet  ! d(ec_lda)/drs, d(ec_lda)/dzet
142      real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn
143      real*8 fnxup,fdnxup          ! d(n*ex)/dnup, d(n*ex)/dndn
144      real*8 fnxdn,fdnxdn          ! d(n*ex)/d|grad nup|, d(n*ex)/d|grad ndn|
145      real*8 fncup,fncdn           ! d(n*ec)/dnup, d(n*ec)/dndn
146      real*8 fdnx_const
147
148      pi = 4.0d0*datan(1.0d0)
149      rs_scale = (0.75d0/pi)**onethird
150      fdnx_const = -3.0d0/(8.0d0*pi)
151
152
153!$OMP DO
154      do i=1,n2ft3d
155         nup     = dn_in(i,1)+ETA
156         agrup   = agr_in(i,1)
157
158         ndn     = dn_in(i,2)+ETA
159         agrdn   = agr_in(i,2)
160
161c        ****************************************************************
162c        ***** calculate polarized Exchange energies and potentials *****
163c        ****************************************************************
164
165c        ************
166c        **** up ****
167c        ************
168         n     = 2.0d0*nup
169         agr   = 2.0d0*agrup
170
171         n_onethird = (3.0d0*n/pi)**onethird
172         ex_lda     = -0.75d0*n_onethird
173
174         kf = (3.0d0*pi*pi*n)**onethird
175         s  = agr/(2.0d0*kf*n)
176         P0 = 1.0d0 + (MU/KAPPA)*s*s
177
178         F   = (1.0d0 + KAPPA - KAPPA/P0)
179         Fs  = 2.0d0*MU/(P0*P0)*s
180
181         exup = ex_lda*F
182         fnxup = fourthird*(exup - ex_lda*Fs*s)
183         fdnxup = fdnx_const*Fs
184
185c        **************
186c        **** down ****
187c        **************
188         n     = 2.0d0*ndn
189         agr   = 2.0d0*agrdn
190
191         n_onethird = (3.0d0*n/pi)**onethird
192         ex_lda     = -0.75d0*n_onethird
193
194         kf = (3.0d0*pi*pi*n)**onethird
195         s  = agr/(2.0d0*kf*n)
196         P0 = 1.0d0 + (MU/KAPPA)*s*s
197
198         F   = (1.0d0 + KAPPA - KAPPA/P0)
199         Fs  = 2.0d0*MU/(P0*P0)*s
200
201         exdn   = ex_lda*F
202         fnxdn  = fourthird*(exdn - ex_lda*Fs*s)
203         fdnxdn = fdnx_const*Fs
204
205         n = nup+ndn
206
207         ex = (exup*nup+ exdn*ndn)/ n
208
209c        *******************************************************************
210c        ***** calculate polarized correlation energies and potentials *****
211c        *******************************************************************
212         agr   = agr_in(i,3)
213
214         zet = (nup-ndn)/n
215c         if (zet.gt.0.0d0) zet = zet - ETA2
216c         if (zet.lt.0.0d0) zet = zet + ETA2
217c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n**2
218c        if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2
219c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n
220c        zet_nup =  2*ndn/n
221c        zet_ndn = -2*nup/n
222         zet_nup = -(zet - 1.0d0)
223         zet_ndn = -(zet + 1.0d0)
224         zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm
225         zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm
226         zetp_1_3  = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2
227         zetm_1_3  = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2
228
229
230         phi = 0.5d0*( zetp_1_3**2 + zetm_1_3**2)
231         phi_zet = alpha_zeta*( zetpm_1_3 - zetmm_1_3)/3.0d0
232         F =(  (1.0d0+zet*alpha_zeta)*zetp_1_3
233     >       + (1.0d0-zet*alpha_zeta)*zetm_1_3
234     >       - 2.0d0)*iGAM
235
236         FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM)
237
238
239
240*        **** calculate Wigner radius ****
241         rs    = rs_scale/(n**onethird)
242         rss   = dsqrt(rs)
243
244*        **** calculate n*drs/dn ****
245c        rs_n = onethirdm*rs/n
246         rs_n = onethirdm*rs
247
248
249
250c        **** calculate t ****
251         kf = (3.0d0*pi*pi*n)**onethird
252         ks = dsqrt(4.0d0*kf/pi)
253
254         twoksg = 2.0d0*ks*phi
255
256         t  = agr/(twoksg*n)
257
258*        *** calculate n*dt/dnup, n*dt/dndn, n*dt/d|grad n| ****
259         t_nup = sevensixthm*t - (phi_zet)*(zet_nup)*t/phi
260         t_ndn = sevensixthm*t - (phi_zet)*(zet_ndn)*t/phi
261         t_agr  = 1.0d0/(twoksg)
262
263
264
265
266c        **************************************************
267c        ***** compute LSDA correlation energy density ****
268c        **************************************************
269         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs)
270         call LSDT(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs)
271         call LSDT(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs)
272
273         z4 = zet**4
274
275         ec_lda = ecu*(1.0d0-F*z4)
276     >          + ecp*F*z4
277     >          - eca*F*(1.0d0-z4)/FZZ
278
279         ec_lda_rs = ecu_rs*(1.0d0-F*z4)
280     >             + ecp_rs*F*z4
281     >             - eca_rs*F*(1.0d0-z4)/FZZ
282
283         ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ)
284     >              - FZ*eca*iFZZ
285
286
287
288
289c        ********************************************
290c        **** calculate PBEsol correlation energy ****
291c        ********************************************
292         phi3 = phi**3
293         phi4 = phi3*phi
294         PON  = -ec_lda/(phi3*GAMMA)
295         tau  = DEXP(PON)
296
297         A = BOG/(tau-1.0d0+ETA)
298         A2 = A*A
299         t2 = t*t
300         t4 = t2*t2
301         t6 = t4*t2
302         Q4 = 1.0d0 + A*t2
303         Q5 = 1.0d0 + 2.0d0*A*t2
304         Q6 = 2.0d0 + A*t2
305         Q7 = 1.0d0+A*t2+A2*t4
306         Q8 = Q7*Q7
307
308         Ipbe = 1.0d0 + BOG*t2*Q4/Q7
309         Hpbe = GAMMA*phi3*DLOG(Ipbe)
310
311         Ipbe_t =  BOG*(2.0d0*t)*Q5/Q8
312         Ipbe_A = -BOG*(A*t6)   *Q6/Q8
313
314         A_ec_lda  = tau/(BETA*phi3)*A2
315         A_phi     = -3.0d0*ec_lda*tau/(BETA*phi4)*A2
316
317
318         Hpbe_ec_lda = (GAMMA*phi3/Ipbe)*Ipbe_A*A_ec_lda
319
320         Hpbe_phi    = 3.0d0*Hpbe/phi
321     >               + (GAMMA*phi3/Ipbe)*Ipbe_A*A_phi
322
323         Hpbe_t      = (GAMMA*phi3/Ipbe)*Ipbe_t
324
325         ec_lda_nup = ec_lda_zet
326     >              - zet * ec_lda_zet
327     >              + rs_n * ec_lda_rs
328         ec_lda_ndn = -ec_lda_zet
329     >              - zet  * ec_lda_zet
330     >              + rs_n * ec_lda_rs
331
332
333
334         Hpbe_nup  = ec_lda_nup   * Hpbe_ec_lda
335     >          + phi_zet*zet_nup * Hpbe_phi
336     >          + t_nup           * Hpbe_t
337
338         Hpbe_ndn  = ec_lda_ndn   * Hpbe_ec_lda
339     >          + phi_zet*zet_ndn * Hpbe_phi
340     >          + t_ndn           * Hpbe_t
341
342
343
344         ec = ec_lda + Hpbe
345
346         fncup  = ec + (ec_lda_nup + Hpbe_nup)
347         fncdn  = ec + (ec_lda_ndn + Hpbe_ndn)
348
349         xce(i)   = x_parameter*ex     + c_parameter*ec
350         fn(i,1)  = x_parameter*fnxup  + c_parameter*fncup
351         fn(i,2)  = x_parameter*fnxdn  + c_parameter*fncdn
352
353         fdn(i,1) = x_parameter*fdnxup
354         fdn(i,2) = x_parameter*fdnxdn
355         fdn(i,3) = c_parameter*t_agr*Hpbe_t
356
357      end do
358!$OMP END DO
359
360      return
361      end
362
363
364
365
366
367*    ************************************
368*    *					*
369*    *	    gen_PBEsol_BW_restricted	*
370*    *					*
371*    ************************************
372*
373*   This routine calculates the PBEsol exchange-correlation
374*   potential(xcp) and energy density(xce).
375*
376*
377*   Entry - n2ft3d     : number of grid points
378*           rho_in(*) :  density (nup+ndn)
379*           agr_in(*): |grad rho_in|
380*           x_parameter: scale parameter for exchange
381*           c_parameter: scale parameter for correlation
382*
383*     Exit  - xce(n2ft3d) : PBEsol exchange correlation energy density
384*             fn(n2ft3d)  : d(n*xce)/dn
385*             fdn(n2ft3d) : d(n*xce/d|grad n|
386*
387      subroutine gen_PBEsol_BW_restricted(n2ft3d,rho_in,agr_in,
388     >                                x_parameter,c_parameter,
389     >                                xce,fn,fdn)
390      implicit none
391
392      integer    n2ft3d
393      real*8     rho_in(n2ft3d)
394      real*8     agr_in(n2ft3d)
395      real*8     x_parameter,c_parameter
396      real*8     xce(n2ft3d)
397      real*8     fn(n2ft3d)
398      real*8     fdn(n2ft3d)
399
400
401*     **** Density cutoff parameter ****
402      real*8 DNS_CUT,ETA
403      parameter (DNS_CUT = 1.0d-20)
404      parameter (ETA     = 1.0d-20)
405
406c     ***** PBEsol GGA exchange constants ******
407      real*8 MU,KAPPA
408      parameter (MU    = (10.0d0/81.0d0))
409      parameter (KAPPA = 0.8040000000000000d0)
410
411c     ****** PBEsol GGA correlation constants ******
412      real*8 GAMMA,BETA,BOG
413      parameter (GAMMA	= 0.031090690869655d0)
414      !parameter (BETA	= 0.066724550603149d0)
415      parameter (BETA	= 0.046d0)
416      parameter (BOG    = BETA/GAMMA)
417
418
419c     ****** Perdew-Wang92 LDA correlation coefficients *******
420      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
421      parameter (A_1  = 0.0310907d0)
422      parameter (A1_1 =	0.2137000d0)
423      parameter (B1_1 =	7.5957000d0)
424      parameter (B2_1 =	3.5876000d0)
425      parameter (B3_1 =	1.6382000d0)
426      parameter (B4_1 =	0.4929400d0)
427
428      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
429      parameter (A_2  =  0.01554535d0)
430      parameter (A1_2 =	 0.20548000d0)
431      parameter (B1_2 =	14.11890000d0)
432      parameter (B2_2 =	 6.19770000d0)
433      parameter (B3_2 =	 3.36620000d0)
434      parameter (B4_2 =	 0.62517000d0)
435
436      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
437      parameter (A_3  =  0.0168869d0)
438      parameter (A1_3 =	 0.1112500d0)
439      parameter (B1_3 =	10.3570000d0)
440      parameter (B2_3 =	 3.6231000d0)
441      parameter (B3_3 =	 0.8802600d0)
442      parameter (B4_3 =	 0.4967100d0)
443
444c     **** other constants ****
445      real*8 onethird,fourthird,sevensixths
446      parameter (onethird=1.0d0/3.0d0)
447      parameter (fourthird=4.0d0/3.0d0)
448      parameter (sevensixths=7.0d0/6.0d0)
449
450c     **** local variables ****
451      integer i
452      real*8 n,agr
453      real*8 kf,ks,s,P0,n_onethird,pi,rs_scale
454      real*8 fdnx_const
455      real*8 rs,rss,t,t2,t4,t6
456      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q8,Q9,B
457      real*8 Ht
458      real*8 B_ec,Hrs,H_B
459      real*8 F,Fs
460
461      real*8 ex_lda,ec_lda
462      real*8 ec_lda_rs
463      real*8 ex,ec,H
464      real*8 fnx,fdnx,fnc,fdnc
465
466
467      pi         = 4.0d0*datan(1.0d0)
468      rs_scale   = (0.75d0/pi)**onethird
469      fdnx_const = -3.0d0/(8.0d0*pi)
470
471!$OMP DO
472      do i=1,n2ft3d
473         n     = rho_in(i)+ETA
474         agr   = agr_in(i)
475
476c        ***** calculate unpolarized Exchange energies and potentials *****
477         n_onethird = (3.0d0*n/pi)**onethird
478         ex_lda     = -0.75d0*n_onethird
479
480         kf = (3.0d0*pi*pi*n)**onethird
481         s  = agr/(2.0d0*kf*n)
482         P0 = 1.0d0 + (MU/KAPPA)*s*s
483
484c        if (n.gt.DNS_CUT) then
485c           F   = (1.0d0 + KAPPA - KAPPA/P0)
486c           Fs  = 2.0d0*MU/(P0*P0)*s
487c        else
488c           F   = 1.0d0
489c           Fs  = 0.0d0
490c        end if
491         F   = (1.0d0 + KAPPA - KAPPA/P0)
492         Fs  = 2.0d0*MU/(P0*P0)*s
493
494         ex   = ex_lda*F
495         fnx  = fourthird*(ex - ex_lda*Fs*s)
496         fdnx = fdnx_const*Fs
497
498
499
500*        *********************************************************************
501c        ***** calculate unpolarized correlation energies and potentials *****
502*        *********************************************************************
503
504c        **** calculate rs and t ****
505         rs    = rs_scale/(n**onethird)
506         rss   = dsqrt(rs)
507
508         kf = (3.0d0*pi*pi*n)**onethird
509         ks = dsqrt(4.0d0*kf/pi)
510         t  = agr/(2.0*ks*n)
511
512
513c        **** unpolarized LDA correlation energy ****
514c        **** ec_p = correlation energy          ****
515c        ****   ec_p_rs = dec_p/drs              ****
516c        ****   uc_p    = dec_p/dn               ****
517         call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs)
518c        **** PBEsol correlation energy  corrections ****
519         t2 = t*t
520         t4 = t2*t2
521         B = -ec_lda/GAMMA
522         B = BOG/(exp(B)-1.0d0+ETA)
523         Q4 = 1.0d0 + B*t2
524         Q5 = 1.0d0 + B*t2 + B*B*t4
525         H = GAMMA*dlog(1.0d0 + BOG*Q4*t2/Q5)
526
527
528c        **** PBEsol correlation fdn and fdnc derivatives ****
529         t6   = t4*t2
530
531         B_ec = (B/BETA)*(BOG+B)
532
533         Q8  = Q5*Q5+BOG*Q4*Q5*t2
534         Q9  = 1.0d0+2*B*t2
535         H_B  = -BETA*B*t6*(2.0d0+B*t2)/Q8
536         Hrs  = H_B*B_ec*ec_lda_rs
537
538         Ht  = 2.0d0*BETA*Q9/Q8*t
539
540         ec   = ec_lda + H
541         fnc = ec  - (onethird*rs*ec_lda_rs)
542     >             - (onethird*rs*Hrs)
543     >             - (sevensixths*t*Ht)
544         fdnc = 0.5d0* Ht/ks
545
546         xce(i) = x_parameter*ex   + c_parameter*ec
547         fn(i)  = x_parameter*fnx  + c_parameter*fnc
548         fdn(i) = x_parameter*fdnx + c_parameter*fdnc
549
550
551c       write(*,*) "pbe96:",i,ec,fnc,fdnc
552
553
554      end do
555!$OMP END DO
556
557      return
558      end
559c $Id$
560