TABLE OF CONTENTS


ABINIT/m_gwls_ComputePoles [ Modules ]

[ Top ] [ 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