TABLE OF CONTENTS
- ABINIT/gwls_ComputePoles
- gwls_ComputePoles/clean_degeneracy_table_for_poles
- gwls_ComputePoles/compute_pole_contribution
- gwls_ComputePoles/compute_Poles
- gwls_ComputePoles/generate_degeneracy_table_for_poles
ABINIT/gwls_ComputePoles [ Modules ]
NAME
gwls_ComputePoles
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2018 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 .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 28 29 module gwls_ComputePoles 30 31 ! local modules 32 use m_gwls_utility 33 use gwls_wf 34 use gwls_hamiltonian 35 use gwls_lineqsolver 36 use gwls_polarisability 37 !use gwls_ComplementSpacePolarizability 38 use gwls_GWlanczos 39 use gwls_GenerateEpsilon 40 use gwls_GWanalyticPart 41 use gwls_TimingLog 42 use gwls_LanczosBasis 43 44 45 ! abinit modules 46 use defs_basis 47 use defs_datatypes 48 use defs_abitypes 49 use defs_wvltypes 50 use m_profiling_abi 51 use m_xmpi 52 use m_pawang 53 use m_errors 54 55 use m_io_tools, only : get_unit 56 use m_paw_dmft, only: paw_dmft_type 57 use m_ebands, only : ebands_init, ebands_free 58 59 use m_gaussian_quadrature, only: gaussian_quadrature_gegenbauer, gaussian_quadrature_legendre 60 61 62 implicit none 63 save 64 private 65 66 integer :: number_of_denerate_sets 67 integer :: largest_degeneracy 68 69 integer, allocatable :: degeneracy_table(:,:) 70 integer, allocatable :: number_of_degenerate_states(:) 71 72 real(dp) :: En_m_omega_2 73 74 public :: compute_Poles 75 public :: generate_degeneracy_table_for_poles 76 public :: clean_degeneracy_table_for_poles 77 78 CONTAINS
gwls_ComputePoles/clean_degeneracy_table_for_poles [ Functions ]
[ Top ] [ gwls_ComputePoles ] [ Functions ]
NAME
clean_degeneracy_table_for_poles
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_ComputeCorrelationEnergy
CHILDREN
block_lanczos_algorithm,diagonalize_lanczos_banded ritz_analysis_general,xmpi_sum
SOURCE
297 subroutine clean_degeneracy_table_for_poles() 298 299 300 !This section has been created automatically by the script Abilint (TD). 301 !Do not modify the following lines by hand. 302 #undef ABI_FUNC 303 #define ABI_FUNC 'clean_degeneracy_table_for_poles' 304 !End of the abilint section 305 306 implicit none 307 ! ************************************************************************* 308 309 if(allocated(degeneracy_table)) then 310 ABI_DEALLOCATE(degeneracy_table) 311 end if 312 if(allocated(number_of_degenerate_states)) then 313 ABI_DEALLOCATE(number_of_degenerate_states) 314 end if 315 316 end subroutine clean_degeneracy_table_for_poles
gwls_ComputePoles/compute_pole_contribution [ Functions ]
[ Top ] [ gwls_ComputePoles ] [ Functions ]
NAME
compute_pole_contribution
FUNCTION
.
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
660 function compute_pole_contribution(epsilon_matrix_function,nseeds,kmax,seeds,debug) 661 !---------------------------------------------------------------------- 662 ! This routine computes the contribution to the poles energy 663 ! coming from the states in the seeds. 664 !---------------------------------------------------------------------- 665 666 !This section has been created automatically by the script Abilint (TD). 667 !Do not modify the following lines by hand. 668 #undef ABI_FUNC 669 #define ABI_FUNC 'compute_pole_contribution' 670 !End of the abilint section 671 672 implicit none 673 interface 674 subroutine epsilon_matrix_function(v_out,v_in,l) 675 676 use defs_basis 677 678 integer, intent(in) :: l 679 complex(dp), intent(out) :: v_out(l) 680 complex(dp), intent(in) :: v_in(l) 681 682 end subroutine epsilon_matrix_function 683 end interface 684 685 real(dp) :: compute_pole_contribution 686 687 integer, intent(in) :: nseeds, kmax 688 complex(dpc), intent(in) :: seeds(npw_k,nseeds) 689 logical, intent(in) :: debug 690 691 ! local variables 692 693 integer :: mpi_communicator 694 695 complex(dpc),allocatable :: local_seeds(:,:) 696 697 complex(dpc),allocatable :: Lbasis(:,:) ! array containing the Lanczos basis 698 complex(dpc),allocatable :: alpha(:,:,:) 699 complex(dpc),allocatable :: beta (:,:,:) 700 real(dp), allocatable :: epsilon_eigenvalues(:) 701 702 real(dp):: matrix_elements 703 704 complex(dpc) :: cmplx_value 705 integer :: l, s 706 integer :: ierr 707 708 ! ************************************************************************* 709 710 ! compute the Lanczos basis 711 ABI_ALLOCATE(alpha,(nseeds,nseeds,kmax)) 712 ABI_ALLOCATE(beta ,(nseeds,nseeds,kmax)) 713 ABI_ALLOCATE(Lbasis,(npw_k,nseeds*kmax)) 714 ABI_ALLOCATE(local_seeds,(npw_k,nseeds)) 715 ABI_ALLOCATE(epsilon_eigenvalues, (nseeds*kmax)) 716 717 718 mpi_communicator = mpi_enreg%comm_bandfft !Missing maybe something for easy access of LA and FFT comms? 719 720 local_seeds(:,:) = seeds(:,:) 721 722 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, & 723 & local_seeds,alpha,beta,Lbasis) 724 725 write(std_out,*) "alpha:" 726 do l=1,kmax 727 do s=1,nseeds 728 write(std_out,*) alpha(:,s,l) 729 end do 730 write(std_out,*) " " 731 end do 732 733 write(std_out,*) "beta:" 734 do l=1,kmax 735 do s=1,nseeds 736 write(std_out,*) beta(:,s,l) 737 end do 738 write(std_out,*) " " 739 end do 740 741 ABI_DEALLOCATE(local_seeds) 742 743 ! Diagonalize the epsilon matrix, which is banded 744 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug) 745 746 if (debug) then 747 call ritz_analysis_general(mpi_communicator ,epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues) 748 end if 749 750 compute_pole_contribution = zero 751 752 do l = 1, nseeds*kmax 753 754 matrix_elements = zero 755 756 do s = 1, nseeds 757 758 cmplx_value = complex_vector_product(seeds(:,s),Lbasis(:,l),npw_k) 759 760 call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors working on FFT! 761 762 matrix_elements = matrix_elements + abs(cmplx_value)**2 763 764 765 end do 766 767 768 compute_pole_contribution = compute_pole_contribution + & 769 matrix_elements *(one/epsilon_eigenvalues(l)-one) 770 end do 771 772 773 774 ABI_DEALLOCATE(alpha) 775 ABI_DEALLOCATE(beta) 776 ABI_DEALLOCATE(Lbasis) 777 ABI_DEALLOCATE(epsilon_eigenvalues) 778 779 end function compute_pole_contribution 780 781 end module gwls_ComputePoles
gwls_ComputePoles/compute_Poles [ Functions ]
[ Top ] [ gwls_ComputePoles ] [ Functions ]
NAME
compute_Poles
FUNCTION
.
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
338 function compute_Poles(external_omega,kmax_poles,debug) 339 !---------------------------------------------------------------------- 340 ! This function extract the Pole contributions to the correlation 341 ! energy, as a function of the external frequency. 342 ! 343 ! The algorithm builds a Lanczos chain for each contributing subspace; 344 ! the number of steps is controlled by kmax_poles. 345 ! 346 ! This function will take in explicit arguments, as it is simpler 347 ! to do this than to define global arrays. 348 !---------------------------------------------------------------------- 349 350 !This section has been created automatically by the script Abilint (TD). 351 !Do not modify the following lines by hand. 352 #undef ABI_FUNC 353 #define ABI_FUNC 'compute_Poles' 354 !End of the abilint section 355 356 implicit none 357 real(dp) :: compute_Poles 358 359 real(dp), intent(in) :: external_omega 360 integer, intent(in) :: kmax_poles 361 logical, intent(in) :: debug 362 363 real(dp) :: energy_tolerance 364 real(dp) :: pole_contribution 365 366 367 368 integer :: number_of_seeds 369 integer :: i_set 370 integer :: n 371 372 real(dp) :: prefactor 373 real(dp) :: En_m_omega 374 375 logical :: pole_is_valence 376 logical :: pole_is_conduction 377 logical :: pole_is_in_gap 378 379 integer :: io_unit, i 380 character(128) :: filename 381 logical :: file_exists 382 383 384 385 complex(dpc), allocatable :: seeds(:,:) 386 387 ! ************************************************************************* 388 389 compute_Poles = zero 390 391 energy_tolerance = 1.0D-8 392 393 394 if (debug .and. mpi_enreg%me == 0) then 395 io_unit = get_unit() 396 397 i = 0 398 399 file_exists = .true. 400 do while (file_exists) 401 i = i+1 402 write (filename,'(A,I0.4,A)') "ComputePoles_",i,".log" 403 inquire(file=filename,exist=file_exists) 404 end do 405 406 407 open(io_unit,file=filename,status=files_status_new) 408 409 write(io_unit,10) " " 410 write(io_unit,10) "#===============================================================================================" 411 write(io_unit,10) "# ComputePoles: debug information for the pole computation " 412 write(io_unit,10) "# -------------------------------------------------------- " 413 write(io_unit,10) "# " 414 write(io_unit,10) "# This file contains data describing the computation of the pole contribution to the " 415 write(io_unit,10) "# correlation self energy. " 416 write(io_unit,10) "# " 417 write(io_unit,10) "#===============================================================================================" 418 write(io_unit,10) " " 419 write(io_unit,10) "#===============================================================================================" 420 write(io_unit,10) "# " 421 write(io_unit,10) "# parameters: " 422 write(io_unit,10) "# " 423 write(io_unit,11) "# external_omega = ",external_omega," Ha " 424 write(io_unit,12) "# kmax_poles = ",kmax_poles 425 write(io_unit,12) "# nbandv = ",nbandv 426 write(io_unit,10) "# " 427 write(io_unit,10) "#===============================================================================================" 428 write(io_unit,10) "# " 429 write(io_unit,10) "# DFT Eigenvalues (Ha) " 430 write(io_unit,10) "#===============================================================================================" 431 write(io_unit,13) eig(:) 432 flush(io_unit) 433 434 end if 435 436 437 438 !-------------------------------------------------------------------------------- 439 ! 440 ! Determine if external frequency corresponds to the valence or the conduction 441 ! manifold. 442 ! 443 !-------------------------------------------------------------------------------- 444 445 compute_Poles = zero 446 447 pole_is_conduction = .false. 448 pole_is_valence = .false. 449 pole_is_in_gap = .false. 450 451 452 453 ! Careful here! there may be only nbandv states in memory; nbandv+1 causes segfaults! 454 if ( external_omega <= eig(nbandv)) then 455 pole_is_valence = .true. 456 else if ( external_omega > eig(nbandv)) then 457 pole_is_conduction = .true. 458 else 459 pole_is_in_gap = .true. 460 end if 461 462 463 464 465 if (debug .and. mpi_enreg%me == 0 ) then 466 write(io_unit,10) "#====================================================================================================" 467 write(io_unit,10) "# " 468 write(io_unit,10) "# Determine where the external energy is: " 469 write(io_unit,10) "# " 470 write(io_unit,14) "# pole_is_valence = ",pole_is_valence 471 write(io_unit,14) "# pole_is_conduction = ",pole_is_conduction 472 write(io_unit,14) "# pole_is_in_gap = ",pole_is_in_gap 473 write(io_unit,10) "# " 474 write(io_unit,10) "#====================================================================================================" 475 flush(io_unit) 476 477 end if 478 479 if ( pole_is_in_gap) return 480 481 !-------------------------------------------------------------------------------- 482 ! 483 ! Loop on all degenerate sets 484 ! 485 !-------------------------------------------------------------------------------- 486 487 if (debug .and. mpi_enreg%me == 0 ) then 488 write(io_unit,10) "#====================================================================================================" 489 write(io_unit,10) "# " 490 write(io_unit,10) "# Iterating over all degenerate sets of eigenvalues: " 491 write(io_unit,10) "# " 492 write(io_unit,10) "#====================================================================================================" 493 flush(io_unit) 494 495 end if 496 497 498 do i_set =1, number_of_denerate_sets 499 500 n = degeneracy_table(i_set,1) 501 502 En_m_omega = eig(n)-external_omega 503 504 if (debug .and. mpi_enreg%me == 0) then 505 write(io_unit,12) "# i_set = ", i_set 506 write(io_unit,12) "# n = ",n 507 write(io_unit,16) "# eig(n)-omega = ",En_m_omega," Ha" 508 flush(io_unit) 509 end if 510 511 !------------------------------------------ 512 ! Test if we need to exit the loop 513 !------------------------------------------ 514 if (pole_is_valence ) then 515 516 ! If the pole is valence, get out when we enter conduction states 517 if ( n > nbandv ) then 518 if (debug.and. mpi_enreg%me == 0) then 519 write(io_unit,10) "#" 520 write(io_unit,10) "# n > nbandv : exit loop!" 521 flush(io_unit) 522 end if 523 524 exit 525 end if 526 527 ! if the valence energy is smaller than the external frequency, 528 ! then there is no contribution 529 530 if (En_m_omega < zero .and. abs(En_m_omega) > energy_tolerance) then 531 ! careful close to zero! 532 if (debug .and. mpi_enreg%me == 0) then 533 write(io_unit,10) "# " 534 write(io_unit,10) "# eig(n) < omega : cycle!" 535 flush(io_unit) 536 end if 537 cycle 538 end if 539 540 541 ! if we are still here, there is a valence contribution 542 prefactor = -one 543 544 else if ( pole_is_conduction ) then 545 546 ! If the pole is conduction, get out when the conduction state is 547 ! larger than the frequency (careful close to zero!) 548 if ( En_m_omega > energy_tolerance ) then 549 if (debug .and. mpi_enreg%me == 0) then 550 write(io_unit,10) "#" 551 write(io_unit,10) "# eig(n) > omega : exit!" 552 flush(io_unit) 553 end if 554 exit 555 end if 556 557 ! If the pole is conduction, there is no contribution while 558 ! we are in the valence states 559 if ( n <= nbandv ) then 560 if (debug .and. mpi_enreg%me == 0) then 561 write(io_unit,10) "#" 562 write(io_unit,10) "# n <= nbandv : cycle!" 563 flush(io_unit) 564 end if 565 566 cycle 567 end if 568 569 ! if we are still here, there is a conduction contribution 570 prefactor = one 571 end if 572 573 574 !------------------------------------------------- 575 ! If we made it this far, we have a contribution! 576 !------------------------------------------------- 577 578 if (abs(En_m_omega) < energy_tolerance ) then 579 580 if (debug .and. mpi_enreg%me == 0) then 581 write(io_unit,10) "# " 582 write(io_unit,10) "# En - omega ~ 0: pole at the origin, multiply by 1/2!" 583 flush(io_unit) 584 end if 585 586 ! The factor of 1/2 accounts for the fact that 587 ! the pole is at the origin! 588 prefactor = 0.5_dp*prefactor 589 end if 590 591 592 593 number_of_seeds = number_of_degenerate_states(i_set) 594 595 ABI_ALLOCATE(seeds, (npw_k,number_of_seeds)) 596 597 call get_seeds(n, number_of_seeds, seeds) !Missing wrappers 598 599 call set_dielectric_function_frequency([En_m_omega,zero]) 600 if (debug .and. mpi_enreg%me == 0) then 601 write(io_unit,10) "# Compute pole contribution:" 602 write(io_unit,12) "# number of seeds = ",number_of_seeds 603 write(io_unit,16) "# || seeds || = ",sqrt(sum(abs(seeds(:,:))**2)) !Missing xmpi_sum 604 write(io_unit,16) "# eig(n)-omega = ",En_m_omega, " Ha" 605 write(io_unit,17) "# prefactor = ",prefactor 606 flush(io_unit) 607 end if 608 if(dtset%zcut > tol12) activate_inf_shift_poles = .true. 609 En_m_omega_2 = En_m_omega 610 pole_contribution = & 611 compute_pole_contribution(matrix_function_epsilon_k, & 612 number_of_seeds, kmax_poles, & 613 seeds,debug) 614 if(dtset%zcut > tol12) activate_inf_shift_poles = .false. 615 if (debug .and. mpi_enreg%me == 0) then 616 write(io_unit,16) "# pole contribution = ",prefactor*pole_contribution, " Ha" 617 flush(io_unit) 618 end if 619 620 621 compute_Poles = compute_Poles + prefactor*pole_contribution 622 ABI_DEALLOCATE(seeds) 623 end do 624 625 if (debug .and. mpi_enreg%me == 0) then 626 close(io_unit) 627 end if 628 629 630 10 format(A) 631 11 format(A,F8.4,A) 632 12 format(A,I5) 633 13 format(1000F16.8) 634 14 format(A,L10) 635 16 format(A,ES12.4,A) 636 17 format(A,F8.4) 637 638 end function compute_Poles
gwls_ComputePoles/generate_degeneracy_table_for_poles [ Functions ]
[ Top ] [ gwls_ComputePoles ] [ Functions ]
NAME
generate_degeneracy_table_for_poles
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_ComputeCorrelationEnergy
CHILDREN
block_lanczos_algorithm,diagonalize_lanczos_banded ritz_analysis_general,xmpi_sum
SOURCE
101 subroutine generate_degeneracy_table_for_poles(debug) 102 !---------------------------------------------------------------------- 103 ! This subroutine groups, once and for all, the indices of 104 ! degenerate eigenstates. This will be useful to compute the poles 105 ! 106 !---------------------------------------------------------------------- 107 108 !This section has been created automatically by the script Abilint (TD). 109 !Do not modify the following lines by hand. 110 #undef ABI_FUNC 111 #define ABI_FUNC 'generate_degeneracy_table_for_poles' 112 !End of the abilint section 113 114 implicit none 115 116 117 logical, intent(in) :: debug 118 119 real(dp) :: degeneracy_tolerance 120 121 integer :: nbands 122 integer :: n, i 123 integer :: i_set, j_deg 124 integer :: n_degeneracy 125 126 real(dp) :: energy 127 128 integer :: io_unit 129 character(128) :: filename 130 logical :: file_exists 131 132 ! ************************************************************************* 133 134 135 degeneracy_tolerance = 1.0D-8 136 !-------------------------------------------------------------------------------- 137 ! 138 ! First, find the largest degeneracy in the eigenvalue spectrum 139 ! 140 !-------------------------------------------------------------------------------- 141 142 nbands = size(eig) 143 144 ! initialize 145 largest_degeneracy = 1 146 number_of_denerate_sets = 1 147 148 n_degeneracy = 1 149 energy = eig(1) 150 151 !============================================================ 152 ! Notes for the code block below: 153 ! 154 ! The electronic eigenvalues can be grouped in degenerate 155 ! sets. The code below counts how many such sets there are, 156 ! and how many eigenvalues belong to each set. 157 ! 158 ! It is important to have a robust algorithm to do this 159 ! properly. In particular, the edge case where the LAST SET 160 ! is degenerate must be treated with care . 161 162 ! The algorithm.I'm looking at at the time of this writing 163 ! has a bug in it and cannot handle a degenerate last set... 164 ! Let's fix that! 165 !============================================================ 166 167 do n = 2, nbands 168 if (abs(eig(n) - energy) < degeneracy_tolerance) then 169 ! A degenerate state! add one 170 n_degeneracy = n_degeneracy + 1 171 else 172 ! We are no longer degenerate. Update 173 if ( n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy 174 175 n_degeneracy = 1 176 energy = eig(n) 177 number_of_denerate_sets = number_of_denerate_sets + 1 178 end if 179 180 ! If this is the last index, update the largest_degeneracy if necessary 181 if ( n == nbands .and. n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy 182 183 end do 184 185 !-------------------------------------------------------------------------------- 186 ! 187 ! Allocate the array which will contain the indices of the degenerate 188 ! states, and populate it. 189 !-------------------------------------------------------------------------------- 190 ABI_ALLOCATE(degeneracy_table, (number_of_denerate_sets,largest_degeneracy)) 191 ABI_ALLOCATE(number_of_degenerate_states, (number_of_denerate_sets)) 192 193 degeneracy_table(:,:) = 0 194 195 i_set = 1 196 j_deg = 1 197 198 ! initialize 199 energy = eig(1) 200 201 degeneracy_table(i_set,j_deg) = 1 202 number_of_degenerate_states(i_set) = 1 203 204 do n = 2, nbands 205 206 if (abs(eig(n) - energy) < degeneracy_tolerance) then 207 ! A degenerate state! add one 208 j_deg = j_deg + 1 209 210 else 211 ! We are no longer degenerate. Update 212 j_deg = 1 213 i_set = i_set+1 214 energy = eig(n) 215 216 end if 217 218 219 number_of_degenerate_states(i_set) = j_deg 220 degeneracy_table(i_set,j_deg) = n 221 222 end do 223 224 if (debug .and. mpi_enreg%me == 0) then 225 io_unit = get_unit() 226 filename = "degeneracy_table.log" 227 228 i = 0 229 inquire(file=filename,exist=file_exists) 230 do while (file_exists) 231 i = i+1 232 write (filename,'(A,I0,A)') "degeneracy_table_",i,".log" 233 inquire(file=filename,exist=file_exists) 234 end do 235 236 io_unit = get_unit() 237 238 open(io_unit,file=filename,status=files_status_new) 239 240 write(io_unit,10) " " 241 write(io_unit,10) "#===============================================================================================" 242 write(io_unit,10) "# Degeneracy table : tabulate the degenerate states " 243 write(io_unit,10) "# ------------------------------------------------------- " 244 write(io_unit,10) "# " 245 write(io_unit,14) "# number_of_denerate_sets = ", number_of_denerate_sets 246 write(io_unit,10) "# " 247 write(io_unit,14) "# largest_degeneracy = ", largest_degeneracy 248 write(io_unit,10) "# " 249 write(io_unit,10) "# Eigenvalues (Ha) " 250 write(io_unit,10) "#===============================================================================================" 251 write(io_unit,16) eig(:) 252 253 254 255 write(io_unit,10) "#===============================================================================================" 256 write(io_unit,10) "# i_set number of states States " 257 write(io_unit,10) "#===============================================================================================" 258 flush(io_unit) 259 260 do i_set = 1, number_of_denerate_sets 261 write(io_unit,12) i_set, number_of_degenerate_states(i_set), degeneracy_table(i_set,:) 262 end do 263 264 flush(io_unit) 265 266 close(io_unit) 267 end if 268 269 10 format(A) 270 12 format(I5,10X,I5,15X,1000I5) 271 14 format(A,I5) 272 16 format(1000F12.8,2X) 273 274 end subroutine generate_degeneracy_table_for_poles