1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief calculates a functional from libxc and its derivatives 8!> \note 9!> LibXC: 10!> (Marques, Oliveira, Burnus, CPC 183, 2272 (2012)). 11!> 12!> WARNING: In the subroutine libxc_lsd_calc, it could be that the 13!> ordering for the 1st index of v2lapltau, v2rholapl, v2rhotau, 14!> v2sigmalapl and v2sigmatau is not correct. For the moment it does not 15!> matter since the calculation of the 2nd derivatives for meta-GGA 16!> functionals is not implemented in CP2K. 17!> 18!> \par History 19!> 01.2013 created [F. Tran] 20!> 07.2014 updates to versions 2.1 [JGH] 21!> 08.2015 refactoring [A. Gloess (agloess)] 22!> 01.2018 refactoring [A. Gloess (agloess)] 23!> 10.2018/04.2019 added hyb_mgga [S. Simko, included by F. Stein] 24!> \author F. Tran 25! ************************************************************************************************** 26MODULE xc_libxc 27 USE bibliography, ONLY: Lehtola2018, & 28 Marques2012, & 29 cite_reference 30 31 USE input_section_types, ONLY: section_vals_type, & 32 section_vals_val_get 33 USE kinds, ONLY: default_string_length, & 34 dp 35 USE xc_derivative_set_types, ONLY: xc_derivative_set_type, & 36 xc_dset_get_derivative 37 USE xc_derivative_types, ONLY: xc_derivative_get, & 38 xc_derivative_type 39 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 40 USE xc_rho_set_types, ONLY: xc_rho_set_get, & 41 xc_rho_set_type 42#if defined (__LIBXC) 43 USE xc_libxc_wrap, ONLY: xc_f03_func_t, & 44 xc_f03_func_init, & 45 xc_f03_func_end, & 46 xc_f03_func_info_t, & 47 xc_f03_func_get_info, & 48 xc_f03_func_info_get_family, & 49 xc_f03_func_info_get_kind, & 50 xc_f03_func_info_get_name, & 51 xc_f03_gga_exc, & 52 xc_f03_gga_exc_vxc, & 53 xc_f03_gga_fxc, & 54 xc_f03_gga_vxc, & 55 xc_f03_lda, & 56 xc_f03_lda_exc, & 57 xc_f03_lda_exc_vxc, & 58 xc_f03_lda_fxc, & 59 xc_f03_lda_kxc, & 60 xc_f03_lda_vxc, & 61 xc_f03_mgga, & 62 xc_f03_mgga_exc, & 63 xc_f03_mgga_exc_vxc, & 64 xc_f03_mgga_fxc, & 65 xc_f03_mgga_vxc, & 66 XC_POLARIZED, & 67 XC_UNPOLARIZED, & 68 XC_FAMILY_LDA, & 69 XC_FAMILY_GGA, & 70 XC_FAMILY_MGGA, & 71 XC_FAMILY_HYB_GGA, & 72 XC_FAMILY_HYB_MGGA, & 73 XC_CORRELATION, & 74 XC_EXCHANGE, & 75 XC_EXCHANGE_CORRELATION, & 76 XC_KINETIC, & 77 xc_libxc_wrap_info_refs, & 78 xc_libxc_wrap_version, & 79 xc_libxc_wrap_functional_get_number, & 80 xc_libxc_wrap_needs_laplace, & 81 xc_libxc_wrap_functional_set_params, & 82 xc_libxc_wrap_is_under_development 83#endif 84 85#include "../base/base_uses.f90" 86 87 IMPLICIT NONE 88 PRIVATE 89 90 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_libxc' 91 92 PUBLIC :: libxc_lda_info, libxc_lda_eval, libxc_lsd_info, libxc_lsd_eval, & 93 libxc_version_info 94 95CONTAINS 96 97! ************************************************************************************************** 98!> \brief info about the functional from libxc 99!> \param libxc_params input parameter (functional name, scaling and parameters) 100!> \param reference string with the reference of the actual functional 101!> \param shortform string with the shortform of the functional name 102!> \param needs the components needed by this functional are set to 103!> true (does not set the unneeded components to false) 104!> \param max_deriv maximum implemented derivative of the xc functional 105!> \param print_warn whether to print warning about development status of a functional 106!> \author F. Tran 107! ************************************************************************************************** 108 SUBROUTINE libxc_lda_info(libxc_params, reference, shortform, needs, max_deriv, print_warn) 109 110 TYPE(section_vals_type), POINTER :: libxc_params 111 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 112 TYPE(xc_rho_cflags_type), & 113 INTENT(inout), OPTIONAL :: needs 114 INTEGER, INTENT(out), OPTIONAL :: max_deriv 115 LOGICAL, INTENT(IN), OPTIONAL :: print_warn 116 117 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lda_info', & 118 routineP = moduleN//':'//routineN 119 120#if defined (__LIBXC) 121 CHARACTER(LEN=128) :: s1, s2 122 CHARACTER(LEN=default_string_length) :: func_name 123 INTEGER :: func_id 124 REAL(KIND=dp) :: func_scale 125 TYPE(xc_f03_func_t) :: xc_func 126 TYPE(xc_f03_func_info_t) :: xc_info 127 128 CALL section_vals_val_get(libxc_params, "functional", c_val=func_name) 129 CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale) 130 131 CALL cite_reference(Marques2012) 132 CALL cite_reference(Lehtola2018) 133 134 IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp 135 136 func_id = xc_libxc_wrap_functional_get_number(func_name) 137!$OMP CRITICAL(libxc_init) 138 CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED) 139 xc_info = xc_f03_func_get_info(xc_func) 140!$OMP END CRITICAL(libxc_init) 141!$OMP BARRIER 142 143 s1 = xc_f03_func_info_get_name(xc_info) 144 SELECT CASE (xc_f03_func_info_get_kind(xc_info)) 145 CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange" 146 CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation" 147 CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation" 148 CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic" 149 CASE default 150 CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.") 151 END SELECT 152 IF (PRESENT(shortform)) THEN 153 shortform = TRIM(s1)//' ('//TRIM(s2)//')' 154 END IF 155 IF (PRESENT(reference)) THEN 156 CALL xc_libxc_wrap_info_refs(xc_info, XC_UNPOLARIZED, func_scale, reference) 157 END IF 158 IF (PRESENT(needs)) THEN 159 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 160 CASE (XC_FAMILY_LDA) 161 needs%rho = .TRUE. 162 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 163 needs%rho = .TRUE. 164 needs%norm_drho = .TRUE. 165 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 166 needs%rho = .TRUE. 167 needs%norm_drho = .TRUE. 168 needs%tau = .TRUE. 169 needs%laplace_rho = xc_libxc_wrap_needs_laplace(func_id) 170 CASE default 171 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 172 END SELECT 173 END IF 174 IF (PRESENT(max_deriv)) THEN 175 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 176 CASE (XC_FAMILY_LDA) 177 max_deriv = 3 178 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 179 max_deriv = 2 180 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 181 max_deriv = 1 182 CASE default 183 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 184 END SELECT 185 END IF 186 IF (PRESENT(print_warn)) THEN 187 IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN 188 CPWARN(TRIM(func_name)//" is under development. Use with caution.") 189 END IF 190 END IF 191 192 CALL xc_f03_func_end(xc_func) 193#else 194 MARK_USED(libxc_params) 195 MARK_USED(reference) 196 MARK_USED(shortform) 197 MARK_USED(needs) 198 MARK_USED(max_deriv) 199 MARK_USED(print_warn) 200 CPABORT("In order to use libxc you need to download and install it") 201#endif 202 203 END SUBROUTINE libxc_lda_info 204 205! ************************************************************************************************** 206!> \brief info about the functional from libxc 207!> \param libxc_params input parameter (functional name, scaling and parameters) 208!> \param reference string with the reference of the actual functional 209!> \param shortform string with the shortform of the functional name 210!> \param needs the components needed by this functional are set to 211!> true (does not set the unneeded components to false) 212!> \param max_deriv maximum implemented derivative of the xc functional 213!> \param print_warn whether to print warning about development status of a functional 214!> \author F. Tran 215! ************************************************************************************************** 216 SUBROUTINE libxc_lsd_info(libxc_params, reference, shortform, needs, max_deriv, print_warn) 217 218 TYPE(section_vals_type), POINTER :: libxc_params 219 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 220 TYPE(xc_rho_cflags_type), & 221 INTENT(inout), OPTIONAL :: needs 222 INTEGER, INTENT(out), OPTIONAL :: max_deriv 223 LOGICAL, INTENT(IN), OPTIONAL :: print_warn 224 225 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lsd_info', & 226 routineP = moduleN//':'//routineN 227 228#if defined (__LIBXC) 229 CHARACTER(LEN=128) :: s1, s2 230 CHARACTER(LEN=default_string_length) :: func_name 231 INTEGER :: func_id 232 REAL(KIND=dp) :: func_scale 233 TYPE(xc_f03_func_t) :: xc_func 234 TYPE(xc_f03_func_info_t) :: xc_info 235 236 CALL section_vals_val_get(libxc_params, "functional", c_val=func_name) 237 CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale) 238 239 CALL cite_reference(Marques2012) 240 CALL cite_reference(Lehtola2018) 241 242 IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp 243 244 func_id = xc_libxc_wrap_functional_get_number(func_name) 245!$OMP CRITICAL(libxc_init) 246 CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED) 247 xc_info = xc_f03_func_get_info(xc_func) 248!$OMP END CRITICAL(libxc_init) 249!$OMP BARRIER 250 251 s1 = xc_f03_func_info_get_name(xc_info) 252 SELECT CASE (xc_f03_func_info_get_kind(xc_info)) 253 CASE (XC_EXCHANGE); WRITE (s2, '(a)') "exchange" 254 CASE (XC_CORRELATION); WRITE (s2, '(a)') "correlation" 255 CASE (XC_EXCHANGE_CORRELATION); WRITE (s2, '(a)') "exchange-correlation" 256 CASE (XC_KINETIC); WRITE (s2, '(a)') "kinetic" 257 CASE default 258 CPABORT(TRIM(func_name)//": this XC_KIND is currently not supported.") 259 END SELECT 260 IF (PRESENT(shortform)) THEN 261 shortform = TRIM(s1)//' ('//TRIM(s2)//')' 262 END IF 263 IF (PRESENT(reference)) THEN 264 CALL xc_libxc_wrap_info_refs(xc_info, XC_POLARIZED, func_scale, reference) 265 END IF 266 IF (PRESENT(needs)) THEN 267 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 268 CASE (XC_FAMILY_LDA) 269 needs%rho_spin = .TRUE. 270 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 271 needs%rho_spin = .TRUE. 272 needs%norm_drho = .TRUE. 273 needs%norm_drho_spin = .TRUE. 274 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 275 needs%rho_spin = .TRUE. 276 needs%norm_drho = .TRUE. 277 needs%norm_drho_spin = .TRUE. 278 needs%tau_spin = .TRUE. 279 needs%laplace_rho_spin = xc_libxc_wrap_needs_laplace(func_id) 280 CASE default 281 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 282 END SELECT 283 END IF 284 IF (PRESENT(max_deriv)) THEN 285 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 286 CASE (XC_FAMILY_LDA) 287 max_deriv = 3 288 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 289 max_deriv = 2 290 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 291 max_deriv = 1 292 CASE default 293 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 294 END SELECT 295 END IF 296 IF (PRESENT(print_warn)) THEN 297 IF (print_warn .AND. xc_libxc_wrap_is_under_development(xc_info)) THEN 298 CPWARN(TRIM(func_name)//" is under development. Use with caution.") 299 END IF 300 END IF 301 302 CALL xc_f03_func_end(xc_func) 303#else 304 MARK_USED(libxc_params) 305 MARK_USED(reference) 306 MARK_USED(shortform) 307 MARK_USED(needs) 308 MARK_USED(max_deriv) 309 MARK_USED(print_warn) 310 CPABORT("In order to use libxc you need to download and install it") 311#endif 312 313 END SUBROUTINE libxc_lsd_info 314 315! ************************************************************************************************** 316!> \brief info about the LibXC version 317!> \param version ... 318!> \author A. Gloess (agloess) 319! ************************************************************************************************** 320 SUBROUTINE libxc_version_info(version) 321 CHARACTER(LEN=*), INTENT(OUT) :: version ! the string that is output 322 323 CHARACTER(LEN=*), PARAMETER :: routineN = 'libxc_version_info', & 324 routineP = moduleN//':'//routineN 325 326#if defined (__LIBXC) 327 CALL xc_libxc_wrap_version(version) 328#else 329 version = "none" 330 CPABORT("In order to use libxc you need to download and install it") 331#endif 332 333 END SUBROUTINE libxc_version_info 334 335! ************************************************************************************************** 336!> \brief evaluates the functional from libxc 337!> \param rho_set the density where you want to evaluate the functional 338!> \param deriv_set place where to store the functional derivatives (they are 339!> added to the derivatives) 340!> \param grad_deriv degree of the derivative that should be evaluated, 341!> if positive all the derivatives up to the given degree are evaluated, 342!> if negative only the given degree is calculated 343!> \param libxc_params input parameter (functional name, scaling and parameters) 344!> \author F. Tran 345! ************************************************************************************************** 346 SUBROUTINE libxc_lda_eval(rho_set, deriv_set, grad_deriv, libxc_params) 347 348 TYPE(xc_rho_set_type), POINTER :: rho_set 349 TYPE(xc_derivative_set_type), POINTER :: deriv_set 350 INTEGER, INTENT(in) :: grad_deriv 351 TYPE(section_vals_type), POINTER :: libxc_params 352 353 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lda_eval', & 354 routineP = moduleN//':'//routineN 355 356#if defined (__LIBXC) 357 CHARACTER(LEN=default_string_length) :: func_name 358 INTEGER :: func_id, handle, npoints 359 INTEGER, DIMENSION(:, :), POINTER :: bo 360 LOGICAL :: has_laplace 361 REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale 362 REAL(KIND=dp), DIMENSION(:), POINTER :: params 363 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rho, & 364 e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_ndrho, & 365 e_ndrho_laplace_rho, e_ndrho_ndrho, e_ndrho_rho, e_ndrho_tau, e_rho, & 366 e_rho_laplace_rho, e_rho_rho, e_rho_rho_rho, e_rho_tau, e_tau, & 367 e_tau_tau, laplace_rho, norm_drho, rho, tau 368 TYPE(xc_derivative_type), POINTER :: deriv 369 TYPE(xc_f03_func_t) :: xc_func 370 TYPE(xc_f03_func_info_t) :: xc_info 371 372 CALL timeset(routineN, handle) 373 374 has_laplace = .FALSE. 375 NULLIFY (bo, dummy) 376 NULLIFY (rho, norm_drho, laplace_rho, tau) 377 378 CPASSERT(ASSOCIATED(rho_set)) 379 CPASSERT(rho_set%ref_count > 0) 380 CPASSERT(ASSOCIATED(deriv_set)) 381 CPASSERT(deriv_set%ref_count > 0) 382 383 CALL section_vals_val_get(libxc_params, "functional", c_val=func_name) 384 CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale) 385 CALL section_vals_val_get(libxc_params, "parameters", r_vals=params) 386 387 IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp 388 389 func_id = xc_libxc_wrap_functional_get_number(func_name) 390!$OMP CRITICAL(libxc_init) 391 CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED) 392 xc_info = xc_f03_func_get_info(xc_func) 393!$OMP END CRITICAL(libxc_init) 394!$OMP BARRIER 395 396 CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., & 397 rho=rho, norm_drho=norm_drho, laplace_rho=laplace_rho, & 398 rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, & 399 tau=tau, local_bounds=bo) 400 401 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 402 403 dummy => rho 404 405 ! due to assumed shape array usage in next routine 406 IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => dummy 407 IF (.NOT. ASSOCIATED(tau)) tau => dummy 408 409 ! only some MGGA functionals really need the Laplacian, 410 ! all others can work with rho (read-only) as dummy 411 IF (ASSOCIATED(laplace_rho)) has_laplace = .TRUE. 412 IF (.NOT. has_laplace) laplace_rho => dummy 413 414 e_0 => dummy 415 e_rho => dummy 416 e_ndrho => dummy 417 e_laplace_rho => dummy 418 e_tau => dummy 419 e_rho_rho => dummy 420 e_ndrho_rho => dummy 421 e_ndrho_ndrho => dummy 422 e_rho_laplace_rho => dummy 423 e_rho_tau => dummy 424 e_ndrho_laplace_rho => dummy 425 e_ndrho_tau => dummy 426 e_laplace_rho_laplace_rho => dummy 427 e_laplace_rho_tau => dummy 428 e_tau_tau => dummy 429 e_rho_rho_rho => dummy 430 431 IF (grad_deriv >= 0) THEN 432 deriv => xc_dset_get_derivative(deriv_set, "", & 433 allocate_deriv=.TRUE.) 434 CALL xc_derivative_get(deriv, deriv_data=e_0) 435 END IF 436 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 437 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 438 CASE (XC_FAMILY_LDA) 439 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 440 allocate_deriv=.TRUE.) 441 CALL xc_derivative_get(deriv, deriv_data=e_rho) 442 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 443 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 444 allocate_deriv=.TRUE.) 445 CALL xc_derivative_get(deriv, deriv_data=e_rho) 446 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 447 allocate_deriv=.TRUE.) 448 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 449 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 450 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 451 allocate_deriv=.TRUE.) 452 CALL xc_derivative_get(deriv, deriv_data=e_rho) 453 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 454 allocate_deriv=.TRUE.) 455 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 456 deriv => xc_dset_get_derivative(deriv_set, "(tau)", & 457 allocate_deriv=.TRUE.) 458 CALL xc_derivative_get(deriv, deriv_data=e_tau) 459 IF (has_laplace) THEN 460 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rho)", & 461 allocate_deriv=.TRUE.) 462 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho) 463 END IF 464 CASE default 465 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 466 END SELECT 467 END IF 468 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 469 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 470 CASE (XC_FAMILY_LDA) 471 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)", & 472 allocate_deriv=.TRUE.) 473 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho) 474 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 475 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)", & 476 allocate_deriv=.TRUE.) 477 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho) 478 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rho)", & 479 allocate_deriv=.TRUE.) 480 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho) 481 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drho)", & 482 allocate_deriv=.TRUE.) 483 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho) 484 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 485 ! not implemented ... 486 CPABORT("derivatives larger than 1 not implemented or checked") 487 488 ! deriv => xc_dset_get_derivative(deriv_set,"(rho)(rho)",& 489 ! allocate_deriv=.TRUE.) 490 ! CALL xc_derivative_get(deriv,deriv_data=e_rho_rho) 491 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rho)",& 492 ! allocate_deriv=.TRUE.) 493 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rho) 494 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(norm_drho)",& 495 ! allocate_deriv=.TRUE.) 496 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho) 497 ! deriv => xc_dset_get_derivative(deriv_set,"(rho)(tau)",& 498 ! allocate_deriv=.TRUE.) 499 ! CALL xc_derivative_get(deriv,deriv_data=e_rho_tau) 500 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(tau)",& 501 ! allocate_deriv=.TRUE.) 502 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_tau) 503 ! deriv => xc_dset_get_derivative(deriv_set,"(tau)(tau)",& 504 ! allocate_deriv=.TRUE.) 505 ! CALL xc_derivative_get(deriv,deriv_data=e_tau_tau) 506 ! IF (has_laplace) THEN 507 ! deriv => xc_dset_get_derivative(deriv_set,"(rho)(laplace_rho)",& 508 ! allocate_deriv=.TRUE.) 509 ! CALL xc_derivative_get(deriv,deriv_data=e_rho_laplace_rho) 510 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_rho)(laplace_rho)",& 511 ! allocate_deriv=.TRUE.) 512 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_laplace_rho) 513 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rho)(laplace_rho)",& 514 ! allocate_deriv=.TRUE.) 515 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rho_laplace_rho) 516 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rho)(tau)",& 517 ! allocate_deriv=.TRUE.) 518 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rho_tau) 519 ! END IF 520 CASE default 521 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 522 END SELECT 523 END IF 524 IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN 525 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 526 CASE (XC_FAMILY_LDA) 527 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)(rho)", & 528 allocate_deriv=.TRUE.) 529 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho) 530 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 531 CPABORT("derivatives larger than 2 not implemented") 532 CASE default 533 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 534 END SELECT 535 END IF 536 IF (grad_deriv >= 4 .OR. grad_deriv <= -4) THEN 537 CPABORT("derivatives larger than 3 not implemented") 538 END IF 539 540!$OMP PARALLEL DEFAULT(NONE), & 541!$OMP SHARED(rho,norm_drho,laplace_rho,tau,e_0,e_rho,e_ndrho,e_laplace_rho),& 542!$OMP SHARED(e_tau,e_rho_rho,e_ndrho_rho,e_ndrho_ndrho,e_rho_laplace_rho),& 543!$OMP SHARED(e_rho_tau,e_ndrho_laplace_rho,e_ndrho_tau,e_laplace_rho_laplace_rho),& 544!$OMP SHARED(e_laplace_rho_tau,e_tau_tau,e_rho_rho_rho),& 545!$OMP SHARED(grad_deriv,npoints),& 546!$OMP SHARED(epsilon_rho,epsilon_tau),& 547!$OMP SHARED(func_name,func_scale,params) 548 549 CALL libxc_lda_calc(rho=rho, norm_drho=norm_drho, & 550 laplace_rho=laplace_rho, tau=tau, & 551 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_laplace_rho=e_laplace_rho, & 552 e_tau=e_tau, e_rho_rho=e_rho_rho, e_ndrho_rho=e_ndrho_rho, & 553 e_ndrho_ndrho=e_ndrho_ndrho, e_rho_laplace_rho=e_rho_laplace_rho, & 554 e_rho_tau=e_rho_tau, e_ndrho_laplace_rho=e_ndrho_laplace_rho, & 555 e_ndrho_tau=e_ndrho_tau, e_laplace_rho_laplace_rho=e_laplace_rho_laplace_rho, & 556 e_laplace_rho_tau=e_laplace_rho_tau, e_tau_tau=e_tau_tau, & 557 e_rho_rho_rho=e_rho_rho_rho, & 558 grad_deriv=grad_deriv, npoints=npoints, & 559 epsilon_rho=epsilon_rho, & 560 epsilon_tau=epsilon_tau, func_name=func_name, & 561 sc=func_scale, params=params) 562 563!$OMP END PARALLEL 564 565 NULLIFY (dummy) 566 567 CALL xc_f03_func_end(xc_func) 568 569 CALL timestop(handle) 570#else 571 MARK_USED(rho_set) 572 MARK_USED(deriv_set) 573 MARK_USED(grad_deriv) 574 MARK_USED(libxc_params) 575 CPABORT("In order to use libxc you need to download and install it") 576#endif 577 END SUBROUTINE libxc_lda_eval 578 579! ************************************************************************************************** 580!> \brief evaluates the functional from libxc 581!> \param rho_set the density where you want to evaluate the functional 582!> \param deriv_set place where to store the functional derivatives (they are 583!> added to the derivatives) 584!> \param grad_deriv degree of the derivative that should be evaluated, 585!> if positive all the derivatives up to the given degree are evaluated, 586!> if negative only the given degree is calculated 587!> \param libxc_params input parameter (functional name, scaling and parameters) 588!> \author F. Tran 589! ************************************************************************************************** 590 SUBROUTINE libxc_lsd_eval(rho_set, deriv_set, grad_deriv, libxc_params) 591 592 TYPE(xc_rho_set_type), POINTER :: rho_set 593 TYPE(xc_derivative_set_type), POINTER :: deriv_set 594 INTEGER, INTENT(in) :: grad_deriv 595 TYPE(section_vals_type), POINTER :: libxc_params 596 597 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lsd_eval', & 598 routineP = moduleN//':'//routineN 599 600#if defined (__LIBXC) 601 CHARACTER(LEN=default_string_length) :: func_name 602 INTEGER :: func_id, handle, npoints 603 INTEGER, DIMENSION(:, :), POINTER :: bo 604 LOGICAL :: has_laplace 605 REAL(KIND=dp) :: epsilon_rho, epsilon_tau, func_scale 606 REAL(KIND=dp), DIMENSION(:), POINTER :: params 607 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, & 608 e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, & 609 e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, e_laplace_rhob, & 610 e_laplace_rhob_laplace_rhob, e_laplace_rhob_tau_a, & 611 e_laplace_rhob_tau_b, e_ndrho, e_ndrho_laplace_rhoa, & 612 e_ndrho_laplace_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, & 613 e_ndrho_rhoa, e_ndrho_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa, & 614 e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, e_ndrhoa_ndrhoa, & 615 e_ndrhoa_ndrhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhoa_tau_a, & 616 e_ndrhoa_tau_b, e_ndrhob 617 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_ndrhob_laplace_rhoa, & 618 e_ndrhob_laplace_rhob, e_ndrhob_ndrhob, e_ndrhob_rhoa, e_ndrhob_rhob, & 619 e_ndrhob_tau_a, e_ndrhob_tau_b, e_rhoa, e_rhoa_laplace_rhoa, & 620 e_rhoa_laplace_rhob, e_rhoa_rhoa, e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, & 621 e_rhoa_rhob, e_rhoa_rhob_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob, & 622 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhob_rhob, & 623 e_rhob_rhob_rhob, e_rhob_tau_a, e_rhob_tau_b, e_tau_a, e_tau_a_tau_a, & 624 e_tau_a_tau_b, e_tau_b, e_tau_b_tau_b, laplace_rhoa, laplace_rhob, & 625 norm_drho, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b 626 TYPE(xc_derivative_type), POINTER :: deriv 627 TYPE(xc_f03_func_t) :: xc_func 628 TYPE(xc_f03_func_info_t) :: xc_info 629 630 CALL timeset(routineN, handle) 631 632 has_laplace = .FALSE. 633 NULLIFY (bo, dummy) 634 NULLIFY (rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, laplace_rhoa, & 635 laplace_rhob, tau_a, tau_b) 636 637 CPASSERT(ASSOCIATED(rho_set)) 638 CPASSERT(rho_set%ref_count > 0) 639 CPASSERT(ASSOCIATED(deriv_set)) 640 CPASSERT(deriv_set%ref_count > 0) 641 642 CALL section_vals_val_get(libxc_params, "functional", c_val=func_name) 643 CALL section_vals_val_get(libxc_params, "scale", r_val=func_scale) 644 CALL section_vals_val_get(libxc_params, "parameters", r_vals=params) 645 646 IF (ABS(func_scale - 1.0_dp) < 1.0e-10_dp) func_scale = 1.0_dp 647 648 func_id = xc_libxc_wrap_functional_get_number(func_name) 649!$OMP CRITICAL(libxc_init) 650 CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED) 651 xc_info = xc_f03_func_get_info(xc_func) 652!$OMP END CRITICAL(libxc_init) 653!$OMP BARRIER 654 655 CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., & 656 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, & 657 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, & 658 laplace_rhoa=laplace_rhoa, laplace_rhob=laplace_rhob, & 659 rho_cutoff=epsilon_rho, tau_cutoff=epsilon_tau, & 660 tau_a=tau_a, tau_b=tau_b, local_bounds=bo) 661 662 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 663 664 dummy => rhoa 665 666 ! due to assumed shape array usage in next routine 667 IF (.NOT. ASSOCIATED(norm_drho)) norm_drho => dummy 668 IF (.NOT. ASSOCIATED(norm_drhoa)) norm_drhoa => dummy 669 IF (.NOT. ASSOCIATED(norm_drhob)) norm_drhob => dummy 670 IF (.NOT. ASSOCIATED(tau_a)) tau_a => dummy 671 IF (.NOT. ASSOCIATED(tau_b)) tau_b => dummy 672 673 ! only some MGGA functionals really need the Laplacian, 674 ! all others can work with rhoa (read-only) as dummy 675 IF (ASSOCIATED(laplace_rhoa) .AND. ASSOCIATED(laplace_rhob)) has_laplace = .TRUE. 676 IF (.NOT. has_laplace) laplace_rhoa => dummy 677 IF (.NOT. has_laplace) laplace_rhob => dummy 678 679 e_0 => dummy 680 e_rhoa => dummy 681 e_rhob => dummy 682 e_ndrho => dummy 683 e_ndrhoa => dummy 684 e_ndrhob => dummy 685 e_laplace_rhoa => dummy 686 e_laplace_rhob => dummy 687 e_tau_a => dummy 688 e_tau_b => dummy 689 e_rhoa_rhoa => dummy 690 e_rhoa_rhob => dummy 691 e_rhob_rhob => dummy 692 e_ndrho_rhoa => dummy 693 e_ndrho_rhob => dummy 694 e_ndrhoa_rhoa => dummy 695 e_ndrhoa_rhob => dummy 696 e_ndrhob_rhoa => dummy 697 e_ndrhob_rhob => dummy 698 e_ndrho_ndrho => dummy 699 e_ndrho_ndrhoa => dummy 700 e_ndrho_ndrhob => dummy 701 e_ndrhoa_ndrhoa => dummy 702 e_ndrhoa_ndrhob => dummy 703 e_ndrhob_ndrhob => dummy 704 e_rhoa_laplace_rhoa => dummy 705 e_rhoa_laplace_rhob => dummy 706 e_rhob_laplace_rhoa => dummy 707 e_rhob_laplace_rhob => dummy 708 e_rhoa_tau_a => dummy 709 e_rhoa_tau_b => dummy 710 e_rhob_tau_a => dummy 711 e_rhob_tau_b => dummy 712 e_ndrho_laplace_rhoa => dummy 713 e_ndrho_laplace_rhob => dummy 714 e_ndrhoa_laplace_rhoa => dummy 715 e_ndrhoa_laplace_rhob => dummy 716 e_ndrhob_laplace_rhoa => dummy 717 e_ndrhob_laplace_rhob => dummy 718 e_ndrho_tau_a => dummy 719 e_ndrho_tau_b => dummy 720 e_ndrhoa_tau_a => dummy 721 e_ndrhoa_tau_b => dummy 722 e_ndrhob_tau_a => dummy 723 e_ndrhob_tau_b => dummy 724 e_laplace_rhoa_laplace_rhoa => dummy 725 e_laplace_rhoa_laplace_rhob => dummy 726 e_laplace_rhob_laplace_rhob => dummy 727 e_laplace_rhoa_tau_a => dummy 728 e_laplace_rhoa_tau_b => dummy 729 e_laplace_rhob_tau_a => dummy 730 e_laplace_rhob_tau_b => dummy 731 e_tau_a_tau_a => dummy 732 e_tau_a_tau_b => dummy 733 e_tau_b_tau_b => dummy 734 e_rhoa_rhoa_rhoa => dummy 735 e_rhoa_rhoa_rhob => dummy 736 e_rhoa_rhob_rhob => dummy 737 e_rhob_rhob_rhob => dummy 738 739 IF (grad_deriv >= 0) THEN 740 deriv => xc_dset_get_derivative(deriv_set, "", & 741 allocate_deriv=.TRUE.) 742 CALL xc_derivative_get(deriv, deriv_data=e_0) 743 END IF 744 IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN 745 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 746 CASE (XC_FAMILY_LDA) 747 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 748 allocate_deriv=.TRUE.) 749 CALL xc_derivative_get(deriv, deriv_data=e_rhoa) 750 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 751 allocate_deriv=.TRUE.) 752 CALL xc_derivative_get(deriv, deriv_data=e_rhob) 753 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 754 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 755 allocate_deriv=.TRUE.) 756 CALL xc_derivative_get(deriv, deriv_data=e_rhoa) 757 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 758 allocate_deriv=.TRUE.) 759 CALL xc_derivative_get(deriv, deriv_data=e_rhob) 760 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 761 allocate_deriv=.TRUE.) 762 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 763 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 764 allocate_deriv=.TRUE.) 765 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa) 766 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 767 allocate_deriv=.TRUE.) 768 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob) 769 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 770 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 771 allocate_deriv=.TRUE.) 772 CALL xc_derivative_get(deriv, deriv_data=e_rhoa) 773 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 774 allocate_deriv=.TRUE.) 775 CALL xc_derivative_get(deriv, deriv_data=e_rhob) 776 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 777 allocate_deriv=.TRUE.) 778 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 779 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 780 allocate_deriv=.TRUE.) 781 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa) 782 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 783 allocate_deriv=.TRUE.) 784 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob) 785 deriv => xc_dset_get_derivative(deriv_set, "(tau_a)", & 786 allocate_deriv=.TRUE.) 787 CALL xc_derivative_get(deriv, deriv_data=e_tau_a) 788 deriv => xc_dset_get_derivative(deriv_set, "(tau_b)", & 789 allocate_deriv=.TRUE.) 790 CALL xc_derivative_get(deriv, deriv_data=e_tau_b) 791 IF (has_laplace) THEN 792 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rhoa)", & 793 allocate_deriv=.TRUE.) 794 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa) 795 deriv => xc_dset_get_derivative(deriv_set, "(laplace_rhob)", & 796 allocate_deriv=.TRUE.) 797 CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob) 798 END IF 799 CASE default 800 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 801 END SELECT 802 END IF 803 IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN 804 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 805 CASE (XC_FAMILY_LDA) 806 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)", & 807 allocate_deriv=.TRUE.) 808 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa) 809 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhob)", & 810 allocate_deriv=.TRUE.) 811 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob) 812 deriv => xc_dset_get_derivative(deriv_set, "(rhob)(rhob)", & 813 allocate_deriv=.TRUE.) 814 CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob) 815 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 816 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)", & 817 allocate_deriv=.TRUE.) 818 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa) 819 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhob)", & 820 allocate_deriv=.TRUE.) 821 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob) 822 deriv => xc_dset_get_derivative(deriv_set, "(rhob)(rhob)", & 823 allocate_deriv=.TRUE.) 824 CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob) 825 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rhoa)", & 826 allocate_deriv=.TRUE.) 827 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhoa) 828 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(rhob)", & 829 allocate_deriv=.TRUE.) 830 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rhob) 831 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)(rhoa)", & 832 allocate_deriv=.TRUE.) 833 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhoa) 834 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)(rhob)", & 835 allocate_deriv=.TRUE.) 836 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_rhob) 837 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)(rhoa)", & 838 allocate_deriv=.TRUE.) 839 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhoa) 840 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)(rhob)", & 841 allocate_deriv=.TRUE.) 842 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_rhob) 843 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drho)", & 844 allocate_deriv=.TRUE.) 845 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho) 846 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drhoa)", & 847 allocate_deriv=.TRUE.) 848 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhoa) 849 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drhob)", & 850 allocate_deriv=.TRUE.) 851 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrhob) 852 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)(norm_drhoa)", & 853 allocate_deriv=.TRUE.) 854 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhoa) 855 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)(norm_drhob)", & 856 allocate_deriv=.TRUE.) 857 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa_ndrhob) 858 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)(norm_drhob)", & 859 allocate_deriv=.TRUE.) 860 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob_ndrhob) 861 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 862 ! not implemented ... 863 CPABORT("derivatives larger than 1 not implemented or checked") 864 865 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhoa)",& 866 ! allocate_deriv=.TRUE.) 867 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_rhoa) 868 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(rhob)",& 869 ! allocate_deriv=.TRUE.) 870 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_rhob) 871 ! deriv => xc_dset_get_derivative(deriv_set,"(rhob)(rhob)",& 872 ! allocate_deriv=.TRUE.) 873 ! CALL xc_derivative_get(deriv,deriv_data=e_rhob_rhob) 874 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhoa)",& 875 ! allocate_deriv=.TRUE.) 876 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rhoa) 877 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(rhob)",& 878 ! allocate_deriv=.TRUE.) 879 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_rhob) 880 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(rhoa)",& 881 ! allocate_deriv=.TRUE.) 882 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_rhoa) 883 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(rhob)",& 884 ! allocate_deriv=.TRUE.) 885 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_rhob) 886 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(rhoa)",& 887 ! allocate_deriv=.TRUE.) 888 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_rhoa) 889 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(rhob)",& 890 ! allocate_deriv=.TRUE.) 891 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_rhob) 892 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(norm_drho)",& 893 ! allocate_deriv=.TRUE.) 894 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho) 895 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(norm_drhoa)",& 896 ! allocate_deriv=.TRUE.) 897 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrhoa) 898 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(norm_drhob)",& 899 ! allocate_deriv=.TRUE.) 900 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_ndrhob) 901 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(norm_drhoa)",& 902 ! allocate_deriv=.TRUE.) 903 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_ndrhoa) 904 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(norm_drhob)",& 905 ! allocate_deriv=.TRUE.) 906 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_ndrhob) 907 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(norm_drhob)",& 908 ! allocate_deriv=.TRUE.) 909 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_ndrhob) 910 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(tau_a)",& 911 ! allocate_deriv=.TRUE.) 912 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_tau_a) 913 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(tau_b)",& 914 ! allocate_deriv=.TRUE.) 915 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_tau_b) 916 ! deriv => xc_dset_get_derivative(deriv_set,"(rhob)(tau_a)",& 917 ! allocate_deriv=.TRUE.) 918 ! CALL xc_derivative_get(deriv,deriv_data=e_rhob_tau_a) 919 ! deriv => xc_dset_get_derivative(deriv_set,"(rhob)(tau_b)",& 920 ! allocate_deriv=.TRUE.) 921 ! CALL xc_derivative_get(deriv,deriv_data=e_rhob_tau_b) 922 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(tau_a)",& 923 ! allocate_deriv=.TRUE.) 924 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_tau_a) 925 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(tau_b)",& 926 ! allocate_deriv=.TRUE.) 927 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_tau_b) 928 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(tau_a)",& 929 ! allocate_deriv=.TRUE.) 930 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_tau_a) 931 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(tau_b)",& 932 ! allocate_deriv=.TRUE.) 933 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_tau_b) 934 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(tau_a)",& 935 ! allocate_deriv=.TRUE.) 936 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_tau_a) 937 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(tau_b)",& 938 ! allocate_deriv=.TRUE.) 939 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_tau_b) 940 ! deriv => xc_dset_get_derivative(deriv_set,"(tau_a)(tau_a)",& 941 ! allocate_deriv=.TRUE.) 942 ! CALL xc_derivative_get(deriv,deriv_data=e_tau_a_tau_a) 943 ! deriv => xc_dset_get_derivative(deriv_set,"(tau_a)(tau_b)",& 944 ! allocate_deriv=.TRUE.) 945 ! CALL xc_derivative_get(deriv,deriv_data=e_tau_a_tau_b) 946 ! deriv => xc_dset_get_derivative(deriv_set,"(tau_b)(tau_b)",& 947 ! allocate_deriv=.TRUE.) 948 ! CALL xc_derivative_get(deriv,deriv_data=e_tau_b_tau_b) 949 ! IF (has_laplace) THEN 950 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(laplace_rhoa)",& 951 ! allocate_deriv=.TRUE.) 952 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_laplace_rhoa) 953 ! deriv => xc_dset_get_derivative(deriv_set,"(rhoa)(laplace_rhob)",& 954 ! allocate_deriv=.TRUE.) 955 ! CALL xc_derivative_get(deriv,deriv_data=e_rhoa_laplace_rhob) 956 ! deriv => xc_dset_get_derivative(deriv_set,"(rhob)(laplace_rhoa)",& 957 ! allocate_deriv=.TRUE.) 958 ! CALL xc_derivative_get(deriv,deriv_data=e_rhob_laplace_rhoa) 959 ! deriv => xc_dset_get_derivative(deriv_set,"(rhob)(laplace_rhob)",& 960 ! allocate_deriv=.TRUE.) 961 ! CALL xc_derivative_get(deriv,deriv_data=e_rhob_laplace_rhob) 962 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(laplace_rhoa)",& 963 ! allocate_deriv=.TRUE.) 964 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_laplace_rhoa) 965 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)(laplace_rhob)",& 966 ! allocate_deriv=.TRUE.) 967 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrho_laplace_rhob) 968 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(laplace_rhoa)",& 969 ! allocate_deriv=.TRUE.) 970 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_laplace_rhoa) 971 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)(laplace_rhob)",& 972 ! allocate_deriv=.TRUE.) 973 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhoa_laplace_rhob) 974 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(laplace_rhoa)",& 975 ! allocate_deriv=.TRUE.) 976 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_laplace_rhoa) 977 ! deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)(laplace_rhob)",& 978 ! allocate_deriv=.TRUE.) 979 ! CALL xc_derivative_get(deriv,deriv_data=e_ndrhob_laplace_rhob) 980 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhoa)(laplace_rhoa)",& 981 ! allocate_deriv=.TRUE.) 982 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhoa_laplace_rhoa) 983 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhoa)(laplace_rhob)",& 984 ! allocate_deriv=.TRUE.) 985 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhoa_laplace_rhob) 986 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhob)(laplace_rhob)",& 987 ! allocate_deriv=.TRUE.) 988 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhob_laplace_rhob) 989 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhoa)(tau_a)",& 990 ! allocate_deriv=.TRUE.) 991 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhoa_tau_a) 992 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhoa)(tau_b)",& 993 ! allocate_deriv=.TRUE.) 994 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhoa_tau_b) 995 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhob)(tau_a)",& 996 ! allocate_deriv=.TRUE.) 997 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhob_tau_a) 998 ! deriv => xc_dset_get_derivative(deriv_set,"(laplace_rhob)(tau_b)",& 999 ! allocate_deriv=.TRUE.) 1000 ! CALL xc_derivative_get(deriv,deriv_data=e_laplace_rhob_tau_b) 1001 ! END IF 1002 CASE default 1003 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 1004 END SELECT 1005 END IF 1006 IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN 1007 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 1008 CASE (XC_FAMILY_LDA) 1009 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)(rhoa)", & 1010 allocate_deriv=.TRUE.) 1011 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhoa) 1012 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)(rhob)", & 1013 allocate_deriv=.TRUE.) 1014 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhoa_rhob) 1015 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)(rhob)(rhob)", & 1016 allocate_deriv=.TRUE.) 1017 CALL xc_derivative_get(deriv, deriv_data=e_rhoa_rhob_rhob) 1018 deriv => xc_dset_get_derivative(deriv_set, "(rhob)(rhob)(rhob)", & 1019 allocate_deriv=.TRUE.) 1020 CALL xc_derivative_get(deriv, deriv_data=e_rhob_rhob_rhob) 1021 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 1022 CPABORT("derivatives larger than 2 not implemented") 1023 CASE default 1024 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 1025 END SELECT 1026 END IF 1027 IF (grad_deriv >= 4 .OR. grad_deriv <= -4) THEN 1028 CPABORT("derivatives larger than 3 not implemented") 1029 END IF 1030 1031!$OMP PARALLEL DEFAULT(NONE), & 1032!$OMP SHARED(rhoa,rhob,norm_drho,norm_drhoa,norm_drhob),& 1033!$OMP SHARED(laplace_rhoa,laplace_rhob,tau_a,tau_b),& 1034!$OMP SHARED(e_0,e_rhoa,e_rhob,e_ndrho,e_ndrhoa,e_ndrhob),& 1035!$OMP SHARED(e_laplace_rhoa,e_laplace_rhob,e_tau_a,e_tau_b),& 1036!$OMP SHARED(e_rhoa_rhoa,e_rhoa_rhob,e_rhob_rhob),& 1037!$OMP SHARED(e_ndrho_rhoa,e_ndrho_rhob),& 1038!$OMP SHARED(e_ndrhoa_rhoa,e_ndrhoa_rhob,e_ndrhob_rhoa,e_ndrhob_rhob),& 1039!$OMP SHARED(e_ndrho_ndrho,e_ndrho_ndrhoa,e_ndrho_ndrhob),& 1040!$OMP SHARED(e_ndrhoa_ndrhoa,e_ndrhoa_ndrhob,e_ndrhob_ndrhob),& 1041!$OMP SHARED(e_rhoa_laplace_rhoa,e_rhoa_laplace_rhob,e_rhob_laplace_rhoa,e_rhob_laplace_rhob),& 1042!$OMP SHARED(e_rhoa_tau_a,e_rhoa_tau_b,e_rhob_tau_a,e_rhob_tau_b),& 1043!$OMP SHARED(e_ndrho_laplace_rhoa,e_ndrho_laplace_rhob),& 1044!$OMP SHARED(e_ndrhoa_laplace_rhoa,e_ndrhoa_laplace_rhob,e_ndrhob_laplace_rhoa,e_ndrhob_laplace_rhob),& 1045!$OMP SHARED(e_ndrho_tau_a,e_ndrho_tau_b),& 1046!$OMP SHARED(e_ndrhoa_tau_a,e_ndrhoa_tau_b,e_ndrhob_tau_a,e_ndrhob_tau_b),& 1047!$OMP SHARED(e_laplace_rhoa_laplace_rhoa,e_laplace_rhoa_laplace_rhob,e_laplace_rhob_laplace_rhob),& 1048!$OMP SHARED(e_laplace_rhoa_tau_a,e_laplace_rhoa_tau_b,e_laplace_rhob_tau_a,e_laplace_rhob_tau_b),& 1049!$OMP SHARED(e_tau_a_tau_a,e_tau_a_tau_b,e_tau_b_tau_b),& 1050!$OMP SHARED(e_rhoa_rhoa_rhoa,e_rhoa_rhoa_rhob,e_rhoa_rhob_rhob,e_rhob_rhob_rhob),& 1051!$OMP SHARED(grad_deriv,npoints),& 1052!$OMP SHARED(epsilon_rho,epsilon_tau),& 1053!$OMP SHARED(func_name,func_scale,params) 1054 1055 CALL libxc_lsd_calc(rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, & 1056 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, laplace_rhoa=laplace_rhoa, & 1057 laplace_rhob=laplace_rhob, tau_a=tau_a, tau_b=tau_b, & 1058 e_0=e_0, e_rhoa=e_rhoa, e_rhob=e_rhob, e_ndrho=e_ndrho, & 1059 e_ndrhoa=e_ndrhoa, e_ndrhob=e_ndrhob, e_laplace_rhoa=e_laplace_rhoa, & 1060 e_laplace_rhob=e_laplace_rhob, e_tau_a=e_tau_a, e_tau_b=e_tau_b, & 1061 e_rhoa_rhoa=e_rhoa_rhoa, e_rhoa_rhob=e_rhoa_rhob, e_rhob_rhob=e_rhob_rhob, & 1062 e_ndrho_rhoa=e_ndrho_rhoa, e_ndrho_rhob=e_ndrho_rhob, & 1063 e_ndrhoa_rhoa=e_ndrhoa_rhoa, e_ndrhoa_rhob=e_ndrhoa_rhob, & 1064 e_ndrhob_rhoa=e_ndrhob_rhoa, e_ndrhob_rhob=e_ndrhob_rhob, & 1065 e_ndrho_ndrho=e_ndrho_ndrho, e_ndrho_ndrhoa=e_ndrho_ndrhoa, & 1066 e_ndrho_ndrhob=e_ndrho_ndrhob, e_ndrhoa_ndrhoa=e_ndrhoa_ndrhoa, & 1067 e_ndrhoa_ndrhob=e_ndrhoa_ndrhob, e_ndrhob_ndrhob=e_ndrhob_ndrhob, & 1068 e_rhoa_laplace_rhoa=e_rhoa_laplace_rhoa, & 1069 e_rhoa_laplace_rhob=e_rhoa_laplace_rhob, & 1070 e_rhob_laplace_rhoa=e_rhob_laplace_rhoa, & 1071 e_rhob_laplace_rhob=e_rhob_laplace_rhob, & 1072 e_rhoa_tau_a=e_rhoa_tau_a, e_rhoa_tau_b=e_rhoa_tau_b, & 1073 e_rhob_tau_a=e_rhob_tau_a, e_rhob_tau_b=e_rhob_tau_b, & 1074 e_ndrho_laplace_rhoa=e_ndrho_laplace_rhoa, & 1075 e_ndrho_laplace_rhob=e_ndrho_laplace_rhob, & 1076 e_ndrhoa_laplace_rhoa=e_ndrhoa_laplace_rhoa, & 1077 e_ndrhoa_laplace_rhob=e_ndrhoa_laplace_rhob, & 1078 e_ndrhob_laplace_rhoa=e_ndrhob_laplace_rhoa, & 1079 e_ndrhob_laplace_rhob=e_ndrhob_laplace_rhob, & 1080 e_ndrho_tau_a=e_ndrho_tau_a, e_ndrho_tau_b=e_ndrho_tau_b, & 1081 e_ndrhoa_tau_a=e_ndrhoa_tau_a, e_ndrhoa_tau_b=e_ndrhoa_tau_b, & 1082 e_ndrhob_tau_a=e_ndrhob_tau_a, e_ndrhob_tau_b=e_ndrhob_tau_b, & 1083 e_laplace_rhoa_laplace_rhoa=e_laplace_rhoa_laplace_rhoa, & 1084 e_laplace_rhoa_laplace_rhob=e_laplace_rhoa_laplace_rhob, & 1085 e_laplace_rhob_laplace_rhob=e_laplace_rhob_laplace_rhob, & 1086 e_laplace_rhoa_tau_a=e_laplace_rhoa_tau_a, & 1087 e_laplace_rhoa_tau_b=e_laplace_rhoa_tau_b, & 1088 e_laplace_rhob_tau_a=e_laplace_rhob_tau_a, & 1089 e_laplace_rhob_tau_b=e_laplace_rhob_tau_b, & 1090 e_tau_a_tau_a=e_tau_a_tau_a, & 1091 e_tau_a_tau_b=e_tau_a_tau_b, & 1092 e_tau_b_tau_b=e_tau_b_tau_b, & 1093 e_rhoa_rhoa_rhoa=e_rhoa_rhoa_rhoa, & 1094 e_rhoa_rhoa_rhob=e_rhoa_rhoa_rhob, & 1095 e_rhoa_rhob_rhob=e_rhoa_rhob_rhob, & 1096 e_rhob_rhob_rhob=e_rhob_rhob_rhob, & 1097 grad_deriv=grad_deriv, npoints=npoints, & 1098 epsilon_rho=epsilon_rho, & 1099 epsilon_tau=epsilon_tau, func_name=func_name, & 1100 sc=func_scale, params=params) 1101 1102!$OMP END PARALLEL 1103 1104 NULLIFY (dummy) 1105 1106 CALL xc_f03_func_end(xc_func) 1107 1108 CALL timestop(handle) 1109#else 1110 MARK_USED(rho_set) 1111 MARK_USED(deriv_set) 1112 MARK_USED(grad_deriv) 1113 MARK_USED(libxc_params) 1114 CPABORT("In order to use libxc you need to download and install it") 1115#endif 1116 END SUBROUTINE libxc_lsd_eval 1117 1118! ************************************************************************************************** 1119!> \brief libxc exchange-correlation functionals 1120!> \param rho density 1121!> \param norm_drho norm of the gradient of the density 1122!> \param laplace_rho laplacian of the density 1123!> \param tau kinetic-energy density 1124!> \param e_0 energy density 1125!> \param e_rho derivative of the energy density with respect to rho 1126!> \param e_ndrho derivative of the energy density with respect to ndrho 1127!> \param e_laplace_rho derivative of the energy density with respect to laplace_rho 1128!> \param e_tau derivative of the energy density with respect to tau 1129!> \param e_rho_rho derivative of the energy density with respect to rho_rho 1130!> \param e_ndrho_rho derivative of the energy density with respect to ndrho_rho 1131!> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho 1132!> \param e_rho_laplace_rho derivative of the energy density with respect to rho_laplace_rho 1133!> \param e_rho_tau derivative of the energy density with respect to rho_tau 1134!> \param e_ndrho_laplace_rho derivative of the energy density with respect to ndrho_laplace_rho 1135!> \param e_ndrho_tau derivative of the energy density with respect to ndrho_tau 1136!> \param e_laplace_rho_laplace_rho derivative of the energy density with respect to laplace_rho_laplace_rho 1137!> \param e_laplace_rho_tau derivative of the energy density with respect to laplace_rho_tau 1138!> \param e_tau_tau derivative of the energy density with respect to tau_tau 1139!> \param e_rho_rho_rho derivative of the energy density with respect to rho_rho_rho 1140!> \param grad_deriv degree of the derivative that should be evaluated, 1141!> if positive all the derivatives up to the given degree are evaluated, 1142!> if negative only the given degree is calculated 1143!> \param npoints number of points on the grid 1144!> \param epsilon_rho ... 1145!> \param epsilon_tau ... 1146!> \param func_name name of the functional 1147!> \param sc scaling factor 1148!> \param params parameters of the functional 1149!> \author F. Tran 1150! ************************************************************************************************** 1151#if defined (__LIBXC) 1152 SUBROUTINE libxc_lda_calc(rho, norm_drho, laplace_rho, tau, & 1153 e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, e_rho_rho, e_ndrho_rho, & 1154 e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, & 1155 e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, & 1156 e_tau_tau, e_rho_rho_rho, & 1157 grad_deriv, npoints, epsilon_rho, & 1158 epsilon_tau, func_name, sc, params) 1159 1160 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho, laplace_rho, tau 1161 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho, e_laplace_rho, e_tau, & 1162 e_rho_rho, e_ndrho_rho, e_ndrho_ndrho, e_rho_laplace_rho, e_rho_tau, e_ndrho_laplace_rho, & 1163 e_ndrho_tau, e_laplace_rho_laplace_rho, e_laplace_rho_tau, e_tau_tau, e_rho_rho_rho 1164 INTEGER, INTENT(in) :: grad_deriv, npoints 1165 REAL(KIND=dp), INTENT(in) :: epsilon_rho, epsilon_tau 1166 CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name 1167 REAL(KIND=dp), INTENT(in) :: sc 1168 REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: params 1169 1170 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lda_calc', routineP = moduleN//':'//routineN 1171 1172 INTEGER :: func_id, ii 1173 LOGICAL :: no_exc 1174 REAL(KIND=dp), DIMENSION(1) :: exc, my_tau, sigma, v2lapl2, v2lapltau, v2rho2, v2rholapl, & 1175 v2rhosigma, v2rhotau, v2sigma2, v2sigmalapl, v2sigmatau, v2tau2, v3rho3, vlapl, vrho, & 1176 vsigma, vtau 1177 TYPE(xc_f03_func_info_t) :: xc_info 1178 TYPE(xc_f03_func_t) :: xc_func 1179 1180 ! init vlapl (prevent libxc-4.0.x bug) 1181 vlapl = 0.0_dp 1182 1183 func_id = xc_libxc_wrap_functional_get_number(func_name) 1184!$OMP CRITICAL(libxc_init) 1185 CALL xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED) 1186 xc_info = xc_f03_func_get_info(xc_func) 1187 CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, params, no_exc) 1188!$OMP END CRITICAL(libxc_init) 1189!$OMP BARRIER 1190 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 1191 CASE (XC_FAMILY_LDA) 1192 IF (grad_deriv == 0) THEN 1193!$OMP DO 1194 DO ii = 1, npoints 1195 IF (rho(ii) > epsilon_rho) THEN 1196 CALL xc_f03_lda_exc(xc_func, 1, rho(ii), exc) 1197 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1198 END IF 1199 END DO 1200!$OMP END DO 1201 ELSE IF (grad_deriv == -1) THEN 1202!$OMP DO 1203 DO ii = 1, npoints 1204 IF (rho(ii) > epsilon_rho) THEN 1205 CALL xc_f03_lda_vxc(xc_func, 1, rho(ii), vrho) 1206 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1207 END IF 1208 END DO 1209!$OMP END DO 1210 ELSE IF (grad_deriv == 1) THEN 1211!$OMP DO 1212 DO ii = 1, npoints 1213 IF (rho(ii) > epsilon_rho) THEN 1214 CALL xc_f03_lda_exc_vxc(xc_func, 1, rho(ii), exc, vrho) 1215 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1216 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1217 END IF 1218 END DO 1219!$OMP END DO 1220 ELSE IF (grad_deriv == -2) THEN 1221!$OMP DO 1222 DO ii = 1, npoints 1223 IF (rho(ii) > epsilon_rho) THEN 1224 CALL xc_f03_lda_fxc(xc_func, 1, rho(ii), v2rho2) 1225 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1226 END IF 1227 END DO 1228!$OMP END DO 1229 ELSE IF (grad_deriv == 2) THEN 1230!$OMP DO 1231 DO ii = 1, npoints 1232 IF (rho(ii) > epsilon_rho) THEN 1233 CALL xc_f03_lda_exc_vxc(xc_func, 1, rho(ii), exc, vrho) 1234 CALL xc_f03_lda_fxc(xc_func, 1, rho(ii), v2rho2) 1235 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1236 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1237 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1238 END IF 1239 END DO 1240!$OMP END DO 1241 ELSE IF (grad_deriv == -3) THEN 1242!$OMP DO 1243 DO ii = 1, npoints 1244 IF (rho(ii) > epsilon_rho) THEN 1245 CALL xc_f03_lda_kxc(xc_func, 1, rho(ii), v3rho3) 1246 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1) 1247 END IF 1248 END DO 1249!$OMP END DO 1250 ELSE IF (grad_deriv == 3) THEN 1251!$OMP DO 1252 DO ii = 1, npoints 1253 IF (rho(ii) > epsilon_rho) THEN 1254 CALL xc_f03_lda(xc_func, 1, rho(ii), exc, vrho, v2rho2, v3rho3) 1255 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1256 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1257 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1258 e_rho_rho_rho(ii) = e_rho_rho_rho(ii) + sc*v3rho3(1) 1259 END IF 1260 END DO 1261!$OMP END DO 1262 END IF 1263 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 1264 IF (grad_deriv == 0) THEN 1265!$OMP DO 1266 DO ii = 1, npoints 1267 IF (rho(ii) > epsilon_rho) THEN 1268 sigma = norm_drho(ii)**2 1269 CALL xc_f03_gga_exc(xc_func, 1, rho(ii), sigma, exc) 1270 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1271 END IF 1272 END DO 1273!$OMP END DO 1274 ELSE IF (grad_deriv == -1) THEN 1275!$OMP DO 1276 DO ii = 1, npoints 1277 IF (rho(ii) > epsilon_rho) THEN 1278 sigma = norm_drho(ii)**2 1279 CALL xc_f03_gga_vxc(xc_func, 1, rho(ii), sigma, vrho, vsigma) 1280 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1281 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1282 END IF 1283 END DO 1284!$OMP END DO 1285 ELSE IF (grad_deriv == 1) THEN 1286!$OMP DO 1287 DO ii = 1, npoints 1288 IF (rho(ii) > epsilon_rho) THEN 1289 sigma = norm_drho(ii)**2 1290 IF (no_exc) THEN 1291 CALL xc_f03_gga_vxc(xc_func, 1, rho(ii), sigma, vrho, vsigma) 1292 exc = 0.0_dp 1293 ELSE 1294 CALL xc_f03_gga_exc_vxc(xc_func, 1, rho(ii), sigma, & 1295 exc, vrho, vsigma) 1296 END IF 1297 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1298 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1299 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1300 END IF 1301 END DO 1302!$OMP END DO 1303 ELSE IF (grad_deriv == -2) THEN 1304!$OMP DO 1305 DO ii = 1, npoints 1306 IF (rho(ii) > epsilon_rho) THEN 1307 sigma = norm_drho(ii)**2 1308 IF (no_exc) THEN 1309 CALL xc_f03_gga_vxc(xc_func, 1, rho(ii), sigma, vrho, vsigma) 1310 CALL xc_f03_gga_fxc(xc_func, 1, rho(ii), sigma, & 1311 v2rho2, v2rhosigma, v2sigma2) 1312 ELSE 1313 CALL xc_f03_gga_exc_vxc(xc_func, 1, rho(ii), sigma, & 1314 exc, vrho, vsigma) 1315 CALL xc_f03_gga_fxc(xc_func, 1, rho(ii), sigma, & 1316 v2rho2, v2rhosigma, v2sigma2) 1317 END IF 1318 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1319 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii) 1320 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1321 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1)) 1322 END IF 1323 END DO 1324!$OMP END DO 1325 ELSE IF (grad_deriv == 2) THEN 1326!$OMP DO 1327 DO ii = 1, npoints 1328 IF (rho(ii) > epsilon_rho) THEN 1329 sigma = norm_drho(ii)**2 1330 IF (no_exc) THEN 1331 CALL xc_f03_gga_vxc(xc_func, 1, rho(ii), sigma, vrho, vsigma) 1332 CALL xc_f03_gga_fxc(xc_func, 1, rho(ii), sigma, & 1333 v2rho2, v2rhosigma, v2sigma2) 1334 exc = 0.0_dp 1335 ELSE 1336 CALL xc_f03_gga_exc_vxc(xc_func, 1, rho(ii), sigma, & 1337 exc, vrho, vsigma) 1338 CALL xc_f03_gga_fxc(xc_func, 1, rho(ii), sigma, & 1339 v2rho2, v2rhosigma, v2sigma2) 1340 END IF 1341 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1342 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1343 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1344 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1345 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii) 1346 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1347 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1)) 1348 END IF 1349 END DO 1350!$OMP END DO 1351 END IF 1352 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 1353 IF (grad_deriv == 0) THEN 1354!$OMP DO 1355 DO ii = 1, npoints 1356 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN 1357 sigma = norm_drho(ii)**2 1358 my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) 1359 CALL xc_f03_mgga_exc(xc_func, 1, rho(ii), sigma, & 1360 laplace_rho(ii), my_tau, exc) 1361 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1362 END IF 1363 END DO 1364!$OMP END DO 1365 ELSE IF (grad_deriv == -1) THEN 1366!$OMP DO 1367 DO ii = 1, npoints 1368 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN 1369 sigma = norm_drho(ii)**2 1370 my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) 1371 CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & 1372 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) 1373 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1374 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1375 e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1) 1376 e_tau(ii) = e_tau(ii) + sc*vtau(1) 1377 END IF 1378 END DO 1379!$OMP END DO 1380 ELSE IF (grad_deriv == 1) THEN 1381!$OMP DO 1382 DO ii = 1, npoints 1383 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN 1384 sigma(1) = norm_drho(ii)**2 1385 my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) 1386 IF (no_exc) THEN 1387 CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & 1388 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) 1389 exc = 0.0_dp 1390 ELSE 1391 CALL xc_f03_mgga_exc_vxc(xc_func, 1, rho(ii), sigma, & 1392 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau) 1393 END IF 1394 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1395 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1396 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1397 e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1) 1398 e_tau(ii) = e_tau(ii) + sc*vtau(1) 1399 END IF 1400 END DO 1401!$OMP END DO 1402 ELSE IF (grad_deriv == -2) THEN 1403!$OMP DO 1404 DO ii = 1, npoints 1405 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN 1406 sigma = norm_drho(ii)**2 1407 my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) 1408 IF (no_exc) THEN 1409 CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & 1410 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) 1411 CALL xc_f03_mgga_fxc(xc_func, 1, rho(ii), sigma, & 1412 laplace_rho(ii), my_tau, & 1413 v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, & 1414 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau) 1415 ELSE 1416 CALL xc_f03_mgga(xc_func, 1, rho(ii), sigma, & 1417 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, & 1418 v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, & 1419 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau) 1420 END IF 1421 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1422 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii) 1423 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1424 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1)) 1425 e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1) 1426 e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1) 1427 e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + & 1428 sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii) 1429 e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii) 1430 e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1) 1431 e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1) 1432 e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1) 1433 END IF 1434 END DO 1435!$OMP END DO 1436 ELSE IF (grad_deriv == 2) THEN 1437!$OMP DO 1438 DO ii = 1, npoints 1439 IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN 1440 sigma = norm_drho(ii)**2 1441 my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) 1442 IF (no_exc) THEN 1443 CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & 1444 laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) 1445 CALL xc_f03_mgga_fxc(xc_func, 1, rho(ii), sigma, & 1446 laplace_rho(ii), my_tau, & 1447 v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, & 1448 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau) 1449 exc = 0.0_dp 1450 ELSE 1451 CALL xc_f03_mgga(xc_func, 1, rho(ii), sigma, & 1452 laplace_rho(ii), my_tau, exc, vrho, vsigma, vlapl, vtau, & 1453 v2rho2, v2sigma2, v2lapl2, v2tau2, v2rhosigma, v2rholapl, & 1454 v2rhotau, v2sigmalapl, v2sigmatau, v2lapltau) 1455 END IF 1456 e_0(ii) = e_0(ii) + sc*exc(1)*rho(ii) 1457 e_rho(ii) = e_rho(ii) + sc*vrho(1) 1458 e_ndrho(ii) = e_ndrho(ii) + sc*2.0_dp*vsigma(1)*norm_drho(ii) 1459 e_laplace_rho(ii) = e_laplace_rho(ii) + sc*vlapl(1) 1460 e_tau(ii) = e_tau(ii) + sc*vtau(1) 1461 e_rho_rho(ii) = e_rho_rho(ii) + sc*v2rho2(1) 1462 e_ndrho_rho(ii) = e_ndrho_rho(ii) + sc*2.0_dp*v2rhosigma(1)*norm_drho(ii) 1463 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1464 sc*2.0_dp*(2.0_dp*sigma(1)*v2sigma2(1) + vsigma(1)) 1465 e_rho_laplace_rho(ii) = e_rho_laplace_rho(ii) + sc*v2rholapl(1) 1466 e_rho_tau(ii) = e_rho_tau(ii) + sc*v2rhotau(1) 1467 e_ndrho_laplace_rho(ii) = e_ndrho_laplace_rho(ii) + & 1468 sc*2.0_dp*v2sigmalapl(1)*norm_drho(ii) 1469 e_ndrho_tau(ii) = e_ndrho_tau(ii) + sc*2.0_dp*v2sigmatau(1)*norm_drho(ii) 1470 e_laplace_rho_laplace_rho(ii) = e_laplace_rho_laplace_rho(ii) + sc*v2lapl2(1) 1471 e_laplace_rho_tau(ii) = e_laplace_rho_tau(ii) + sc*v2lapltau(1) 1472 e_tau_tau(ii) = e_tau_tau(ii) + sc*v2tau2(1) 1473 END IF 1474 END DO 1475!$OMP END DO 1476 END IF 1477 CASE default 1478 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 1479 END SELECT 1480 1481 CALL xc_f03_func_end(xc_func) 1482 1483 END SUBROUTINE libxc_lda_calc 1484#endif 1485 1486! ************************************************************************************************** 1487!> \brief libxc exchange-correlation functionals 1488!> \param rhoa alpha density 1489!> \param rhob beta density 1490!> \param norm_drho ... 1491!> \param norm_drhoa norm of the gradient of the alpha density 1492!> \param norm_drhob norm of the gradient of the beta density 1493!> \param laplace_rhoa laplacian of the alpha density 1494!> \param laplace_rhob laplacian of the beta density 1495!> \param tau_a alpha kinetic-energy density 1496!> \param tau_b beta kinetic-energy density 1497!> \param e_0 energy density 1498!> \param e_rhoa derivative of the energy density with respect to rhoa 1499!> \param e_rhob derivative of the energy density with respect to rhob 1500!> \param e_ndrho derivative of the energy density with respect to ndrho 1501!> \param e_ndrhoa derivative of the energy density with respect to ndrhoa 1502!> \param e_ndrhob derivative of the energy density with respect to ndrhob 1503!> \param e_laplace_rhoa derivative of the energy density with respect to laplace_rhoa 1504!> \param e_laplace_rhob derivative of the energy density with respect to laplace_rhob 1505!> \param e_tau_a derivative of the energy density with respect to tau_a 1506!> \param e_tau_b derivative of the energy density with respect to tau_b 1507!> \param e_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa 1508!> \param e_rhoa_rhob derivative of the energy density with respect to rhoa_rhob 1509!> \param e_rhob_rhob derivative of the energy density with respect to rhob_rhob 1510!> \param e_ndrho_rhoa derivative of the energy density with respect to ndrho_rhoa 1511!> \param e_ndrho_rhob derivative of the energy density with respect to ndrho_rhob 1512!> \param e_ndrhoa_rhoa derivative of the energy density with respect to ndrhoa_rhoa 1513!> \param e_ndrhoa_rhob derivative of the energy density with respect to ndrhoa_rhob 1514!> \param e_ndrhob_rhoa derivative of the energy density with respect to ndrhob_rhoa 1515!> \param e_ndrhob_rhob derivative of the energy density with respect to ndrhob_rhob 1516!> \param e_ndrho_ndrho derivative of the energy density with respect to ndrho_ndrho 1517!> \param e_ndrho_ndrhoa derivative of the energy density with respect to ndrho_ndrhoa 1518!> \param e_ndrho_ndrhob derivative of the energy density with respect to ndrho_ndrhob 1519!> \param e_ndrhoa_ndrhoa derivative of the energy density with respect to ndrhoa_ndrhoa 1520!> \param e_ndrhoa_ndrhob derivative of the energy density with respect to ndrhoa_ndrhob 1521!> \param e_ndrhob_ndrhob derivative of the energy density with respect to ndrhob_ndrhob 1522!> \param e_rhoa_laplace_rhoa derivative of the energy density with respect to rhoa_laplace_rhoa 1523!> \param e_rhoa_laplace_rhob derivative of the energy density with respect to rhoa_laplace_rhob 1524!> \param e_rhob_laplace_rhoa derivative of the energy density with respect to rhob_laplace_rhoa 1525!> \param e_rhob_laplace_rhob derivative of the energy density with respect to rhob_laplace_rhob 1526!> \param e_rhoa_tau_a derivative of the energy density with respect to rhoa_tau_a 1527!> \param e_rhoa_tau_b derivative of the energy density with respect to rhoa_tau_b 1528!> \param e_rhob_tau_a derivative of the energy density with respect to rhob_tau_a 1529!> \param e_rhob_tau_b derivative of the energy density with respect to rhob_tau_b 1530!> \param e_ndrho_laplace_rhoa derivative of the energy density with respect to ndrho_laplace_rhoa 1531!> \param e_ndrho_laplace_rhob derivative of the energy density with respect to ndrho_laplace_rhob 1532!> \param e_ndrhoa_laplace_rhoa derivative of the energy density with respect to ndrhoa_laplace_rhoa 1533!> \param e_ndrhoa_laplace_rhob derivative of the energy density with respect to ndrhoa_laplace_rhob 1534!> \param e_ndrhob_laplace_rhoa derivative of the energy density with respect to ndrhob_laplace_rhoa 1535!> \param e_ndrhob_laplace_rhob derivative of the energy density with respect to ndrhob_laplace_rhob 1536!> \param e_ndrho_tau_a derivative of the energy density with respect to ndrho_tau_a 1537!> \param e_ndrho_tau_b derivative of the energy density with respect to ndrho_tau_b 1538!> \param e_ndrhoa_tau_a derivative of the energy density with respect to ndrhoa_tau_a 1539!> \param e_ndrhoa_tau_b derivative of the energy density with respect to ndrhoa_tau_b 1540!> \param e_ndrhob_tau_a derivative of the energy density with respect to ndrhob_tau_a 1541!> \param e_ndrhob_tau_b derivative of the energy density with respect to ndrhob_tau_b 1542!> \param e_laplace_rhoa_laplace_rhoa derivative of the energy density with respect to laplace_rhoa_laplace_rhoa 1543!> \param e_laplace_rhoa_laplace_rhob derivative of the energy density with respect to laplace_rhoa_laplace_rhob 1544!> \param e_laplace_rhob_laplace_rhob derivative of the energy density with respect to laplace_rhob_laplace_rhob 1545!> \param e_laplace_rhoa_tau_a derivative of the energy density with respect to laplace_rhoa_tau_a 1546!> \param e_laplace_rhoa_tau_b derivative of the energy density with respect to laplace_rhoa_tau_b 1547!> \param e_laplace_rhob_tau_a derivative of the energy density with respect to laplace_rhob_tau_a 1548!> \param e_laplace_rhob_tau_b derivative of the energy density with respect to laplace_rhob_tau_b 1549!> \param e_tau_a_tau_a derivative of the energy density with respect to tau_a_tau_a 1550!> \param e_tau_a_tau_b derivative of the energy density with respect to tau_a_tau_b 1551!> \param e_tau_b_tau_b derivative of the energy density with respect to tau_b_tau_b 1552!> \param e_rhoa_rhoa_rhoa derivative of the energy density with respect to rhoa_rhoa_rhoa 1553!> \param e_rhoa_rhoa_rhob derivative of the energy density with respect to rhoa_rhoa_rhob 1554!> \param e_rhoa_rhob_rhob derivative of the energy density with respect to rhoa_rhob_rhob 1555!> \param e_rhob_rhob_rhob derivative of the energy density with respect to rhob_rhob_rhob 1556!> \param grad_deriv degree of the derivative that should be evaluated, 1557!> if positive all the derivatives up to the given degree are evaluated, 1558!> if negative only the given degree is calculated 1559!> \param npoints number of points on the grid 1560!> \param epsilon_rho ... 1561!> \param epsilon_tau ... 1562!> \param func_name name of the functional 1563!> \param sc scaling factor 1564!> \param params parameters of the functional 1565!> \author F. Tran 1566! ************************************************************************************************** 1567#if defined (__LIBXC) 1568 SUBROUTINE libxc_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, & 1569 norm_drhob, laplace_rhoa, laplace_rhob, tau_a, tau_b, & 1570 e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, e_ndrhob, & 1571 e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, & 1572 e_rhoa_rhoa, e_rhoa_rhob, e_rhob_rhob, & 1573 e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, & 1574 e_ndrhoa_rhob, e_ndrhob_rhoa, e_ndrhob_rhob, & 1575 e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, & 1576 e_ndrhoa_ndrhoa, e_ndrhoa_ndrhob, e_ndrhob_ndrhob, & 1577 e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, & 1578 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, & 1579 e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, e_rhob_tau_b, & 1580 e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, & 1581 e_ndrhoa_laplace_rhoa, e_ndrhoa_laplace_rhob, & 1582 e_ndrhob_laplace_rhoa, e_ndrhob_laplace_rhob, & 1583 e_ndrho_tau_a, e_ndrho_tau_b, & 1584 e_ndrhoa_tau_a, e_ndrhoa_tau_b, & 1585 e_ndrhob_tau_a, e_ndrhob_tau_b, & 1586 e_laplace_rhoa_laplace_rhoa, & 1587 e_laplace_rhoa_laplace_rhob, & 1588 e_laplace_rhob_laplace_rhob, & 1589 e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, & 1590 e_laplace_rhob_tau_a, e_laplace_rhob_tau_b, & 1591 e_tau_a_tau_a, e_tau_a_tau_b, e_tau_b_tau_b, & 1592 e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, & 1593 e_rhoa_rhob_rhob, e_rhob_rhob_rhob, & 1594 grad_deriv, npoints, epsilon_rho, & 1595 epsilon_tau, func_name, sc, params) 1596 1597 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, norm_drho, norm_drhoa, & 1598 norm_drhob, laplace_rhoa, & 1599 laplace_rhob, tau_a, tau_b 1600 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rhoa, e_rhob, e_ndrho, e_ndrhoa, & 1601 e_ndrhob, e_laplace_rhoa, e_laplace_rhob, e_tau_a, e_tau_b, e_rhoa_rhoa, e_rhoa_rhob, & 1602 e_rhob_rhob, e_ndrho_rhoa, e_ndrho_rhob, e_ndrhoa_rhoa, e_ndrhoa_rhob, e_ndrhob_rhoa, & 1603 e_ndrhob_rhob, e_ndrho_ndrho, e_ndrho_ndrhoa, e_ndrho_ndrhob, e_ndrhoa_ndrhoa, & 1604 e_ndrhoa_ndrhob, e_ndrhob_ndrhob, e_rhoa_laplace_rhoa, e_rhoa_laplace_rhob, & 1605 e_rhob_laplace_rhoa, e_rhob_laplace_rhob, e_rhoa_tau_a, e_rhoa_tau_b, e_rhob_tau_a, & 1606 e_rhob_tau_b, e_ndrho_laplace_rhoa, e_ndrho_laplace_rhob, e_ndrhoa_laplace_rhoa 1607 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_ndrhoa_laplace_rhob, e_ndrhob_laplace_rhoa, & 1608 e_ndrhob_laplace_rhob, e_ndrho_tau_a, e_ndrho_tau_b, e_ndrhoa_tau_a, e_ndrhoa_tau_b, & 1609 e_ndrhob_tau_a, e_ndrhob_tau_b, e_laplace_rhoa_laplace_rhoa, e_laplace_rhoa_laplace_rhob, & 1610 e_laplace_rhob_laplace_rhob, e_laplace_rhoa_tau_a, e_laplace_rhoa_tau_b, & 1611 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, & 1612 e_rhoa_rhoa_rhoa, e_rhoa_rhoa_rhob, e_rhoa_rhob_rhob, e_rhob_rhob_rhob 1613 INTEGER, INTENT(in) :: grad_deriv, npoints 1614 REAL(KIND=dp), INTENT(in) :: epsilon_rho, epsilon_tau 1615 CHARACTER(LEN=default_string_length), INTENT(IN) :: func_name 1616 REAL(KIND=dp), INTENT(in) :: sc 1617 REAL(KIND=dp), DIMENSION(:), INTENT(IN), POINTER :: params 1618 1619 CHARACTER(len=*), PARAMETER :: routineN = 'libxc_lsd_calc', routineP = moduleN//':'//routineN 1620 1621 INTEGER :: func_id, ii 1622 LOGICAL :: no_exc 1623 REAL(KIND=dp) :: my_norm_drho, my_norm_drhoa, & 1624 my_norm_drhob, my_rhoa, my_rhob, & 1625 my_tau_a, my_tau_b 1626 REAL(KIND=dp), DIMENSION(1) :: exc 1627 REAL(KIND=dp), DIMENSION(2, 1) :: laplace_rhov, rhov, tauv, vlapl, vrho, & 1628 vtau 1629 REAL(KIND=dp), DIMENSION(3, 1) :: sigmav, v2lapl2, v2rho2, v2tau2, vsigma 1630 REAL(KIND=dp), DIMENSION(4, 1) :: v2lapltau, v2rholapl, v2rhotau, v3rho3 1631 REAL(KIND=dp), DIMENSION(6, 1) :: v2rhosigma, v2sigma2, v2sigmalapl, & 1632 v2sigmatau 1633 TYPE(xc_f03_func_info_t) :: xc_info 1634 TYPE(xc_f03_func_t) :: xc_func 1635 1636 vlapl(1, 1) = 0.0_dp 1637 vlapl(2, 1) = 0.0_dp 1638 1639 func_id = xc_libxc_wrap_functional_get_number(func_name) 1640!$OMP CRITICAL(libxc_init) 1641 CALL xc_f03_func_init(xc_func, func_id, XC_POLARIZED) 1642 xc_info = xc_f03_func_get_info(xc_func) 1643 CALL xc_libxc_wrap_functional_set_params(xc_func, xc_info, params, no_exc) 1644!$OMP END CRITICAL(libxc_init) 1645!$OMP BARRIER 1646 1647 SELECT CASE (xc_f03_func_info_get_family(xc_info)) 1648 CASE (XC_FAMILY_LDA) 1649 IF (grad_deriv == 0) THEN 1650!$OMP DO 1651 DO ii = 1, npoints 1652 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1653 my_rhob = MAX(rhob(ii), 0.0_dp) 1654 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1655 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1656 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1657 CALL xc_f03_lda_exc(xc_func, 1, rhov(1, 1), exc) 1658 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1659 END IF 1660 END DO 1661!$OMP END DO 1662 ELSE IF (grad_deriv == -1) THEN 1663!$OMP DO 1664 DO ii = 1, npoints 1665 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1666 my_rhob = MAX(rhob(ii), 0.0_dp) 1667 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1668 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1669 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1670 CALL xc_f03_lda_vxc(xc_func, 1, rhov(1, 1), vrho(1, 1)) 1671 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1672 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1673 END IF 1674 END DO 1675!$OMP END DO 1676 ELSE IF (grad_deriv == 1) THEN 1677!$OMP DO 1678 DO ii = 1, npoints 1679 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1680 my_rhob = MAX(rhob(ii), 0.0_dp) 1681 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1682 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1683 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1684 CALL xc_f03_lda_exc_vxc(xc_func, 1, rhov(1, 1), exc, vrho(1, 1)) 1685 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1686 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1687 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1688 END IF 1689 END DO 1690!$OMP END DO 1691 ELSE IF (grad_deriv == -2) THEN 1692!$OMP DO 1693 DO ii = 1, npoints 1694 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1695 my_rhob = MAX(rhob(ii), 0.0_dp) 1696 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1697 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1698 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1699 CALL xc_f03_lda_fxc(xc_func, 1, rhov(1, 1), v2rho2(1, 1)) 1700 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 1701 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 1702 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 1703 END IF 1704 END DO 1705!$OMP END DO 1706 ELSE IF (grad_deriv == 2) THEN 1707!$OMP DO 1708 DO ii = 1, npoints 1709 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1710 my_rhob = MAX(rhob(ii), 0.0_dp) 1711 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1712 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1713 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1714 CALL xc_f03_lda_exc_vxc(xc_func, 1, rhov(1, 1), exc, vrho(1, 1)) 1715 CALL xc_f03_lda_fxc(xc_func, 1, rhov(1, 1), v2rho2(1, 1)) 1716 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1717 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1718 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1719 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 1720 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 1721 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 1722 END IF 1723 END DO 1724!$OMP END DO 1725 ELSE IF (grad_deriv == -3) THEN 1726!$OMP DO 1727 DO ii = 1, npoints 1728 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1729 my_rhob = MAX(rhob(ii), 0.0_dp) 1730 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1731 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1732 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1733 CALL xc_f03_lda_kxc(xc_func, 1, rhov(1, 1), v3rho3(1, 1)) 1734 e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1) 1735 e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1) 1736 e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1) 1737 e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1) 1738 END IF 1739 END DO 1740!$OMP END DO 1741 ELSE IF (grad_deriv == 3) THEN 1742!$OMP DO 1743 DO ii = 1, npoints 1744 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1745 my_rhob = MAX(rhob(ii), 0.0_dp) 1746 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1747 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1748 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1749 CALL xc_f03_lda(xc_func, 1, rhov(1, 1), exc, vrho(1, 1), v2rho2(1, 1), v3rho3(1, 1)) 1750 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1751 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1752 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1753 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 1754 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 1755 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 1756 e_rhoa_rhoa_rhoa(ii) = e_rhoa_rhoa_rhoa(ii) + sc*v3rho3(1, 1) 1757 e_rhoa_rhoa_rhob(ii) = e_rhoa_rhoa_rhob(ii) + sc*v3rho3(2, 1) 1758 e_rhoa_rhob_rhob(ii) = e_rhoa_rhob_rhob(ii) + sc*v3rho3(3, 1) 1759 e_rhob_rhob_rhob(ii) = e_rhob_rhob_rhob(ii) + sc*v3rho3(4, 1) 1760 END IF 1761 END DO 1762!$OMP END DO 1763 END IF 1764 CASE (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 1765 IF (grad_deriv == 0) THEN 1766!$OMP DO 1767 DO ii = 1, npoints 1768 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1769 my_rhob = MAX(rhob(ii), 0.0_dp) 1770 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1771 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1772 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1773 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1774 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1775 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1776 sigmav(1, 1) = my_norm_drhoa**2 1777 sigmav(3, 1) = my_norm_drhob**2 1778 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1779 CALL xc_f03_gga_exc(xc_func, 1, rhov(1, 1), sigmav(1, 1), exc) 1780 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1781 END IF 1782 END DO 1783!$OMP END DO 1784 ELSE IF (grad_deriv == -1) THEN 1785!$OMP DO 1786 DO ii = 1, npoints 1787 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1788 my_rhob = MAX(rhob(ii), 0.0_dp) 1789 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1790 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1791 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1792 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1793 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1794 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1795 sigmav(1, 1) = my_norm_drhoa**2 1796 sigmav(3, 1) = my_norm_drhob**2 1797 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1798 CALL xc_f03_gga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1)) 1799 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1800 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1801 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 1802 e_ndrhoa(ii) = e_ndrhoa(ii) + & 1803 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 1804 e_ndrhob(ii) = e_ndrhob(ii) + & 1805 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 1806 END IF 1807 END DO 1808!$OMP END DO 1809 ELSE IF (grad_deriv == 1) THEN 1810!$OMP DO 1811 DO ii = 1, npoints 1812 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1813 my_rhob = MAX(rhob(ii), 0.0_dp) 1814 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1815 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1816 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1817 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1818 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1819 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1820 sigmav(1, 1) = my_norm_drhoa**2 1821 sigmav(3, 1) = my_norm_drhob**2 1822 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1823 IF (no_exc) THEN 1824 CALL xc_f03_gga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1)) 1825 exc = 0.0_dp 1826 ELSE 1827 CALL xc_f03_gga_exc_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1)) 1828 END IF 1829 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1830 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1831 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1832 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 1833 e_ndrhoa(ii) = e_ndrhoa(ii) + & 1834 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 1835 e_ndrhob(ii) = e_ndrhob(ii) + & 1836 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 1837 END IF 1838 END DO 1839!$OMP END DO 1840 ELSE IF (grad_deriv == -2) THEN 1841!$OMP DO 1842 DO ii = 1, npoints 1843 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1844 my_rhob = MAX(rhob(ii), 0.0_dp) 1845 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1846 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1847 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1848 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1849 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1850 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1851 sigmav(1, 1) = my_norm_drhoa**2 1852 sigmav(3, 1) = my_norm_drhob**2 1853 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1854 IF (no_exc) THEN 1855 CALL xc_f03_gga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1)) 1856 CALL xc_f03_gga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 1857 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1)) 1858 ELSE 1859 CALL xc_f03_gga_exc_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1)) 1860 CALL xc_f03_gga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 1861 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1)) 1862 END IF 1863 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 1864 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 1865 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 1866 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho 1867 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho 1868 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + & 1869 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa 1870 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + & 1871 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa 1872 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + & 1873 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob 1874 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + & 1875 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob 1876 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1877 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1)) 1878 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + & 1879 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa 1880 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + & 1881 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob 1882 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + & 1883 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( & 1884 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1))) 1885 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + & 1886 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - & 1887 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob 1888 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + & 1889 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( & 1890 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))) 1891 END IF 1892 END DO 1893!$OMP END DO 1894 ELSE IF (grad_deriv == 2) THEN 1895!$OMP DO 1896 DO ii = 1, npoints 1897 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1898 my_rhob = MAX(rhob(ii), 0.0_dp) 1899 IF ((my_rhoa + my_rhob) > epsilon_rho) THEN 1900 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1901 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1902 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1903 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1904 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1905 sigmav(1, 1) = my_norm_drhoa**2 1906 sigmav(3, 1) = my_norm_drhob**2 1907 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1908 IF (no_exc) THEN 1909 CALL xc_f03_gga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), vrho(1, 1), vsigma(1, 1)) 1910 CALL xc_f03_gga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 1911 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1)) 1912 exc = 0.0_dp 1913 ELSE 1914 CALL xc_f03_gga_exc_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), exc, vrho(1, 1), vsigma(1, 1)) 1915 CALL xc_f03_gga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 1916 v2rho2(1, 1), v2rhosigma(1, 1), v2sigma2(1, 1)) 1917 END IF 1918 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1919 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 1920 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 1921 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 1922 e_ndrhoa(ii) = e_ndrhoa(ii) + & 1923 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 1924 e_ndrhob(ii) = e_ndrhob(ii) + & 1925 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 1926 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 1927 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 1928 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 1929 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho 1930 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho 1931 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + & 1932 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa 1933 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + & 1934 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa 1935 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + & 1936 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob 1937 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + & 1938 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob 1939 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 1940 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1)) 1941 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + & 1942 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa 1943 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + & 1944 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob 1945 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + & 1946 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( & 1947 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1))) 1948 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + & 1949 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - & 1950 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob 1951 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + & 1952 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( & 1953 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))) 1954 END IF 1955 END DO 1956!$OMP END DO 1957 END IF 1958 CASE (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 1959 IF (grad_deriv == 0) THEN 1960!$OMP DO 1961 DO ii = 1, npoints 1962 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1963 my_rhob = MAX(rhob(ii), 0.0_dp) 1964 my_tau_a = MAX(tau_a(ii), 0.0_dp) 1965 my_tau_b = MAX(tau_b(ii), 0.0_dp) 1966 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN 1967 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1968 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1969 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1970 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1971 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 1972 sigmav(1, 1) = my_norm_drhoa**2 1973 sigmav(3, 1) = my_norm_drhob**2 1974 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 1975 laplace_rhov(1, 1) = laplace_rhoa(ii) 1976 laplace_rhov(2, 1) = laplace_rhob(ii) 1977 tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp) 1978 tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp) 1979 tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1))) 1980 tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1))) 1981 CALL xc_f03_mgga_exc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 1982 laplace_rhov(1, 1), tauv(1, 1), exc) 1983 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 1984 END IF 1985 END DO 1986!$OMP END DO 1987 ELSE IF (grad_deriv == -1) THEN 1988!$OMP DO 1989 DO ii = 1, npoints 1990 my_rhoa = MAX(rhoa(ii), 0.0_dp) 1991 my_rhob = MAX(rhob(ii), 0.0_dp) 1992 my_tau_a = MAX(tau_a(ii), 0.0_dp) 1993 my_tau_b = MAX(tau_b(ii), 0.0_dp) 1994 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN 1995 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 1996 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 1997 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 1998 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 1999 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 2000 sigmav(1, 1) = my_norm_drhoa**2 2001 sigmav(3, 1) = my_norm_drhob**2 2002 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 2003 laplace_rhov(1, 1) = laplace_rhoa(ii) 2004 laplace_rhov(2, 1) = laplace_rhob(ii) 2005 tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp) 2006 tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp) 2007 tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1))) 2008 tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1))) 2009 CALL xc_f03_mgga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2010 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1)) 2011 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 2012 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 2013 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 2014 e_ndrhoa(ii) = e_ndrhoa(ii) + & 2015 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 2016 e_ndrhob(ii) = e_ndrhob(ii) + & 2017 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 2018 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1) 2019 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1) 2020 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1) 2021 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1) 2022 END IF 2023 END DO 2024!$OMP END DO 2025 ELSE IF (grad_deriv == 1) THEN 2026!$OMP DO 2027 DO ii = 1, npoints 2028 my_rhoa = MAX(rhoa(ii), 0.0_dp) 2029 my_rhob = MAX(rhob(ii), 0.0_dp) 2030 my_tau_a = MAX(tau_a(ii), 0.0_dp) 2031 my_tau_b = MAX(tau_b(ii), 0.0_dp) 2032 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN 2033 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 2034 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 2035 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 2036 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 2037 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 2038 sigmav(1, 1) = my_norm_drhoa**2 2039 sigmav(3, 1) = my_norm_drhob**2 2040 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 2041 laplace_rhov(1, 1) = laplace_rhoa(ii) 2042 laplace_rhov(2, 1) = laplace_rhob(ii) 2043 tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp) 2044 tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp) 2045 tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1))) 2046 tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1))) 2047 IF (no_exc) THEN 2048 CALL xc_f03_mgga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2049 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), & 2050 vlapl(1, 1), vtau(1, 1)) 2051 exc = 0.0_dp 2052 ELSE 2053 CALL xc_f03_mgga_exc_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2054 laplace_rhov(1, 1), tauv(1, 1), exc, & 2055 vrho(1, 1), vsigma(1, 1), vlapl(1, 1), vtau(1, 1)) 2056 END IF 2057 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 2058 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 2059 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 2060 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 2061 e_ndrhoa(ii) = e_ndrhoa(ii) + & 2062 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 2063 e_ndrhob(ii) = e_ndrhob(ii) + & 2064 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 2065 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1) 2066 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1) 2067 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1) 2068 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1) 2069 END IF 2070 END DO 2071!$OMP END DO 2072 ELSE IF (grad_deriv == -2) THEN 2073!$OMP DO 2074 DO ii = 1, npoints 2075 my_rhoa = MAX(rhoa(ii), 0.0_dp) 2076 my_rhob = MAX(rhob(ii), 0.0_dp) 2077 my_tau_a = MAX(tau_a(ii), 0.0_dp) 2078 my_tau_b = MAX(tau_b(ii), 0.0_dp) 2079 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN 2080 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 2081 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 2082 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 2083 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 2084 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 2085 sigmav(1, 1) = my_norm_drhoa**2 2086 sigmav(3, 1) = my_norm_drhob**2 2087 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 2088 laplace_rhov(1, 1) = laplace_rhoa(ii) 2089 laplace_rhov(2, 1) = laplace_rhob(ii) 2090 tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp) 2091 tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp) 2092 tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1))) 2093 tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1))) 2094 IF (no_exc) THEN 2095 CALL xc_f03_mgga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2096 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), & 2097 vlapl(1, 1), vtau(1, 1)) 2098 CALL xc_f03_mgga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2099 laplace_rhov(1, 1), tauv(1, 1), & 2100 v2rho2(1, 1), v2sigma2(1, 1), v2lapl2(1, 1), v2tau2(1, 1), & 2101 v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), & 2102 v2sigmalapl(1, 1), v2sigmatau(1, 1), v2lapltau(1, 1)) 2103 ELSE 2104 CALL xc_f03_mgga(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2105 laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), & 2106 vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2sigma2(1, 1), & 2107 v2lapl2(1, 1), v2tau2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), & 2108 v2rhotau(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), v2lapltau(1, 1)) 2109 END IF 2110 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 2111 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 2112 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 2113 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho 2114 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho 2115 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + & 2116 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa 2117 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + & 2118 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa 2119 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + & 2120 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob 2121 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + & 2122 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob 2123 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 2124 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1)) 2125 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + & 2126 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa 2127 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + & 2128 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob 2129 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + & 2130 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( & 2131 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1))) 2132 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + & 2133 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - & 2134 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob 2135 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + & 2136 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( & 2137 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))) 2138 e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1) 2139 e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1) 2140 e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1) 2141 e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1) 2142 e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1) 2143 e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1) 2144 e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1) 2145 e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1) 2146 e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho 2147 e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho 2148 e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + & 2149 sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa 2150 e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + & 2151 sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa 2152 e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + & 2153 sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob 2154 e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + & 2155 sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob 2156 e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho 2157 e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho 2158 e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + & 2159 sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa 2160 e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + & 2161 sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa 2162 e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + & 2163 sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob 2164 e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + & 2165 sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob 2166 e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1) 2167 e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1) 2168 e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1) 2169 e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1) 2170 e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1) 2171 e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1) 2172 e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1) 2173 e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1) 2174 e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1) 2175 e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1) 2176 END IF 2177 END DO 2178!$OMP END DO 2179 ELSE IF (grad_deriv == 2) THEN 2180!$OMP DO 2181 DO ii = 1, npoints 2182 my_rhoa = MAX(rhoa(ii), 0.0_dp) 2183 my_rhob = MAX(rhob(ii), 0.0_dp) 2184 my_tau_a = MAX(tau_a(ii), 0.0_dp) 2185 my_tau_b = MAX(tau_b(ii), 0.0_dp) 2186 IF (((my_rhoa + my_rhob) > epsilon_rho) .AND. ((my_tau_a + my_tau_b) > epsilon_tau)) THEN 2187 rhov(1, 1) = MAX(my_rhoa, EPSILON(0.0_dp)*1.e4_dp) 2188 rhov(2, 1) = MAX(my_rhob, EPSILON(0.0_dp)*1.e4_dp) 2189 my_norm_drhoa = MAX(norm_drhoa(ii), EPSILON(0.0_dp)*1.e4_dp) 2190 my_norm_drhob = MAX(norm_drhob(ii), EPSILON(0.0_dp)*1.e4_dp) 2191 my_norm_drho = MAX(norm_drho(ii), EPSILON(0.0_dp)*1.e4_dp) 2192 sigmav(1, 1) = my_norm_drhoa**2 2193 sigmav(3, 1) = my_norm_drhob**2 2194 sigmav(2, 1) = 0.5_dp*(my_norm_drho**2 - sigmav(1, 1) - sigmav(3, 1)) 2195 laplace_rhov(1, 1) = laplace_rhoa(ii) 2196 laplace_rhov(2, 1) = laplace_rhob(ii) 2197 tauv(1, 1) = MAX(my_tau_a, EPSILON(0.0_dp)*1.e4_dp) 2198 tauv(2, 1) = MAX(my_tau_b, EPSILON(0.0_dp)*1.e4_dp) 2199 tauv(1, 1) = MAX(tauv(1, 1), sigmav(1, 1)/(8.0_dp*rhov(1, 1))) 2200 tauv(2, 1) = MAX(tauv(2, 1), sigmav(3, 1)/(8.0_dp*rhov(2, 1))) 2201 IF (no_exc) THEN 2202 CALL xc_f03_mgga_vxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2203 laplace_rhov(1, 1), tauv(1, 1), vrho(1, 1), vsigma(1, 1), & 2204 vlapl(1, 1), vtau(1, 1)) 2205 CALL xc_f03_mgga_fxc(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2206 laplace_rhov(1, 1), tauv(1, 1), & 2207 v2rho2(1, 1), v2sigma2(1, 1), v2lapl2(1, 1), v2tau2(1, 1), & 2208 v2rhosigma(1, 1), v2rholapl(1, 1), v2rhotau(1, 1), & 2209 v2sigmalapl(1, 1), v2sigmatau(1, 1), v2lapltau(1, 1)) 2210 exc = 0.0_dp 2211 ELSE 2212 CALL xc_f03_mgga(xc_func, 1, rhov(1, 1), sigmav(1, 1), & 2213 laplace_rhov(1, 1), tauv(1, 1), exc, vrho(1, 1), vsigma(1, 1), & 2214 vlapl(1, 1), vtau(1, 1), v2rho2(1, 1), v2sigma2(1, 1), & 2215 v2lapl2(1, 1), v2tau2(1, 1), v2rhosigma(1, 1), v2rholapl(1, 1), & 2216 v2rhotau(1, 1), v2sigmalapl(1, 1), v2sigmatau(1, 1), v2lapltau(1, 1)) 2217 END IF 2218 e_0(ii) = e_0(ii) + sc*exc(1)*(rhov(1, 1) + rhov(2, 1)) 2219 e_rhoa(ii) = e_rhoa(ii) + sc*vrho(1, 1) 2220 e_rhob(ii) = e_rhob(ii) + sc*vrho(2, 1) 2221 e_ndrho(ii) = e_ndrho(ii) + sc*vsigma(2, 1)*my_norm_drho 2222 e_ndrhoa(ii) = e_ndrhoa(ii) + & 2223 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1))*my_norm_drhoa 2224 e_ndrhob(ii) = e_ndrhob(ii) + & 2225 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1))*my_norm_drhob 2226 e_laplace_rhoa(ii) = e_laplace_rhoa(ii) + sc*vlapl(1, 1) 2227 e_laplace_rhob(ii) = e_laplace_rhob(ii) + sc*vlapl(2, 1) 2228 e_tau_a(ii) = e_tau_a(ii) + sc*vtau(1, 1) 2229 e_tau_b(ii) = e_tau_b(ii) + sc*vtau(2, 1) 2230 e_rhoa_rhoa(ii) = e_rhoa_rhoa(ii) + sc*v2rho2(1, 1) 2231 e_rhoa_rhob(ii) = e_rhoa_rhob(ii) + sc*v2rho2(2, 1) 2232 e_rhob_rhob(ii) = e_rhob_rhob(ii) + sc*v2rho2(3, 1) 2233 e_ndrho_rhoa(ii) = e_ndrho_rhoa(ii) + sc*v2rhosigma(2, 1)*my_norm_drho 2234 e_ndrho_rhob(ii) = e_ndrho_rhob(ii) + sc*v2rhosigma(5, 1)*my_norm_drho 2235 e_ndrhoa_rhoa(ii) = e_ndrhoa_rhoa(ii) + & 2236 sc*(2.0_dp*v2rhosigma(1, 1) - v2rhosigma(2, 1))*my_norm_drhoa 2237 e_ndrhoa_rhob(ii) = e_ndrhoa_rhob(ii) + & 2238 sc*(2.0_dp*v2rhosigma(4, 1) - v2rhosigma(5, 1))*my_norm_drhoa 2239 e_ndrhob_rhoa(ii) = e_ndrhob_rhoa(ii) + & 2240 sc*(2.0_dp*v2rhosigma(3, 1) - v2rhosigma(2, 1))*my_norm_drhob 2241 e_ndrhob_rhob(ii) = e_ndrhob_rhob(ii) + & 2242 sc*(2.0_dp*v2rhosigma(6, 1) - v2rhosigma(5, 1))*my_norm_drhob 2243 e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) + & 2244 sc*(vsigma(2, 1) + my_norm_drho**2*v2sigma2(4, 1)) 2245 e_ndrho_ndrhoa(ii) = e_ndrho_ndrhoa(ii) + & 2246 sc*(2.0_dp*v2sigma2(2, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhoa 2247 e_ndrho_ndrhob(ii) = e_ndrho_ndrhob(ii) + & 2248 sc*(2.0_dp*v2sigma2(5, 1) - v2sigma2(4, 1))*my_norm_drho*my_norm_drhob 2249 e_ndrhoa_ndrhoa(ii) = e_ndrhoa_ndrhoa(ii) + & 2250 sc*(2.0_dp*vsigma(1, 1) - vsigma(2, 1) + my_norm_drhoa**2*( & 2251 4.0_dp*v2sigma2(1, 1) - 4.0_dp*v2sigma2(2, 1) + v2sigma2(4, 1))) 2252 e_ndrhoa_ndrhob(ii) = e_ndrhoa_ndrhob(ii) + & 2253 sc*(4.0_dp*v2sigma2(3, 1) - 2.0_dp*v2sigma2(2, 1) - & 2254 2.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))*my_norm_drhoa*my_norm_drhob 2255 e_ndrhob_ndrhob(ii) = e_ndrhob_ndrhob(ii) + & 2256 sc*(2.0_dp*vsigma(3, 1) - vsigma(2, 1) + my_norm_drhob**2*( & 2257 4.0_dp*v2sigma2(6, 1) - 4.0_dp*v2sigma2(5, 1) + v2sigma2(4, 1))) 2258 e_rhoa_laplace_rhoa(ii) = e_rhoa_laplace_rhoa(ii) + sc*v2rholapl(1, 1) 2259 e_rhoa_laplace_rhob(ii) = e_rhoa_laplace_rhob(ii) + sc*v2rholapl(2, 1) 2260 e_rhob_laplace_rhoa(ii) = e_rhob_laplace_rhoa(ii) + sc*v2rholapl(3, 1) 2261 e_rhob_laplace_rhob(ii) = e_rhob_laplace_rhob(ii) + sc*v2rholapl(4, 1) 2262 e_rhoa_tau_a(ii) = e_rhoa_tau_a(ii) + sc*v2rhotau(1, 1) 2263 e_rhoa_tau_b(ii) = e_rhoa_tau_b(ii) + sc*v2rhotau(2, 1) 2264 e_rhob_tau_a(ii) = e_rhob_tau_a(ii) + sc*v2rhotau(3, 1) 2265 e_rhob_tau_b(ii) = e_rhob_tau_b(ii) + sc*v2rhotau(4, 1) 2266 e_ndrho_laplace_rhoa(ii) = e_ndrho_laplace_rhoa(ii) + sc*v2sigmalapl(3, 1)*my_norm_drho 2267 e_ndrho_laplace_rhob(ii) = e_ndrho_laplace_rhob(ii) + sc*v2sigmalapl(4, 1)*my_norm_drho 2268 e_ndrhoa_laplace_rhoa(ii) = e_ndrhoa_laplace_rhoa(ii) + & 2269 sc*(2.0_dp*v2sigmalapl(1, 1) - v2sigmalapl(3, 1))*my_norm_drhoa 2270 e_ndrhoa_laplace_rhob(ii) = e_ndrhoa_laplace_rhob(ii) + & 2271 sc*(2.0_dp*v2sigmalapl(2, 1) - v2sigmalapl(4, 1))*my_norm_drhoa 2272 e_ndrhob_laplace_rhoa(ii) = e_ndrhob_laplace_rhoa(ii) + & 2273 sc*(2.0_dp*v2sigmalapl(5, 1) - v2sigmalapl(3, 1))*my_norm_drhob 2274 e_ndrhob_laplace_rhob(ii) = e_ndrhob_laplace_rhob(ii) + & 2275 sc*(2.0_dp*v2sigmalapl(6, 1) - v2sigmalapl(4, 1))*my_norm_drhob 2276 e_ndrho_tau_a(ii) = e_ndrho_tau_a(ii) + sc*v2sigmatau(3, 1)*my_norm_drho 2277 e_ndrho_tau_b(ii) = e_ndrho_tau_b(ii) + sc*v2sigmatau(4, 1)*my_norm_drho 2278 e_ndrhoa_tau_a(ii) = e_ndrhoa_tau_a(ii) + & 2279 sc*(2.0_dp*v2sigmatau(1, 1) - v2sigmatau(3, 1))*my_norm_drhoa 2280 e_ndrhoa_tau_b(ii) = e_ndrhoa_tau_b(ii) + & 2281 sc*(2.0_dp*v2sigmatau(2, 1) - v2sigmatau(4, 1))*my_norm_drhoa 2282 e_ndrhob_tau_a(ii) = e_ndrhob_tau_a(ii) + & 2283 sc*(2.0_dp*v2sigmatau(5, 1) - v2sigmatau(3, 1))*my_norm_drhob 2284 e_ndrhob_tau_b(ii) = e_ndrhob_tau_b(ii) + & 2285 sc*(2.0_dp*v2sigmatau(6, 1) - v2sigmatau(4, 1))*my_norm_drhob 2286 e_laplace_rhoa_laplace_rhoa(ii) = e_laplace_rhoa_laplace_rhoa(ii) + sc*v2lapl2(1, 1) 2287 e_laplace_rhoa_laplace_rhob(ii) = e_laplace_rhoa_laplace_rhob(ii) + sc*v2lapl2(2, 1) 2288 e_laplace_rhob_laplace_rhob(ii) = e_laplace_rhob_laplace_rhob(ii) + sc*v2lapl2(3, 1) 2289 e_laplace_rhoa_tau_a(ii) = e_laplace_rhoa_tau_a(ii) + sc*v2lapltau(1, 1) 2290 e_laplace_rhoa_tau_b(ii) = e_laplace_rhoa_tau_b(ii) + sc*v2lapltau(2, 1) 2291 e_laplace_rhob_tau_a(ii) = e_laplace_rhob_tau_a(ii) + sc*v2lapltau(3, 1) 2292 e_laplace_rhob_tau_b(ii) = e_laplace_rhob_tau_b(ii) + sc*v2lapltau(4, 1) 2293 e_tau_a_tau_a(ii) = e_tau_a_tau_a(ii) + sc*v2tau2(1, 1) 2294 e_tau_a_tau_b(ii) = e_tau_a_tau_b(ii) + sc*v2tau2(2, 1) 2295 e_tau_b_tau_b(ii) = e_tau_b_tau_b(ii) + sc*v2tau2(3, 1) 2296 END IF 2297 END DO 2298!$OMP END DO 2299 END IF 2300 CASE default 2301 CPABORT(TRIM(func_name)//": this XC_FAMILY is currently not supported.") 2302 END SELECT 2303 2304 CALL xc_f03_func_end(xc_func) 2305 2306 END SUBROUTINE libxc_lsd_calc 2307#endif 2308 2309END MODULE xc_libxc 2310