52 USE iso_c_binding,
ONLY: c_size_t, c_int, c_double
61 xc_f03_functional_get_name, &
62 xc_f03_func_get_info, &
63 xc_f03_func_info_get_family, &
64 xc_f03_func_info_get_kind, &
65 xc_f03_func_info_get_n_ext_params, &
66 xc_f03_func_info_get_name, &
67 xc_f03_available_functional_numbers, &
68 xc_f03_available_functional_names, &
69 xc_f03_maximum_name_length, &
70 xc_f03_number_of_functionals, &
71 xc_f03_func_info_get_ext_params_name, &
72 xc_f03_func_info_get_ext_params_description, &
73 xc_f03_func_info_get_ext_params_default_value, &
76 xc_f03_gga_exc_vxc_fxc, &
83 xc_f03_lda_exc_vxc_fxc, &
89 xc_f03_mgga_exc_vxc, &
92 xc_f03_mgga_vxc_fxc, &
100 xc_family_hyb_mgga, &
103 xc_exchange_correlation, &
105 xc_libxc_wrap_info_refs, &
106 xc_libxc_wrap_version, &
107 xc_libxc_wrap_functional_get_number, &
108 xc_libxc_wrap_needs_laplace, &
109 xc_libxc_wrap_functional_set_params, &
110 xc_libxc_wrap_is_under_development, &
111 xc_libxc_get_reference_length, &
112 xc_libxc_check_functional
115#include "../base/base_uses.f90"
120 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'xc_libxc'
127 INTEGER(C_SIZE_T),
PARAMETER,
PRIVATE :: one = 1
144 exists = xc_libxc_check_functional(libxc_params%section%name)
146 mark_used(libxc_params)
161 LOGICAL,
INTENT(IN) :: lsd
165 CHARACTER(len=*),
PARAMETER :: routinen =
'libxc_get_reference_length'
167 CHARACTER(LEN=default_string_length) :: func_name
168 INTEGER :: func_id, handle
169 TYPE(xc_f03_func_t) :: xc_func
170 TYPE(xc_f03_func_info_t) :: xc_info
172 CALL timeset(routinen, handle)
174 func_name = libxc_params%section%name
176 func_id = xc_libxc_wrap_functional_get_number(func_name)
179 CALL xc_f03_func_init(xc_func, func_id, xc_polarized)
181 CALL xc_f03_func_init(xc_func, func_id, xc_unpolarized)
183 xc_info = xc_f03_func_get_info(xc_func)
187 length = xc_libxc_get_reference_length(xc_info)
189 CALL xc_f03_func_end(xc_func)
191 CALL timestop(handle)
193 mark_used(libxc_params)
196 cpabort(
"In order to use LibXC you have to download and install it!")
210 CHARACTER(len=*),
PARAMETER :: routinen =
'libxc_add_sections'
214 INTEGER :: handle, no_func, len_name, ii, func_id, n_param, iparam
215 REAL(kind=c_double) :: default_val
216 CHARACTER(LEN=128) :: func_name, param_name, param_descr, description
217 CHARACTER(LEN=2*default_string_length) :: warning
218 INTEGER(KIND=C_INT),
DIMENSION(:),
ALLOCATABLE :: func_ids
219 TYPE(xc_f03_func_t) :: xc_func
220 TYPE(xc_f03_func_info_t) :: xc_info
222 CALL timeset(routinen, handle)
224 cpassert(
ASSOCIATED(section))
225 NULLIFY (subsection, keyword)
227 no_func = xc_f03_number_of_functionals()
228 len_name = xc_f03_maximum_name_length()
230 ALLOCATE (func_ids(no_func))
232 CALL xc_f03_available_functional_numbers(func_ids)
236 func_id = func_ids(ii)
238 CALL xc_f03_func_init(xc_func, func_id, xc_unpolarized)
239 xc_info = xc_f03_func_get_info(xc_func)
243 func_name = xc_f03_functional_get_name(func_id)
244 description = xc_f03_func_info_get_name(xc_info)
245 n_param = xc_f03_func_info_get_n_ext_params(xc_info)
248 CALL section_create(subsection, __location__, name=trim(func_name), description=trim(description), &
249 n_keywords=2 + n_param, n_subsections=0, repeats=.false.)
251 IF (description(1:1) ==
"_")
THEN
252 warning =
" This parameter is an internal parameter of the functional. Changing this "// &
253 "parameter effectively changes the functional."
259 CALL keyword_create(keyword, __location__, name=
"_SECTION_PARAMETERS_", &
260 description=
"Activates the functional."//trim(warning), &
261 lone_keyword_l_val=.true., default_l_val=.false.)
265 CALL keyword_create(keyword, __location__, name=
"SCALE", description=
"Scales this functional", &
266 default_r_val=1.0_dp)
270 DO iparam = 1, n_param
271 param_name = xc_f03_func_info_get_ext_params_name(xc_info, iparam - 1)
272 param_descr = xc_f03_func_info_get_ext_params_description(xc_info, iparam - 1)
273 default_val = xc_f03_func_info_get_ext_params_default_value(xc_info, iparam - 1)
275 CALL keyword_create(keyword, __location__, name=trim(param_name), &
276 description=trim(param_descr), default_r_val=default_val)
284 CALL xc_f03_func_end(xc_func)
288 CALL timestop(handle)
307 SUBROUTINE libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
310 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
312 INTENT(inout),
OPTIONAL :: needs
313 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
314 LOGICAL,
INTENT(IN),
OPTIONAL :: print_warn
317 CHARACTER(LEN=128) :: s1, s2
318 CHARACTER(LEN=default_string_length) :: func_name
320 REAL(kind=
dp) :: func_scale
321 TYPE(xc_f03_func_t) :: xc_func
322 TYPE(xc_f03_func_info_t) :: xc_info
324 func_name = libxc_params%section%name
330 IF (abs(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
332 func_id = xc_libxc_wrap_functional_get_number(func_name)
334 CALL xc_f03_func_init(xc_func, func_id, xc_unpolarized)
335 xc_info = xc_f03_func_get_info(xc_func)
339 s1 = xc_f03_func_info_get_name(xc_info)
340 SELECT CASE (xc_f03_func_info_get_kind(xc_info))
341 CASE (xc_exchange);
WRITE (s2,
'(a)')
"exchange"
342 CASE (xc_correlation);
WRITE (s2,
'(a)')
"correlation"
343 CASE (xc_exchange_correlation);
WRITE (s2,
'(a)')
"exchange-correlation"
344 CASE (xc_kinetic);
WRITE (s2,
'(a)')
"kinetic"
346 cpabort(trim(func_name)//
": this XC_KIND is currently not supported.")
348 IF (
PRESENT(shortform))
THEN
349 shortform = trim(s1)//
' ('//trim(s2)//
')'
351 IF (
PRESENT(reference))
THEN
352 CALL xc_libxc_wrap_info_refs(xc_info, xc_unpolarized, func_scale, reference)
354 IF (
PRESENT(needs))
THEN
355 SELECT CASE (xc_f03_func_info_get_family(xc_info))
356 CASE (xc_family_lda, xc_family_hyb_lda)
358 CASE (xc_family_gga, xc_family_hyb_gga)
360 needs%norm_drho = .true.
361 CASE (xc_family_mgga, xc_family_hyb_mgga)
363 needs%norm_drho = .true.
365 needs%laplace_rho = xc_libxc_wrap_needs_laplace(func_id)
367 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
370 IF (
PRESENT(max_deriv))
THEN
371 SELECT CASE (xc_f03_func_info_get_family(xc_info))
372 CASE (xc_family_lda, xc_family_hyb_lda)
374 CASE (xc_family_gga, xc_family_hyb_gga)
376 CASE (xc_family_mgga, xc_family_hyb_mgga)
379 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
382 IF (
PRESENT(print_warn))
THEN
383 IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info))
THEN
384 cpwarn(trim(func_name)//
" is under development. Use with caution.")
388 CALL xc_f03_func_end(xc_func)
390 mark_used(libxc_params)
395 mark_used(print_warn)
397 CALL cp_abort(__location__,
"Unknown functional! If you are asking "// &
398 "for a functional of the LibXC library, "// &
399 "you have to download and install the library!")
307 SUBROUTINE libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
…
415 SUBROUTINE libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
418 CHARACTER(LEN=*),
INTENT(OUT),
OPTIONAL :: reference, shortform
420 INTENT(inout),
OPTIONAL :: needs
421 INTEGER,
INTENT(out),
OPTIONAL :: max_deriv
422 LOGICAL,
INTENT(IN),
OPTIONAL :: print_warn
425 CHARACTER(LEN=128) :: s1, s2
426 CHARACTER(LEN=default_string_length) :: func_name
428 REAL(kind=
dp) :: func_scale
429 TYPE(xc_f03_func_t) :: xc_func
430 TYPE(xc_f03_func_info_t) :: xc_info
432 func_name = libxc_params%section%name
438 IF (abs(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
440 func_id = xc_libxc_wrap_functional_get_number(func_name)
442 CALL xc_f03_func_init(xc_func, func_id, xc_polarized)
443 xc_info = xc_f03_func_get_info(xc_func)
447 s1 = xc_f03_func_info_get_name(xc_info)
448 SELECT CASE (xc_f03_func_info_get_kind(xc_info))
449 CASE (xc_exchange);
WRITE (s2,
'(a)')
"exchange"
450 CASE (xc_correlation);
WRITE (s2,
'(a)')
"correlation"
451 CASE (xc_exchange_correlation);
WRITE (s2,
'(a)')
"exchange-correlation"
452 CASE (xc_kinetic);
WRITE (s2,
'(a)')
"kinetic"
454 cpabort(trim(func_name)//
": this XC_KIND is currently not supported.")
456 IF (
PRESENT(shortform))
THEN
457 shortform = trim(s1)//
' ('//trim(s2)//
')'
459 IF (
PRESENT(reference))
THEN
460 CALL xc_libxc_wrap_info_refs(xc_info, xc_polarized, func_scale, reference)
462 IF (
PRESENT(needs))
THEN
463 SELECT CASE (xc_f03_func_info_get_family(xc_info))
464 CASE (xc_family_lda, xc_family_hyb_lda)
465 needs%rho_spin = .true.
466 CASE (xc_family_gga, xc_family_hyb_gga)
467 needs%rho_spin = .true.
468 needs%norm_drho = .true.
469 needs%norm_drho_spin = .true.
470 CASE (xc_family_mgga, xc_family_hyb_mgga)
471 needs%rho_spin = .true.
472 needs%norm_drho = .true.
473 needs%norm_drho_spin = .true.
474 needs%tau_spin = .true.
475 needs%laplace_rho_spin = xc_libxc_wrap_needs_laplace(func_id)
477 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
480 IF (
PRESENT(max_deriv))
THEN
481 SELECT CASE (xc_f03_func_info_get_family(xc_info))
482 CASE (xc_family_lda, xc_family_hyb_lda)
484 CASE (xc_family_gga, xc_family_hyb_gga)
486 CASE (xc_family_mgga, xc_family_hyb_mgga)
489 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
492 IF (
PRESENT(print_warn))
THEN
493 IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info))
THEN
494 cpwarn(trim(func_name)//
" is under development. Use with caution.")
498 CALL xc_f03_func_end(xc_func)
500 mark_used(libxc_params)
505 mark_used(print_warn)
507 CALL cp_abort(__location__,
"Unknown functional! If you are "// &
508 "asking for a functional of the LibXC library, "// &
509 "you have to download and install the library!")
415 SUBROUTINE libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
…
520 CHARACTER(LEN=*),
INTENT(OUT) :: version
523 CALL xc_libxc_wrap_version(version)
526 cpabort(
"In order to use libxc you need to download and install it")
546 INTEGER,
INTENT(in) :: grad_deriv
550 CHARACTER(len=*),
PARAMETER :: routinen =
'libxc_lda_eval'
552 CHARACTER(LEN=default_string_length) :: func_name
553 INTEGER :: func_id, handle, npoints
554 INTEGER,
DIMENSION(2, 3) :: bo
555 LOGICAL :: has_laplace, no_exc
556 REAL(kind=
dp) :: epsilon_rho, epsilon_tau, func_scale
557 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: dummy, e_0, e_laplace_rho, &
558 e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_ndrho, &
559 e_ndrho_laplace_rho, e_ndrho_ndrho, e_ndrho_rho, e_ndrho_tau, e_rho, &
560 e_rho_laplace_rho, e_rho_rho, e_rho_rho_rho, e_rho_tau, e_tau, &
561 e_tau_tau, laplace_rho, norm_drho, rho, tau
563 TYPE(xc_f03_func_t) :: xc_func
564 TYPE(xc_f03_func_info_t) :: xc_info
566 CALL timeset(routinen, handle)
568 has_laplace = .false.
570 NULLIFY (rho, norm_drho, laplace_rho, tau)
572 func_name = libxc_params%section%name
575 IF (abs(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
577 func_id = xc_libxc_wrap_functional_get_number(func_name)
578 CALL xc_f03_func_init(xc_func, func_id, xc_unpolarized)
579 xc_info = xc_f03_func_get_info(xc_func)
580 CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
583 rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, &
584 rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
585 tau=tau, local_bounds=bo)
587 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
592 IF (.NOT.
ASSOCIATED(norm_drho)) norm_drho => dummy
593 IF (.NOT.
ASSOCIATED(tau)) tau => dummy
597 has_laplace = xc_libxc_wrap_needs_laplace(func_id)
598 IF (.NOT. has_laplace) laplace_rho => dummy
603 e_laplace_rho => dummy
607 e_ndrho_ndrho => dummy
608 e_rho_laplace_rho => dummy
610 e_ndrho_laplace_rho => dummy
612 e_laplace_rho_laplace_rho => dummy
613 e_laplace_rho_tau => dummy
615 e_rho_rho_rho => dummy
617 IF (grad_deriv >= 0)
THEN
619 allocate_deriv=.true.)
622 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
623 SELECT CASE (xc_f03_func_info_get_family(xc_info))
624 CASE (xc_family_lda, xc_family_hyb_lda)
626 allocate_deriv=.true.)
628 CASE (xc_family_gga, xc_family_hyb_gga)
630 allocate_deriv=.true.)
633 allocate_deriv=.true.)
635 CASE (xc_family_mgga, xc_family_hyb_mgga)
637 allocate_deriv=.true.)
640 allocate_deriv=.true.)
643 allocate_deriv=.true.)
645 IF (has_laplace)
THEN
647 allocate_deriv=.true.)
651 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
654 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
655 SELECT CASE (xc_f03_func_info_get_family(xc_info))
656 CASE (xc_family_lda, xc_family_hyb_lda)
658 allocate_deriv=.true.)
660 CASE (xc_family_gga, xc_family_hyb_gga)
662 allocate_deriv=.true.)
665 allocate_deriv=.true.)
668 allocate_deriv=.true.)
670 CASE (xc_family_mgga, xc_family_hyb_mgga)
674 allocate_deriv=.true.)
677 allocate_deriv=.true.)
680 allocate_deriv=.true.)
683 allocate_deriv=.true.)
686 allocate_deriv=.true.)
689 allocate_deriv=.true.)
691 IF (has_laplace)
THEN
693 allocate_deriv=.true.)
696 allocate_deriv=.true.)
699 allocate_deriv=.true.)
702 allocate_deriv=.true.)
706 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
709 IF (grad_deriv >= 3 .OR. grad_deriv == -3)
THEN
710 SELECT CASE (xc_f03_func_info_get_family(xc_info))
711 CASE (xc_family_lda, xc_family_hyb_lda)
713 allocate_deriv=.true.)
715 CASE (xc_family_gga, xc_family_hyb_gga, xc_family_mgga, xc_family_hyb_mgga)
716 cpabort(
"derivatives larger than 2 not implemented")
718 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
721 IF (grad_deriv >= 4 .OR. grad_deriv <= -4)
THEN
722 cpabort(
"derivatives larger than 3 not implemented")
734 CALL libxc_lda_calc(rho=rho, norm_drho=norm_drho, &
735 laplace_rho=laplace_rho, tau=tau, &
736 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_laplace_rho=e_laplace_rho, &
737 e_tau=e_tau, e_rho_rho=e_rho_rho, e_ndrho_rho=e_ndrho_rho, &
738 e_ndrho_ndrho=e_ndrho_ndrho, e_rho_laplace_rho=e_rho_laplace_rho, &
739 e_rho_tau=e_rho_tau, e_ndrho_laplace_rho=e_ndrho_laplace_rho, &
740 e_ndrho_tau=e_ndrho_tau, e_laplace_rho_laplace_rho=e_laplace_rho_laplace_rho, &
741 e_laplace_rho_tau=e_laplace_rho_tau, e_tau_tau=e_tau_tau, &
742 e_rho_rho_rho=e_rho_rho_rho, &
743 grad_deriv=grad_deriv, npoints=npoints, &
744 epsilon_rho=epsilon_rho, &
745 epsilon_tau=epsilon_tau, func_name=func_name, &
746 sc=func_scale, xc_func=xc_func, xc_info=xc_info, no_exc=no_exc, has_laplace=has_laplace)
752 CALL xc_f03_func_end(xc_func)
754 CALL timestop(handle)
758 mark_used(grad_deriv)
759 mark_used(libxc_params)
760 CALL cp_abort(__location__,
"Unknown functional! If you are asking "// &
761 "for a functional of the LibXC library, "// &
762 "you have to download and install the library!")
781 INTEGER,
INTENT(in) :: grad_deriv
785 CHARACTER(len=*),
PARAMETER :: routinen =
'libxc_lsd_eval'
787 CHARACTER(LEN=default_string_length) :: func_name
788 INTEGER :: func_id, handle, npoints
789 INTEGER,
DIMENSION(2, 3) :: bo
790 LOGICAL :: has_laplace, no_exc
791 REAL(kind=
dp) :: epsilon_rho, epsilon_tau, func_scale
792 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: dummy, e_0, e_laplace_rhoa, &
793 e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
794 e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, e_laplace_rhob, &
795 e_laplace_rhob_laplace_rhob, e_laplace_rhob_tau_a, &
796 e_laplace_rhob_tau_b, e_ndrho, e_ndrho_laplace_rhoa, &
797 e_ndrho_laplace_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
798 e_ndrho_rhoa, e_ndrho_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa, &
799 e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, e_ndrhoa_ndrhoa, &
800 e_ndrhoa_ndrhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhoa_tau_a, &
801 e_ndrhoa_tau_b, e_ndrhob
802 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :, :),
POINTER :: e_ndrhob_laplace_rhoa, &
803 e_ndrhob_laplace_rhob, e_ndrhob_ndrhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
804 e_ndrhob_tau_a, e_ndrhob_tau_b, e_rhoa, e_rhoa_laplace_rhoa, &
805 e_rhoa_laplace_rhob, e_rhoa_rhoa, e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
806 e_rhoa_rhob, e_rhoa_rhob_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob, &
807 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhob_rhob, &
808 e_rhob_rhob_rhob, e_rhob_tau_a, e_rhob_tau_b, e_tau_a, e_tau_a_tau_a, &
809 e_tau_a_tau_b, e_tau_b, e_tau_b_tau_b, laplace_rhoa, laplace_rhob, &
810 norm_drho, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
812 TYPE(xc_f03_func_t) :: xc_func
813 TYPE(xc_f03_func_info_t) :: xc_info
815 CALL timeset(routinen, handle)
818 NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, &
819 laplace_rhob, tau_a, tau_b)
821 func_name = libxc_params%section%name
824 IF (abs(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp
826 func_id = xc_libxc_wrap_functional_get_number(func_name)
827 CALL xc_f03_func_init(xc_func, func_id, xc_polarized)
828 xc_info = xc_f03_func_get_info(xc_func)
829 CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, libxc_params, no_exc)
832 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
833 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
834 laplace_rhoa=laplace_rhoa, laplace_rhob=laplace_rhob, &
835 rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, &
836 tau_a=tau_a, tau_b=tau_b, local_bounds=bo)
838 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
843 IF (.NOT.
ASSOCIATED(norm_drho)) norm_drho => dummy
844 IF (.NOT.
ASSOCIATED(norm_drhoa)) norm_drhoa => dummy
845 IF (.NOT.
ASSOCIATED(norm_drhob)) norm_drhob => dummy
846 IF (.NOT.
ASSOCIATED(tau_a)) tau_a => dummy
847 IF (.NOT.
ASSOCIATED(tau_b)) tau_b => dummy
851 has_laplace = xc_libxc_wrap_needs_laplace(func_id)
852 IF (.NOT. has_laplace) laplace_rhoa => dummy
853 IF (.NOT. has_laplace) laplace_rhob => dummy
861 e_laplace_rhoa => dummy
862 e_laplace_rhob => dummy
868 e_ndrho_rhoa => dummy
869 e_ndrho_rhob => dummy
870 e_ndrhoa_rhoa => dummy
871 e_ndrhoa_rhob => dummy
872 e_ndrhob_rhoa => dummy
873 e_ndrhob_rhob => dummy
874 e_ndrho_ndrho => dummy
875 e_ndrho_ndrhoa => dummy
876 e_ndrho_ndrhob => dummy
877 e_ndrhoa_ndrhoa => dummy
878 e_ndrhoa_ndrhob => dummy
879 e_ndrhob_ndrhob => dummy
880 e_rhoa_laplace_rhoa => dummy
881 e_rhoa_laplace_rhob => dummy
882 e_rhob_laplace_rhoa => dummy
883 e_rhob_laplace_rhob => dummy
884 e_rhoa_tau_a => dummy
885 e_rhoa_tau_b => dummy
886 e_rhob_tau_a => dummy
887 e_rhob_tau_b => dummy
888 e_ndrho_laplace_rhoa => dummy
889 e_ndrho_laplace_rhob => dummy
890 e_ndrhoa_laplace_rhoa => dummy
891 e_ndrhoa_laplace_rhob => dummy
892 e_ndrhob_laplace_rhoa => dummy
893 e_ndrhob_laplace_rhob => dummy
894 e_ndrho_tau_a => dummy
895 e_ndrho_tau_b => dummy
896 e_ndrhoa_tau_a => dummy
897 e_ndrhoa_tau_b => dummy
898 e_ndrhob_tau_a => dummy
899 e_ndrhob_tau_b => dummy
900 e_laplace_rhoa_laplace_rhoa => dummy
901 e_laplace_rhoa_laplace_rhob => dummy
902 e_laplace_rhob_laplace_rhob => dummy
903 e_laplace_rhoa_tau_a => dummy
904 e_laplace_rhoa_tau_b => dummy
905 e_laplace_rhob_tau_a => dummy
906 e_laplace_rhob_tau_b => dummy
907 e_tau_a_tau_a => dummy
908 e_tau_a_tau_b => dummy
909 e_tau_b_tau_b => dummy
910 e_rhoa_rhoa_rhoa => dummy
911 e_rhoa_rhoa_rhob => dummy
912 e_rhoa_rhob_rhob => dummy
913 e_rhob_rhob_rhob => dummy
915 IF (grad_deriv >= 0)
THEN
917 allocate_deriv=.true.)
920 IF (grad_deriv >= 1 .OR. grad_deriv == -1)
THEN
921 SELECT CASE (xc_f03_func_info_get_family(xc_info))
922 CASE (xc_family_lda, xc_family_hyb_lda)
924 allocate_deriv=.true.)
927 allocate_deriv=.true.)
929 CASE (xc_family_gga, xc_family_hyb_gga)
931 allocate_deriv=.true.)
934 allocate_deriv=.true.)
937 allocate_deriv=.true.)
940 allocate_deriv=.true.)
943 allocate_deriv=.true.)
945 CASE (xc_family_mgga, xc_family_hyb_mgga)
947 allocate_deriv=.true.)
950 allocate_deriv=.true.)
953 allocate_deriv=.true.)
956 allocate_deriv=.true.)
959 allocate_deriv=.true.)
962 allocate_deriv=.true.)
965 allocate_deriv=.true.)
967 IF (has_laplace)
THEN
969 allocate_deriv=.true.)
972 allocate_deriv=.true.)
976 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
979 IF (grad_deriv >= 2 .OR. grad_deriv == -2)
THEN
980 SELECT CASE (xc_f03_func_info_get_family(xc_info))
981 CASE (xc_family_lda, xc_family_hyb_lda)
983 allocate_deriv=.true.)
986 allocate_deriv=.true.)
989 allocate_deriv=.true.)
991 CASE (xc_family_gga, xc_family_hyb_gga)
993 allocate_deriv=.true.)
996 allocate_deriv=.true.)
999 allocate_deriv=.true.)
1002 allocate_deriv=.true.)
1005 allocate_deriv=.true.)
1008 allocate_deriv=.true.)
1011 allocate_deriv=.true.)
1014 allocate_deriv=.true.)
1017 allocate_deriv=.true.)
1020 allocate_deriv=.true.)
1023 allocate_deriv=.true.)
1026 allocate_deriv=.true.)
1029 allocate_deriv=.true.)
1032 allocate_deriv=.true.)
1035 allocate_deriv=.true.)
1037 CASE (xc_family_mgga, xc_family_hyb_mgga)
1040 allocate_deriv=.true.)
1043 allocate_deriv=.true.)
1046 allocate_deriv=.true.)
1049 allocate_deriv=.true.)
1052 allocate_deriv=.true.)
1055 allocate_deriv=.true.)
1058 allocate_deriv=.true.)
1061 allocate_deriv=.true.)
1064 allocate_deriv=.true.)
1067 allocate_deriv=.true.)
1070 allocate_deriv=.true.)
1073 allocate_deriv=.true.)
1076 allocate_deriv=.true.)
1079 allocate_deriv=.true.)
1082 allocate_deriv=.true.)
1085 allocate_deriv=.true.)
1088 allocate_deriv=.true.)
1091 allocate_deriv=.true.)
1094 allocate_deriv=.true.)
1097 allocate_deriv=.true.)
1100 allocate_deriv=.true.)
1103 allocate_deriv=.true.)
1106 allocate_deriv=.true.)
1109 allocate_deriv=.true.)
1112 allocate_deriv=.true.)
1115 allocate_deriv=.true.)
1118 allocate_deriv=.true.)
1121 allocate_deriv=.true.)
1123 IF (has_laplace)
THEN
1125 allocate_deriv=.true.)
1128 allocate_deriv=.true.)
1131 allocate_deriv=.true.)
1134 allocate_deriv=.true.)
1137 allocate_deriv=.true.)
1140 allocate_deriv=.true.)
1143 allocate_deriv=.true.)
1146 allocate_deriv=.true.)
1149 allocate_deriv=.true.)
1152 allocate_deriv=.true.)
1155 allocate_deriv=.true.)
1158 allocate_deriv=.true.)
1161 allocate_deriv=.true.)
1164 allocate_deriv=.true.)
1167 allocate_deriv=.true.)
1170 allocate_deriv=.true.)
1173 allocate_deriv=.true.)
1177 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
1180 IF (grad_deriv >= 3 .OR. grad_deriv == -3)
THEN
1181 SELECT CASE (xc_f03_func_info_get_family(xc_info))
1182 CASE (xc_family_lda, xc_family_hyb_lda)
1184 allocate_deriv=.true.)
1187 allocate_deriv=.true.)
1190 allocate_deriv=.true.)
1193 allocate_deriv=.true.)
1195 CASE (xc_family_gga, xc_family_hyb_gga, xc_family_mgga, xc_family_hyb_mgga)
1196 cpabort(
"derivatives larger than 2 not implemented")
1198 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
1201 IF (grad_deriv >= 4 .OR. grad_deriv <= -4)
THEN
1202 cpabort(
"derivatives larger than 3 not implemented")
1229 CALL libxc_lsd_calc(rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, &
1230 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, laplace_rhoa=laplace_rhoa, &
1231 laplace_rhob=laplace_rhob, tau_a=tau_a, tau_b=tau_b, &
1232 e_0=e_0, e_rhoa=e_rhoa, e_rhob=e_rhob, e_ndrho=e_ndrho, &
1233 e_ndrhoa=e_ndrhoa, e_ndrhob=e_ndrhob, e_laplace_rhoa=e_laplace_rhoa, &
1234 e_laplace_rhob=e_laplace_rhob, e_tau_a=e_tau_a, e_tau_b=e_tau_b, &
1235 e_rhoa_rhoa=e_rhoa_rhoa, e_rhoa_rhob=e_rhoa_rhob, e_rhob_rhob=e_rhob_rhob, &
1236 e_ndrho_rhoa=e_ndrho_rhoa, e_ndrho_rhob=e_ndrho_rhob, &
1237 e_ndrhoa_rhoa=e_ndrhoa_rhoa, e_ndrhoa_rhob=e_ndrhoa_rhob, &
1238 e_ndrhob_rhoa=e_ndrhob_rhoa, e_ndrhob_rhob=e_ndrhob_rhob, &
1239 e_ndrho_ndrho=e_ndrho_ndrho, e_ndrho_ndrhoa=e_ndrho_ndrhoa, &
1240 e_ndrho_ndrhob=e_ndrho_ndrhob, e_ndrhoa_ndrhoa=e_ndrhoa_ndrhoa, &
1241 e_ndrhoa_ndrhob=e_ndrhoa_ndrhob, e_ndrhob_ndrhob=e_ndrhob_ndrhob, &
1242 e_rhoa_laplace_rhoa=e_rhoa_laplace_rhoa, &
1243 e_rhoa_laplace_rhob=e_rhoa_laplace_rhob, &
1244 e_rhob_laplace_rhoa=e_rhob_laplace_rhoa, &
1245 e_rhob_laplace_rhob=e_rhob_laplace_rhob, &
1246 e_rhoa_tau_a=e_rhoa_tau_a, e_rhoa_tau_b=e_rhoa_tau_b, &
1247 e_rhob_tau_a=e_rhob_tau_a, e_rhob_tau_b=e_rhob_tau_b, &
1248 e_ndrho_laplace_rhoa=e_ndrho_laplace_rhoa, &
1249 e_ndrho_laplace_rhob=e_ndrho_laplace_rhob, &
1250 e_ndrhoa_laplace_rhoa=e_ndrhoa_laplace_rhoa, &
1251 e_ndrhoa_laplace_rhob=e_ndrhoa_laplace_rhob, &
1252 e_ndrhob_laplace_rhoa=e_ndrhob_laplace_rhoa, &
1253 e_ndrhob_laplace_rhob=e_ndrhob_laplace_rhob, &
1254 e_ndrho_tau_a=e_ndrho_tau_a, e_ndrho_tau_b=e_ndrho_tau_b, &
1255 e_ndrhoa_tau_a=e_ndrhoa_tau_a, e_ndrhoa_tau_b=e_ndrhoa_tau_b, &
1256 e_ndrhob_tau_a=e_ndrhob_tau_a, e_ndrhob_tau_b=e_ndrhob_tau_b, &
1257 e_laplace_rhoa_laplace_rhoa=e_laplace_rhoa_laplace_rhoa, &
1258 e_laplace_rhoa_laplace_rhob=e_laplace_rhoa_laplace_rhob, &
1259 e_laplace_rhob_laplace_rhob=e_laplace_rhob_laplace_rhob, &
1260 e_laplace_rhoa_tau_a=e_laplace_rhoa_tau_a, &
1261 e_laplace_rhoa_tau_b=e_laplace_rhoa_tau_b, &
1262 e_laplace_rhob_tau_a=e_laplace_rhob_tau_a, &
1263 e_laplace_rhob_tau_b=e_laplace_rhob_tau_b, &
1264 e_tau_a_tau_a=e_tau_a_tau_a, &
1265 e_tau_a_tau_b=e_tau_a_tau_b, &
1266 e_tau_b_tau_b=e_tau_b_tau_b, &
1267 e_rhoa_rhoa_rhoa=e_rhoa_rhoa_rhoa, &
1268 e_rhoa_rhoa_rhob=e_rhoa_rhoa_rhob, &
1269 e_rhoa_rhob_rhob=e_rhoa_rhob_rhob, &
1270 e_rhob_rhob_rhob=e_rhob_rhob_rhob, &
1271 grad_deriv=grad_deriv, npoints=npoints, &
1272 epsilon_rho=epsilon_rho, &
1273 epsilon_tau=epsilon_tau, func_name=func_name, &
1274 sc=func_scale, xc_func=xc_func, xc_info=xc_info, no_exc=no_exc, has_laplace=has_laplace)
1280 CALL xc_f03_func_end(xc_func)
1282 CALL timestop(handle)
1285 mark_used(deriv_set)
1286 mark_used(grad_deriv)
1287 mark_used(libxc_params)
1289 CALL cp_abort(__location__,
"Unknown functional! If you are asking "// &
1290 "for a functional of the LibXC library, "// &
1291 "you have to download and install the library!")
1331#if defined (__LIBXC)
1332 SUBROUTINE libxc_lda_calc(rho, norm_drho, laplace_rho, tau, &
1333 e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, e_rho_rho, e_ndrho_rho, &
1334 e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1335 e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, &
1336 e_tau_tau, e_rho_rho_rho, &
1337 grad_deriv, npoints, epsilon_rho, &
1338 epsilon_tau, func_name, sc, xc_func, xc_info, no_exc, has_laplace)
1340 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rho, norm_drho, laplace_rho, tau
1341 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, &
1342 e_rho_rho, e_ndrho_rho, e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, &
1343 e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_tau_tau, e_rho_rho_rho
1344 INTEGER,
INTENT(in) :: grad_deriv, npoints
1345 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, epsilon_tau
1346 CHARACTER(LEN=default_string_length),
INTENT(IN) :: func_name
1347 REAL(kind=
dp),
INTENT(in) :: sc
1348 TYPE(xc_f03_func_t),
INTENT(IN) :: xc_func
1349 TYPE(xc_f03_func_info_t),
INTENT(IN) :: xc_info
1350 LOGICAL,
INTENT(IN) :: no_exc, has_laplace
1353 REAL(kind=
dp),
DIMENSION(1) :: exc, my_tau, sigma, v2lapl2, v2lapltau, v2rho2, v2rholapl, &
1354 v2rhosigma, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2tau2, v3rho3, vlapl, vrho, &
1360 SELECT CASE (xc_f03_func_info_get_family(xc_info))
1361 CASE (xc_family_lda, xc_family_hyb_lda)
1362 IF (grad_deriv == 0)
THEN
1365 IF (rho(ii) > epsilon_rho)
THEN
1366 CALL xc_f03_lda_exc(xc_func, one, rho(ii), exc)
1367 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1371 ELSE IF (grad_deriv == -1)
THEN
1374 IF (rho(ii) > epsilon_rho)
THEN
1375 CALL xc_f03_lda_vxc(xc_func, one, rho(ii), vrho)
1376 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1380 ELSE IF (grad_deriv == 1)
THEN
1383 IF (rho(ii) > epsilon_rho)
THEN
1384 CALL xc_f03_lda_exc_vxc(xc_func, one, rho(ii), exc, vrho)
1385 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1386 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1390 ELSE IF (grad_deriv == -2)
THEN
1393 IF (rho(ii) > epsilon_rho)
THEN
1394 CALL xc_f03_lda_fxc(xc_func, one, rho(ii), v2rho2)
1395 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1399 ELSE IF (grad_deriv == 2)
THEN
1402 IF (rho(ii) > epsilon_rho)
THEN
1403 CALL xc_f03_lda_exc_vxc_fxc(xc_func, one, rho(ii), exc, vrho, v2rho2)
1404 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1405 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1406 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1410 ELSE IF (grad_deriv == -3)
THEN
1413 IF (rho(ii) > epsilon_rho)
THEN
1414 CALL xc_f03_lda_kxc(xc_func, one, rho(ii), v3rho3)
1415 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1)
1419 ELSE IF (grad_deriv == 3)
THEN
1422 IF (rho(ii) > epsilon_rho)
THEN
1423 CALL xc_f03_lda(xc_func, one, rho(ii), exc, vrho, v2rho2, v3rho3)
1424 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1425 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1426 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1427 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1)
1432 CASE (xc_family_gga, xc_family_hyb_gga)
1433 IF (grad_deriv == 0)
THEN
1436 IF (rho(ii) > epsilon_rho)
THEN
1437 sigma = norm_drho(ii)**2
1438 CALL xc_f03_gga_exc(xc_func, one, rho(ii), sigma, exc)
1439 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1443 ELSE IF (grad_deriv == -1)
THEN
1446 IF (rho(ii) > epsilon_rho)
THEN
1447 sigma = norm_drho(ii)**2
1448 CALL xc_f03_gga_vxc(xc_func, one, rho(ii), sigma, vrho, vsigma)
1449 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1450 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1454 ELSE IF (grad_deriv == 1)
THEN
1457 IF (rho(ii) > epsilon_rho)
THEN
1458 sigma = norm_drho(ii)**2
1460 CALL xc_f03_gga_vxc(xc_func, one, rho(ii), sigma, vrho, vsigma)
1463 CALL xc_f03_gga_exc_vxc(xc_func, one, rho(ii), sigma, &
1466 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1467 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1468 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1472 ELSE IF (grad_deriv == -2)
THEN
1475 IF (rho(ii) > epsilon_rho)
THEN
1476 sigma = norm_drho(ii)**2
1478 CALL xc_f03_gga_vxc_fxc(xc_func, one, rho(ii), sigma, vrho, vsigma, &
1479 v2rho2, v2rhosigma, v2sigma2)
1481 CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rho(ii), sigma, &
1482 exc, vrho, vsigma, v2rho2, &
1483 v2rhosigma, v2sigma2)
1485 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1486 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1487 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1488 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1492 ELSE IF (grad_deriv == 2)
THEN
1495 IF (rho(ii) > epsilon_rho)
THEN
1496 sigma = norm_drho(ii)**2
1498 CALL xc_f03_gga_vxc_fxc(xc_func, one, rho(ii), sigma, vrho, vsigma, &
1499 v2rho2, v2rhosigma, v2sigma2)
1502 CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rho(ii), sigma, &
1503 exc, vrho, vsigma, &
1504 v2rho2, v2rhosigma, v2sigma2)
1506 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1507 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1508 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1509 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1510 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1511 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1512 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1517 CASE (xc_family_mgga, xc_family_hyb_mgga)
1518 IF (grad_deriv == 0)
THEN
1521 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau))
THEN
1522 sigma = norm_drho(ii)**2
1523 my_tau(1) = max(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1524 CALL xc_f03_mgga_exc(xc_func, one, rho(ii), sigma, &
1525 laplace_rho(ii), my_tau, exc)
1526 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1530 ELSE IF (grad_deriv == -1)
THEN
1533 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau))
THEN
1534 sigma = norm_drho(ii)**2
1535 my_tau(1) = max(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1536 CALL xc_f03_mgga_vxc(xc_func, one, rho(ii), sigma, &
1537 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau)
1538 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1539 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1540 IF (has_laplace) e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1541 e_tau(ii) = e_tau(ii) + sc*vtau(1)
1545 ELSE IF (grad_deriv == 1)
THEN
1548 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau))
THEN
1549 sigma(1) = norm_drho(ii)**2
1550 my_tau(1) = max(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1552 CALL xc_f03_mgga_vxc(xc_func, one, rho(ii), sigma, &
1553 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau)
1556 CALL xc_f03_mgga_exc_vxc(xc_func, one, rho(ii), sigma, &
1557 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau)
1559 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1560 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1561 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1562 IF (has_laplace) e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1563 e_tau(ii) = e_tau(ii) + sc*vtau(1)
1567 ELSE IF (grad_deriv == -2)
THEN
1570 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau))
THEN
1571 sigma = norm_drho(ii)**2
1572 my_tau(1) = max(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1574 CALL xc_f03_mgga_vxc_fxc(xc_func, one, rho(ii), sigma, &
1575 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau, &
1576 v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1577 v2lapl2, v2lapltau, v2tau2)
1579 CALL xc_f03_mgga(xc_func, one, rho(ii), sigma, &
1580 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, &
1581 v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1582 v2lapl2, v2lapltau, v2tau2)
1584 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1585 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1586 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1587 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1588 e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1)
1589 e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii)
1590 e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1)
1591 IF (has_laplace)
THEN
1592 e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1)
1593 e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
1594 sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii)
1595 e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1)
1596 e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1)
1601 ELSE IF (grad_deriv == 2)
THEN
1604 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau))
THEN
1605 sigma = norm_drho(ii)**2
1606 my_tau(1) = max(tau(ii), sigma(1)/(8.0_dp*rho(ii)))
1608 CALL xc_f03_mgga_vxc_fxc(xc_func, one, rho(ii), sigma, &
1609 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau, &
1610 v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1611 v2lapl2, v2lapltau, v2tau2)
1614 CALL xc_f03_mgga(xc_func, one, rho(ii), sigma, &
1615 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, &
1616 v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, &
1617 v2lapl2, v2lapltau, v2tau2)
1619 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii)
1620 e_rho(ii) = e_rho(ii) + sc*vrho(1)
1621 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii)
1622 e_tau(ii) = e_tau(ii) + sc*vtau(1)
1623 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1)
1624 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii)
1625 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
1626 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1))
1627 e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1)
1628 e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii)
1629 e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1)
1630 IF (has_laplace)
THEN
1631 e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1)
1632 e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1)
1633 e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + &
1634 sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii)
1635 e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1)
1636 e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1)
1643 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
1646 END SUBROUTINE libxc_lda_calc
1733#if defined (__LIBXC)
1734 SUBROUTINE libxc_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, &
1735 norm_drhob, laplace_rhoa, laplace_rhob, tau_a, tau_b, &
1736 e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, e_ndrhob, &
1737 e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, &
1738 e_rhoa_rhoa, e_rhoa_rhob, e_rhob_rhob, &
1739 e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, &
1740 e_ndrhoa_rhob, e_ndrhob_rhoa, e_ndrhob_rhob, &
1741 e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, &
1742 e_ndrhoa_ndrhoa, e_ndrhoa_ndrhob, e_ndrhob_ndrhob, &
1743 e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
1744 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, &
1745 e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, e_rhob_tau_b, &
1746 e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, &
1747 e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, &
1748 e_ndrhob_laplace_rhoa, e_ndrhob_laplace_rhob, &
1749 e_ndrho_tau_a, e_ndrho_tau_b, &
1750 e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
1751 e_ndrhob_tau_a, e_ndrhob_tau_b, &
1752 e_laplace_rhoa_laplace_rhoa, &
1753 e_laplace_rhoa_laplace_rhob, &
1754 e_laplace_rhob_laplace_rhob, &
1755 e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
1756 e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, &
1757 e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
1758 e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, &
1759 e_rhoa_rhob_rhob, e_rhob_rhob_rhob, &
1760 grad_deriv, npoints, epsilon_rho, &
1761 epsilon_tau, func_name, sc, xc_func, xc_info, no_exc, has_laplace)
1763 REAL(kind=
dp),
DIMENSION(*),
INTENT(IN) :: rhoa, rhob, norm_drho, norm_drhoa, &
1764 norm_drhob, laplace_rhoa, &
1765 laplace_rhob, tau_a, tau_b
1766 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, &
1767 e_ndrhob, e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, e_rhoa_rhoa, e_rhoa_rhob, &
1768 e_rhob_rhob, e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhob_rhoa, &
1769 e_ndrhob_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, e_ndrhoa_ndrhoa, &
1770 e_ndrhoa_ndrhob, e_ndrhob_ndrhob, e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, &
1771 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, &
1772 e_rhob_tau_b, e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, e_ndrhoa_laplace_rhoa
1773 REAL(kind=
dp),
DIMENSION(*),
INTENT(INOUT) :: e_ndrhoa_laplace_rhob, e_ndrhob_laplace_rhoa, &
1774 e_ndrhob_laplace_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa_tau_a, e_ndrhoa_tau_b, &
1775 e_ndrhob_tau_a, e_ndrhob_tau_b, e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, &
1776 e_laplace_rhob_laplace_rhob, e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, &
1777 e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, &
1778 e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, e_rhoa_rhob_rhob, e_rhob_rhob_rhob
1779 INTEGER,
INTENT(in) :: grad_deriv, npoints
1780 REAL(kind=
dp),
INTENT(in) :: epsilon_rho, epsilon_tau
1781 CHARACTER(LEN=default_string_length),
INTENT(IN) :: func_name
1782 REAL(kind=
dp),
INTENT(in) :: sc
1783 TYPE(xc_f03_func_t),
INTENT(IN) :: xc_func
1784 TYPE(xc_f03_func_info_t),
INTENT(IN) :: xc_info
1785 LOGICAL,
INTENT(IN) :: no_exc, has_laplace
1788 REAL(kind=
dp) :: my_norm_drho, my_norm_drhoa, &
1789 my_norm_drhob, my_rhoa, my_rhob, &
1791 REAL(kind=
dp),
DIMENSION(1) :: exc
1792 REAL(kind=
dp),
DIMENSION(2, 1) :: laplace_rhov, rhov, tauv, vlapl, vrho, &
1794 REAL(kind=
dp),
DIMENSION(3, 1) :: sigmav, v2lapl2, v2rho2, v2tau2, vsigma
1795 REAL(kind=
dp),
DIMENSION(4, 1) :: v2lapltau, v2rholapl, v2rhotau, v3rho3
1796 REAL(kind=
dp),
DIMENSION(6, 1) :: v2rhosigma, v2sigma2, v2sigmalapl, &
1799 vlapl(1, 1) = 0.0_dp
1800 vlapl(2, 1) = 0.0_dp
1802 SELECT CASE (xc_f03_func_info_get_family(xc_info))
1803 CASE (xc_family_lda, xc_family_hyb_lda)
1804 IF (grad_deriv == 0)
THEN
1807 my_rhoa = max(rhoa(ii), 0.0_dp)
1808 my_rhob = max(rhob(ii), 0.0_dp)
1809 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1810 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1811 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1812 CALL xc_f03_lda_exc(xc_func, one, rhov(1, 1), exc)
1813 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1817 ELSE IF (grad_deriv == -1)
THEN
1820 my_rhoa = max(rhoa(ii), 0.0_dp)
1821 my_rhob = max(rhob(ii), 0.0_dp)
1822 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1823 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1824 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1825 CALL xc_f03_lda_vxc(xc_func, one, rhov(1, 1), vrho(1, 1))
1826 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1827 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1831 ELSE IF (grad_deriv == 1)
THEN
1834 my_rhoa = max(rhoa(ii), 0.0_dp)
1835 my_rhob = max(rhob(ii), 0.0_dp)
1836 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1837 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1838 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1839 CALL xc_f03_lda_exc_vxc(xc_func, one, rhov(1, 1), exc, vrho(1, 1))
1840 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1841 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1842 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1846 ELSE IF (grad_deriv == -2)
THEN
1849 my_rhoa = max(rhoa(ii), 0.0_dp)
1850 my_rhob = max(rhob(ii), 0.0_dp)
1851 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1852 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1853 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1854 CALL xc_f03_lda_fxc(xc_func, one, rhov(1, 1), v2rho2(1, 1))
1855 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1856 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1857 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1861 ELSE IF (grad_deriv == 2)
THEN
1864 my_rhoa = max(rhoa(ii), 0.0_dp)
1865 my_rhob = max(rhob(ii), 0.0_dp)
1866 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1867 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1868 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1869 CALL xc_f03_lda_exc_vxc_fxc(xc_func, one, rhov(1, 1), exc, vrho(1, 1), v2rho2(1, 1))
1870 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1871 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1872 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1873 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1874 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1875 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1879 ELSE IF (grad_deriv == -3)
THEN
1882 my_rhoa = max(rhoa(ii), 0.0_dp)
1883 my_rhob = max(rhob(ii), 0.0_dp)
1884 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1885 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1886 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1887 CALL xc_f03_lda_kxc(xc_func, one, rhov(1, 1), v3rho3(1, 1))
1888 e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1)
1889 e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1)
1890 e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1)
1891 e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1)
1895 ELSE IF (grad_deriv == 3)
THEN
1898 my_rhoa = max(rhoa(ii), 0.0_dp)
1899 my_rhob = max(rhob(ii), 0.0_dp)
1900 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1901 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1902 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1903 CALL xc_f03_lda(xc_func, one, rhov(1, 1), exc, vrho(1, 1), v2rho2(1, 1), v3rho3(1, 1))
1904 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1905 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1906 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1907 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
1908 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
1909 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
1910 e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1)
1911 e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1)
1912 e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1)
1913 e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1)
1918 CASE (xc_family_gga, xc_family_hyb_gga)
1919 IF (grad_deriv == 0)
THEN
1922 my_rhoa = max(rhoa(ii), 0.0_dp)
1923 my_rhob = max(rhob(ii), 0.0_dp)
1924 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1925 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1926 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1927 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
1928 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
1929 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
1930 sigmav(1, 1) = my_norm_drhoa**2
1931 sigmav(3, 1) = my_norm_drhob**2
1932 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
1933 CALL xc_f03_gga_exc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc)
1934 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1938 ELSE IF (grad_deriv == -1)
THEN
1941 my_rhoa = max(rhoa(ii), 0.0_dp)
1942 my_rhob = max(rhob(ii), 0.0_dp)
1943 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1944 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1945 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1946 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
1947 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
1948 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
1949 sigmav(1, 1) = my_norm_drhoa**2
1950 sigmav(3, 1) = my_norm_drhob**2
1951 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
1952 CALL xc_f03_gga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1))
1953 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1954 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1955 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
1956 e_ndrhoa(ii) = e_ndrhoa(ii) + &
1957 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
1958 e_ndrhob(ii) = e_ndrhob(ii) + &
1959 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
1963 ELSE IF (grad_deriv == 1)
THEN
1966 my_rhoa = max(rhoa(ii), 0.0_dp)
1967 my_rhob = max(rhob(ii), 0.0_dp)
1968 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
1969 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
1970 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
1971 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
1972 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
1973 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
1974 sigmav(1, 1) = my_norm_drhoa**2
1975 sigmav(3, 1) = my_norm_drhob**2
1976 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
1978 CALL xc_f03_gga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1))
1981 CALL xc_f03_gga_exc_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1))
1983 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
1984 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
1985 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
1986 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
1987 e_ndrhoa(ii) = e_ndrhoa(ii) + &
1988 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
1989 e_ndrhob(ii) = e_ndrhob(ii) + &
1990 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
1994 ELSE IF (grad_deriv == -2)
THEN
1997 my_rhoa = max(rhoa(ii), 0.0_dp)
1998 my_rhob = max(rhob(ii), 0.0_dp)
1999 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
2000 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2001 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2002 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2003 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2004 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2005 sigmav(1, 1) = my_norm_drhoa**2
2006 sigmav(3, 1) = my_norm_drhob**2
2007 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2009 CALL xc_f03_gga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1), &
2010 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2012 CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2013 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2015 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2016 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2017 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2018 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2019 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2020 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2021 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2022 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2023 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2024 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2025 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2026 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2027 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2028 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2029 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2030 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2031 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2032 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2033 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2034 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2035 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2036 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2037 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2038 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2039 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2040 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2041 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2042 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2046 ELSE IF (grad_deriv == 2)
THEN
2049 my_rhoa = max(rhoa(ii), 0.0_dp)
2050 my_rhob = max(rhob(ii), 0.0_dp)
2051 IF ((my_rhoa + my_rhob) > epsilon_rho)
THEN
2052 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2053 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2054 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2055 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2056 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2057 sigmav(1, 1) = my_norm_drhoa**2
2058 sigmav(3, 1) = my_norm_drhob**2
2059 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2061 CALL xc_f03_gga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1), &
2062 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2065 CALL xc_f03_gga_exc_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2066 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1))
2068 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2069 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2070 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2071 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2072 e_ndrhoa(ii) = e_ndrhoa(ii) + &
2073 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2074 e_ndrhob(ii) = e_ndrhob(ii) + &
2075 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2076 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2077 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2078 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2079 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2080 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2081 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2082 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2083 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2084 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2085 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2086 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2087 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2088 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2089 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2090 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2091 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2092 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2093 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2094 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2095 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2096 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2097 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2098 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2099 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2100 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2101 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2102 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2103 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2108 CASE (xc_family_mgga, xc_family_hyb_mgga)
2109 IF (grad_deriv == 0)
THEN
2112 my_rhoa = max(rhoa(ii), 0.0_dp)
2113 my_rhob = max(rhob(ii), 0.0_dp)
2114 my_tau_a = max(tau_a(ii), 0.0_dp)
2115 my_tau_b = max(tau_b(ii), 0.0_dp)
2116 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau))
THEN
2117 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2118 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2119 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2120 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2121 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2122 sigmav(1, 1) = my_norm_drhoa**2
2123 sigmav(3, 1) = my_norm_drhob**2
2124 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2125 tauv(1, 1) = max(my_tau_a, epsilon(0.0_dp)*1.e4_dp)
2126 tauv(2, 1) = max(my_tau_b, epsilon(0.0_dp)*1.e4_dp)
2127 tauv(1, 1) = max(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2128 tauv(2, 1) = max(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2129 laplace_rhov(1, 1) = laplace_rhoa(ii)
2130 laplace_rhov(2, 1) = laplace_rhob(ii)
2131 CALL xc_f03_mgga_exc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2132 laplace_rhov(1, 1), tauv(1, 1), exc)
2133 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2137 ELSE IF (grad_deriv == -1)
THEN
2140 my_rhoa = max(rhoa(ii), 0.0_dp)
2141 my_rhob = max(rhob(ii), 0.0_dp)
2142 my_tau_a = max(tau_a(ii), 0.0_dp)
2143 my_tau_b = max(tau_b(ii), 0.0_dp)
2144 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau))
THEN
2145 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2146 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2147 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2148 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2149 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2150 sigmav(1, 1) = my_norm_drhoa**2
2151 sigmav(3, 1) = my_norm_drhob**2
2152 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2153 laplace_rhov(1, 1) = laplace_rhoa(ii)
2154 laplace_rhov(2, 1) = laplace_rhob(ii)
2155 tauv(1, 1) = max(my_tau_a, epsilon(0.0_dp)*1.e4_dp)
2156 tauv(2, 1) = max(my_tau_b, epsilon(0.0_dp)*1.e4_dp)
2157 tauv(1, 1) = max(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2158 tauv(2, 1) = max(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2159 CALL xc_f03_mgga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2160 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1))
2161 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2162 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2163 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2164 e_ndrhoa(ii) = e_ndrhoa(ii) + &
2165 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2166 e_ndrhob(ii) = e_ndrhob(ii) + &
2167 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2168 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2169 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2170 IF (has_laplace)
THEN
2171 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2172 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2177 ELSE IF (grad_deriv == 1)
THEN
2180 my_rhoa = max(rhoa(ii), 0.0_dp)
2181 my_rhob = max(rhob(ii), 0.0_dp)
2182 my_tau_a = max(tau_a(ii), 0.0_dp)
2183 my_tau_b = max(tau_b(ii), 0.0_dp)
2184 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau))
THEN
2185 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2186 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2187 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2188 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2189 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2190 sigmav(1, 1) = my_norm_drhoa**2
2191 sigmav(3, 1) = my_norm_drhob**2
2192 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2193 laplace_rhov(1, 1) = laplace_rhoa(ii)
2194 laplace_rhov(2, 1) = laplace_rhob(ii)
2195 tauv(1, 1) = max(my_tau_a, epsilon(0.0_dp)*1.e4_dp)
2196 tauv(2, 1) = max(my_tau_b, epsilon(0.0_dp)*1.e4_dp)
2197 tauv(1, 1) = max(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2198 tauv(2, 1) = max(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2200 CALL xc_f03_mgga_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2201 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2202 vlapl(1, 1), vtau(1, 1))
2205 CALL xc_f03_mgga_exc_vxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2206 laplace_rhov(1, 1), tauv(1, 1), exc, &
2207 vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1))
2209 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2210 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2211 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2212 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2213 e_ndrhoa(ii) = e_ndrhoa(ii) + &
2214 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2215 e_ndrhob(ii) = e_ndrhob(ii) + &
2216 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2217 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2218 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2219 IF (has_laplace)
THEN
2220 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2221 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2226 ELSE IF (grad_deriv == -2)
THEN
2229 my_rhoa = max(rhoa(ii), 0.0_dp)
2230 my_rhob = max(rhob(ii), 0.0_dp)
2231 my_tau_a = max(tau_a(ii), 0.0_dp)
2232 my_tau_b = max(tau_b(ii), 0.0_dp)
2233 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau))
THEN
2234 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2235 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2236 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2237 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2238 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2239 sigmav(1, 1) = my_norm_drhoa**2
2240 sigmav(3, 1) = my_norm_drhob**2
2241 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2242 laplace_rhov(1, 1) = laplace_rhoa(ii)
2243 laplace_rhov(2, 1) = laplace_rhob(ii)
2244 tauv(1, 1) = max(my_tau_a, epsilon(0.0_dp)*1.e4_dp)
2245 tauv(2, 1) = max(my_tau_b, epsilon(0.0_dp)*1.e4_dp)
2246 tauv(1, 1) = max(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2247 tauv(2, 1) = max(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2249 CALL xc_f03_mgga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2250 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2251 vlapl(1, 1), vtau(1, 1), &
2252 v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), &
2253 v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2254 v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2256 CALL xc_f03_mgga(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2257 laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2258 vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), &
2259 v2rhotau(1, 1), v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2260 v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2262 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2263 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2264 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2265 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2266 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2267 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2268 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2269 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2270 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2271 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2272 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2273 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2274 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2275 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2276 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2277 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2278 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2279 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2280 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2281 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2282 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2283 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2284 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2285 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2286 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2287 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2288 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2289 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2290 e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1)
2291 e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1)
2292 e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1)
2293 e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1)
2294 e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho
2295 e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho
2296 e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2297 sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa
2298 e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2299 sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa
2300 e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2301 sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob
2302 e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2303 sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob
2304 e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1)
2305 e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1)
2306 e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1)
2307 IF (has_laplace)
THEN
2308 e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1)
2309 e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1)
2310 e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1)
2311 e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1)
2312 e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho
2313 e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho
2314 e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2315 sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa
2316 e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2317 sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa
2318 e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2319 sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob
2320 e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2321 sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob
2322 e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1)
2323 e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1)
2324 e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1)
2325 e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1)
2326 e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1)
2327 e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1)
2328 e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1)
2333 ELSE IF (grad_deriv == 2)
THEN
2336 my_rhoa = max(rhoa(ii), 0.0_dp)
2337 my_rhob = max(rhob(ii), 0.0_dp)
2338 my_tau_a = max(tau_a(ii), 0.0_dp)
2339 my_tau_b = max(tau_b(ii), 0.0_dp)
2340 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau))
THEN
2341 rhov(1, 1) = max(my_rhoa, epsilon(0.0_dp)*1.e4_dp)
2342 rhov(2, 1) = max(my_rhob, epsilon(0.0_dp)*1.e4_dp)
2343 my_norm_drhoa = max(norm_drhoa(ii), epsilon(0.0_dp)*1.e4_dp)
2344 my_norm_drhob = max(norm_drhob(ii), epsilon(0.0_dp)*1.e4_dp)
2345 my_norm_drho = max(norm_drho(ii), epsilon(0.0_dp)*1.e4_dp)
2346 sigmav(1, 1) = my_norm_drhoa**2
2347 sigmav(3, 1) = my_norm_drhob**2
2348 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1))
2349 laplace_rhov(1, 1) = laplace_rhoa(ii)
2350 laplace_rhov(2, 1) = laplace_rhob(ii)
2351 tauv(1, 1) = max(my_tau_a, epsilon(0.0_dp)*1.e4_dp)
2352 tauv(2, 1) = max(my_tau_b, epsilon(0.0_dp)*1.e4_dp)
2353 tauv(1, 1) = max(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1)))
2354 tauv(2, 1) = max(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1)))
2356 CALL xc_f03_mgga_vxc_fxc(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2357 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), &
2358 vlapl(1, 1), vtau(1, 1), &
2359 v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), &
2360 v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2361 v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2364 CALL xc_f03_mgga(xc_func, one, rhov(1, 1), sigmav(1, 1), &
2365 laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), &
2366 vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), &
2367 v2rhotau(1, 1), v2sigma2(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), &
2368 v2lapl2(1, 1), v2lapltau(1, 1), v2tau2(1, 1))
2370 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1))
2371 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1)
2372 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1)
2373 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho
2374 e_ndrhoa(ii) = e_ndrhoa(ii) + &
2375 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa
2376 e_ndrhob(ii) = e_ndrhob(ii) + &
2377 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob
2378 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1)
2379 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1)
2380 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1)
2381 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1)
2382 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1)
2383 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho
2384 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho
2385 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + &
2386 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa
2387 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + &
2388 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa
2389 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + &
2390 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob
2391 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + &
2392 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob
2393 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + &
2394 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1))
2395 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + &
2396 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa
2397 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + &
2398 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob
2399 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + &
2400 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( &
2401 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1)))
2402 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + &
2403 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - &
2404 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob
2405 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + &
2406 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( &
2407 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1)))
2408 e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1)
2409 e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1)
2410 e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1)
2411 e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1)
2412 e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho
2413 e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho
2414 e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + &
2415 sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa
2416 e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + &
2417 sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa
2418 e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + &
2419 sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob
2420 e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + &
2421 sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob
2422 e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1)
2423 e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1)
2424 e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1)
2425 IF (has_laplace)
THEN
2426 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1)
2427 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1)
2428 e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1)
2429 e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1)
2430 e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1)
2431 e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1)
2432 e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho
2433 e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho
2434 e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + &
2435 sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa
2436 e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + &
2437 sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa
2438 e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + &
2439 sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob
2440 e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + &
2441 sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob
2442 e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1)
2443 e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1)
2444 e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1)
2445 e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1)
2446 e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1)
2447 e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1)
2448 e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1)
2455 cpabort(trim(func_name)//
": this XC_FAMILY is currently not supported.")
2458 END SUBROUTINE libxc_lsd_calc
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public marques2012
integer, save, public lehtola2018
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_norm_drho
integer, parameter, public deriv_laplace_rhob
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_tau
integer, parameter, public deriv_tau_b
integer, parameter, public deriv_tau_a
integer, parameter, public deriv_laplace_rhoa
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
integer, parameter, public deriv_laplace_rho
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
Includes all necessary routines, functions and parameters from libxc. Provides CP2K routines/function...
calculates a functional from libxc and its derivatives
subroutine, public libxc_lda_eval(rho_set, deriv_set, grad_deriv, libxc_params)
evaluates the functional from libxc
subroutine, public libxc_lsd_eval(rho_set, deriv_set, grad_deriv, libxc_params)
evaluates the functional from libxc
subroutine, public libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
info about the functional from libxc
subroutine, public libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn)
info about the functional from libxc
subroutine, public libxc_add_sections(section)
...
logical function, public libxc_check_existence_in_libxc(libxc_params)
This function checks whether a functional name belongs to LibXC.
integer function, public libxc_get_reference_length(libxc_params, lsd)
This function returns the maximum length of the reference string for a given LibXC functional.
subroutine, public libxc_version_info(version)
info about the LibXC version
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
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation