1!
2MODULE corr_gga                       !<GPU:corr_gga=>corr_gga_gpu>
3  !
4  USE corr_lda,  ONLY : pw, pw_spin   !<GPU:pw_spin=>pw_spin_d,pw=>pw_d,corr_lda=>corr_lda_gpu>
5  !
6CONTAINS
7!
8!-----------------------------------------------------------------------
9SUBROUTINE perdew86( rho, grho, sc, v1c, v2c )                    !<GPU:DEVICE>
10  !-----------------------------------------------------------------------
11  !! Perdew gradient correction on correlation: PRB 33, 8822 (1986).
12  !
13  USE kinds, ONLY : DP
14  !
15  IMPLICIT NONE
16  !
17  REAL(DP), INTENT(IN) :: rho, grho
18  REAL(DP), INTENT(OUT) :: sc, v1c, v2c
19  !
20  ! ... local variables
21  !
22  REAL(DP), PARAMETER :: p1=0.023266_DP, p2=7.389d-6, p3=8.723_DP, &
23                         p4=0.472_DP
24  REAL(DP), PARAMETER :: pc1=0.001667_DP, pc2=0.002568_DP, pci=pc1 + pc2
25  REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP
26                                          ! pi34=(3/4pi)^(1/3)
27  REAL(DP) :: rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs
28  REAL(DP) :: dcna, dcnb, dcn, phi, ephi
29  !
30  rho13 = rho**third
31  rho43 = rho13**4
32  rs  = pi34 / rho13
33  rs2 = rs * rs
34  rs3 = rs * rs2
35  !
36  cna = pc2 + p1 * rs + p2 * rs2
37  cnb = 1._DP + p3 * rs + p4 * rs2 + 1.d4 * p2 * rs3
38  cn = pc1 + cna / cnb
39  !
40  drs  = - third * pi34 / rho43
41  dcna = (p1 + 2._DP * p2 * rs) * drs
42  dcnb = (p3 + 2._DP * p4 * rs + 3.d4 * p2 * rs2) * drs
43  dcn  = dcna / cnb - cna / (cnb * cnb) * dcnb
44  !
45  phi = 0.192_DP * pci / cn * SQRT(grho) * rho**(-7._DP/6._DP)
46  ! SdG: in the original paper 1.745*0.11=0.19195 is used
47  ephi = EXP( -phi )
48  !
49  sc  = grho / rho43 * cn * ephi
50  v1c = sc * ( (1._DP+phi) * dcn / cn - ((4._DP/3._DP)-(7._DP/ &
51                  6._DP)*phi) / rho )
52  v2c = cn * ephi / rho43 * (2._DP - phi)
53  !
54  RETURN
55  !
56END SUBROUTINE perdew86
57!
58!
59!-----------------------------------------------------------------------
60SUBROUTINE ggac( rho, grho, sc, v1c, v2c )                    !<GPU:DEVICE>
61  !-----------------------------------------------------------------------
62  !! Perdew-Wang GGA (PW91) correlation part
63  !
64  USE kinds,    ONLY: DP
65  !
66  IMPLICIT NONE
67  !
68  REAL(DP), INTENT(IN) :: rho, grho
69  REAL(DP), INTENT(OUT) :: sc, v1c, v2c
70  !
71  ! ... local variables
72  !
73  REAL(DP) :: rs, ec, vc
74  !
75  REAL(DP), PARAMETER :: al=0.09_DP,  pa=0.023266_DP, pb=7.389d-6, &
76                         pc=8.723_DP, pd=0.472_DP,                 &
77                         cx=-0.001667_DP, cxc0=0.002568_DP, cc0=-cx+cxc0
78  !
79  REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP, &
80                         nu=15.755920349483144_DP, be=nu*cc0,        &
81                         xkf=1.919158292677513_DP, xks=1.128379167095513_DP
82                         ! pi34=(3/4pi)^(1/3),  nu=(16/pi)*(3 pi^2)^(1/3)
83                         ! xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
84  REAL(DP) :: kf, ks, rs1, rs2, rs3, t, expe, af, bf, y, xy, qy, s1
85  REAL(DP) :: h0, dh0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, &
86              dh1, ddh1
87  !
88  rs = pi34 / rho**third
89  !
90  CALL pw( rs, 1, ec, vc )                                  !<GPU:pw=>pw_d>
91  !
92  rs1 = rs
93  rs2 = rs1 * rs1
94  rs3 = rs1 * rs2
95  !
96  kf = xkf / rs1
97  ks = xks * SQRT(kf)
98  t = SQRT(grho) / (2._DP * ks * rho)
99  !
100  expe = EXP( - 2._DP * al * ec / (be * be) )
101  af = 2._DP * al / be * (1._DP / (expe-1._DP) )
102  bf = expe * (vc - ec)
103  !
104  y = af * t * t
105  xy = (1._DP + y) / (1._DP + y + y * y)
106  qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2
107  s1 = 1._DP + 2._DP * al / be * t * t * xy
108  !
109  h0 = be * be / (2._DP * al) * LOG(s1)
110  dh0 = be * t * t / s1 * ( - 7._DP / 3._DP * xy - qy * (af * bf / &
111                            be-7._DP / 3._DP) )
112  ddh0 = be / (2._DP * ks * ks * rho) * (xy - qy) / s1
113  !
114  ee = - 100._DP * (ks / kf * t)**2
115  !
116  cna = cxc0 + pa * rs1 + pb * rs2
117  dcna = pa * rs1 + 2._DP * pb * rs2
118  cnb = 1._DP + pc * rs1 + pd * rs2 + 1.d4 * pb * rs3
119  dcnb = pc * rs1 + 2._DP * pd * rs2 + 3.d4 * pb * rs3
120  cn = cna / cnb - cx
121  dcn = dcna / cnb - cna * dcnb / (cnb * cnb)
122  !
123  h1 = nu * (cn - cc0 - 3._DP / 7._DP * cx) * t * t * EXP(ee)
124  dh1 = - third * ( h1 * (7._DP + 8._DP * ee) + nu * t * t * EXP(ee) &
125                    * dcn )
126  ddh1 = 2._DP * h1 * (1._DP + ee) * rho / grho
127  !
128  sc = rho * (h0 + h1)
129  v1c = h0 + h1 + dh0 + dh1
130  v2c = ddh0 + ddh1
131  !
132  RETURN
133  !
134END SUBROUTINE ggac
135!
136!
137!-----------------------------------------------------------------------
138SUBROUTINE glyp( rho, grho, sc, v1c, v2c )                    !<GPU:DEVICE>
139  !-----------------------------------------------------------------------
140  !! Lee Yang Parr: gradient correction part.
141  !
142  USE kinds, ONLY: DP
143  !
144  IMPLICIT NONE
145  !
146  REAL(DP), INTENT(IN) :: rho, grho
147  REAL(DP), INTENT(OUT) :: sc, v1c, v2c
148  !
149  ! ... local varibles
150  !
151  REAL(DP), PARAMETER :: a=0.04918_DP, b=0.132_DP, c=0.2533_DP, &
152                         d=0.349_DP
153  REAL(DP) :: rhom13, rhom43, rhom53, om, xl, ff, dom, dxl
154  !
155  rhom13 = rho**(-1._DP/3._DP)
156  om = EXP(-c*rhom13) / (1._DP+d*rhom13)
157  xl = 1._DP + (7._DP/3._DP) * ( c*rhom13 + d * rhom13 / (1._DP + &
158                                 d * rhom13) )
159  ff = a * b * grho / 24._DP
160  rhom53 = rhom13**5
161  !
162  sc = ff * rhom53 * om * xl
163  !
164  dom = - om * (c + d+c * d * rhom13) / (1.d0 + d * rhom13)
165  dxl = (7.d0 / 3.d0) * (c + d+2.d0 * c * d * rhom13 + c * d * d * &
166       rhom13**2) / (1.d0 + d * rhom13) **2
167  rhom43 = rhom13**4
168  !
169  v1c = - ff * rhom43 / 3.d0 * ( 5.d0 * rhom43 * om * xl + rhom53 * &
170                                    dom * xl + rhom53 * om * dxl )
171  v2c = 2.d0 * sc / grho
172  !
173  RETURN
174  !
175END SUBROUTINE glyp
176!
177!
178!---------------------------------------------------------------
179SUBROUTINE pbec( rho, grho, iflag, sc, v1c, v2c )                    !<GPU:DEVICE>
180  !---------------------------------------------------------------
181  !! PBE correlation (without LDA part)
182  !
183  !! * iflag=1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
184  !! * iflag=2: J.P.Perdew et al., PRL 100, 136406 (2008).
185  !! * iflag=3: L. Chiodo et al, PRL 108, 126402 (2012)  (PBEQ2D)
186  !
187  USE kinds,    ONLY: DP
188  !
189  IMPLICIT NONE
190  !
191  INTEGER,  INTENT(IN) :: iflag              !<GPU:VALUE>
192  REAL(DP), INTENT(IN) :: rho, grho
193  ! input: charge and squared gradient
194  REAL(DP), INTENT(OUT) :: sc, v1c, v2c
195  ! output: energy, potential
196  REAL(DP), PARAMETER :: ga = 0.0310906908696548950_DP
197  REAL(DP) :: be(3)
198  !             pbe           pbesol   pbeq2d
199  DATA be / 0.06672455060314922_DP, 0.046_DP, 0.06672455060314922_DP/
200  REAL(DP), PARAMETER :: third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0
201  REAL(DP), PARAMETER :: xkf = 1.919158292677513d0, xks = 1.128379167095513d0
202  ! pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
203  !
204  REAL(DP) :: rs, ec, vc
205  !
206  REAL(DP) :: kf, ks, t, expe, af, bf, y, xy, qy
207  REAL(DP) :: s1, h0, dh0, ddh0, sc2D, v1c2D, v2c2D
208  !
209  rs = pi34 / rho**third
210  !
211  CALL pw( rs, 1, ec, vc )                      !<GPU:pw=>pw_d>
212  !
213  kf = xkf / rs
214  ks = xks * SQRT(kf)
215  t = SQRT(grho) / (2._DP * ks * rho)
216  expe = EXP( - ec / ga )
217  af = be(iflag) / ga * (1._DP / (expe-1._DP))
218  bf = expe * (vc - ec)
219  y = af * t * t
220  xy = (1._DP + y) / (1._DP + y + y * y)
221  qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2
222  s1 = 1._DP + be(iflag) / ga * t * t * xy
223  h0 = ga * LOG(s1)
224  dh0 = be(iflag) * t * t / s1 * ( - 7._DP / 3._DP * xy - qy * (af * bf / &
225        be(iflag)-7._DP / 3._DP) )
226  ddh0 = be(iflag) / (2._DP * ks * ks * rho) * (xy - qy) / s1
227  !
228  sc  = rho * h0
229  v1c = h0 + dh0
230  v2c = ddh0
231  ! q2D
232  IF (iflag == 3) THEN
233     CALL cpbe2d( rho, grho, sc2D, v1c2D, v2c2D )       !<GPU:cpbe2d=>cpbe2d_d>
234     sc  = sc  + sc2D
235     v1c = v1c + v1c2D
236     v2c = v2c + v2c2D
237  ENDIF
238  !
239  RETURN
240  !
241END SUBROUTINE pbec
242!
243!
244! ===========> SPIN <===========
245!
246!-----------------------------------------------------------------------
247SUBROUTINE perdew86_spin( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c )                    !<GPU:DEVICE>
248  !---------------------------------------------------------------------
249  !! Perdew gradient correction on correlation: PRB 33, 8822 (1986)
250  !! spin-polarized case.
251  !
252  USE kinds,    ONLY: DP
253  !
254  IMPLICIT NONE
255  !
256  REAL(DP), INTENT(IN) :: rho
257  !! the total charge density
258  REAL(DP), INTENT(IN) :: zeta
259  !! the magnetization
260  REAL(DP), INTENT(IN) :: grho
261  !! the gradient of the charge squared
262  REAL(DP), INTENT(OUT) :: sc
263  !! correlation energies
264  REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw !v1c(2)
265  !! derivative of correlation wr. rho
266  REAL(DP), INTENT(OUT) :: v2c
267  !! derivatives of correlation wr. grho
268  !
269  ! ... local variables
270  !
271  REAL(DP), PARAMETER :: p1=0.023266_DP, p2=7.389D-6, p3=8.723_DP, &
272                         p4=0.472_DP
273  REAL(DP), PARAMETER :: pc1=0.001667_DP, pc2 = 0.002568_DP, pci=pc1+pc2
274  REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP
275  !                                         pi34=(3/4pi)^(1/3)
276  !
277  REAL(DP) :: rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs
278  REAL(DP) :: dcna, dcnb, dcn, phi, ephi, dd, ddd
279  !
280  rho13 = rho**third
281  rho43 = rho13**4
282  rs  = pi34 / rho13
283  rs2 = rs * rs
284  rs3 = rs * rs2
285  !
286  cna = pc2 + p1 * rs + p2 * rs2
287  cnb = 1._DP + p3 * rs + p4 * rs2 + 1.D4 * p2 * rs3
288  cn  = pc1 + cna / cnb
289  !
290  drs  = -third * pi34 / rho43
291  dcna = (p1 + 2._DP * p2 * rs) * drs
292  dcnb = (p3 + 2._DP * p4 * rs + 3.D4 * p2 * rs2) * drs
293  dcn  = dcna / cnb - cna / (cnb * cnb) * dcnb
294  phi  = 0.192_DP * pci / cn * SQRT(grho) * rho**(-7._DP/6._DP)
295  !SdG: in the original paper 1.745*0.11=0.19195 is used
296  !
297  dd  = (2._DP)**third * SQRT( ( (1.0_DP + zeta) * 0.5_DP)**(5.0_DP/ &
298       3._DP) + ( (1.0_DP - zeta) * 0.5d0)**(5._DP/3._DP) )
299  !
300  ddd = (2._DP)**(-4.0_DP/3._DP) * 5._DP * ( ((1._DP + zeta) &
301       * 0.5_DP)**(2._DP/3._DP) - ((1._DP - zeta)*0.5d0)**(2._DP/ &
302       3._DP))/(3._DP*dd)
303  !
304  ephi = EXP(-phi)
305  sc = grho / rho43 * cn * ephi / dd
306  !
307  v1c_up = sc * ( (1._DP + phi) * dcn / cn - ( (4._DP / 3._DP) - &
308       (7._DP/6._DP)*phi)/rho) - sc * ddd / dd * (1._DP - zeta)/rho
309  v1c_dw = sc * ( (1._DP + phi) * dcn / cn - ( (4._DP / 3._DP) - &
310       (7._DP/6._DP)*phi)/rho) + sc * ddd / dd * (1._DP + zeta)/rho
311  v2c = cn * ephi / rho43 * (2._DP - phi) / dd
312  !
313  RETURN
314  !
315END SUBROUTINE perdew86_spin
316!
317!
318!-----------------------------------------------------------------------
319SUBROUTINE ggac_spin( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c )                    !<GPU:DEVICE>
320  !---------------------------------------------------------------------
321  !! Perdew-Wang GGA (PW91) correlation part - spin-polarized
322  !
323  USE kinds, ONLY: DP
324  !
325  IMPLICIT NONE
326  !
327  REAL(DP), INTENT(IN) :: rho
328  !! the total charge density
329  REAL(DP), INTENT(IN) :: zeta
330  !! the magnetization
331  REAL(DP), INTENT(IN) :: grho
332  !! the gradient of the charge squared
333  REAL(DP), INTENT(OUT) :: sc
334  !! correlation energies
335  REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw
336  !! derivative of correlation wr. rho
337  REAL(DP), INTENT(OUT) :: v2c
338  !! derivatives of correlation wr. grho
339  !
340  ! ... local variables
341  !
342  REAL(DP), PARAMETER :: al=0.09_DP,  pa=0.023266_DP, pb=7.389D-6, &
343                         pc=8.723_DP, pd=0.472_DP
344  REAL(DP), PARAMETER :: cx=-0.001667_DP, cxc0=0.002568_DP, cc0=-cx+cxc0
345  REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP
346  !                                         pi34=(3/4pi)^(1/3)
347  REAL(DP), PARAMETER :: nu=15.755920349483144_DP, be=nu*cc0
348  !                      nu=(16/pi)*(3 pi^2)^(1/3)
349  REAL(DP), PARAMETER :: xkf=1.919158292677513_DP, xks=1.128379167095513_DP
350  !                      xkf=(9 pi/4)^(1/3)      , xks=sqrt(4/pi)
351  !
352  REAL(DP) :: rs, ec, vc_up, vc_dn
353  REAL(DP) :: kf, ks, rs2, rs3, t, expe, af, y, xy, qy, s1, h0, ddh0, ee, &
354              cn, dcn, cna, dcna, cnb, dcnb, h1, dh1, ddh1, fz, fz2, fz3, &
355              fz4, dfz, bfup, bfdw, dh0up, dh0dw, dh0zup, dh0zdw, dh1zup, &
356              dh1zdw
357  !
358  rs = pi34 / rho**third
359  !
360  CALL pw_spin( rs, zeta, ec, vc_up, vc_dn )                                 !<GPU:pw_spin=>pw_spin_d>
361  !
362  rs2 = rs * rs
363  rs3 = rs * rs2
364  kf = xkf / rs
365  ks = xks * SQRT(kf)
366  !
367  fz = 0.5_DP * ( (1._DP+zeta)**(2._DP/3._DP) + (1._DP-zeta)**(2._DP/3._DP) )
368  fz2 = fz  * fz
369  fz3 = fz2 * fz
370  fz4 = fz3 * fz
371  dfz = ( (1._DP+zeta)**(-1._DP/3._DP) - (1._DP-zeta)**(-1._DP/3._DP) )/3._DP
372  !
373  t  = SQRT(grho) / (2._DP * fz * ks * rho)
374  expe = EXP( - 2._DP * al * ec / (fz3 * be * be) )
375  af = 2._DP * al / be * (1._DP / (expe-1._DP) )
376  bfup = expe * (vc_up - ec) / fz3
377  bfdw = expe * (vc_dn - ec) / fz3
378  !
379  y = af * t * t
380  xy = (1._DP + y) / (1._DP + y + y * y)
381  qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2
382  s1 = 1._DP + 2._DP * al / be * t * t * xy
383  !
384  h0 = fz3 * be * be / (2._DP * al) * LOG(s1)
385  dh0up  = be * t * t * fz3 / s1 * ( - 7._DP / 3._DP * xy - qy *   &
386           (af * bfup / be-7._DP / 3._DP) )
387  dh0dw  = be * t * t * fz3 / s1 * ( - 7._DP / 3._DP * xy - qy *   &
388           (af * bfdw / be-7._DP / 3._DP) )
389  dh0zup = (3._DP * h0 / fz - be * t * t * fz2 / s1 * (2._DP *xy - &
390           qy * (3._DP * af * expe * ec / fz3 / be+2._DP) ) ) *    &
391           dfz * (1._DP - zeta )
392  dh0zdw = -(3._DP * h0 / fz - be * t * t * fz3 / s1 * (2._DP*xy - &
393           qy * (3._DP * af * expe * ec / fz3 / be + 2._DP) ) ) *  &
394           dfz * (1._DP + zeta )
395  ddh0   = be * fz / (2._DP * ks * ks * rho) * (xy - qy) / s1
396  !
397  ee = -100._DP * fz4 * (ks / kf * t)**2
398  cna  = cxc0 + pa * rs + pb * rs2
399  dcna = pa * rs + 2._DP * pb * rs2
400  cnb  = 1._DP + pc * rs + pd * rs2 + 1.D4 * pb * rs3
401  dcnb = pc * rs + 2._DP * pd * rs2 + 3.D4 * pb * rs3
402  cn  = cna / cnb - cx
403  dcn = dcna / cnb - cna * dcnb / (cnb * cnb)
404  h1  = nu * (cn - cc0 - 3._DP / 7._DP * cx) * fz3 * t * t * EXP(ee)
405  dh1 = - third * (h1 * (7._DP + 8._DP * ee) + fz3 * nu * t * t * &
406        EXP(ee) * dcn)
407  !
408  ddh1 = 2._DP * h1 * (1._DP + ee) * rho / grho
409  dh1zup =  (1._DP - zeta) * dfz * h1 * (1._DP + 2._DP * ee / fz)
410  dh1zdw = -(1._DP + zeta) * dfz * h1 * (1._DP + 2._DP * ee / fz)
411  !
412  sc     = rho * (h0 + h1)
413  v1c_up = h0 + h1 + dh0up + dh1 + dh0zup + dh1zup
414  v1c_dw = h0 + h1 + dh0up + dh1 + dh0zdw + dh1zdw
415  v2c    = ddh0 + ddh1
416  !
417  !
418  RETURN
419  !
420END SUBROUTINE ggac_spin
421!
422!
423!-------------------------------------------------------------------
424SUBROUTINE pbec_spin( rho, zeta, grho, iflag, sc, v1c_up, v1c_dw, v2c )                    !<GPU:DEVICE>
425  !-----------------------------------------------------------------
426  !! PBE correlation (without LDA part) - spin-polarized.
427  !
428  !! * iflag = 1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996);
429  !! * iflag = 2: J.P.Perdew et al., PRL 100, 136406 (2008).
430  !
431  USE kinds, ONLY : DP
432  !
433  IMPLICIT NONE
434  !
435  INTEGER, INTENT(IN) :: iflag        !<GPU:VALUE>
436  !! see main comments
437  REAL(DP), INTENT(IN) :: rho
438  !! the total charge density
439  REAL(DP), INTENT(IN) :: zeta
440  !! the magnetization
441  REAL(DP), INTENT(IN) :: grho
442  !! the gradient of the charge squared
443  REAL(DP), INTENT(OUT) :: sc
444  !! correlation energies
445  REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw
446  !! derivative of correlation wr. rho
447  REAL(DP), INTENT(OUT) :: v2c
448  !! derivatives of correlation wr. grho
449  !
450  ! ... local variables
451  !
452  REAL(DP) :: rs, ec, vc_up, vc_dn
453  !
454  REAL(DP), PARAMETER :: ga=0.031091_DP
455  REAL(DP) :: be(2)
456  DATA be / 0.06672455060314922_DP, 0.046_DP /
457  REAL(DP), PARAMETER :: third=1.D0/3.D0, pi34=0.6203504908994_DP
458  !                                         pi34=(3/4pi)^(1/3)
459  REAL(DP), PARAMETER :: xkf=1.919158292677513_DP, xks=1.128379167095513_DP
460  !                      xkf=(9 pi/4)^(1/3)      , xks=sqrt(4/pi)
461  REAL(DP) :: kf, ks, t, expe, af, y, xy, qy, s1, h0, ddh0
462  REAL(DP) :: fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, &
463              dh0zup, dh0zdw
464  !
465  !
466  rs = pi34 / rho**third
467  !
468  CALL pw_spin( rs, zeta, ec, vc_up, vc_dn )                                 !<GPU:pw_spin=>pw_spin_d>
469  !
470  kf = xkf / rs
471  ks = xks * SQRT(kf)
472  !
473  fz = 0.5_DP*( (1._DP+zeta)**(2._DP/3._DP) + (1._DP-zeta)**(2._DP/3._DP) )
474  fz2 = fz * fz
475  fz3 = fz2 * fz
476  fz4 = fz3 * fz
477  dfz = ( (1._DP+zeta)**(-1._DP/3._DP) - (1._DP - zeta)**(-1._DP/3._DP) ) &
478         / 3._DP
479  !
480  t  = SQRT(grho) / (2._DP * fz * ks * rho)
481  expe = EXP( - ec / (fz3 * ga) )
482  af   = be(iflag) / ga * (1._DP / (expe-1._DP) )
483  bfup = expe * (vc_up - ec) / fz3
484  bfdw = expe * (vc_dn - ec) / fz3
485  !
486  y  = af * t * t
487  xy = (1._DP + y) / (1._DP + y + y * y)
488  qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2
489  s1 = 1._DP + be(iflag) / ga * t * t * xy
490  !
491  h0 = fz3 * ga * LOG(s1)
492  !
493  dh0up = be(iflag) * t * t * fz3 / s1 * ( -7._DP/3._DP * xy - qy * &
494          (af * bfup / be(iflag)-7._DP/3._DP) )
495  !
496  dh0dw = be(iflag) * t * t * fz3 / s1 * ( -7._DP/3._DP * xy - qy * &
497          (af * bfdw / be(iflag)-7._DP/3._DP) )
498  !
499  dh0zup =   (3._DP * h0 / fz - be(iflag) * t * t * fz2 / s1 *  &
500             (2._DP * xy - qy * (3._DP * af * expe * ec / fz3 / &
501             be(iflag)+2._DP) ) ) * dfz * (1._DP - zeta)
502  !
503  dh0zdw = - (3._DP * h0 / fz - be(iflag) * t * t * fz2 / s1 *  &
504             (2._DP * xy - qy * (3._DP * af * expe * ec / fz3 / &
505             be(iflag)+2._DP) ) ) * dfz * (1._DP + zeta)
506  !
507  ddh0 = be(iflag) * fz / (2._DP * ks * ks * rho) * (xy - qy) / s1
508  !
509  sc     = rho * h0
510  v1c_up = h0 + dh0up + dh0zup
511  v1c_dw = h0 + dh0dw + dh0zdw
512  v2c    = ddh0
513  !
514  !
515  RETURN
516  !
517END SUBROUTINE pbec_spin
518!
519!
520!------------------------------------------------------------------------
521SUBROUTINE lsd_glyp( rho_in_up, rho_in_dw, grho_up, grho_dw, grho_ud, sc, v1c_up, v1c_dw, v2c_up, v2c_dw, v2c_ud )                     !<GPU:DEVICE>
522  !----------------------------------------------------------------------
523  !! Lee, Yang, Parr: gradient correction part.
524  !
525  USE kinds, ONLY: DP
526  !
527  IMPLICIT NONE
528  !
529  REAL(DP), INTENT(IN) :: rho_in_up, rho_in_dw
530  !! the total charge density
531  REAL(DP), INTENT(IN) :: grho_up, grho_dw
532  !! the gradient of the charge squared
533  REAL(DP), INTENT(IN) :: grho_ud
534  !! gradient off-diagonal term up-down
535  REAL(DP), INTENT(OUT) :: sc
536  !! correlation energy
537  REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw
538  !! derivative of correlation wr. rho
539  REAL(DP), INTENT(OUT) :: v2c_up, v2c_dw
540  !! derivatives of correlation wr. grho
541  REAL(DP), INTENT(OUT) :: v2c_ud
542  !! derivative of correlation wr. grho, off-diag. term
543  !
544  ! ... local variables
545  !
546  REAL(DP) :: ra, rb, rho, grhoaa, grhoab, grhobb
547  REAL(DP) :: rm3, dr, or, dor, der, dder
548  REAL(DP) :: dlaa, dlab, dlbb, dlaaa, dlaab, dlaba, dlabb, dlbba, dlbbb
549  REAL(DP), PARAMETER :: a=0.04918_DP, b=0.132_DP, c=0.2533_DP, d=0.349_DP
550  !
551  ra = rho_in_up
552  rb = rho_in_dw
553  rho = ra + rb
554  rm3 = rho**(-1._DP/3._DP)
555  !
556  dr = ( 1._DP + d*rm3 )
557  or = EXP(-c*rm3) / dr * rm3**11._DP
558  dor = -1._DP/3._DP * rm3**4 * or * (11._DP/rm3-c-d/dr)
559  !
560  der  = c*rm3 + d*rm3/dr
561  dder = 1._DP/3._DP * (d*d*rm3**5/dr/dr - der/rho)
562  !
563  dlaa = -a*b*or*( ra*rb/9._DP*(1._DP-3*der-(der-11._DP)*ra/rho)-rb*rb )
564  dlab = -a*b*or*( ra*rb/9._DP*(47._DP-7._DP*der)-4._DP/3._DP*rho*rho  )
565  dlbb = -a*b*or*( ra*rb/9._DP*(1._DP-3*der-(der-11._DP)*rb/rho)-ra*ra )
566  !
567  dlaaa = dor/or*dlaa-a*b*or*(rb/9._DP*(1._DP-3*der-(der-11._DP)*ra/rho)-     &
568          ra*rb/9._DP*((3._DP+ra/rho)*dder+(der-11._DP)*rb/rho/rho))
569  dlaab = dor/or*dlaa-a*b*or*(ra/9._DP*(1._DP-3._DP*der-(der-11._DP)*ra/rho)- &
570          ra*rb/9._DP*((3._DP+ra/rho)*dder-(der-11._DP)*ra/rho/rho)-2._DP*rb)
571  dlaba = dor/or*dlab-a*b*or*(rb/9._DP*(47._DP-7._DP*der)-7._DP/9._DP*ra*rb*dder- &
572          8._DP/3._DP*rho)
573  dlabb = dor/or*dlab-a*b*or*(ra/9._DP*(47._DP-7._DP*der)-7._DP/9._DP*ra*rb*dder- &
574          8._DP/3._DP*rho)
575  dlbba = dor/or*dlbb-a*b*or*(rb/9._DP*(1._DP-3._DP*der-(der-11._DP)*rb/rho)- &
576          ra*rb/9._DP*((3._DP+rb/rho)*dder-(der-11._DP)*rb/rho/rho)-2._DP*ra)
577  dlbbb = dor/or*dlbb-a*b*or*(ra/9._DP*(1._DP-3*der-(der-11._DP)*rb/rho)-     &
578          ra*rb/9._DP*((3._DP+rb/rho)*dder+(der-11._DP)*ra/rho/rho))
579  !
580  grhoaa = grho_up
581  grhobb = grho_dw
582  grhoab = grho_ud
583  !
584  sc     = dlaa *grhoaa + dlab *grhoab + dlbb *grhobb
585  v1c_up = dlaaa*grhoaa + dlaba*grhoab + dlbba*grhobb
586  v1c_dw = dlaab*grhoaa + dlabb*grhoab + dlbbb*grhobb
587  v2c_up = 2._DP*dlaa
588  v2c_dw = 2._DP*dlbb
589  v2c_ud = dlab
590  !
591  !
592  RETURN
593  !
594END SUBROUTINE lsd_glyp
595!
596!
597!---------------------------------------------------------------
598SUBROUTINE cpbe2d( rho, grho, sc, v1c, v2c )                    !<GPU:DEVICE>
599  !---------------------------------------------------------------
600  !! 2D correction (last term of Eq. 5, PRL 108, 126402 (2012))
601  !
602  USE kinds,      ONLY: DP
603  !
604  IMPLICIT NONE
605  !
606  REAL(DP), INTENT(IN)  :: rho, grho
607  REAL(DP), INTENT(OUT) :: sc, v1c, v2c
608  !
609  ! ... local variables
610  !
611  REAL(8), PARAMETER :: pi=3.14159265358979323846d0
612  REAL(DP), PARAMETER :: ex1=0.333333333333333333_DP, ex2=1.166666666666667_DP
613  REAL(DP), PARAMETER :: ex3=ex2+1.0_DP
614  REAL(DP) :: fac1, fac2, zeta, phi, gr, rs, drsdn, akf, aks, t, dtdn, dtdgr
615  REAL(DP) :: p, a, g, alpha1, beta1,beta2,beta3,beta4, dgdrs, epsc, depscdrs
616  REAL(DP) :: c, gamma1, beta, aa, cg, adddepsc, h, dhdaa, dhdt, dhdrs
617  REAL(DP) :: epscpbe, depscpbedrs, depscpbedt, a0,a1,a2, b0,b1,b2, c0,c1,c2
618  REAL(DP) :: e0,e1,e2, f0,f1,f2, g0,g1,g2, h0,h1,h2, d0,d1,d2, ff, dffdt
619  REAL(DP) :: rs3d, rs2d, drs2ddrs3d, eps2d, deps2ddrs2, depsGGAdrs, depsGGAdt
620  REAL(DP) :: drs2ddt, rs2, ec, decdn, decdgr, daadepsc
621  !
622  fac1 = (3.d0*pi*pi)**ex1
623  fac2 = SQRT(4.d0*fac1/pi)
624  !
625  zeta = 0.d0
626  phi  = 1.d0
627  !
628  gr = SQRT(grho)
629  !
630  rs = (3.d0/4.d0/pi/rho)**ex1
631  drsdn = -DBLE(3 ** (0.1D1 / 0.3D1)) * DBLE(2 ** (0.1D1 / 0.3D1)) * &
632  0.3141592654D1 ** (-0.1D1 / 0.3D1) * (0.1D1 / rho) ** (-0.2D1 /    &
633  0.3D1) / rho ** 2 / 0.6D1
634  !
635  akf  = (3.d0*pi*pi*rho)**(1.d0/3.d0)
636  aks  = DSQRT(4.d0*akf/pi)
637  t    = gr/2.d0 / phi / aks / rho
638  dtdn = -7.d0/6.d0 * gr/2.d0 / phi/DSQRT(4.d0/pi)/   &
639         ((3.d0*pi*pi)**(1.d0/6.d0)) / (rho**(13.d0/6.d0))
640  dtdgr = 1.d0/2.d0/phi/aks/rho
641  !
642  ! for the LDA correlation
643  p = 1.d0
644  A = 0.0310906908696548950_DP
645  alpha1 = 0.21370d0
646  beta1  = 7.5957d0
647  beta2  = 3.5876d0
648  beta3  = 1.6382d0
649  beta4  = 0.49294d0
650  G = -0.2D1 * A * DBLE(1 + alpha1 * rs) * LOG(0.1D1 + 0.1D1 / A / ( &
651  beta1 * SQRT(DBLE(rs)) + DBLE(beta2 * rs) + DBLE(beta3 * rs ** (   &
652  0.3D1 / 0.2D1)) + DBLE(beta4 * rs ** (p + 1))) / 0.2D1)
653  !
654  dGdrs = -0.2D1 * A * alpha1 * LOG(0.1D1 + 0.1D1 / A / (beta1 * SQRT(rs) &
655   + beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs **           &
656  (p + 1)) / 0.2D1) + (0.1D1 + alpha1 * rs) / (beta1 * SQRT(rs) +         &
657  beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1))     &
658  ** 2 * (beta1 * rs ** (-0.1D1 / 0.2D1) / 0.2D1 + beta2 + 0.3D1 /        &
659  0.2D1 * beta3 * SQRT(rs) + beta4 * rs ** (p + 1) * DBLE(p + 1) /        &
660  rs) / (0.1D1 + 0.1D1 / A / (beta1 * SQRT(rs) + beta2 * rs + beta3 *     &
661  rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1)) / 0.2D1)
662  !
663  epsc = G
664  depscdrs = dGdrs
665  !
666  ! PBE
667  c = 1.d0
668  gamma1 = 0.0310906908696548950_dp
669  beta = 0.06672455060314922_dp
670  !
671  AA = beta / gamma1 / (EXP(-epsc / gamma1 / phi ** 3) - 0.1D1)
672  cg = beta / gamma1 ** 2 / (EXP(-epsc/ gamma1 / phi ** 3) - 0.1D1) &
673       ** 2 / phi ** 3 * EXP(-epsc / gamma1 / phi ** 3)
674  dAAdepsc = cg
675  !
676  IF (t <= 10.d0) THEN
677     !
678     H = DBLE(gamma1) * phi ** 3 * LOG(DBLE(1 + beta / gamma1 * t ** 2   &
679         * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4)))
680     !
681     dHdAA = gamma1 * phi ** 3 * (beta / gamma1 * t ** 4 / (1 + c * AA   &
682      * t ** 2 + AA ** 2 * t ** 4) - beta / gamma1 * t ** 2 * (1 + AA *  &
683      t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) ** 2 * (c * t **&
684      2 + 2 * AA * t ** 4)) / (1 + beta / gamma1 * t ** 2 * (1 + AA *    &
685      t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4))
686     !
687     dHdt = gamma1 * phi ** 3 * (2 * beta / gamma1 * t * (1 + AA * t **  &
688      2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) + 2 * beta / gamma1  &
689      * t ** 3 * AA / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) - beta /  &
690      gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** &
691      2 * t ** 4) ** 2 * (2 * c * AA * t + 4 * AA ** 2 * t ** 3)) / (1   &
692      + beta / gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t **  &
693      2 + AA ** 2 * t ** 4))
694     !
695  ELSE
696     !
697     H = gamma1 * (phi**3) * DLOG(1.d0+(beta/gamma1)*(1.d0/AA))
698     !
699     dHdAA = gamma1 * (phi**3) * 1.d0/(1.d0+(beta/gamma1)*(1.d0/AA))* &
700             (beta/gamma1) * (-1.d0/AA/AA)
701     !
702     dHdt = 0.d0
703     !
704  ENDIF
705  !
706  dHdrs = dHdAA*dAAdepsc*depscdrs
707  !
708  epscPBE     = epsc+H
709  depscPBEdrs = depscdrs+dHdrs
710  depscPBEdt  = dHdt
711  !
712  ! START THE 2D CORRECTION
713  !
714  beta = 1.3386d0
715  a0 = -0.1925d0
716  a1 =  0.117331d0
717  a2 =  0.0234188d0
718  b0 =  0.0863136d0
719  b1 = -0.03394d0
720  b2 = -0.037093d0
721  c0 =  0.057234d0
722  c1 = -0.00766765d0
723  c2 =  0.0163618d0
724  e0 =  1.0022d0
725  e1 =  0.4133d0
726  e2 =  1.424301d0
727  f0 = -0.02069d0
728  f1 =  0.d0
729  f2 =  0.d0
730  g0 =  0.340d0
731  g1 =  0.0668467d0
732  g2 =  0.d0
733  h0 =  0.01747d0
734  h1 =  0.0007799d0
735  h2 =  1.163099d0
736  d0 = -a0*h0
737  d1 = -a1*h1
738  d2 = -a2*h2
739  !
740  ff = t ** 4 * (1 + t ** 2) / (1000000 + t ** 6)
741  dffdt = 4 * t ** 3 * (1 + t ** 2) / (1000000 + t ** 6) + 2 * t ** &
742  5 / (1000000 + t ** 6) - 6 * t ** 9 * (1 + t ** 2) / (1000000 + t &
743  ** 6) ** 2
744!
745  rs3d=rs
746  rs2d = 0.4552100000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** ( &
747  0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** (  &
748  0.5D1 / 0.4D1) * SQRT(t)
749
750  cg = 0.5690125000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** (       &
751  0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** (0.1D1 &
752   / 0.4D1) * SQRT(t)
753  drs2ddrs3d=cg
754
755  cg = 0.2276050000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** (       &
756  0.5D1 / 0.12D2)) * DBLE((1 / pi) ** (-0.5D1 / 0.12D2)) * DBLE(rs3d ** &
757   (0.5D1 / 0.4D1)) * DBLE(t ** (-0.1D1 / 0.2D1))
758  drs2ddt=cg
759  rs2=rs2d
760  !
761  eps2d = (EXP(-beta * rs2) - 0.1D1) * (-0.2D1 / 0.3D1 * SQRT(0.2D1)   &
762   * DBLE((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 /       &
763  0.2D1)) / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE(  &
764  zeta ** 2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi / &
765   rs2) + a0 + (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 ** 3) * LOG(0.1D1  &
766   + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 **     &
767  2 + h0 * rs2 ** 3)) + (a1 + (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 **  &
768   3) * LOG(0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1)    &
769   + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * DBLE(zeta ** 2) + (a2 + (b2    &
770  * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) * LOG(0.1D1 + 0.1D1 / (e2 *   &
771   rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3   &
772  ))) * DBLE(zeta ** 4)
773  !
774  cg = -beta * EXP(-beta * rs2) * (-0.2D1 / 0.3D1 * SQRT(0.2D1) *      &
775  DBLE((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1))  &
776  / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE(zeta **   &
777   2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi / rs2)   &
778   + (EXP(-beta * rs2) - 0.1D1) * (0.2D1 / 0.3D1 * SQRT(0.2D1) * DBLE  &
779  ((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1)) /    &
780   pi / rs2 ** 2 - 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE(zeta  &
781   ** 2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi /     &
782  rs2 ** 2) + (b0 + 0.2D1 * c0 * rs2 + 0.3D1 * d0 * rs2 ** 2) * LOG(   &
783  0.1D1 + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2   &
784   ** 2 + h0 * rs2 ** 3)) - (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 **    &
785  3) / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0    &
786  * rs2 ** 3) ** 2 * (e0 + 0.3D1 / 0.2D1 * f0 * SQRT(rs2) + 0.2D1 *    &
787  g0 * rs2 + 0.3D1 * h0 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e0 * rs2 + f0  &
788   * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0 * rs2 ** 3)) + ((     &
789  b1 + 0.2D1 * c1 * rs2 + 0.3D1 * d1 * rs2 ** 2) * LOG(0.1D1 + 0.1D1   &
790  / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 *     &
791  rs2 ** 3)) - (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 ** 3) / (e1 * rs2  &
792   + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3) **   &
793   2 * (e1 + 0.3D1 / 0.2D1 * f1 * SQRT(rs2) + 0.2D1 * g1 * rs2 +       &
794  0.3D1 * h1 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** (  &
795  0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * DBLE(zeta ** 2)  &
796  + ((b2 + 0.2D1 * c2 * rs2 + 0.3D1 * d2 * rs2 ** 2) * LOG(0.1D1 +     &
797  0.1D1 / (e2 * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 &
798   * rs2 ** 3)) - (b2 * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) / (e2     &
799  * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 **    &
800  3) ** 2 * (e2 + 0.3D1 / 0.2D1 * f2 * SQRT(rs2) + 0.2D1 * g2 * rs2    &
801  + 0.3D1 * h2 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e2 * rs2 + f2 * rs2 **  &
802   (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3))) * DBLE(zeta **   &
803   4)
804  !
805  deps2ddrs2=cg
806  !
807  ! GGA-2D
808  !
809  depsGGAdrs = ff*(-depscPBEdrs+deps2ddrs2*drs2ddrs3d)
810  depsGGAdt  = dffdt*(-epscPBE+eps2d)+ff*      &
811               (-depscPBEdt+deps2ddrs2*drs2ddt)
812  !
813  ec = rho*(ff*(-epscPBE+eps2d))
814  !
815  decdn = ff*(-epscPBE+eps2d)+rho*depsGGAdrs*drsdn+ &
816          rho*depsGGAdt*dtdn
817  !
818  decdgr = rho*depsGGAdt*dtdgr
819  !
820  sc = ec
821  v1c = decdn
822  v2c = decdgr/gr
823  !
824  RETURN
825  !
826END SUBROUTINE cpbe2d
827! !
828END MODULE
829
830