34 #include "../base/base_uses.f90"
39 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
40 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_lyp'
41 REAL(kind=
dp),
PARAMETER,
PRIVATE :: a = 0.04918_dp, b = 0.132_dp, &
42 c = 0.2533_dp, d = 0.349_dp
59 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
60 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
61 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
63 IF (
PRESENT(reference))
THEN
64 reference =
"C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}"
66 IF (
PRESENT(shortform))
THEN
67 shortform =
"Lee-Yang-Parr correlation energy functional (LDA)"
69 IF (
PRESENT(needs))
THEN
71 needs%rho_1_3 = .true.
72 needs%norm_drho = .true.
74 IF (
PRESENT(max_deriv)) max_deriv = 3
90 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
91 TYPE(xc_rho_cflags_type),
INTENT(inout),
OPTIONAL :: needs
92 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
94 IF (
PRESENT(reference))
THEN
95 reference =
"C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}"
97 IF (
PRESENT(shortform))
THEN
98 shortform =
"Lee-Yang-Parr correlation energy functional (LSD)"
100 IF (
PRESENT(needs))
THEN
101 needs%rho_spin = .true.
102 needs%norm_drho_spin = .true.
103 needs%norm_drho = .true.
105 IF (
PRESENT(max_deriv)) max_deriv = 2
124 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
125 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
126 INTEGER,
INTENT(in) :: grad_deriv
127 TYPE(section_vals_type),
POINTER :: lyp_params
129 CHARACTER(len=*),
PARAMETER :: routinen =
'lyp_lda_eval'
131 INTEGER :: handle, npoints
132 INTEGER,
DIMENSION(2, 3) :: bo
133 REAL(kind=
dp) :: epsilon_norm_drho, epsilon_rho, sc
134 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: dummy, e_0, e_ndrho, &
135 e_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, &
136 e_rho_rho_rho, norm_drho, rho, rho_1_3
137 TYPE(xc_derivative_type),
POINTER :: deriv
139 CALL timeset(routinen, handle)
145 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
146 drho_cutoff=epsilon_norm_drho)
147 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
156 e_ndrho_ndrho => dummy
157 e_rho_rho_rho => dummy
158 e_ndrho_rho_rho => dummy
159 e_ndrho_ndrho_rho => dummy
161 IF (grad_deriv >= 0)
THEN
163 allocate_deriv=.true.)
166 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
168 allocate_deriv=.true.)
171 allocate_deriv=.true.)
174 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
176 allocate_deriv=.true.)
179 allocate_deriv=.true.)
185 IF (grad_deriv >= 3 .OR. grad_deriv == -3)
THEN
187 allocate_deriv=.true.)
200 IF (grad_deriv > 3 .OR. grad_deriv < -3)
THEN
201 cpabort(
"derivatives bigger than 3 not implemented")
211 CALL lyp_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
212 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
213 e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
214 e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
215 e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
216 grad_deriv=grad_deriv, &
217 npoints=npoints, epsilon_rho=epsilon_rho, sc=sc)
221 CALL timestop(handle)
249 SUBROUTINE lyp_lda_calc(rho, rho_1_3, norm_drho, &
250 e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
251 e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
252 grad_deriv, npoints, epsilon_rho, sc)
253 INTEGER,
INTENT(in) :: npoints, grad_deriv
254 REAL(kind=
dp),
DIMENSION(1:npoints),
INTENT(inout) :: e_ndrho_ndrho_rho, e_ndrho_rho_rho, &
255 e_rho_rho_rho, e_ndrho_ndrho, &
256 e_ndrho_rho, e_rho_rho, e_ndrho, &
258 REAL(kind=
dp),
DIMENSION(1:npoints),
INTENT(in) :: norm_drho, rho_1_3, rho
259 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, sc
262 REAL(kind=
dp) :: cf, epsilon_rho13, my_ndrho, my_rho, my_rho_1_3, t1, t102, t103, t105, &
263 t106, t11, t112, t123, t124, t127, t128, t13, t133, t134, t14, t161, t165, t166, t173, &
264 t176, t184, t185, t189, t192, t193, t196, t199, t2, t200, t201, t202, t203, t215, t22, &
265 t227, t228, t231, t245, t246, t248, t26, t265, t278, t279, t3, t332, t37, t38, t39, t4, &
266 t40, t41, t44, t45, t48, t5, t52, t6, t61, t62, t69, t7, t70, t77, t78, t80, t87, t88, &
267 t89, t93, t94, t95, t98, t99
269 epsilon_rho13 = epsilon_rho**(1.0_dp/3.0_dp)
270 cf = 0.3_dp*(3._dp*
pi*
pi)**(2._dp/3._dp)
271 SELECT CASE (grad_deriv)
276 IF (my_rho > epsilon_rho)
THEN
277 my_rho_1_3 = rho_1_3(ii)
278 my_ndrho = norm_drho(ii)
286 t11 = 0.1e1_dp/my_rho_1_3
291 t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
292 t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
293 & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
294 &*t26*c + 0.7e1_dp*t14*t22*c*d
301 + (t4*t41/0.72e2_dp)*sc
310 t77 = 0.1e1_dp/my_rho
316 t93 = my_rho_1_3*my_rho
319 t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
320 & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
321 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
322 &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
323 t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
324 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
325 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
330 t105 = 0.1e1_dp/t39/t38
333 e_rho(ii) = e_rho(ii) &
334 - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
335 & + t103*t106/0.108e3_dp)*sc
337 t112 = my_rho_1_3*my_ndrho
338 t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
339 &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
343 e_ndrho(ii) = e_ndrho(ii) &
344 + (t4*t124/0.72e2_dp)*sc
352 IF (my_rho > epsilon_rho)
THEN
353 my_rho_1_3 = rho_1_3(ii)
354 my_ndrho = norm_drho(ii)
362 t11 = 0.1e1_dp/my_rho_1_3
367 t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
368 t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
369 & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
370 &*t26*c + 0.7e1_dp*t14*t22*c*d
376 IF (grad_deriv >= 0)
THEN
378 + (t4*t41/0.72e2_dp)*sc
389 t77 = 0.1e1_dp/my_rho
395 t93 = my_rho_1_3*my_rho
398 t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
399 & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
400 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
401 &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
402 t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
403 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
404 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
409 t105 = 0.1e1_dp/t39/t38
411 t112 = my_rho_1_3*my_ndrho
412 t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
413 &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
417 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
418 e_rho(ii) = e_rho(ii) &
419 - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
420 & + t103*t106/0.108e3_dp)*sc
421 e_ndrho(ii) = e_ndrho(ii) &
422 + (t4*t124/0.72e2_dp)*sc
425 t127 = 0.1e1_dp/t1/t6
437 t192 = -0.560e3_dp*t93 - 0.432e3_dp*my_rho*d - 0.128e3_dp&
438 & *t52*my_rho*t13*cf - 0.8e1_dp*t88*t1*t13* &
439 cf - 0.560e3_dp*t14*t93*cf - 0.112e3_dp*t52*t1&
440 & *t62 - 0.8e1_dp*t88*my_rho_1_3*t62 - 0.432e3_dp* &
441 t14*my_rho*cf*d - 0.14e2_dp/0.9e1_dp*t52* &
442 t161*t22 - 0.11e2_dp/0.9e1_dp*t88*t166*t22 - &
443 0.2e1_dp/0.3e1_dp*t14*t94*t22 - 0.20e2_dp/ &
444 0.9e1_dp*t173*t80 - 0.2e1_dp*t176*t80 - &
445 0.20e2_dp/0.9e1_dp*t14*t3*t22*d + 0.7e1_dp/ &
446 0.9e1_dp*t184*t185*t22 + 0.7e1_dp/0.9e1_dp* &
455 t215 = t13*my_ndrho*d
456 t227 = 0.20e2_dp/0.3e1_dp*t52*t70*my_ndrho + 0.4e1_dp* &
457 t14*t11*my_ndrho + 0.20e2_dp/0.3e1_dp*t78*t215&
458 & + 0.20e2_dp/0.3e1_dp*t14*t69*my_ndrho*d + &
459 0.14e2_dp/0.3e1_dp*t88*t89*my_ndrho + 0.14e2_dp &
463 t245 = 0.6e1_dp*t14*t1 + 0.20e2_dp*t14*my_rho_1_3*d + &
464 0.14e2_dp*t14*my_rho_1_3*c + 0.14e2_dp*t14*c* &
468 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
469 e_rho_rho(ii) = e_rho_rho(ii) &
470 + (0.5e1_dp/0.81e2_dp*t128*t41 - 0.5e1_dp/0.108e3_dp&
471 & *t45*t99 + t134*t106/0.27e2_dp + t4*t193/ &
472 0.72e2_dp - t103*t196/0.54e2_dp + t200*t203/ &
474 e_ndrho_rho(ii) = e_ndrho_rho(ii) &
475 - (0.5e1_dp/0.216e3_dp*t45*t124 - t4*t228/ &
476 0.72e2_dp + t103*t231/0.108e3_dp)*sc
477 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
478 + (t4*t246/0.72e2_dp)*sc
485 t332 = -0.432e3_dp*d - 0.2240e4_dp/0.3e1_dp*my_rho_1_3 - &
486 0.74e2_dp/0.27e2_dp*t184*t127*t80 + 0.100e3_dp/ &
487 0.27e2_dp*t14*t44*t22*d + 0.7e1_dp/0.27e2_dp* &
488 t279*t127*t13*t22 - 0.8e1_dp/0.3e1_dp*t184* &
489 t70*cf - 0.2240e4_dp/0.3e1_dp*t14*my_rho_1_3* &
490 cf - 0.48e2_dp*t88*t11*t13*cf - 0.944e3_dp/ &
491 0.3e1_dp*t52*t61 - 0.40e2_dp*t88*t69*t62 - &
492 0.656e3_dp/0.3e1_dp*t52*t11*t62 + 0.7e1_dp/ &
493 0.27e2_dp*t279*t265*t80 + 0.64e2_dp/0.27e2_dp* &
494 t52*t44*t13*t22 - 0.8e1_dp/0.3e1_dp*t184*t77&
495 & *t62 - 0.432e3_dp*t14*cf*d + 0.52e2_dp/ &
496 0.27e2_dp*t88*t199*t13*t22 - 0.20e2_dp/ &
497 0.9e1_dp*t184*t133*t13*t22 + 0.8e1_dp/0.9e1_dp&
498 & *t14*t102*t22 + 0.100e3_dp/0.27e2_dp*t52*t199&
499 & *t80 + 0.106e3_dp/0.27e2_dp*t88*t133*t80
501 IF (grad_deriv >= 3 .OR. grad_deriv == -3)
THEN
502 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
503 - (0.55e2_dp/0.243e3_dp*a/t1/t248*t41 - 0.5e1_dp &
504 &/0.27e2_dp*t128*t99 + 0.40e2_dp/0.243e3_dp*a/ &
505 my_rho_1_3/t248*t106 + 0.5e1_dp/0.72e2_dp*t45* &
506 t193 - t134*t196/0.9e1_dp + 0.7e1_dp/0.108e3_dp* &
507 a*t265*t203 - t4*t332*t40/0.72e2_dp + t103* &
508 t192*t105/0.36e2_dp - t200*t98*t202/0.36e2_dp &
509 &+ t128*t37/t201/t38/0.81e2_dp)*sc
510 e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
511 + (0.5e1_dp/0.81e2_dp*t128*t124 - 0.5e1_dp/ &
512 0.108e3_dp*t45*t228 + t134*t231/0.27e2_dp + t4*&
513 & (-0.28e2_dp/0.9e1_dp*t52*t161*my_ndrho - &
514 0.22e2_dp/0.9e1_dp*t88*t166*my_ndrho - 0.4e1_dp &
515 &/0.3e1_dp*t14*t94*my_ndrho - 0.40e2_dp/0.9e1_dp &
516 &*t173*t215 - 0.4e1_dp*t176*t215 - 0.40e2_dp/ &
517 0.9e1_dp*t14*t3*my_ndrho*d + 0.14e2_dp/ &
518 0.9e1_dp*t184*t185*my_ndrho + 0.14e2_dp/0.9e1_dp&
519 & *t189*t215)*t40/0.72e2_dp - t103*t227*t105/ &
520 0.54e2_dp + t200*t123*t202/0.108e3_dp)*sc
521 e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
522 - (0.5e1_dp/0.216e3_dp*t45*t246 - t4*(0.20e2_dp/ &
523 0.3e1_dp*t52*t70 + 0.4e1_dp*t14*t11 + 0.20e2_dp &
524 &/0.3e1_dp*t52*t89*d + 0.20e2_dp/0.3e1_dp*t14* &
525 t69*d + 0.14e2_dp/0.3e1_dp*t88*t89 + 0.14e2_dp/ &
526 0.3e1_dp*t88*t94*t13*d)*t40/0.72e2_dp + t103&
527 & *t245*t105/0.108e3_dp)*sc
713 END SUBROUTINE lyp_lda_calc
727 TYPE(xc_rho_set_type),
INTENT(IN) :: rho_set
728 TYPE(xc_derivative_set_type),
INTENT(IN) :: deriv_set
729 INTEGER,
INTENT(in) :: grad_deriv
730 TYPE(section_vals_type),
POINTER :: lyp_params
732 CHARACTER(len=*),
PARAMETER :: routinen =
'lyp_lsd_eval'
734 INTEGER :: handle, npoints
735 INTEGER,
DIMENSION(2, 3) :: bo
736 REAL(kind=
dp) :: epsilon_drho, epsilon_rho, sc
737 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
738 e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, &
739 e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, &
740 norm_drhob, rhoa, rhob
741 TYPE(xc_derivative_type),
POINTER :: deriv
743 CALL timeset(routinen, handle)
750 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
751 norm_drhob=norm_drhob, norm_drho=norm_drho, &
752 rho_cutoff=epsilon_rho, &
753 drho_cutoff=epsilon_drho, local_bounds=bo)
754 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
776 IF (grad_deriv >= 0)
THEN
778 allocate_deriv=.true.)
781 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
783 allocate_deriv=.true.)
786 allocate_deriv=.true.)
789 allocate_deriv=.true.)
792 allocate_deriv=.true.)
795 allocate_deriv=.true.)
798 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
800 allocate_deriv=.true.)
803 allocate_deriv=.true.)
806 allocate_deriv=.true.)
809 allocate_deriv=.true.)
812 allocate_deriv=.true.)
815 allocate_deriv=.true.)
818 allocate_deriv=.true.)
821 allocate_deriv=.true.)
843 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
844 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
845 e_ndra_ra=e_ndra_ra, e_ndra_rb=e_ndra_rb, e_ndrb_ra&
846 &=e_ndrb_ra, e_ndrb_rb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, &
847 e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, &
848 e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, &
849 e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ndr_ra=e_ndr_ra, &
851 grad_deriv=grad_deriv, npoints=npoints, &
852 epsilon_rho=epsilon_rho, sc=sc)
856 CALL timestop(handle)
893 SUBROUTINE lyp_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
894 e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, &
895 e_ndra_ndra, e_ndrb_ndrb, e_ndr, &
896 e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb, &
897 grad_deriv, npoints, epsilon_rho, sc)
898 REAL(kind=
dp),
DIMENSION(*),
INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
900 REAL(kind=
dp),
DIMENSION(*),
INTENT(inout) :: e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, &
901 e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, &
902 e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb
903 INTEGER,
INTENT(in) :: grad_deriv, npoints
904 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, sc
907 REAL(kind=
dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rho, my_rhoa, my_rhob, t1, t100, &
908 t101, t102, t103, t104, t107, t108, t109, t112, t115, t118, t12, t120, t124, t129, t13, &
909 t130, t132, t135, t138, t14, t140, t142, t145, t15, t153, t155, t159, t16, t164, t168, &
910 t17, t171, t176, t18, t181, t182, t185, t189, t194, t195, t198, t2, t20, t202, t205, &
911 t206, t21, t210, t214, t215, t218, t22, t222, t228, t23, t232, t236, t238, t24, t243, &
912 t249, t25, t252, t253, t255, t257, t26, t260, t265, t268, t27, t270, t29, t3, t30, t304, &
913 t31, t310, t313, t316, t319, t322, t326, t329, t332, t334, t341, t35
914 REAL(kind=
dp) :: t354, t356, t360, t363, t37, t373, t376, t381, t39, t391, t4, t40, t408, &
915 t41, t415, t419, t422, t430, t44, t445, t449, t45, t452, t46, t467, t47, t48, t483, t487, &
916 t49, t490, t5, t50, t505, t515, t519, t52, t520, t54, t56, t57, t6, t61, t62, t64, t66, &
917 t7, t72, t75, t78, t8, t82, t85, t86, t87, t88, t90, t92, t94, t95, t98, t99
919 cf = 0.3_dp*(3._dp*
pi*
pi)**(2._dp/3._dp)
924 my_rhoa = max(rhoa(ii), 0.0_dp)
925 my_rhob = max(rhob(ii), 0.0_dp)
926 my_rho = my_rhoa + my_rhob
927 IF (my_rho > epsilon_rho)
THEN
928 my_ndrho = norm_drho(ii)
929 my_ndrhoa = norm_drhoa(ii)
930 my_ndrhob = norm_drhob(ii)
932 t1 = my_rho**(0.1e1_dp/0.3e1_dp)
947 t20 = 0.1e1_dp/t18/t17
949 t22 = 2**(0.1e1_dp/0.3e1_dp)
953 t26 = my_rhoa**(0.1e1_dp/0.3e1_dp)
956 t30 = my_rhob**(0.1e1_dp/0.3e1_dp)
958 t35 = 0.8e1_dp*t24*(t27*t25 + t31*t29)
960 t39 = 0.47e2_dp/0.18e2_dp - 0.7e1_dp/0.18e2_dp*t13 - &
961 0.7e1_dp/0.18e2_dp*t37
964 t44 = 0.5e1_dp/0.2e1_dp - t13/0.18e2_dp - t37/0.18e2_dp
969 t49 = t13 + t37 - 0.11e2_dp
972 t54 = t50*t45 + t52*t46
973 t56 = t49*t54/0.9e1_dp
974 t57 = t35 + t41 - t48 - t56
975 t61 = 0.2e1_dp/0.3e1_dp*t16
978 t66 = t7*t57 - 0.2e1_dp/0.3e1_dp*t16*t40 + t62*t46 + &
981 IF (grad_deriv >= 0 .AND. my_rho > epsilon_rho)
THEN
983 - (0.4e1_dp*t6*t7*t8 + t15*t21*t66)*sc
989 t78 = 0.64e2_dp/0.3e1_dp*t24*t72 - t75*t45/0.9e1_dp
990 t82 = my_rhob*t57 + t7*t78 - 0.2e1_dp*my_rhoa*t46
995 t90 = 0.1e1_dp/t1/t16
997 t94 = 0.4e1_dp/0.3e1_dp*t88*t92
999 t98 = 0.4e1_dp*t6*t7*t95
1003 t102 = 0.1e1_dp/t101
1006 t107 = t99*t103*t104/0.3e1_dp
1009 t112 = t15*t108*t109/0.3e1_dp
1011 t118 = 0.11e2_dp/0.3e1_dp*t15*t115*t66
1012 t120 = 0.1e1_dp/t1/my_rho
1014 t129 = c*t120 + d*t120*t5 - t124/t18/my_rho*t86
1015 t130 = 0.7e1_dp/0.54e2_dp*t129
1016 t132 = t129/0.54e2_dp
1017 t135 = -t129/0.3e1_dp
1020 t142 = -t138*t45 - t140*t46
1021 t145 = t130*t40 - t132*t47 - t135*t54/0.9e1_dp - t49* &
1023 t153 = t7*t145 - 0.4e1_dp/0.3e1_dp*my_rho*t40 + &
1024 0.4e1_dp/0.3e1_dp*my_rho*t46 + 0.4e1_dp/0.3e1_dp&
1028 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
1029 e_ra(ii) = e_ra(ii) &
1030 - (0.4e1_dp*t6*t52 + t15*t21*t82 + t94 - t98 + &
1031 t107 + t112 - t118 + t155)*sc
1035 t164 = 0.64e2_dp/0.3e1_dp*t24*t159 - t75*t46/0.9e1_dp
1036 t168 = my_rhoa*t57 + t7*t164 - 0.2e1_dp*my_rhob*t45
1038 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
1039 e_rb(ii) = e_rb(ii) &
1040 - (0.4e1_dp*t6*t50 + t15*t21*t168 + t94 - t98 + &
1041 t107 + t112 - t118 + t155)*sc
1045 t176 = 0.2e1_dp*t7*t171 - 0.4e1_dp/0.3e1_dp*t16*my_ndrho
1047 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
1048 e_ndr(ii) = e_ndr(ii) &
1054 t185 = -0.2e1_dp*t44*my_ndrhoa - 0.2e1_dp/0.9e1_dp*t181*t182
1055 t189 = t7*t185 + 0.2e1_dp*t64*my_ndrhoa
1057 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
1058 e_ndra(ii) = e_ndra(ii) &
1064 t198 = -0.2e1_dp*t44*my_ndrhob - 0.2e1_dp/0.9e1_dp*t194*t195
1065 t202 = t7*t198 + 0.2e1_dp*t62*my_ndrhob
1067 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
1068 e_ndrb(ii) = e_ndrb(ii) &
1074 t206 = 0.1e1_dp/t205
1075 t210 = 0.26e2_dp/0.9e1_dp*t99*t206*t14*t104
1076 t214 = 0.2e1_dp/0.3e1_dp*t99*t103*t5*t153
1078 t218 = 0.1e1_dp/t1/t205
1079 t222 = t12*t215*t218*t14*t104/0.9e1_dp
1080 t228 = 0.2e1_dp/0.9e1_dp*t12*c*t218*t14*t86*t109
1081 t232 = 0.26e2_dp/0.9e1_dp*t15*t86*t206*t109
1082 t236 = 0.2e1_dp/0.3e1_dp*t15*t108*t153*d
1083 t238 = 0.1e1_dp/t85/t4
1084 t243 = 0.2e1_dp/0.9e1_dp*t15*t238*t218*t66*t124
1085 t249 = 0.154e3_dp/0.9e1_dp*t15*t5/t18/t101*t66
1086 t252 = 0.22e2_dp/0.3e1_dp*t15*t115*t153
1091 t265 = t124/t18/t16*t86
1093 t270 = t124*d*t268*t238
1094 t304 = t15*t21*(t7*((-0.14e2_dp/0.81e2_dp*t257 - &
1095 0.14e2_dp/0.81e2_dp*t260 + 0.7e1_dp/0.27e2_dp* &
1096 t265 - 0.7e1_dp/0.81e2_dp*t270)*t40 - (-0.2e1_dp/ &
1097 0.81e2_dp*t257 - 0.2e1_dp/0.81e2_dp*t260 + t265/ &
1098 0.27e2_dp - t270/0.81e2_dp)*t47 - (0.4e1_dp/ &
1099 0.9e1_dp*t257 + 0.4e1_dp/0.9e1_dp*t260 - 0.2e1_dp &
1100 &/0.3e1_dp*t265 + 0.2e1_dp/0.9e1_dp*t270)*t54/ &
1101 0.9e1_dp - 0.2e1_dp/0.9e1_dp*t135*t142 - t49*&
1102 & (0.2e1_dp*my_rhoa*t268*t45 + 0.2e1_dp*my_rhob* &
1103 t268*t46)/0.9e1_dp) - 0.4e1_dp/0.3e1_dp*t40 + &
1104 0.4e1_dp/0.3e1_dp*t46 + 0.4e1_dp/0.3e1_dp*t45)
1105 t310 = 0.40e2_dp/0.9e1_dp*t88*my_rhob/t1/t17*d
1107 t316 = 0.8e1_dp/0.9e1_dp*a*t238*my_rhoa*t313*t124
1108 t319 = 0.8e1_dp*t6*t7*t268
1109 t322 = t99*t103*t5*t82
1110 t326 = t15*t108*t82*d
1114 t341 = t15*t21*(my_rhob*t145 + t7*(-t332*t45/ &
1115 0.9e1_dp + t334*t45/0.9e1_dp))
1117 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1118 e_ra_ra(ii) = e_ra_ra(ii) &
1119 + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
1120 t252 + 0.8e1_dp*t253 - 0.8e1_dp/0.3e1_dp*t255 - &
1121 t304 + t310 - t316 - t319 - 0.2e1_dp/0.3e1_dp*t322 - &
1122 0.2e1_dp/0.3e1_dp*t326 + 0.22e2_dp/0.3e1_dp*t329&
1123 & - 0.2e1_dp*t341 - t15*t21*(0.2e1_dp*my_rhob* &
1124 t78 + 0.320e3_dp/0.9e1_dp*t72*my_rhob*t24 - &
1128 t354 = t99*t103*t5*t168
1130 t360 = t87*my_rhoa*t90*d
1131 t363 = t15*t115*t168
1132 t373 = t15*t21*(my_rhoa*t145 + t7*(-t332*t46/ &
1133 0.9e1_dp + t334*t46/0.9e1_dp))
1134 t376 = t15*t108*t168*d
1136 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1137 e_rb_rb(ii) = e_rb_rb(ii) &
1138 + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
1139 t252 - t304 + t310 - t316 - t319 + 0.8e1_dp*t356 - &
1140 0.8e1_dp/0.3e1_dp*t360 + 0.22e2_dp/0.3e1_dp* &
1141 t363 - 0.2e1_dp*t373 - 0.2e1_dp/0.3e1_dp*t354 - &
1142 0.2e1_dp/0.3e1_dp*t376 - t15*t21*(0.2e1_dp* &
1143 my_rhoa*t164 + 0.320e3_dp/0.9e1_dp*my_rhoa* &
1144 t159*t24 - 0.2e1_dp*t45))*sc
1147 t381 = -t354/0.3e1_dp + 0.4e1_dp*t356 - 0.4e1_dp/0.3e1_dp&
1148 & *t360 + 0.11e2_dp/0.3e1_dp*t363 - t373 - t341 - &
1149 t376/0.3e1_dp + t310 - 0.4e1_dp*t6*t8 + 0.4e1_dp* &
1150 t253 + t210 - t214 - t222
1151 t391 = -t228 + t232 - t236 - t243 - t249 + t252 - 0.4e1_dp/ &
1152 0.3e1_dp*t255 - t304 - t316 - t319 - t322/0.3e1_dp - &
1153 t326/0.3e1_dp + 0.11e2_dp/0.3e1_dp*t329 - t15* &
1154 t21*(t35 + t41 - t48 - t56 + my_rhob*t164 + my_rhoa &
1157 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1158 e_ra_rb(ii) = e_ra_rb(ii) &
1163 t415 = t99*t103*t5*t176/0.3e1_dp
1164 t419 = t15*t108*t176*d/0.3e1_dp
1165 t422 = 0.11e2_dp/0.3e1_dp*t15*t115*t176
1166 t430 = t15*t21*(0.2e1_dp*t7*t130*my_ndrho - 0.8e1_dp &
1167 &/0.3e1_dp*my_rho*my_ndrho)
1169 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1170 e_ndr_ra(ii) = e_ndr_ra(ii) &
1171 - (0.2e1_dp*t408*t313*t171 + t415 + t419 - t422 + &
1173 e_ndr_rb(ii) = e_ndr_rb(ii) &
1174 - (0.2e1_dp*t408*t20*my_rhoa*t171 + t415 + t419 - &
1178 t445 = t99*t103*t5*t189/0.3e1_dp
1179 t449 = t15*t108*t189*d/0.3e1_dp
1180 t452 = 0.11e2_dp/0.3e1_dp*t15*t115*t189
1181 t467 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhoa - &
1182 0.2e1_dp/0.9e1_dp*t135*my_rhoa*t182 + 0.2e1_dp/ &
1183 0.9e1_dp*t181*t95*my_ndrhoa) + 0.8e1_dp/0.3e1_dp&
1184 & *my_rho*my_ndrhoa)
1186 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1187 e_ndra_ra(ii) = e_ndra_ra(ii) &
1188 - (t15*t21*(my_rhob*t185 - 0.2e1_dp/0.9e1_dp*t7&
1189 & *t75*my_ndrhoa) + t445 + t449 - t452 + t467)*sc
1190 e_ndra_rb(ii) = e_ndra_rb(ii) &
1191 - (t15*t21*(my_rhoa*t185 - 0.4e1_dp*my_rhob* &
1192 my_ndrhoa) + t445 + t449 - t452 + t467)*sc
1195 t483 = t99*t103*t5*t202/0.3e1_dp
1196 t487 = t15*t108*t202*d/0.3e1_dp
1197 t490 = 0.11e2_dp/0.3e1_dp*t15*t115*t202
1198 t505 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhob - &
1199 0.2e1_dp/0.9e1_dp*t135*my_rhob*t195 + 0.2e1_dp/ &
1200 0.9e1_dp*t194*t95*my_ndrhob) + 0.8e1_dp/0.3e1_dp&
1201 & *my_rho*my_ndrhob)
1203 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1204 e_ndrb_ra(ii) = e_ndrb_ra(ii) &
1205 - (t15*t21*(my_rhob*t198 - 0.4e1_dp*my_rhoa* &
1206 my_ndrhob) + t483 + t487 - t490 + t505)*sc
1207 e_ndrb_rb(ii) = e_ndrb_rb(ii) &
1208 - (t15*t21*(my_rhoa*t198 - 0.2e1_dp/0.9e1_dp*t7&
1209 & *t75*my_ndrhob) + t483 + t487 - t490 + t505)*sc
1212 t515 = 0.4e1_dp/0.3e1_dp*t16
1214 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1215 e_ndr_ndr(ii) = e_ndr_ndr(ii) &
1216 - (t15*t21*(0.2e1_dp*t7*t39 - t515))*sc
1222 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
1223 e_ndra_ndra(ii) = e_ndra_ndra(ii) &
1224 - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1225 &/0.9e1_dp*t181*t8) + t515 - 0.2e1_dp*t29))*sc
1227 e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) &
1228 - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1229 &/0.9e1_dp*t194*t8) + t515 - 0.2e1_dp*t25))*sc
1235 END SUBROUTINE lyp_lsd_calc
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public lee1988
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
calculates the lyp correlation functional
subroutine, public lyp_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
subroutine, public lyp_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_params)
evaluates the becke 88 exchange functional for lsd
subroutine, public lyp_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
subroutine, public lyp_lda_eval(rho_set, deriv_set, grad_deriv, lyp_params)
evaluates the lyp correlation functional for lda
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