exchange_val_hf.F90 24.7 KB
 Gregor Michalicek committed Jun 18, 2018 1 2 3 4 5 6 ``````!-------------------------------------------------------------------------------- ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany ! This file is part of FLEUR and available as free software under the conditions ! of the MIT license as expressed in the LICENSE file in more detail. !-------------------------------------------------------------------------------- `````` Matthias Redies committed Jul 15, 2019 7 ``````! Calculates the HF exchange term `````` Daniel Wortmann committed Jun 21, 2017 8 9 10 11 12 13 14 15 16 17 18 ``````! ! s s* s s* ! phi (r) phi (r) phi (r') phi (r') ! occ. n_1k n'k+q n'k+q n_2k ! exchange(n,q) = - SUM INT INT ------------------------------------------- dr dr' ! k,n' | r - r' | ! ! occ s s ~ ~ s s ! = - SUM SUM v < phi | phi M > < M phi | phi > ! k,n' I,J k,IJ n'k+q n_1k q,I q,J n_2k n'k+q ! `````` Matthias Redies committed Jul 15, 2019 19 ``````! for the different combinations of n_1 and n_2 and where n' runs only over the valence states. `````` Daniel Wortmann committed Jun 21, 2017 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 ``````! ( n_1,n_2: valence-valence, core-core,core-valence ) ! ! ! At the Gamma point (k=0) v diverges. After diagonalization of v at k=0 the divergence is ! restricted to the head element I=1. Furthermore, we expand <...> with kp perturbation theory. ! As a result, the total I=1 element is given by a sum of a divergent 1/k**2-term and an ! angular dependent term. The former is separated from the numerical k-summation and treated ! analytically while the latter is spherically averaged and added to the k=0 contribution of ! the numerical k-summation. (A better knowledge of the integrand's behavior at the BZ edges ! might further improve the integration.) ! ! The divergence at the Gamma point is integrated with one of the following algorithms: ! (1) Switching-Off Function ! In a sphere of radius k0=radshmin/2 a switching-off function g(k)=1-(k/k0)**n*(n+1-n*k/k0) ! (n=npot) is defined. The 1/k**2 divergence is subtracted from the BZ integral in the form ! g(k)/k**2 and integrated analytically. The non-divergent rest is integrated numerically. ! (2) Periodic Function (similar to the one used by Massidda PRB 48, 5058) ! The function F(k) = SUM(G) exp(-expo*|k+G|**3) / |k+G|**2 is subtracted from the BZ integral ! and integrated analytically. The non-divergent rest is integrated numerically. ! The parameter expo is chosen such that exp(-expo*q**3)=1/2 ! with q = radius of sphere with same volume as BZ. ! (3) Periodic Function (same as Massidda's) with expo->0 ! The function F(k) = lim(expo->0) SUM(G) exp(-expo*|k+G|**2) / |k+G|**2 is subtracted from ! the BZ integral and integrated analytically. The contribution to the BZ integral including ! the "tail" is ! vol/(8*pi**3) INT F(k) d^3k - P SUM(k) F(k) ( P = principal value ) . ! For expo->0 the two terms diverge. Therefore a cutoff radius q0 is introduced and related to `````` Matthias Redies committed Jul 15, 2019 47 ``````! expo by exp(-expo*q0**2)=delta ( delta = small value, e.g., delta = 1*10.0**-10 ) . `````` Daniel Wortmann committed Jun 21, 2017 48 49 50 51 ``````! The resulting formula ! vol/(4*pi**1.5*sqrt(expo)) * erf(sqrt(a)*q0) - sum(q,0 sym%nop) .and. (.not. mat_ex%l_real)) THEN coulomb_mt2_c = conjg(coulomb_mt2_c) coulomb_mtir_c = conjg(coulomb_mtir_c) END IF `````` Daniel Wortmann committed Jun 21, 2017 259 ``````#else `````` Matthias Redies committed Jul 15, 2019 260 261 262 `````` if (.not. mat_ex%l_real) THEN IF (kpts%bksym(ikpt0) > sym%nop) coulomb = conjg(coulomb) endif `````` Daniel Wortmann committed Jun 21, 2017 263 ``````#endif `````` Matthias Redies committed Jul 15, 2019 264 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 265 `````` `````` Matthias Redies committed Jul 15, 2019 266 `````` DO ibando = 1, mnobd, psize `````` Daniel Wortmann committed Jun 21, 2017 267 `````` `````` Matthias Redies committed Jul 15, 2019 268 `````` IF (mat_ex%l_real) THEN `````` Daniel Wortmann committed Jun 21, 2017 269 ``````#ifdef CPP_IRAPPROX `````` Matthias Redies committed Jul 15, 2019 270 271 `````` CALL wavefproducts_inv(1, hybdat, dimension, input, jsp, atoms, lapw, obsolete, kpts, nk, ikpt0, & mnobd, hybrid, parent, cell, sym, noco, nkqpt, cprod_vv) `````` Daniel Wortmann committed Jun 21, 2017 272 ``````#else `````` Matthias Redies committed Jul 15, 2019 273 274 275 `````` CALL wavefproducts_inv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, dimension, input, jsp, atoms, & lapw, kpts, nk, ikpt0, hybdat, mnobd, hybrid, parent, cell, hybrid%nbasp, sym, & noco, nkqpt, cprod_vv_r) `````` Daniel Wortmann committed Jun 21, 2017 276 ``````#endif `````` Matthias Redies committed Jul 15, 2019 277 `````` ELSE `````` Daniel Wortmann committed Jun 21, 2017 278 ``````#ifdef CPP_IRAPPROX `````` Matthias Redies committed Jul 15, 2019 279 280 `````` CALL wavefproducts_noinv(1, hybdat, nk, ikpt0, dimension, input, jsp, cell, atoms, hybrid, kpts, mnobd, lapw, sym, noco, nkqpt, cprod_vv) `````` Daniel Wortmann committed Jun 21, 2017 281 ``````#else `````` Matthias Redies committed Jul 15, 2019 282 283 `````` CALL wavefproducts_noinv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, nk, ikpt0, dimension, input, jsp, &!jsp,& cell, atoms, hybrid, hybdat, kpts, mnobd, lapw, sym, hybrid%nbasp, noco, nkqpt, cprod_vv_c) `````` Daniel Wortmann committed Jun 21, 2017 284 ``````#endif `````` Matthias Redies committed Jul 15, 2019 285 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 286 `````` `````` Matthias Redies committed Jul 15, 2019 287 288 289 290 291 `````` ! The sparse matrix technique is not feasible for the HSE ! functional. Thus, a dynamic adjustment is implemented ! The mixed basis functions and the potential difference ! are Fourier transformed, so that the exchange can be calculated ! in Fourier space `````` Daniel Wortmann committed Jun 21, 2017 292 ``````#ifndef CPP_NOSPMVEC `````` Matthias Redies committed Jul 15, 2019 293 294 295 296 297 298 299 300 301 302 303 304 `````` IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN iband1 = hybrid%nobd(nkqpt) exch_vv = exch_vv + & dynamic_hse_adjustment(atoms%rmsh, atoms%rmt, atoms%dx, atoms%jri, atoms%jmtd, kpts%bkf(:, ikpt0), ikpt0, & kpts%nkptf, cell%bmat, cell%omtil, atoms%ntype, atoms%neq, atoms%nat, atoms%taual, & hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, hybrid%maxindxm1, hybrid%gptm, & hybrid%ngptm(ikpt0), hybrid%pgptm(:, ikpt0), hybrid%gptmd, hybrid%basm1, & hybrid%nbasm(ikpt0), iband1, hybrid%nbands(nk), nsest, ibando, psize, indx_sest, & atoms%invsat, sym%invsatnr, mpi%irank, cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), & cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), mat_ex%l_real, wl_iks(:iband1, nkqpt), n_q(ikpt)) END IF `````` Daniel Wortmann committed Jun 21, 2017 305 306 ``````#endif `````` Matthias Redies committed Jul 15, 2019 307 308 309 310 311 312 313 314 `````` ! the Coulomb matrix is only evaluated at the irrecuible k-points ! bra_trafo transforms cprod instead of rotating the Coulomb matrix ! from IBZ to current k-point IF (kpts%bkp(ikpt0) /= ikpt0) THEN CALL bra_trafo2(mat_ex%l_real, carr3_vv_r(:hybrid%nbasm(ikpt0), :, :), cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), & carr3_vv_c(:hybrid%nbasm(ikpt0), :, :), cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), & hybrid%nbasm(ikpt0), psize, hybrid%nbands(nk), kpts%bkp(ikpt0), ikpt0, kpts%bksym(ikpt0), sym, & hybrid, kpts, cell, atoms, phase_vv) `````` Gregor Michalicek committed Jun 18, 2018 315 `````` IF (mat_ex%l_real) THEN `````` Matthias Redies committed Jul 15, 2019 316 `````` cprod_vv_r(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_r(:hybrid%nbasm(ikpt0), :, :) `````` Gregor Michalicek committed Jun 18, 2018 317 `````` ELSE `````` Matthias Redies committed Jul 15, 2019 318 319 320 321 322 `````` cprod_vv_c(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_c(:hybrid%nbasm(ikpt0), :, :) ENDIF ELSE phase_vv(:, :) = (1.0, 0.0) END IF `````` Daniel Wortmann committed Jun 21, 2017 323 `````` `````` Matthias Redies committed Jul 15, 2019 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 `````` ! calculate exchange matrix at ikpt0 call timestart("exchange matrix") DO n1 = 1, hybrid%nbands(nk) DO iband = 1, psize IF ((ibando + iband - 1) > hybrid%nobd(nkqpt)) CYCLE cdum = wl_iks(ibando + iband - 1, nkqpt)*conjg(phase_vv(iband, n1))/n_q(ikpt) IF (mat_ex%l_real) THEN carr1_v_r(:n) = 0 CALL spmvec_invs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, & coulomb_mtir_r, cprod_vv_r(:n, iband, n1), carr1_v_r(:n)) ELSE carr1_v_c(:n) = 0 CALL spmvec_noinvs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, & coulomb_mtir_c, cprod_vv_c(:n, iband, n1), carr1_v_c(:n)) END IF IF (mat_ex%l_real) THEN DO n2 = 1, nsest(n1)!n1 nn2 = indx_sest(n2, n1) exch_vv(nn2, n1) = exch_vv(nn2, n1) + cdum*phase_vv(iband, nn2)* & dotprod(carr1_v_r(:n), cprod_vv_r(:n, iband, nn2)) END DO !n2 ELSE DO n2 = 1, nsest(n1)!n1 nn2 = indx_sest(n2, n1) exch_vv(nn2, n1) = exch_vv(nn2, n1) + cdum*phase_vv(iband, nn2)* & dotprod(carr1_v_c(:n), cprod_vv_c(:n, iband, nn2)) END DO !n2 END IF END DO END DO !n1 call timestop("exchange matrix") END DO !ibando END DO !ikpt `````` Daniel Wortmann committed Jun 21, 2017 361 `````` `````` Gregor Michalicek committed Aug 22, 2018 362 363 364 365 366 367 368 ``````! WRITE(7001,'(a,i7)') 'nk: ', nk ! DO n1=1,hybrid%nbands(nk) ! DO n2=1,n1 ! WRITE(7001,'(2i7,2f15.8)') n2, n1, exch_vv(n2,n1) ! END DO ! END DO `````` Matthias Redies committed Jul 15, 2019 369 `````` ! add contribution of the gamma point to the different cases (exch_vv,exch_cv,exch_cc) `````` Daniel Wortmann committed Jun 21, 2017 370 `````` `````` Matthias Redies committed Jul 15, 2019 371 `````` ! valence-valence-valence-valence exchange `````` Gregor Michalicek committed Jun 18, 2018 372 `````` `````` Matthias Redies committed Jul 15, 2019 373 374 375 376 377 378 `````` IF ((.not. xcpot%is_name("hse")) .AND. (.not. xcpot%is_name("vhse"))) THEN ! no gamma point correction needed for HSE functional IF (zero_order .and. .not. ibs_corr) THEN WRITE (6, '(A)') ' Take zero order terms into account.' ELSE IF (zero_order .and. ibs_corr) THEN WRITE (6, '(A)') ' Take zero order terms and ibs-correction into account.' END IF `````` Gregor Michalicek committed Jun 18, 2018 379 `````` `````` Matthias Redies committed Jul 15, 2019 380 381 382 `````` IF (zero_order) THEN CALL dwavefproducts(dcprod, nk, 1, hybrid%nbands(nk), 1, hybrid%nbands(nk), .false., atoms, hybrid, & cell, hybdat, kpts, kpts%nkpt, lapw, dimension, jsp, eig_irr) `````` Gregor Michalicek committed Jun 18, 2018 383 `````` `````` Matthias Redies committed Jul 15, 2019 384 385 386 387 388 389 `````` ! make dcprod hermitian DO n1 = 1, hybrid%nbands(nk) DO n2 = 1, n1 dcprod(n1, n2, :) = (dcprod(n1, n2, :) - conjg(dcprod(n2, n1, :)))/2 dcprod(n2, n1, :) = -conjg(dcprod(n1, n2, :)) END DO `````` Daniel Wortmann committed Jun 21, 2017 390 `````` END DO `````` Gregor Michalicek committed Jun 18, 2018 391 `````` `````` Matthias Redies committed Jul 15, 2019 392 393 394 395 `````` IF (ibs_corr) THEN CALL ibs_correction(nk, atoms, dimension, input, jsp, hybdat, hybrid, lapw, kpts, kpts%nkpt, cell, mnobd, & sym, proj_ibsc, olap_ibsc) END IF `````` Gregor Michalicek committed Jun 18, 2018 396 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 397 `````` `````` Matthias Redies committed Jul 15, 2019 398 399 400 401 402 403 404 405 406 `````` !This should be done with w_iks I guess!TODO occup = .false. DO i = 1, hybrid%ne_eig(nk) IF (results%ef >= eig_irr(i, nk)) THEN occup(i) = .true. ELSE IF ((eig_irr(i, nk) - results%ef) <= 1E-06) THEN occup(i) = .true. END IF END DO `````` Gregor Michalicek committed Jun 18, 2018 407 `````` `````` Matthias Redies committed Jul 15, 2019 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 `````` DO n1 = 1, hybrid%nbands(nk) DO n2 = 1, nsest(n1)!n1 nn2 = indx_sest(n2, n1) exchcorrect = 0 exch0 = 0 ! if zero_order = .true. add averaged k-dependent term to the numerical integration at Gamma-point contribution ! if we start with a system with a small DFT band gap (like GaAs), the contribution ! of the highest occupied and lowest unoccupied state in Hessian is typically ! large; a correct numerical integration requires a dense k-point mesh, so ! we don't add the contribution exchcorrect for such materials IF (zero_order) THEN hessian = 0 IF (occup(n1) .and. occup(nn2)) THEN DO i = 1, 3 j = i DO iband = 1, hybrid%nbands(nk) IF (occup(iband)) THEN hessian(i, j) = hessian(i, j) + conjg(dcprod(iband, n1, i))*dcprod(iband, nn2, j) END IF hessian(i, j) = hessian(i, j) - dcprod(iband, nn2, i)*conjg(dcprod(iband, n1, j)) `````` Gregor Michalicek committed Jun 18, 2018 431 `````` END DO `````` Matthias Redies committed Jul 15, 2019 432 433 434 435 436 437 438 `````` ! ibs correction IF (ibs_corr) THEN hessian(i, j) = hessian(i, j) - olap_ibsc(i, j, n1, nn2)/cell%omtil DO iband = 1, hybrid%nbands(nk) hessian(i, j) = hessian(i, j) + conjg(proj_ibsc(i, nn2, iband))*proj_ibsc(j, n1, iband)/cell%omtil END DO `````` Gregor Michalicek committed Jun 18, 2018 439 440 `````` END IF END DO `````` Matthias Redies committed Jul 15, 2019 441 442 443 444 445 446 447 448 449 450 `````` ELSE DO i = 1, 3 j = i DO iband = 1, hybrid%nbands(nk) IF (occup(iband)) THEN hessian(i, j) = hessian(i, j) + conjg(dcprod(iband, n1, i))*dcprod(iband, nn2, j) END IF END DO END DO END IF `````` Daniel Wortmann committed Jun 21, 2017 451 `````` `````` Matthias Redies committed Jul 15, 2019 452 453 454 `````` exchcorrect(1) = fpi_const/3*(hessian(1, 1) + hessian(2, 2) + hessian(3, 3)) exch0 = exchcorrect(1)/kpts%nkptf END IF `````` Daniel Wortmann committed Jun 21, 2017 455 `````` `````` Matthias Redies committed Jul 15, 2019 456 `````` ! tail correction/contribution from all other k-points (it goes into exchcorrect ) `````` Daniel Wortmann committed Jun 21, 2017 457 `````` `````` Matthias Redies committed Jul 15, 2019 458 `````` ! Analytic contribution `````` Daniel Wortmann committed Jun 21, 2017 459 `````` `````` Matthias Redies committed Jul 15, 2019 460 461 462 463 464 465 `````` cdum2 = 0 !multiply divergent contribution with occupation number; !this only affects metals IF (n1 == nn2) THEN cdum2 = fpi_const/cell%omtil*divergence*wl_iks(n1, nk)*kpts%nkptf END IF `````` Daniel Wortmann committed Jun 21, 2017 466 `````` `````` Matthias Redies committed Jul 15, 2019 467 `````` ! due to the symmetrization afterwards the factor 1/n_q(1) must be added `````` Daniel Wortmann committed Jun 21, 2017 468 `````` `````` Matthias Redies committed Jul 15, 2019 469 470 `````` IF (n1 == nn2) hybrid%div_vv(n1, nk, jsp) = REAL(cdum2) exch_vv(nn2, n1) = exch_vv(nn2, n1) + (exch0 + cdum2)/n_q(1) `````` Daniel Wortmann committed Jun 21, 2017 471 `````` `````` Matthias Redies committed Jul 15, 2019 472 473 474 `````` END DO !n2 END DO !n1 END IF ! xcpot%icorr .ne. icorr_hse `````` Daniel Wortmann committed Jun 21, 2017 475 `````` `````` Matthias Redies committed Jul 15, 2019 476 477 478 479 `````` IF (mat_ex%l_real) THEN IF (any(abs(aimag(exch_vv)) > 1E-08)) CALL judft_warn('unusally large imaginary part of exch_vv', & calledby='exchange_val_hf.F90') END IF `````` Daniel Wortmann committed Jun 21, 2017 480 `````` `````` Gregor Michalicek committed Aug 22, 2018 481 482 483 484 485 486 487 ``````! WRITE(7000,'(a,i7)') 'nk: ', nk ! DO n1=1,hybrid%nbands(nk) ! DO n2=1,n1 ! WRITE(7000,'(2i7,2f15.8)') n2, n1, exch_vv(n2,n1) ! END DO ! END DO `````` Matthias Redies committed Jul 15, 2019 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 `````` ! write exch_vv in mat_ex CALL mat_ex%alloc(matsize1=hybrid%nbands(nk)) IF (mat_ex%l_real) THEN mat_ex%data_r = exch_vv ELSE mat_ex%data_c = exch_vv END IF CALL timestop("valence exchange calculation") END SUBROUTINE exchange_valence_hf SUBROUTINE calc_divergence(cell, kpts, divergence) USE m_util, ONLY: cerf USE m_types USE m_constants IMPLICIT NONE TYPE(t_cell), INTENT(IN) :: cell TYPE(t_kpts), INTENT(IN) :: kpts REAL, INTENT(OUT) :: divergence INTEGER :: ix, iy, iz, sign, n logical :: found REAL :: expo, rrad, k(3), kv1(3), kv2(3), kv3(3), knorm2 COMPLEX :: cdum expo = 5*10.0**-3 rrad = sqrt(-log(5*10.0**-3)/expo) cdum = sqrt(expo)*rrad divergence = cell%omtil/(tpi_const**2)*sqrt(pi_const/expo)*cerf(cdum) rrad = rrad**2 kv1 = cell%bmat(1, :)/kpts%nkpt3(1) kv2 = cell%bmat(2, :)/kpts%nkpt3(2) kv3 = cell%bmat(3, :)/kpts%nkpt3(3) n = 1 found = .true. DO WHILE (found) found = .false. DO ix = -n, n DO iy = -(n - abs(ix)), n - abs(ix) iz = n - abs(ix) - abs(iy) DO sign = -1, 1, 2 iz = sign*iz k(1) = ix*kv1(1) + iy*kv2(1) + iz*kv3(1) k(2) = ix*kv1(2) + iy*kv2(2) + iz*kv3(2) k(3) = ix*kv1(3) + iy*kv2(3) + iz*kv3(3) knorm2 = k(1)**2 + k(2)**2 + k(3)**2 IF (knorm2 < rrad) THEN found = .true. divergence = divergence - exp(-expo*knorm2)/knorm2/kpts%nkptf END IF IF (iz == 0) exit END DO END DO `````` Gregor Michalicek committed Jun 18, 2018 545 `````` END DO `````` Matthias Redies committed Jul 15, 2019 546 `````` n = n + 1 `````` Gregor Michalicek committed Jun 18, 2018 547 548 `````` END DO `````` Matthias Redies committed Jul 15, 2019 549 `````` END SUBROUTINE calc_divergence `````` Gregor Michalicek committed Jun 18, 2018 550 551 `````` END MODULE m_exchange_valence_hf``````