1! Copyright (C) 2018 Andre Severo Pereira Gomes, Christoph Jacob, Lucas Visscher and collaborators 2! 3! This file is part of Embed, a program implementing the Frozen Density Embedding (FDE) framework 4! 5! This Source Code Form is subject to the terms of the Mozilla Public 6! License, v. 2.0. If a copy of the MPL was not distributed with this 7! file, You can obtain one at http://mozilla.org/MPL/2.0/. 8! 9 10module fde_nadd_derv 11 12 use fde_types 13 14 use fde_max_block_length 15 use fde_xcfun_interface 16 17 public fde_initialize_nadd_functionals 18 public fde_print_nadd_functionals 19 20 public get_fde_nadd_fun_derv 21 public fde_calc_xc_nadd_fun_derv 22 public fde_calc_ke_nadd_fun_derv 23 24 public fde_set_skip_nadd_xc 25 public fde_set_skip_nadd_ke 26 public fde_set_nadd_kef 27 public fde_set_nadd_xcf 28 public fde_set_nadd_all_alda 29 30 public fde_xc_fun_is_lda 31 public fde_xc_fun_is_gga 32 public fde_ke_fun_is_lda 33 public fde_ke_fun_is_gga 34 public is_fdersp_alda 35 36 logical,save :: skip_kin_contribution = .false. 37 logical,save :: skip_xc_contribution = .false. 38 logical,save :: is_fdersp_alda = .false. 39 40 type(functional), save :: fde_kef, fde_xcf 41 type(functional), save :: fde_kef_alda, fde_xcf_alda 42 type(functional), save :: fde_kef_xalda, fde_xcf_xalda 43 44! the default (orbital-free) kinetic energy functional for now is the thomas-fermi one 45 character(len=80), save :: kefun_fde='kin_tf' 46 character(len=80), save :: xcfun_fde='lda' 47! character(len=80), save :: kefun_fde='kin_pw91' 48! character(len=80), save :: xcfun_fde='pbe' 49 50 character(len=80), save :: kefun_fde_alda='kin_tf' 51 character(len=80), save :: xcfun_fde_alda='lda' 52 53 logical, save :: fde_set_functionals(20) 54 logical, save :: functionals_initialized = .false. 55 56 real(8) :: hfx_out, mu_out, beta_out 57 58! derv_lenght here is set to the appropriate parameter from xc_derv 59! looking at xcint, this would be done differently if using openrsp... 60 integer :: derv_length = nr_nonzero_derv 61 62 public fde_new_derv_result 63 public fde_delete_derv_result 64 65 interface fde_new_derv_result 66 module procedure fde_new_derv_result_partial 67 end interface 68 69 interface fde_delete_derv_result 70 module procedure fde_delete_derv_result_partial 71 end interface 72 73 contains 74 75 76! ---------------------------------------------------------------------------- 77 subroutine fde_new_derv_result_partial(npoints,derv) 78! ---------------------------------------------------------------------------- 79 real(kind=8), allocatable :: derv(:,:) 80 integer :: npoints 81 82 allocate(derv(npoints, 0:derv_length)) 83 derv = 0.0d0 84 end subroutine fde_new_derv_result_partial 85 86 87! ---------------------------------------------------------------------------- 88 subroutine fde_del_derv_result_partial(derv) 89! ---------------------------------------------------------------------------- 90 real(kind=8), allocatable :: derv(:,:) 91 92 deallocate(derv) 93 end subroutine fde_del_derv_result_partial 94 95 96! ---------------------------------------------------------------------------- 97 subroutine calc_xcke_energy_density(gf, derv) 98! ---------------------------------------------------------------------------- 99 type(grid_function), intent(in) :: gf 100 real(kind=8), intent(out) :: derv(gf%npoints,0:derv_length) 101 real(kind=8) :: dervtmp(max_block_length, 0:derv_length) 102 103 integer :: i 104 integer :: order = 0 105 integer :: bllen = 1 106 real(kind=8) :: r_0(1) 107 real(kind=8) :: z_0(1) 108 real(kind=8) :: w(1) 109 110 dervtmp = 0.0d0 111 w(1) = 1.0d0 112 113 do i = 1, gf%npoints 114 r_0(1) = gf%n(i) 115 z_0(1) = gf%gn(i,1)*gf%gn(i,1) & 116 + gf%gn(i,2)*gf%gn(i,2) & 117 + gf%gn(i,3)*gf%gn(i,3) 118 119 call get_xcke_fun_derv(order, & 120 bllen, & 121 w, & 122 r_0, & 123 z_0, & 124 dervtmp) 125 126 derv(i,:) = dervtmp(1,:) 127 enddo 128 end subroutine calc_xcke_energy_density 129 130 131! ---------------------------------------------------------------------------- 132 subroutine fde_delete_derv_result_partial(derv) 133! ---------------------------------------------------------------------------- 134 real(kind=8), allocatable :: derv(:,:) 135 136 deallocate(derv) 137 end subroutine fde_delete_derv_result_partial 138 139 140! ---------------------------------------------------------------------------- 141 function fde_xc_fun_is_lda() 142! ---------------------------------------------------------------------------- 143 logical :: fde_xc_fun_is_lda 144 145 fde_xc_fun_is_lda = (xc_get_type(fde_xcf%id) == 0) 146 end function 147 148 149! ---------------------------------------------------------------------------- 150 function fde_xc_fun_is_gga() 151! ---------------------------------------------------------------------------- 152 logical :: fde_xc_fun_is_gga 153 154 fde_xc_fun_is_gga = (xc_get_type(fde_xcf%id) == 1) 155 end function 156 157 158! ---------------------------------------------------------------------------- 159 function fde_ke_fun_is_lda() 160! ---------------------------------------------------------------------------- 161 logical :: fde_ke_fun_is_lda 162 163 fde_ke_fun_is_lda = (xc_get_type(fde_kef%id) == 0) 164 end function 165 166 167! ---------------------------------------------------------------------------- 168 function fde_ke_fun_is_gga() 169! ---------------------------------------------------------------------------- 170 logical :: fde_ke_fun_is_gga 171 172 fde_ke_fun_is_gga = (xc_get_type(fde_kef%id) == 1) 173 end function 174 175 176! ---------------------------------------------------------------------------- 177 subroutine fde_set_skip_nadd_xc(value) 178! ---------------------------------------------------------------------------- 179 logical :: value 180 skip_xc_contribution = value 181 end subroutine 182 183 184! ---------------------------------------------------------------------------- 185 subroutine fde_set_skip_nadd_ke(value) 186! ---------------------------------------------------------------------------- 187 logical :: value 188 skip_kin_contribution = value 189 end subroutine 190 191 192! ---------------------------------------------------------------------------- 193 subroutine fde_initialize_nadd_functionals 194! ---------------------------------------------------------------------------- 195 integer :: i, order = 5 196 logical :: parallel_fde 197 198 parallel_fde = .false. 199 200 if (.not.functionals_initialized) then 201 functionals_initialized = .true. 202 203 call parse_functional(kefun_fde, fde_kef, hfx_out, mu_out, beta_out, .false.) 204 call parse_functional(xcfun_fde, fde_xcf, hfx_out, mu_out, beta_out, .false.) 205 call parse_functional(kefun_fde_alda, fde_kef_alda, hfx_out, mu_out, beta_out, .false.) 206 call parse_functional(xcfun_fde_alda, fde_xcf_alda, hfx_out, mu_out, beta_out, .false.) 207 208 209! ! we set them now so that we can verify whether or not they 210 call xcfun_set_functional(fde_xcf , order, parallel_fde) 211 call xcfun_set_functional(fde_kef, order, parallel_fde) 212 call xcfun_set_functional(fde_xcf_alda, order, parallel_fde) 213 call xcfun_set_functional(fde_kef_alda, order, parallel_fde) 214 endif 215 216 end subroutine fde_initialize_nadd_functionals 217 218 219! ---------------------------------------------------------------------------- 220 subroutine fde_set_nadd_kef(id) 221! ---------------------------------------------------------------------------- 222 character(len=80) :: id 223 224 write(kefun_fde,'(A)') id 225 end subroutine fde_set_nadd_kef 226 227 228! ---------------------------------------------------------------------------- 229 subroutine fde_set_nadd_xcf(id) 230! ---------------------------------------------------------------------------- 231 character(len=80) :: id 232 233 write(xcfun_fde,'(A)') id 234 end subroutine fde_set_nadd_xcf 235 236 237! ---------------------------------------------------------------------------- 238 subroutine fde_set_nadd_all_alda 239! ---------------------------------------------------------------------------- 240 is_fdersp_alda = .true. 241 end subroutine fde_set_nadd_all_alda 242 243 244! ---------------------------------------------------------------------------- 245 subroutine fde_print_nadd_functionals(unit) 246! ---------------------------------------------------------------------------- 247 integer :: unit 248 249 write(unit,'(A,A20)') ' * Non-additive functional : (KE) ',kefun_fde 250 write(unit,'(A,A20)') ' (XC) ',xcfun_fde 251 write(unit,'(A)') ' ' 252 write(unit,'(A,A20)') ' if .RSPLDA set will use : (XC) ',xcfun_fde_alda 253 write(unit,'(A,A20)') ' (KE) ',kefun_fde_alda 254 write(unit,'(A)') ' ' 255 end subroutine fde_print_nadd_functionals 256 257 258! ---------------------------------------------------------------------------- 259 subroutine get_fde_nadd_fun_derv(order, & 260 block_length, & 261 w, & 262 r_0, & 263 r_0_t, & 264 z_0, & 265 z_0_t, & 266 derv) 267! ---------------------------------------------------------------------------- 268 integer, intent(in) :: order 269 integer, intent(in) :: block_length 270 real(8), intent(in) :: w(*) 271 real(8), intent(in) :: r_0(*) 272 real(8), intent(in) :: z_0(*) 273 real(8), intent(in) :: r_0_t(*) 274 real(8), intent(in) :: z_0_t(*) 275 real(8), intent(inout) :: derv(max_block_length, 0:derv_length) 276 277 derv = 0.0d0 278 279 if (.not.skip_xc_contribution) & 280 call get_nadd_fun_derv(fde_xcf, & 281 fde_xcf_alda, & 282 fde_xcf_alda, & 283 order, & 284 block_length, & 285 w, & 286 r_0, & 287 r_0_t, & 288 z_0, & 289 z_0_t, & 290 derv) 291 292 if (.not.skip_kin_contribution) & 293 call get_nadd_fun_derv(fde_kef, & 294 fde_kef_alda, & 295 fde_kef_alda, & 296 order, & 297 block_length, & 298 w, & 299 r_0, & 300 r_0_t, & 301 z_0, & 302 z_0_t, & 303 derv) 304 305 end subroutine get_fde_nadd_fun_derv 306 307 308! ---------------------------------------------------------------------------- 309 subroutine get_xc_nadd_fun_derv(order, & 310 block_length, & 311 w, & 312 r_0, & 313 r_0_t, & 314 z_0, & 315 z_0_t, & 316 derv) 317! ---------------------------------------------------------------------------- 318 integer, intent(in) :: order 319 integer, intent(in) :: block_length 320 real(8), intent(in) :: w(*) 321 real(8), intent(in) :: r_0(*) 322 real(8), intent(in) :: z_0(*) 323 real(8), intent(in) :: r_0_t(*) 324 real(8), intent(in) :: z_0_t(*) 325 real(8), intent(inout) :: derv(max_block_length, 0:derv_length) 326 327 call get_nadd_fun_derv(fde_xcf, & 328 fde_xcf_alda, & 329 fde_xcf_alda, & 330 order, & 331 block_length, & 332 w, & 333 r_0, & 334 r_0_t, & 335 z_0, & 336 z_0_t, & 337 derv) 338 339 end subroutine get_xc_nadd_fun_derv 340 341 342! ---------------------------------------------------------------------------- 343 subroutine get_ke_nadd_fun_derv(order, & 344 block_length, & 345 w, & 346 r_0, & 347 r_0_t, & 348 z_0, & 349 z_0_t, & 350 derv) 351! ---------------------------------------------------------------------------- 352 integer, intent(in) :: order 353 integer, intent(in) :: block_length 354 real(8), intent(in) :: w(*) 355 real(8), intent(in) :: r_0(*) 356 real(8), intent(in) :: z_0(*) 357 real(8), intent(in) :: r_0_t(*) 358 real(8), intent(in) :: z_0_t(*) 359 real(8), intent(inout) :: derv(max_block_length, 0:derv_length) 360 361 call get_nadd_fun_derv(fde_kef, & 362 fde_kef_alda, & 363 fde_kef_alda, & 364 order, & 365 block_length, & 366 w, & 367 r_0, & 368 r_0_t, & 369 z_0, & 370 z_0_t, & 371 derv) 372 373 end subroutine get_ke_nadd_fun_derv 374 375 376 377! ---------------------------------------------------------------------------- 378 subroutine get_nadd_fun_derv(f, & 379 f_alda, & 380 f_xalda, & 381 order, & 382 block_length, & 383 w, & 384 r_0, & 385 r_0_t, & 386 z_0, & 387 z_0_t, & 388 derv) 389! ---------------------------------------------------------------------------- 390 type(functional) :: f 391 type(functional) :: f_alda, f_xalda 392 integer, intent(in) :: order 393 integer, intent(in) :: block_length 394 real(8), intent(in) :: w(*) 395 real(8), intent(in) :: r_0(*) 396 real(8), intent(in) :: z_0(*) 397 real(8), intent(in) :: r_0_t(*) 398 real(8), intent(in) :: z_0_t(*) 399 real(8), intent(inout) :: derv(max_block_length, 0:derv_length) 400 real(8) :: dervtmp(max_block_length, 0:derv_length) 401 integer :: i, j 402 403 dervtmp = 0.0d0 404 405 call get_functional_derv(f, & 406 f_alda, & 407 f_xalda, & 408 order, & 409 block_length, & 410 w, & 411 r_0_t, & 412 z_0_t, & 413 derv_length, & 414 dervtmp, & 415 alda_real=is_fdersp_alda, & 416 alda_imag=is_fdersp_alda) 417 418 do i = 1, max_block_length 419 do j = 0, derv_length 420 derv(i,j) = derv(i,j) + dervtmp(i,j) 421 enddo 422 enddo 423 424 call get_functional_derv(f, & 425 f_alda, & 426 f_xalda, & 427 order, & 428 block_length, & 429 w, & 430 r_0, & 431 z_0, & 432 derv_length, & 433 dervtmp, & 434 alda_real=is_fdersp_alda, & 435 alda_imag=is_fdersp_alda) 436 437 do i = 1, max_block_length 438 do j = 0, derv_length 439 derv(i,j) = derv(i,j) - dervtmp(i,j) 440 enddo 441 enddo 442 443 end subroutine get_nadd_fun_derv 444 445 446! ---------------------------------------------------------------------------- 447 subroutine get_xcke_fun_derv(order, & 448 block_length, & 449 w, & 450 r_0, & 451 z_0, & 452 derv) 453! ---------------------------------------------------------------------------- 454 integer, intent(in) :: order 455 integer, intent(in) :: block_length 456 real(8), intent(in) :: w(*) 457 real(8), intent(in) :: r_0(*) 458 real(8), intent(in) :: z_0(*) 459 real(8), intent(inout) :: derv(max_block_length, 0:derv_length) 460 integer :: i, j 461 real(8) :: dervtmp(max_block_length, 0:derv_length) 462 463 derv = 0.0d0 464 dervtmp = 0.0d0 465 466 if (.not.skip_xc_contribution) then 467 call get_functional_derv(fde_xcf, & 468 fde_xcf_alda, & 469 fde_xcf_alda, & 470 order, & 471 block_length, & 472 w, & 473 r_0, & 474 z_0, & 475 derv_length, & 476 dervtmp, & 477 alda_real=is_fdersp_alda, & 478 alda_imag=is_fdersp_alda) 479 480 do i = 1, max_block_length 481 do j = 0, derv_length 482 derv(i,j) = derv(i,j) + dervtmp(i,j) 483 enddo 484 enddo 485 endif 486 487 if (.not.skip_kin_contribution) then 488 call get_functional_derv(fde_kef, & 489 fde_kef_alda, & 490 fde_kef_alda, & 491 order, & 492 block_length, & 493 w, & 494 r_0, & 495 z_0, & 496 derv_length, & 497 dervtmp, & 498 alda_real=is_fdersp_alda, & 499 alda_imag=is_fdersp_alda) 500 501 502 do i = 1, max_block_length 503 do j = 0, derv_length 504 derv(i,j) = derv(i,j) + dervtmp(i,j) 505 enddo 506 enddo 507 endif 508 509 end subroutine get_xcke_fun_derv 510 511end module fde_nadd_derv 512 513