TABLE OF CONTENTS
- ABINIT/m_gwls_ComputePoles
- m_gwls_ComputePoles/clean_degeneracy_table_for_poles
- m_gwls_ComputePoles/compute_pole_contribution
- m_gwls_ComputePoles/compute_Poles
- m_gwls_ComputePoles/generate_degeneracy_table_for_poles
ABINIT/m_gwls_ComputePoles [ Modules ]
NAME
m_gwls_ComputePoles
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2022 ABINIT group (JLJ, BR, MC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_gwls_ComputePoles 23 24 use m_gwls_utility 25 use m_gwls_wf 26 use m_gwls_hamiltonian 27 use m_gwls_lineqsolver 28 use m_gwls_polarisability 29 use m_gwls_GWlanczos 30 use m_gwls_GenerateEpsilon 31 use m_gwls_GWanalyticPart 32 use m_gwls_TimingLog 33 use m_gwls_LanczosBasis 34 35 use defs_basis 36 use defs_wvltypes 37 use m_abicore 38 use m_xmpi 39 use m_errors 40 41 use m_io_tools, only : get_unit 42 43 44 implicit none 45 save 46 private 47 48 integer :: number_of_denerate_sets 49 integer :: largest_degeneracy 50 51 integer, allocatable :: degeneracy_table(:,:) 52 integer, allocatable :: number_of_degenerate_states(:) 53 54 real(dp) :: En_m_omega_2 55 56 public :: compute_Poles 57 public :: generate_degeneracy_table_for_poles 58 public :: clean_degeneracy_table_for_poles 59 60 CONTAINS
m_gwls_ComputePoles/clean_degeneracy_table_for_poles [ Functions ]
[ Top ] [ m_gwls_ComputePoles ] [ Functions ]
NAME
clean_degeneracy_table_for_poles
FUNCTION
.
INPUTS
OUTPUT
SOURCE
256 subroutine clean_degeneracy_table_for_poles() 257 258 ! ************************************************************************* 259 260 if(allocated(degeneracy_table)) then 261 ABI_FREE(degeneracy_table) 262 end if 263 if(allocated(number_of_degenerate_states)) then 264 ABI_FREE(number_of_degenerate_states) 265 end if 266 267 end subroutine clean_degeneracy_table_for_poles
m_gwls_ComputePoles/compute_pole_contribution [ Functions ]
[ Top ] [ m_gwls_ComputePoles ] [ Functions ]
NAME
compute_pole_contribution
FUNCTION
.
INPUTS
OUTPUT
SOURCE
595 function compute_pole_contribution(epsilon_matrix_function,nseeds,kmax,seeds,debug) 596 !---------------------------------------------------------------------- 597 ! This routine computes the contribution to the poles energy 598 ! coming from the states in the seeds. 599 !---------------------------------------------------------------------- 600 interface 601 subroutine epsilon_matrix_function(v_out,v_in,l) 602 603 use defs_basis 604 605 integer, intent(in) :: l 606 complex(dp), intent(out) :: v_out(l) 607 complex(dp), intent(in) :: v_in(l) 608 609 end subroutine epsilon_matrix_function 610 end interface 611 612 real(dp) :: compute_pole_contribution 613 614 integer, intent(in) :: nseeds, kmax 615 complex(dpc), intent(in) :: seeds(npw_k,nseeds) 616 logical, intent(in) :: debug 617 618 ! local variables 619 620 integer :: mpi_communicator 621 622 complex(dpc),allocatable :: local_seeds(:,:) 623 624 complex(dpc),allocatable :: Lbasis(:,:) ! array containing the Lanczos basis 625 complex(dpc),allocatable :: alpha(:,:,:) 626 complex(dpc),allocatable :: beta (:,:,:) 627 real(dp), allocatable :: epsilon_eigenvalues(:) 628 629 real(dp):: matrix_elements 630 631 complex(dpc) :: cmplx_value 632 integer :: l, s 633 integer :: ierr 634 635 ! ************************************************************************* 636 637 ! compute the Lanczos basis 638 ABI_MALLOC(alpha,(nseeds,nseeds,kmax)) 639 ABI_MALLOC(beta ,(nseeds,nseeds,kmax)) 640 ABI_MALLOC(Lbasis,(npw_k,nseeds*kmax)) 641 ABI_MALLOC(local_seeds,(npw_k,nseeds)) 642 ABI_MALLOC(epsilon_eigenvalues, (nseeds*kmax)) 643 644 645 mpi_communicator = mpi_enreg%comm_bandfft !Missing maybe something for easy access of LA and FFT comms? 646 647 local_seeds(:,:) = seeds(:,:) 648 649 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, & 650 & local_seeds,alpha,beta,Lbasis) 651 652 write(std_out,*) "alpha:" 653 do l=1,kmax 654 do s=1,nseeds 655 write(std_out,*) alpha(:,s,l) 656 end do 657 write(std_out,*) " " 658 end do 659 660 write(std_out,*) "beta:" 661 do l=1,kmax 662 do s=1,nseeds 663 write(std_out,*) beta(:,s,l) 664 end do 665 write(std_out,*) " " 666 end do 667 668 ABI_FREE(local_seeds) 669 670 ! Diagonalize the epsilon matrix, which is banded 671 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug) 672 673 if (debug) then 674 call ritz_analysis_general(mpi_communicator ,epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues) 675 end if 676 677 compute_pole_contribution = zero 678 679 do l = 1, nseeds*kmax 680 681 matrix_elements = zero 682 683 do s = 1, nseeds 684 685 cmplx_value = complex_vector_product(seeds(:,s),Lbasis(:,l),npw_k) 686 687 call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors working on FFT! 688 689 matrix_elements = matrix_elements + abs(cmplx_value)**2 690 691 692 end do 693 694 695 compute_pole_contribution = compute_pole_contribution + & 696 matrix_elements *(one/epsilon_eigenvalues(l)-one) 697 end do 698 699 700 701 ABI_FREE(alpha) 702 ABI_FREE(beta) 703 ABI_FREE(Lbasis) 704 ABI_FREE(epsilon_eigenvalues) 705 706 end function compute_pole_contribution 707 708 end module m_gwls_ComputePoles
m_gwls_ComputePoles/compute_Poles [ Functions ]
[ Top ] [ m_gwls_ComputePoles ] [ Functions ]
NAME
compute_Poles
FUNCTION
.
INPUTS
OUTPUT
SOURCE
285 function compute_Poles(external_omega,kmax_poles,debug) 286 !---------------------------------------------------------------------- 287 ! This function extract the Pole contributions to the correlation 288 ! energy, as a function of the external frequency. 289 ! 290 ! The algorithm builds a Lanczos chain for each contributing subspace; 291 ! the number of steps is controlled by kmax_poles. 292 ! 293 ! This function will take in explicit arguments, as it is simpler 294 ! to do this than to define global arrays. 295 !---------------------------------------------------------------------- 296 real(dp) :: compute_Poles 297 298 real(dp), intent(in) :: external_omega 299 integer, intent(in) :: kmax_poles 300 logical, intent(in) :: debug 301 302 real(dp) :: energy_tolerance 303 real(dp) :: pole_contribution 304 305 306 307 integer :: number_of_seeds 308 integer :: i_set 309 integer :: n 310 311 real(dp) :: prefactor 312 real(dp) :: En_m_omega 313 314 logical :: pole_is_valence 315 logical :: pole_is_conduction 316 logical :: pole_is_in_gap 317 318 integer :: io_unit, i 319 character(128) :: filename 320 logical :: file_exists 321 322 323 324 complex(dpc), allocatable :: seeds(:,:) 325 326 ! ************************************************************************* 327 328 compute_Poles = zero 329 330 energy_tolerance = 1.0D-8 331 332 333 if (debug .and. mpi_enreg%me == 0) then 334 io_unit = get_unit() 335 336 i = 0 337 338 file_exists = .true. 339 do while (file_exists) 340 i = i+1 341 write (filename,'(A,I0.4,A)') "ComputePoles_",i,".log" 342 inquire(file=filename,exist=file_exists) 343 end do 344 345 346 open(io_unit,file=filename,status=files_status_new) 347 348 write(io_unit,10) " " 349 write(io_unit,10) "#===============================================================================================" 350 write(io_unit,10) "# ComputePoles: debug information for the pole computation " 351 write(io_unit,10) "# -------------------------------------------------------- " 352 write(io_unit,10) "# " 353 write(io_unit,10) "# This file contains data describing the computation of the pole contribution to the " 354 write(io_unit,10) "# correlation self energy. " 355 write(io_unit,10) "# " 356 write(io_unit,10) "#===============================================================================================" 357 write(io_unit,10) " " 358 write(io_unit,10) "#===============================================================================================" 359 write(io_unit,10) "# " 360 write(io_unit,10) "# parameters: " 361 write(io_unit,10) "# " 362 write(io_unit,11) "# external_omega = ",external_omega," Ha " 363 write(io_unit,12) "# kmax_poles = ",kmax_poles 364 write(io_unit,12) "# nbandv = ",nbandv 365 write(io_unit,10) "# " 366 write(io_unit,10) "#===============================================================================================" 367 write(io_unit,10) "# " 368 write(io_unit,10) "# DFT Eigenvalues (Ha) " 369 write(io_unit,10) "#===============================================================================================" 370 write(io_unit,13) eig(:) 371 flush(io_unit) 372 373 end if 374 375 376 377 !-------------------------------------------------------------------------------- 378 ! 379 ! Determine if external frequency corresponds to the valence or the conduction 380 ! manifold. 381 ! 382 !-------------------------------------------------------------------------------- 383 384 compute_Poles = zero 385 386 pole_is_conduction = .false. 387 pole_is_valence = .false. 388 pole_is_in_gap = .false. 389 390 391 392 ! Careful here! there may be only nbandv states in memory; nbandv+1 causes segfaults! 393 if ( external_omega <= eig(nbandv)) then 394 pole_is_valence = .true. 395 else if ( external_omega > eig(nbandv)) then 396 pole_is_conduction = .true. 397 else 398 pole_is_in_gap = .true. 399 end if 400 401 402 403 404 if (debug .and. mpi_enreg%me == 0 ) then 405 write(io_unit,10) "#====================================================================================================" 406 write(io_unit,10) "# " 407 write(io_unit,10) "# Determine where the external energy is: " 408 write(io_unit,10) "# " 409 write(io_unit,14) "# pole_is_valence = ",pole_is_valence 410 write(io_unit,14) "# pole_is_conduction = ",pole_is_conduction 411 write(io_unit,14) "# pole_is_in_gap = ",pole_is_in_gap 412 write(io_unit,10) "# " 413 write(io_unit,10) "#====================================================================================================" 414 flush(io_unit) 415 416 end if 417 418 if ( pole_is_in_gap) return 419 420 !-------------------------------------------------------------------------------- 421 ! 422 ! Loop on all degenerate sets 423 ! 424 !-------------------------------------------------------------------------------- 425 426 if (debug .and. mpi_enreg%me == 0 ) then 427 write(io_unit,10) "#====================================================================================================" 428 write(io_unit,10) "# " 429 write(io_unit,10) "# Iterating over all degenerate sets of eigenvalues: " 430 write(io_unit,10) "# " 431 write(io_unit,10) "#====================================================================================================" 432 flush(io_unit) 433 434 end if 435 436 437 do i_set =1, number_of_denerate_sets 438 439 n = degeneracy_table(i_set,1) 440 441 En_m_omega = eig(n)-external_omega 442 443 if (debug .and. mpi_enreg%me == 0) then 444 write(io_unit,12) "# i_set = ", i_set 445 write(io_unit,12) "# n = ",n 446 write(io_unit,16) "# eig(n)-omega = ",En_m_omega," Ha" 447 flush(io_unit) 448 end if 449 450 !------------------------------------------ 451 ! Test if we need to exit the loop 452 !------------------------------------------ 453 if (pole_is_valence ) then 454 455 ! If the pole is valence, get out when we enter conduction states 456 if ( n > nbandv ) then 457 if (debug.and. mpi_enreg%me == 0) then 458 write(io_unit,10) "#" 459 write(io_unit,10) "# n > nbandv : exit loop!" 460 flush(io_unit) 461 end if 462 463 exit 464 end if 465 466 ! if the valence energy is smaller than the external frequency, 467 ! then there is no contribution 468 469 if (En_m_omega < zero .and. abs(En_m_omega) > energy_tolerance) then 470 ! careful close to zero! 471 if (debug .and. mpi_enreg%me == 0) then 472 write(io_unit,10) "# " 473 write(io_unit,10) "# eig(n) < omega : cycle!" 474 flush(io_unit) 475 end if 476 cycle 477 end if 478 479 480 ! if we are still here, there is a valence contribution 481 prefactor = -one 482 483 else if ( pole_is_conduction ) then 484 485 ! If the pole is conduction, get out when the conduction state is 486 ! larger than the frequency (careful close to zero!) 487 if ( En_m_omega > energy_tolerance ) then 488 if (debug .and. mpi_enreg%me == 0) then 489 write(io_unit,10) "#" 490 write(io_unit,10) "# eig(n) > omega : exit!" 491 flush(io_unit) 492 end if 493 exit 494 end if 495 496 ! If the pole is conduction, there is no contribution while 497 ! we are in the valence states 498 if ( n <= nbandv ) then 499 if (debug .and. mpi_enreg%me == 0) then 500 write(io_unit,10) "#" 501 write(io_unit,10) "# n <= nbandv : cycle!" 502 flush(io_unit) 503 end if 504 505 cycle 506 end if 507 508 ! if we are still here, there is a conduction contribution 509 prefactor = one 510 end if 511 512 513 !------------------------------------------------- 514 ! If we made it this far, we have a contribution! 515 !------------------------------------------------- 516 517 if (abs(En_m_omega) < energy_tolerance ) then 518 519 if (debug .and. mpi_enreg%me == 0) then 520 write(io_unit,10) "# " 521 write(io_unit,10) "# En - omega ~ 0: pole at the origin, multiply by 1/2!" 522 flush(io_unit) 523 end if 524 525 ! The factor of 1/2 accounts for the fact that 526 ! the pole is at the origin! 527 prefactor = 0.5_dp*prefactor 528 end if 529 530 531 532 number_of_seeds = number_of_degenerate_states(i_set) 533 534 ABI_MALLOC(seeds, (npw_k,number_of_seeds)) 535 536 call get_seeds(n, number_of_seeds, seeds) !Missing wrappers 537 538 call set_dielectric_function_frequency([En_m_omega,zero]) 539 if (debug .and. mpi_enreg%me == 0) then 540 write(io_unit,10) "# Compute pole contribution:" 541 write(io_unit,12) "# number of seeds = ",number_of_seeds 542 write(io_unit,16) "# || seeds || = ",sqrt(sum(abs(seeds(:,:))**2)) !Missing xmpi_sum 543 write(io_unit,16) "# eig(n)-omega = ",En_m_omega, " Ha" 544 write(io_unit,17) "# prefactor = ",prefactor 545 flush(io_unit) 546 end if 547 if(dtset%zcut > tol12) activate_inf_shift_poles = .true. 548 En_m_omega_2 = En_m_omega 549 pole_contribution = & 550 compute_pole_contribution(matrix_function_epsilon_k, & 551 number_of_seeds, kmax_poles, & 552 seeds,debug) 553 if(dtset%zcut > tol12) activate_inf_shift_poles = .false. 554 if (debug .and. mpi_enreg%me == 0) then 555 write(io_unit,16) "# pole contribution = ",prefactor*pole_contribution, " Ha" 556 flush(io_unit) 557 end if 558 559 560 compute_Poles = compute_Poles + prefactor*pole_contribution 561 ABI_FREE(seeds) 562 end do 563 564 if (debug .and. mpi_enreg%me == 0) then 565 close(io_unit) 566 end if 567 568 569 10 format(A) 570 11 format(A,F8.4,A) 571 12 format(A,I5) 572 13 format(1000F16.8) 573 14 format(A,L10) 574 16 format(A,ES12.4,A) 575 17 format(A,F8.4) 576 577 end function compute_Poles
m_gwls_ComputePoles/generate_degeneracy_table_for_poles [ Functions ]
[ Top ] [ m_gwls_ComputePoles ] [ Functions ]
NAME
generate_degeneracy_table_for_poles
FUNCTION
.
INPUTS
OUTPUT
SOURCE
76 subroutine generate_degeneracy_table_for_poles(debug) 77 !---------------------------------------------------------------------- 78 ! This subroutine groups, once and for all, the indices of 79 ! degenerate eigenstates. This will be useful to compute the poles 80 ! 81 !---------------------------------------------------------------------- 82 83 logical, intent(in) :: debug 84 85 real(dp) :: degeneracy_tolerance 86 87 integer :: nbands 88 integer :: n, i 89 integer :: i_set, j_deg 90 integer :: n_degeneracy 91 92 real(dp) :: energy 93 94 integer :: io_unit 95 character(128) :: filename 96 logical :: file_exists 97 98 ! ************************************************************************* 99 100 101 degeneracy_tolerance = 1.0D-8 102 !-------------------------------------------------------------------------------- 103 ! 104 ! First, find the largest degeneracy in the eigenvalue spectrum 105 ! 106 !-------------------------------------------------------------------------------- 107 108 nbands = size(eig) 109 110 ! initialize 111 largest_degeneracy = 1 112 number_of_denerate_sets = 1 113 114 n_degeneracy = 1 115 energy = eig(1) 116 117 !============================================================ 118 ! Notes for the code block below: 119 ! 120 ! The electronic eigenvalues can be grouped in degenerate 121 ! sets. The code below counts how many such sets there are, 122 ! and how many eigenvalues belong to each set. 123 ! 124 ! It is important to have a robust algorithm to do this 125 ! properly. In particular, the edge case where the LAST SET 126 ! is degenerate must be treated with care . 127 128 ! The algorithm.I'm looking at at the time of this writing 129 ! has a bug in it and cannot handle a degenerate last set... 130 ! Let's fix that! 131 !============================================================ 132 133 do n = 2, nbands 134 if (abs(eig(n) - energy) < degeneracy_tolerance) then 135 ! A degenerate state! add one 136 n_degeneracy = n_degeneracy + 1 137 else 138 ! We are no longer degenerate. Update 139 if ( n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy 140 141 n_degeneracy = 1 142 energy = eig(n) 143 number_of_denerate_sets = number_of_denerate_sets + 1 144 end if 145 146 ! If this is the last index, update the largest_degeneracy if necessary 147 if ( n == nbands .and. n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy 148 149 end do 150 151 !-------------------------------------------------------------------------------- 152 ! 153 ! Allocate the array which will contain the indices of the degenerate 154 ! states, and populate it. 155 !-------------------------------------------------------------------------------- 156 ABI_MALLOC(degeneracy_table, (number_of_denerate_sets,largest_degeneracy)) 157 ABI_MALLOC(number_of_degenerate_states, (number_of_denerate_sets)) 158 159 degeneracy_table(:,:) = 0 160 161 i_set = 1 162 j_deg = 1 163 164 ! initialize 165 energy = eig(1) 166 167 degeneracy_table(i_set,j_deg) = 1 168 number_of_degenerate_states(i_set) = 1 169 170 do n = 2, nbands 171 172 if (abs(eig(n) - energy) < degeneracy_tolerance) then 173 ! A degenerate state! add one 174 j_deg = j_deg + 1 175 176 else 177 ! We are no longer degenerate. Update 178 j_deg = 1 179 i_set = i_set+1 180 energy = eig(n) 181 182 end if 183 184 185 number_of_degenerate_states(i_set) = j_deg 186 degeneracy_table(i_set,j_deg) = n 187 188 end do 189 190 if (debug .and. mpi_enreg%me == 0) then 191 io_unit = get_unit() 192 filename = "degeneracy_table.log" 193 194 i = 0 195 inquire(file=filename,exist=file_exists) 196 do while (file_exists) 197 i = i+1 198 write (filename,'(A,I0,A)') "degeneracy_table_",i,".log" 199 inquire(file=filename,exist=file_exists) 200 end do 201 202 io_unit = get_unit() 203 204 open(io_unit,file=filename,status=files_status_new) 205 206 write(io_unit,10) " " 207 write(io_unit,10) "#===============================================================================================" 208 write(io_unit,10) "# Degeneracy table : tabulate the degenerate states " 209 write(io_unit,10) "# ------------------------------------------------------- " 210 write(io_unit,10) "# " 211 write(io_unit,14) "# number_of_denerate_sets = ", number_of_denerate_sets 212 write(io_unit,10) "# " 213 write(io_unit,14) "# largest_degeneracy = ", largest_degeneracy 214 write(io_unit,10) "# " 215 write(io_unit,10) "# Eigenvalues (Ha) " 216 write(io_unit,10) "#===============================================================================================" 217 write(io_unit,16) eig(:) 218 219 220 221 write(io_unit,10) "#===============================================================================================" 222 write(io_unit,10) "# i_set number of states States " 223 write(io_unit,10) "#===============================================================================================" 224 flush(io_unit) 225 226 do i_set = 1, number_of_denerate_sets 227 write(io_unit,12) i_set, number_of_degenerate_states(i_set), degeneracy_table(i_set,:) 228 end do 229 230 flush(io_unit) 231 232 close(io_unit) 233 end if 234 235 10 format(A) 236 12 format(I5,10X,I5,15X,1000I5) 237 14 format(A,I5) 238 16 format(1000F12.8,2X) 239 240 end subroutine generate_degeneracy_table_for_poles