44 #include "../base/base_uses.f90"
52 REAL(KIND=
dp),
PARAMETER :: f13 = 1.0_dp/3.0_dp, &
57 REAL(KIND=
dp) :: cx, flda, flsd, sfac, t13
58 REAL(KIND=
dp) :: fact, tact
59 REAL(KIND=
dp) :: eps_rho
60 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_exchange_gga'
77 SUBROUTINE xgga_info(functional, lsd, reference, shortform, needs, max_deriv)
78 INTEGER,
INTENT(in) :: functional
79 LOGICAL,
INTENT(in) :: lsd
80 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
81 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
82 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
84 IF (
PRESENT(reference))
THEN
85 SELECT CASE (functional)
87 reference =
"A. Becke, Phys. Rev. A 38, 3098 (1988)"
89 reference =
"J.P. Perdew and Y. Wang, Phys. Rev. B, 33, 8800 (1986)"
91 reference =
"J.P. Perdew et al., Phys. Rev. B, 46, 6671 (1992)"
93 reference =
"J.P. Perdew, K. Burke, M Ernzerhof, Phys. Rev. Lett, 77, 3865 (1996)"
95 reference =
"Y. Zang et al., PRL, 80, 890 (1998) (Revised PBEX)"
97 reference =
"Wee-Meng Hoe, A.J. Cohen, N.C. Handy, CPL, 341, 319 (2001)"
99 reference =
"E. Engel and S.H. Vosko, Phys. Rev. B, 47, 13164 (1993)"
101 cpabort(
"Invalid functional requested ("//cp_to_string(functional)//
")")
104 IF (len_trim(reference) + 6 < len(reference))
THEN
105 reference(len_trim(reference):len_trim(reference) + 6) =
' {LDA}'
109 IF (
PRESENT(shortform))
THEN
110 SELECT CASE (functional)
112 shortform =
"Becke 1988 Exchange Functional"
114 shortform =
"Perdew-Wang 1986 Functional (exchange energy)"
116 shortform =
"Perdew-Wang 1991 Functional (exchange energy)"
118 shortform =
"PBE exchange energy functional"
120 shortform =
"Revised PBEX by Zang et al."
122 shortform =
"OPTX exchange energy functional"
124 shortform =
"Engel-Vosko exchange energy from virial relation"
126 cpabort(
"Invalid functional requested ("//cp_to_string(functional)//
")")
129 IF (len_trim(shortform) + 6 < len(shortform))
THEN
130 shortform(len_trim(shortform):len_trim(shortform) + 6) =
' {LDA}'
134 IF (
PRESENT(needs))
THEN
136 needs%rho_spin = .true.
137 needs%rho_spin_1_3 = .true.
138 needs%norm_drho_spin = .true.
141 needs%rho_1_3 = .true.
142 needs%norm_drho = .true.
145 IF (
PRESENT(max_deriv)) max_deriv = 3
159 SUBROUTINE xgga_eval(functional, lsd, rho_set, deriv_set, order)
160 INTEGER,
INTENT(in) :: functional
161 LOGICAL,
INTENT(in) :: lsd
162 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
163 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
164 INTEGER,
INTENT(in) :: order
166 CHARACTER(len=*),
PARAMETER :: routinen =
'xgga_eval'
168 INTEGER :: handle, ispin, m, npoints, nspin
169 INTEGER,
DIMENSION(2) :: norm_drho_spin_name, rho_spin_name
170 INTEGER,
DIMENSION(2, 3) :: bo
171 REAL(kind=
dp) :: drho_cutoff, rho_cutoff
172 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: s
173 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fs
174 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
175 e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
177 TYPE(cp_3d_r_cp_type),
DIMENSION(2) :: norm_drho, rho, rho_1_3
178 TYPE(xc_derivative_type),
POINTER :: deriv
180 CALL timeset(routinen, handle)
181 NULLIFY (e_0, e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_rho_ndrho_ndrho, &
182 e_rho_ndrho, e_rho_rho_ndrho, e_rho, e_rho_rho, e_rho_rho_rho)
184 NULLIFY (norm_drho(ispin)%array, rho(ispin)%array, rho_1_3(ispin)%array)
189 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
190 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
191 norm_drhob=norm_drho(2)%array, rho_cutoff=rho_cutoff, &
192 drho_cutoff=drho_cutoff, local_bounds=bo)
197 CALL xc_rho_set_get(rho_set, rho=rho(1)%array, rho_1_3=rho_1_3(1)%array, &
198 norm_drho=norm_drho(1)%array, local_bounds=bo, rho_cutoff=rho_cutoff, &
199 drho_cutoff=drho_cutoff)
204 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
206 CALL xgga_init(rho_cutoff)
208 ALLOCATE (s(npoints))
209 ALLOCATE (fs(npoints, m + 1))
220 norm_drho(ispin)%array, s)
223 SELECT CASE (functional)
225 CALL efactor_b88(s, fs, m)
227 CALL efactor_pw86(s, fs, m)
235 CALL efactor_pw91(s, fs, m)
240 CALL efactor_pbex(s, fs, m, 1)
241 IF (lsd) tact = 1.0_dp
244 CALL efactor_pbex(s, fs, m, 2)
245 IF (lsd) tact = 1.0_dp
247 CALL efactor_optx(s, fs, m)
250 CALL efactor_ev93(s, fs, m)
251 IF (lsd) tact = 1.0_dp
253 cpabort(
"Unsupported functional")
258 allocate_deriv=.true.)
261 CALL x_p_0(rho(ispin)%array, rho_1_3(ispin)%array, fs, e_0, &
264 IF (order >= 1 .OR. order == -1)
THEN
266 allocate_deriv=.true.)
269 allocate_deriv=.true.)
272 CALL x_p_1(rho(ispin)%array, &
273 rho_1_3(ispin)%array, s, fs, e_rho, e_ndrho, npoints)
275 IF (order >= 2 .OR. order == -2)
THEN
277 rho_spin_name(ispin)], allocate_deriv=.true.)
280 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
283 norm_drho_spin_name(ispin)], allocate_deriv=.true.)
286 CALL x_p_2(rho(ispin)%array, &
287 rho_1_3(ispin)%array, s, fs, e_rho_rho, e_rho_ndrho, &
288 e_ndrho_ndrho, npoints)
290 IF (order >= 3 .OR. order == -3)
THEN
292 rho_spin_name(ispin), rho_spin_name(ispin)], &
293 allocate_deriv=.true.)
296 rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
297 allocate_deriv=.true.)
300 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
301 allocate_deriv=.true.)
304 norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
305 allocate_deriv=.true.)
308 CALL x_p_3(rho(ispin)%array, &
309 rho_1_3(ispin)%array, s, fs, e_rho_rho_rho, &
310 e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
313 IF (order > 3 .OR. order < -3)
THEN
314 cpabort(
"derivatives bigger than 3 not implemented")
321 CALL timestop(handle)
328 SUBROUTINE xgga_init(cutoff)
330 REAL(kind=
dp),
INTENT(IN) :: cutoff
335 cx = -0.75_dp*(3.0_dp/
pi)**f13
340 sfac = 1.0_dp/(2.0_dp*(3.0_dp*
pi*
pi)**f13)
342 END SUBROUTINE xgga_init
352 SUBROUTINE x_p_0(rho, r13, fs, e_0, npoints)
354 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13
355 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: fs
356 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_0
357 INTEGER,
INTENT(in) :: npoints
367 IF (rho(ip) > eps_rho)
THEN
368 e_0(ip) = e_0(ip) + fact*r13(ip)*rho(ip)*fs(ip, 1)
387 SUBROUTINE x_p_1(rho, r13, s, fs, e_rho, e_ndrho, npoints)
389 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13, s
390 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: fs
391 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho, e_ndrho
392 INTEGER,
INTENT(in) :: npoints
395 REAL(kind=
dp) :: a0, a1, sx, sy
404 IF (rho(ip) > eps_rho)
THEN
406 a0 = fact*r13(ip)*rho(ip)
407 a1 = f43*fact*r13(ip)
408 sx = -f43*s(ip)/rho(ip)
409 sy = sfac*tact/(r13(ip)*rho(ip))
410 e_rho(ip) = e_rho(ip) + a1*fs(ip, 1) + a0*fs(ip, 2)*sx
411 e_ndrho(ip) = e_ndrho(ip) + a0*fs(ip, 2)*sy
432 SUBROUTINE x_p_2(rho, r13, s, fs, e_rho_rho, e_rho_ndrho, &
433 e_ndrho_ndrho, npoints)
435 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13, s
436 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: fs
437 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
438 INTEGER,
INTENT(in) :: npoints
441 REAL(kind=
dp) :: a0, a1, a2, sx, sxx, sxy, sy
450 IF (rho(ip) > eps_rho)
THEN
452 a0 = fact*r13(ip)*rho(ip)
453 a1 = f43*fact*r13(ip)
454 a2 = f13*f43*fact/(r13(ip)*r13(ip))
455 sx = -f43*s(ip)/rho(ip)
456 sy = sfac*tact/(r13(ip)*rho(ip))
457 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
458 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
459 e_rho_rho(ip) = e_rho_rho(ip) + a2*fs(ip, 1) + 2.0_dp*a1*fs(ip, 2)*sx + &
460 a0*fs(ip, 3)*sx*sx + a0*fs(ip, 2)*sxx
461 e_rho_ndrho(ip) = e_rho_ndrho(ip) &
462 + a1*fs(ip, 2)*sy + a0*fs(ip, 3)*sx*sy + a0*fs(ip, 2)*sxy
463 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + a0*fs(ip, 3)*sy*sy
485 SUBROUTINE x_p_3(rho, r13, s, fs, e_rho_rho_rho, e_rho_rho_ndrho, &
486 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
488 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, r13, s
489 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: fs
490 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
491 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
492 INTEGER,
INTENT(in) :: npoints
495 REAL(kind=
dp) :: a0, a1, a2, a3, sx, sxx, sxxx, sxxy, &
507 IF (rho(ip) > eps_rho)
THEN
509 a0 = fact*r13(ip)*rho(ip)
510 a1 = f43*fact*r13(ip)
511 a2 = f13*f43*fact/(r13(ip)*r13(ip))
512 a3 = -f23*f13*f43*fact/(r13(ip)*r13(ip)*rho(ip))
513 sx = -f43*s(ip)/rho(ip)
514 sy = sfac*tact/(r13(ip)*rho(ip))
515 sxx = 28.0_dp/9.0_dp*s(ip)/(rho(ip)*rho(ip))
516 sxy = -f43*sfac*tact/(r13(ip)*rho(ip)*rho(ip))
517 sxxx = -280.0_dp/27.0_dp*s(ip)/(rho(ip)*rho(ip)*rho(ip))
518 sxxy = 28.0_dp/9.0_dp*sfac*tact/(r13(ip)*rho(ip)*rho(ip)*rho(ip))
519 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
520 + a3*fs(ip, 1) + 3.0_dp*a2*fs(ip, 2)*sx + &
521 3.0_dp*a1*fs(ip, 3)*sx*sx + 3.0_dp*a1*fs(ip, 2)*sxx + &
522 a0*fs(ip, 4)*sx*sx*sx + 3.0_dp*a0*fs(ip, 3)*sx*sxx + &
524 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
525 + a2*fs(ip, 2)*sy + 2.0_dp*a1*fs(ip, 3)*sx*sy + &
526 2.0_dp*a1*fs(ip, 2)*sxy + a0*fs(ip, 4)*sx*sx*sy + &
527 2.0_dp*a0*fs(ip, 3)*sx*sxy + a0*fs(ip, 3)*sxx*sy + &
529 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
530 + a1*fs(ip, 3)*sy*sy + a0*fs(ip, 4)*sx*sy*sy + &
531 2.0_dp*a0*fs(ip, 3)*sxy*sy
532 e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) &
533 + a0*fs(ip, 4)*sy*sy*sy
550 SUBROUTINE efactor_b88(s, fs, m)
551 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
552 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: fs
553 INTEGER,
INTENT(IN) :: m
556 REAL(kind=
dp) :: as, asp, beta, bs, f0, p, q, sas, sbs, sbs3, t1, t10, t13, t15, t16, t19, &
557 t2, t22, t24, t25, t32, t34, t36, t39, t4, t40, t41, t44, t48, t49, t5, t6, t65, t8, t87, &
575 sbs = sqrt(x*x + 1.0_dp)
578 ys = 1.0_dp/(1.0_dp + q*sas)
581 fs(ip, 1) = 1.0_dp + p*x*x*ys
584 fs(ip, 1) = 1.0_dp + p*x*x*ys
585 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
588 sbs3 = 1.0_dp/(sbs*sbs*sbs)
589 fs(ip, 1) = 1.0_dp + p*x*x*ys
590 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
591 fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
592 - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
593 + 3.0_dp*q - sbs) - sbs))
596 sbs3 = 1.0_dp/(sbs*sbs*sbs)
597 fs(ip, 1) = 1.0_dp + p*x*x*ys
598 fs(ip, 2) = (2.0_dp*p*x*ys - p*q*x*x*asp*ys*ys)*f0
599 fs(ip, 3) = -f0*f0*p*ys**3*sbs3*(q*x*x*x*x*(q*sas + 5.0_dp &
600 - 2.0_dp*q*sbs) + 2.0_dp*(x*x*(q*q*sas &
601 + 3.0_dp*q - sbs) - sbs))
613 t19 = q*t6 + t1*t15*t16
623 t44 = 2*q*t15*t16 + t1*t36*t16 - t1*t39*t41
627 t87 = -6*p*t10*t19 + 12*t22*t24*t25 - 6*t22*t10*t44 - 6*t48/t49*t25*t19 + &
628 6*t48*t24*t19*t44 - t48*t10*(3*q*t36*t16 - 3*q*t39*t41 + 3*t1*(1/t65/t4* &
629 t2*x - t34*x)*t16 - 3*t1*t36*t41*t15 + 2*t1*t39*t15/t40/t5)
632 fs(ip, 4) = f0*f0*f0*fs(ip, 4)
635 cpabort(
"Illegal order")
641 END SUBROUTINE efactor_b88
648 SUBROUTINE efactor_pw86(s, fs, m)
649 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
650 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: fs
651 INTEGER,
INTENT(IN) :: m
654 REAL(kind=
dp) :: f15, p0, p1, p15, s2, s4, s6, t1, t10, &
655 t12, t13, t14, t19, t2, t25, t3, t8, t9
673 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
676 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
677 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
680 fs(ip, 2) = f15*p1*p15/p0
682 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
683 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
686 fs(ip, 2) = f15*p1*p15/p0
687 t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
689 fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
690 1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
692 p0 = 1.0_dp + t1*s2 + t2*s4 + t3*s6
693 p1 = s(ip)*(2.0_dp*t1 + 4.0_dp*t2*s2 + 6.0_dp*t3*s4)
696 fs(ip, 2) = f15*p1*p15/p0
697 t9 = p15**2; t10 = t9**2; t12 = t10**2; t13 = t12*t10*t9
699 fs(ip, 3) = -14.0_dp/225.0_dp/t13/p0*t25 + &
700 1.0_dp/t13*(2.0_dp*t1 + 12*t2*s2 + 30.0_dp*t3*s4)/15.0_dp
701 t8 = p0**2; t9 = p0**f15; t14 = p0/t9; t19 = s2*s(ip)
702 fs(ip, 4) = 406.0_dp/3375.0_dp/t14/t8*p1*p1*p1 - 14.0_dp/ &
703 75.0_dp/t14/p0*p1*(2*t1 + 12*t2*s2 + 30*t3*s4) + &
704 1/t14*(24*t2*s(ip) + 120*t3*t19)*f15
706 cpabort(
"Illegal order")
712 END SUBROUTINE efactor_pw86
719 SUBROUTINE efactor_ev93(s, fs, m)
720 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
721 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: fs
722 INTEGER,
INTENT(IN) :: m
725 REAL(kind=
dp) :: a1, a2, a3, b1, b2, b3, d0, d1, d2, d3, &
726 f0, f1, f2, n0, n1, n2, n3, s2, s4, &
751 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
752 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
755 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
756 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
757 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
758 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
761 fs(ip, 2) = (n1 - f0*d1)/d0*scale_s
763 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
764 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
765 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
766 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
767 n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
768 d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
772 fs(ip, 2) = f1*scale_s
773 fs(ip, 3) = (n2 - f0*d2 - 2._dp*f1*d1)/d0*scale_s*scale_s
775 n0 = 1._dp + a1*s2 + a2*s4 + a3*s6
776 d0 = 1._dp + b1*s2 + b2*s4 + b3*s6
777 n1 = ss*(2._dp*a1 + 4._dp*a2*s2 + 6._dp*a3*s4)
778 d1 = ss*(2._dp*b1 + 4._dp*b2*s2 + 6._dp*b3*s4)
779 n2 = 2._dp*a1 + 12._dp*a2*s2 + 30._dp*a3*s4
780 d2 = 2._dp*b1 + 12._dp*b2*s2 + 30._dp*b3*s4
781 n3 = ss*(24._dp*a2 + 120._dp*a3*s2)
782 d3 = ss*(24._dp*b2 + 120._dp*b3*s2)
785 f2 = (n2 - f0*d2 - 2._dp*f1*d1)/d0
787 fs(ip, 2) = f1*scale_s
788 fs(ip, 3) = f2*scale_s*scale_s
789 fs(ip, 4) = (n3 - f0*d3 - 3._dp*f2*d1 - 3._dp*f1*d2)/d0*scale_s*scale_s*scale_s
791 cpabort(
"Illegal order")
797 END SUBROUTINE efactor_ev93
804 SUBROUTINE efactor_optx(s, fs, m)
805 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
806 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: fs
807 INTEGER,
INTENT(IN) :: m
809 REAL(kind=
dp),
PARAMETER :: a1 = 1.05151_dp, a2 = 1.43169_dp, &
813 REAL(kind=
dp) :: a, b, f0, x, y
825 y = 1.0_dp/(1.0_dp + a)
828 fs(ip, 1) = a1 + b*a*a*y*y
830 fs(ip, 1) = a1 + b*a*a*y*y
831 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
833 fs(ip, 1) = a1 + b*a*a*y*y
834 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
835 fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
837 fs(ip, 1) = a1 + b*a*a*y*y
838 fs(ip, 2) = 4.0_dp*b*f0*a*gamma_bo*x*y*y*y
839 fs(ip, 3) = -12.0_dp*b*f0*f0*gamma_bo*a*(a - 1.0_dp)*y*y*y*y
840 fs(ip, 4) = 24.0_dp*b*f0*f0*f0*gamma_bo*gamma_bo*x* &
841 (1.0_dp - 5.0_dp*a + 2.0_dp*a*a)*y*y*y*y*y
843 cpabort(
"Illegal order")
849 END SUBROUTINE efactor_optx
857 SUBROUTINE efactor_pw91(s, fs, m)
858 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
859 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: fs
860 INTEGER,
INTENT(IN) :: m
863 REAL(
dp) :: t1, t10, t101, t106, t109, t111, t113, t119, t12, t123, t124, t13, t14, t15, &
864 t16, t17, t18, t19, t2, t20, t21, t22, t23, t25, t26, t27, t28, t29, t3, t30, t31, t33, &
865 t35, t37, t38, t39, t4, t40, t44, t47, t48, t5, t50, t51, t53, t55, t56, t57, t58, t59, &
866 t6, t60, t64, t65, t69, t7, t70, t71, t73, t77, t78, t8, t80, t82, t9, t90, t93, t94, &
868 REAL(kind=
dp) :: a1, a2, a3, a4, a5, b, o, x
890 fs(ip, 1) = (o + t10 + (a2 - a3*t12)*t4)/(o + t10 + a5*t17)
909 t17 = t10*(b + 1/t6*t2*x)/t7
915 t33 = o + t30 + a5*t31
918 (t9 + t17 + 2._dp*a3*a4*t19*t21 + 2._dp*t26*x)/ &
919 t33 - (o + t30 + t26*t3)/t38*(t9 + t17 + 4._dp*a5*t19)
938 t15 = 2._dp*a1*t9*t13
942 t25 = t16*(-o/t17/t5*t20*t2 + t7)*t13
953 t50 = o + t48 + a5*t39
960 t69 = t53 + t55 + 4._dp*a5*t56
961 t73 = o + t48 + t60*t2
964 (t15 + t25 - t30 + 10._dp*t31*t2*t33 - 4._dp*a3*t37*t39*t33 + &
965 2._dp*a2 - 2._dp*t44)/t50 - 2._dp* &
966 (t53 + t55 + 2._dp*t31*t56*t33 + 2._dp*t60*x)* &
967 t65*t69 + 2._dp*t73/t64/t50*t77 - t73*t65*(t15 + t25 - t30 + 12._dp*a5*t2)
981 t5 = sqrt(0.1e1_dp + t1*t2)
998 t40 = 3*t29*(1/t30/t5*t1*t9*t35 - t10*x)*t18
999 t44 = 3*t29*t14*t26*t22
1000 t50 = 2*t29*t23*t22/t25/t17
1009 t73 = 0.1e1_dp + t71 + a5*t64
1018 t101 = t96 + t98 + 4*a5*t35
1020 t109 = t96 + t98 + 2*t51*t59 + 2*t106*x
1023 t119 = t78 + t80 - t82 + 12*a5*t2
1024 t123 = 0.1e1_dp + t71 + t106*t2
1027 (t20 - t28 + t40 - t44 + t50 + 24*t51*x*t53 - 36._dp*t58*t59 + 8._dp*a3*t57*a4*t64* &
1028 x*t53)/t73 - 3._dp*(t78 + t80 - t82 + 10._dp*t51*t2*t53 - &
1029 4._dp*t58*t64*t53 + 2._dp*a2 - 2._dp*t90)*t94*t101 + &
1030 6._dp*t109*t111*t113 - 3._dp*t109*t94*t119 - 6*t123/t124*t113*t101 + &
1031 6._dp*t123*t111*t101*t119 - t123*t94*(t20 - t28 + t40 - t44 + t50 + 24._dp*a5*x)
1038 END SUBROUTINE efactor_pw91
1047 SUBROUTINE efactor_pbex(s, fs, m, pset)
1048 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: s
1049 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: fs
1050 INTEGER,
INTENT(IN) :: m, pset
1052 REAL(kind=
dp),
PARAMETER :: kappa1 = 0.804_dp, kappa2 = 1.245_dp, &
1053 mu = 0.2195149727645171_dp
1056 REAL(kind=
dp) :: f0, mk, x, x2, y
1058 IF (pset == 1) mk = mu/kappa1
1059 IF (pset == 2) mk = mu/kappa2
1070 y = 1.0_dp/(1.0_dp + mk*x2)
1073 fs(ip, 1) = 1.0_dp + mu*x2*y
1075 fs(ip, 1) = 1.0_dp + mu*x2*y
1076 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1078 fs(ip, 1) = 1.0_dp + mu*x2*y
1079 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1080 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1082 fs(ip, 1) = 1.0_dp + mu*x2*y
1083 fs(ip, 2) = 2.0_dp*mu*x*y*y*f0
1084 fs(ip, 3) = -2.0_dp*mu*(3.0_dp*mk*x2 - 1.0_dp)*y*y*y*f0*f0
1085 fs(ip, 4) = 24.0_dp*mu*mk*x*(mk*x2 - 1.0_dp)*y*y*y*y*f0*f0*f0
1087 cpabort(
"Illegal order")
1093 END SUBROUTINE efactor_pbex
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
various routines to log and control the output. The idea is that decisions about where to log should ...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_norm_drho
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Calculate several different exchange energy functionals with a GGA form.
subroutine, public xgga_eval(functional, lsd, rho_set, deriv_set, order)
evaluates different exchange gga
subroutine, public xgga_info(functional, lsd, reference, shortform, needs, max_deriv)
return various information on the xgga functionals
Utility routines for the functional calculations.
subroutine, public calc_wave_vector(tag, rho, grho, s)
...
subroutine, public set_util(cutoff)
...
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set