exchange_val_hf.F90 24.3 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 247 ``````#else `````` Matthias Redies committed Jul 15, 2019 248 249 250 `````` if (.not. mat_ex%l_real) THEN IF (kpts%bksym(ikpt0) > sym%nop) coulomb = conjg(coulomb) endif `````` Daniel Wortmann committed Jun 21, 2017 251 ``````#endif `````` Matthias Redies committed Jul 15, 2019 252 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 253 `````` `````` Matthias Redies committed Jul 15, 2019 254 `````` DO ibando = 1, mnobd, psize `````` Daniel Wortmann committed Jun 21, 2017 255 `````` `````` Matthias Redies committed Jul 15, 2019 256 `````` IF (mat_ex%l_real) THEN `````` Daniel Wortmann committed Jun 21, 2017 257 ``````#ifdef CPP_IRAPPROX `````` Matthias Redies committed Jul 15, 2019 258 259 `````` 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 260 ``````#else `````` Matthias Redies committed Jul 15, 2019 261 262 263 `````` 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 264 ``````#endif `````` Matthias Redies committed Jul 15, 2019 265 `````` ELSE `````` Daniel Wortmann committed Jun 21, 2017 266 ``````#ifdef CPP_IRAPPROX `````` Matthias Redies committed Jul 15, 2019 267 268 `````` 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 269 ``````#else `````` Matthias Redies committed Jul 15, 2019 270 271 `````` 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 272 ``````#endif `````` Matthias Redies committed Jul 15, 2019 273 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 274 `````` `````` Matthias Redies committed Jul 15, 2019 275 276 277 278 279 `````` ! 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 280 ``````#ifndef CPP_NOSPMVEC `````` Matthias Redies committed Jul 15, 2019 281 282 283 284 285 286 287 288 289 290 291 292 `````` 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 293 294 ``````#endif `````` Matthias Redies committed Jul 15, 2019 295 296 297 298 299 300 301 302 `````` ! 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 303 `````` IF (mat_ex%l_real) THEN `````` Matthias Redies committed Jul 15, 2019 304 `````` cprod_vv_r(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_r(:hybrid%nbasm(ikpt0), :, :) `````` Gregor Michalicek committed Jun 18, 2018 305 `````` ELSE `````` Matthias Redies committed Jul 15, 2019 306 307 308 309 310 `````` 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 311 `````` `````` Matthias Redies committed Jul 15, 2019 312 313 314 315 316 317 318 319 320 321 322 323 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 `````` ! 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 349 `````` `````` Gregor Michalicek committed Aug 22, 2018 350 351 352 353 354 355 356 ``````! 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 357 `````` ! add contribution of the gamma point to the different cases (exch_vv,exch_cv,exch_cc) `````` Daniel Wortmann committed Jun 21, 2017 358 `````` `````` Matthias Redies committed Jul 15, 2019 359 `````` ! valence-valence-valence-valence exchange `````` Gregor Michalicek committed Jun 18, 2018 360 `````` `````` Matthias Redies committed Jul 15, 2019 361 362 363 364 365 366 `````` 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 367 `````` `````` Matthias Redies committed Jul 15, 2019 368 369 370 `````` 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 371 `````` `````` Matthias Redies committed Jul 15, 2019 372 373 374 375 376 377 `````` ! 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 378 `````` END DO `````` Gregor Michalicek committed Jun 18, 2018 379 `````` `````` Matthias Redies committed Jul 15, 2019 380 381 382 383 `````` 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 384 `````` END IF `````` Daniel Wortmann committed Jun 21, 2017 385 `````` `````` Matthias Redies committed Jul 15, 2019 386 387 388 389 390 391 392 393 394 `````` !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 395 `````` `````` Matthias Redies committed Jul 15, 2019 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 `````` 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 419 `````` END DO `````` Matthias Redies committed Jul 15, 2019 420 421 422 423 424 425 426 `````` ! 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 427 428 `````` END IF END DO `````` Matthias Redies committed Jul 15, 2019 429 430 431 432 433 434 435 436 437 438 `````` 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 439 `````` `````` Matthias Redies committed Jul 15, 2019 440 441 442 `````` 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 443 `````` `````` Matthias Redies committed Jul 15, 2019 444 `````` ! tail correction/contribution from all other k-points (it goes into exchcorrect ) `````` Daniel Wortmann committed Jun 21, 2017 445 `````` `````` Matthias Redies committed Jul 15, 2019 446 `````` ! Analytic contribution `````` Daniel Wortmann committed Jun 21, 2017 447 `````` `````` Matthias Redies committed Jul 15, 2019 448 449 450 451 452 453 `````` 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 454 `````` `````` Matthias Redies committed Jul 15, 2019 455 `````` ! due to the symmetrization afterwards the factor 1/n_q(1) must be added `````` Daniel Wortmann committed Jun 21, 2017 456 `````` `````` Matthias Redies committed Jul 15, 2019 457 458 `````` 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 459 `````` `````` Matthias Redies committed Jul 15, 2019 460 461 462 `````` END DO !n2 END DO !n1 END IF ! xcpot%icorr .ne. icorr_hse `````` Daniel Wortmann committed Jun 21, 2017 463 `````` `````` Matthias Redies committed Jul 15, 2019 464 465 466 467 `````` 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 468 `````` `````` Gregor Michalicek committed Aug 22, 2018 469 470 471 472 473 474 475 ``````! 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 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 `````` ! 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 `````` Matthias Redies committed Jul 15, 2019 504 505 `````` expo = 5e-3 rrad = sqrt(-log(5e-3)/expo) `````` Matthias Redies committed Jul 15, 2019 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 `````` 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 533 `````` END DO `````` Matthias Redies committed Jul 15, 2019 534 `````` n = n + 1 `````` Gregor Michalicek committed Jun 18, 2018 535 536 `````` END DO `````` Matthias Redies committed Jul 15, 2019 537 `````` END SUBROUTINE calc_divergence `````` Gregor Michalicek committed Jun 18, 2018 538 539 `````` END MODULE m_exchange_valence_hf``````