TABLE OF CONTENTS


ABINIT/m_gwls_ComputePoles [ Modules ]

[ Top ] [ Modules ]

NAME

 m_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 m_gwls_ComputePoles
30 
31 ! local modules
32 use m_gwls_utility
33 use m_gwls_wf
34 use m_gwls_hamiltonian
35 use m_gwls_lineqsolver
36 use m_gwls_polarisability
37 use m_gwls_GWlanczos
38 use m_gwls_GenerateEpsilon
39 use m_gwls_GWanalyticPart
40 use m_gwls_TimingLog
41 use m_gwls_LanczosBasis
42 
43 
44 ! abinit modules
45 use defs_basis
46 use defs_datatypes
47 use defs_abitypes
48 use defs_wvltypes
49 use m_abicore
50 use m_xmpi
51 use m_pawang
52 use m_errors
53 
54 use m_io_tools,         only : get_unit
55 use m_paw_dmft,         only: paw_dmft_type
56 use m_ebands,           only : ebands_init, ebands_free
57 
58 use m_gaussian_quadrature, only: gaussian_quadrature_gegenbauer, gaussian_quadrature_legendre
59 
60 
61 implicit none
62 save
63 private
64 
65 integer  :: number_of_denerate_sets
66 integer  :: largest_degeneracy
67 
68 integer, allocatable :: degeneracy_table(:,:)
69 integer, allocatable :: number_of_degenerate_states(:)
70 
71 real(dp) :: En_m_omega_2
72 
73 public :: compute_Poles
74 public :: generate_degeneracy_table_for_poles
75 public :: clean_degeneracy_table_for_poles
76 
77 CONTAINS

m_gwls_ComputePoles/clean_degeneracy_table_for_poles [ Functions ]

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

296 subroutine clean_degeneracy_table_for_poles()
297 
298 
299 !This section has been created automatically by the script Abilint (TD).
300 !Do not modify the following lines by hand.
301 #undef ABI_FUNC
302 #define ABI_FUNC 'clean_degeneracy_table_for_poles'
303 !End of the abilint section
304 
305 implicit none
306 ! *************************************************************************
307 
308 if(allocated(degeneracy_table)) then
309   ABI_DEALLOCATE(degeneracy_table)
310 end if
311 if(allocated(number_of_degenerate_states)) then
312   ABI_DEALLOCATE(number_of_degenerate_states)
313 end if
314 
315 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

PARENTS

CHILDREN

SOURCE

659 function compute_pole_contribution(epsilon_matrix_function,nseeds,kmax,seeds,debug)
660 !----------------------------------------------------------------------
661 ! This routine computes the contribution to the  poles energy
662 ! coming from the states in the seeds.
663 !----------------------------------------------------------------------
664 
665 !This section has been created automatically by the script Abilint (TD).
666 !Do not modify the following lines by hand.
667 #undef ABI_FUNC
668 #define ABI_FUNC 'compute_pole_contribution'
669 !End of the abilint section
670 
671 implicit none
672 interface
673   subroutine epsilon_matrix_function(v_out,v_in,l)
674 
675   use defs_basis
676 
677   integer,     intent(in)  :: l
678   complex(dp), intent(out) :: v_out(l)
679   complex(dp), intent(in)  :: v_in(l)
680 
681   end subroutine epsilon_matrix_function
682 end interface
683 
684 real(dp) :: compute_pole_contribution
685 
686 integer,       intent(in) :: nseeds, kmax
687 complex(dpc),  intent(in) :: seeds(npw_k,nseeds)
688 logical,       intent(in) :: debug
689 
690 ! local variables
691 
692 integer  :: mpi_communicator
693 
694 complex(dpc),allocatable :: local_seeds(:,:)
695 
696 complex(dpc),allocatable :: Lbasis(:,:)  ! array containing the Lanczos basis
697 complex(dpc),allocatable :: alpha(:,:,:)
698 complex(dpc),allocatable :: beta (:,:,:)
699 real(dp),    allocatable :: epsilon_eigenvalues(:)
700 
701 real(dp):: matrix_elements
702 
703 complex(dpc) :: cmplx_value
704 integer :: l, s
705 integer :: ierr
706 
707 ! *************************************************************************
708 
709 ! compute the Lanczos basis
710 ABI_ALLOCATE(alpha,(nseeds,nseeds,kmax))
711 ABI_ALLOCATE(beta ,(nseeds,nseeds,kmax))
712 ABI_ALLOCATE(Lbasis,(npw_k,nseeds*kmax))
713 ABI_ALLOCATE(local_seeds,(npw_k,nseeds))
714 ABI_ALLOCATE(epsilon_eigenvalues, (nseeds*kmax))
715 
716 
717 mpi_communicator = mpi_enreg%comm_bandfft !Missing maybe something for easy access of LA and FFT comms?
718 
719 local_seeds(:,:) = seeds(:,:)
720 
721 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k,        &
722 &                                local_seeds,alpha,beta,Lbasis)
723 
724 write(std_out,*) "alpha:"
725 do l=1,kmax
726 do s=1,nseeds
727 write(std_out,*) alpha(:,s,l)
728 end do
729 write(std_out,*) " "
730 end do
731 
732 write(std_out,*) "beta:"
733 do l=1,kmax
734 do s=1,nseeds
735 write(std_out,*) beta(:,s,l)
736 end do
737 write(std_out,*) " "
738 end do
739 
740 ABI_DEALLOCATE(local_seeds)
741 
742 ! Diagonalize the epsilon matrix, which is banded
743 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug)
744 
745 if (debug) then
746   call ritz_analysis_general(mpi_communicator ,epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues)
747 end if
748 
749 compute_pole_contribution = zero
750 
751 do l = 1, nseeds*kmax
752 
753 matrix_elements = zero
754 
755 do s = 1, nseeds
756 
757 cmplx_value = complex_vector_product(seeds(:,s),Lbasis(:,l),npw_k)
758 
759 call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors working on FFT!
760 
761 matrix_elements = matrix_elements + abs(cmplx_value)**2
762 
763 
764 end do
765 
766 
767 compute_pole_contribution = compute_pole_contribution  + &
768 matrix_elements *(one/epsilon_eigenvalues(l)-one)
769 end do
770 
771 
772 
773 ABI_DEALLOCATE(alpha)
774 ABI_DEALLOCATE(beta)
775 ABI_DEALLOCATE(Lbasis)
776 ABI_DEALLOCATE(epsilon_eigenvalues)
777 
778 end function compute_pole_contribution
779 
780 end module m_gwls_ComputePoles

m_gwls_ComputePoles/compute_Poles [ Functions ]

[ Top ] [ m_gwls_ComputePoles ] [ Functions ]

NAME

  compute_Poles

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

337 function compute_Poles(external_omega,kmax_poles,debug)
338 !----------------------------------------------------------------------
339 ! This function extract the Pole contributions to the correlation
340 ! energy, as a function of the external frequency.
341 !
342 ! The algorithm builds a Lanczos chain for each contributing subspace;
343 ! the number of steps is controlled by kmax_poles. 
344 !
345 ! This function will take in explicit arguments, as it is simpler
346 ! to do this than to define global arrays.
347 !----------------------------------------------------------------------
348 
349 !This section has been created automatically by the script Abilint (TD).
350 !Do not modify the following lines by hand.
351 #undef ABI_FUNC
352 #define ABI_FUNC 'compute_Poles'
353 !End of the abilint section
354 
355 implicit none
356 real(dp) :: compute_Poles
357 
358 real(dp),     intent(in) :: external_omega
359 integer,      intent(in) :: kmax_poles
360 logical,      intent(in) :: debug
361 
362 real(dp) :: energy_tolerance
363 real(dp) :: pole_contribution
364 
365 
366 
367 integer :: number_of_seeds
368 integer :: i_set
369 integer :: n
370 
371 real(dp)     :: prefactor
372 real(dp)     :: En_m_omega
373 
374 logical      :: pole_is_valence
375 logical      :: pole_is_conduction
376 logical      :: pole_is_in_gap
377 
378 integer        :: io_unit, i
379 character(128) :: filename
380 logical        :: file_exists
381 
382 
383 
384 complex(dpc), allocatable :: seeds(:,:)
385 
386 ! *************************************************************************
387 
388 compute_Poles    =  zero
389 
390 energy_tolerance = 1.0D-8
391 
392 
393 if (debug .and. mpi_enreg%me == 0) then
394   io_unit = get_unit()
395 
396   i = 0
397 
398   file_exists = .true.
399   do while (file_exists)
400   i = i+1
401   write (filename,'(A,I0.4,A)') "ComputePoles_",i,".log"
402   inquire(file=filename,exist=file_exists)
403   end do
404 
405 
406   open(io_unit,file=filename,status=files_status_new)
407 
408   write(io_unit,10) " "
409   write(io_unit,10) "#==============================================================================================="
410   write(io_unit,10) "#                     ComputePoles: debug information for the pole computation                  "
411   write(io_unit,10) "#                     --------------------------------------------------------                  "
412   write(io_unit,10) "#                                                                                               "
413   write(io_unit,10) "# This file contains data describing the computation of the pole contribution to the            "
414   write(io_unit,10) "# correlation self energy.                                                                      "
415   write(io_unit,10) "#                                                                                               "
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) "#  parameters:                                                                                  "
421   write(io_unit,10) "#                                                                                               "
422   write(io_unit,11) "#          external_omega  = ",external_omega," Ha                                              "
423   write(io_unit,12) "#               kmax_poles = ",kmax_poles
424   write(io_unit,12) "#               nbandv     = ",nbandv
425   write(io_unit,10) "#                                                                                               "
426   write(io_unit,10) "#==============================================================================================="
427   write(io_unit,10) "#                                                                                               "
428   write(io_unit,10) "#   DFT Eigenvalues (Ha)                                                                        "
429   write(io_unit,10) "#==============================================================================================="
430   write(io_unit,13) eig(:)
431   flush(io_unit)
432 
433 end if
434 
435 
436 
437 !--------------------------------------------------------------------------------
438 !
439 ! Determine if external frequency corresponds to the valence or the conduction
440 ! manifold.
441 !
442 !--------------------------------------------------------------------------------
443 
444 compute_Poles = zero
445 
446 pole_is_conduction = .false.
447 pole_is_valence    = .false.
448 pole_is_in_gap     = .false.
449 
450 
451 
452 ! Careful here! there may be only nbandv states in memory; nbandv+1 causes segfaults!
453 if ( external_omega <= eig(nbandv)) then
454   pole_is_valence    = .true.
455 else if ( external_omega > eig(nbandv)) then
456   pole_is_conduction = .true.
457 else
458   pole_is_in_gap = .true.
459 end if
460 
461 
462 
463 
464 if (debug .and. mpi_enreg%me == 0 ) then
465   write(io_unit,10) "#===================================================================================================="
466   write(io_unit,10) "#                                                                                                    "
467   write(io_unit,10) "#  Determine where the external energy is:                                                           "
468   write(io_unit,10) "#                                                                                                    "
469   write(io_unit,14) "#             pole_is_valence    = ",pole_is_valence
470   write(io_unit,14) "#             pole_is_conduction = ",pole_is_conduction 
471   write(io_unit,14) "#             pole_is_in_gap     = ",pole_is_in_gap 
472   write(io_unit,10) "#                                                                                                    "
473   write(io_unit,10) "#===================================================================================================="
474   flush(io_unit)
475 
476 end if
477 
478 if ( pole_is_in_gap) return
479 
480 !--------------------------------------------------------------------------------
481 !
482 ! Loop on all degenerate sets
483 !
484 !--------------------------------------------------------------------------------
485 
486 if (debug .and. mpi_enreg%me == 0 ) then
487   write(io_unit,10) "#===================================================================================================="
488   write(io_unit,10) "#                                                                                                    "
489   write(io_unit,10) "#  Iterating over all degenerate sets of eigenvalues:                                                "
490   write(io_unit,10) "#                                                                                                    "
491   write(io_unit,10) "#===================================================================================================="
492   flush(io_unit)
493 
494 end if
495 
496 
497 do i_set =1, number_of_denerate_sets
498 
499 n = degeneracy_table(i_set,1)
500 
501 En_m_omega = eig(n)-external_omega
502 
503 if (debug .and. mpi_enreg%me == 0) then
504   write(io_unit,12) "# i_set = ", i_set
505   write(io_unit,12) "#                           n = ",n
506   write(io_unit,16) "#                eig(n)-omega = ",En_m_omega," Ha"
507   flush(io_unit)
508 end if
509 
510 !------------------------------------------
511 ! Test if we need to exit the loop
512 !------------------------------------------
513 if (pole_is_valence ) then 
514 
515   ! If the pole is valence, get out when we enter conduction states
516   if (  n > nbandv  ) then
517     if (debug.and. mpi_enreg%me == 0) then
518       write(io_unit,10) "#"
519       write(io_unit,10) "#                n > nbandv : exit loop!"
520       flush(io_unit)
521     end if
522 
523     exit 
524   end if
525 
526   ! if the valence energy is smaller than the external frequency,
527   ! then there is no contribution
528 
529   if (En_m_omega < zero  .and. abs(En_m_omega) > energy_tolerance) then 
530     ! careful close to zero!
531     if (debug .and. mpi_enreg%me == 0) then
532       write(io_unit,10) "# "
533       write(io_unit,10) "#                 eig(n) < omega : cycle!"
534       flush(io_unit)
535     end if
536     cycle
537   end if
538 
539 
540   ! if we are still here, there is a valence contribution
541   prefactor = -one
542 
543 else if ( pole_is_conduction ) then 
544 
545   ! If the pole is conduction, get out when the conduction state is
546   ! larger than the frequency (careful close to zero!)
547   if ( En_m_omega > energy_tolerance ) then
548     if (debug .and. mpi_enreg%me == 0) then
549       write(io_unit,10) "#"
550       write(io_unit,10) "#                eig(n) > omega : exit!"
551       flush(io_unit)
552     end if
553     exit 
554   end if
555 
556   ! If the pole is conduction, there is no contribution while
557   ! we are in the valence states
558   if (  n <= nbandv  ) then
559     if (debug .and. mpi_enreg%me == 0) then
560       write(io_unit,10) "#"
561       write(io_unit,10) "#                n <= nbandv : cycle!"
562       flush(io_unit)
563     end if
564 
565     cycle
566   end if
567 
568   ! if we are still here, there is a conduction contribution
569   prefactor = one
570 end if
571 
572 
573 !-------------------------------------------------
574 ! If we made it this far, we have a contribution!
575 !-------------------------------------------------
576 
577 if (abs(En_m_omega) < energy_tolerance ) then
578 
579   if (debug .and. mpi_enreg%me == 0) then
580     write(io_unit,10) "# "
581     write(io_unit,10) "#                En - omega ~ 0: pole at the origin, multiply by 1/2!"
582     flush(io_unit)
583   end if
584 
585   ! The factor of 1/2 accounts for the fact that
586   ! the pole is at the origin!
587   prefactor = 0.5_dp*prefactor  
588 end if
589 
590 
591 
592 number_of_seeds = number_of_degenerate_states(i_set)
593 
594 ABI_ALLOCATE(seeds, (npw_k,number_of_seeds))
595 
596 call get_seeds(n, number_of_seeds, seeds) !Missing wrappers
597 
598 call set_dielectric_function_frequency([En_m_omega,zero])
599 if (debug .and. mpi_enreg%me == 0) then
600   write(io_unit,10) "#                Compute pole contribution:"
601   write(io_unit,12) "#                        number of seeds = ",number_of_seeds 
602   write(io_unit,16) "#                        ||   seeds   || = ",sqrt(sum(abs(seeds(:,:))**2))  !Missing xmpi_sum
603   write(io_unit,16) "#                        eig(n)-omega    = ",En_m_omega, " Ha"
604   write(io_unit,17) "#                        prefactor       = ",prefactor
605   flush(io_unit)
606 end if
607 if(dtset%zcut > tol12) activate_inf_shift_poles = .true.
608 En_m_omega_2 = En_m_omega
609 pole_contribution =                                             &
610 compute_pole_contribution(matrix_function_epsilon_k,    &
611 number_of_seeds, kmax_poles,    &
612 seeds,debug)
613 if(dtset%zcut > tol12) activate_inf_shift_poles = .false.
614 if (debug .and. mpi_enreg%me == 0) then
615   write(io_unit,16) "#                      pole contribution = ",prefactor*pole_contribution, " Ha"
616   flush(io_unit)
617 end if
618 
619 
620 compute_Poles = compute_Poles + prefactor*pole_contribution 
621 ABI_DEALLOCATE(seeds)
622 end do
623 
624 if (debug .and. mpi_enreg%me == 0) then
625   close(io_unit)
626 end if
627 
628 
629 10 format(A)
630 11 format(A,F8.4,A)
631 12 format(A,I5)
632 13 format(1000F16.8)
633 14 format(A,L10)
634 16 format(A,ES12.4,A)
635 17 format(A,F8.4)
636 
637 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

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      block_lanczos_algorithm,diagonalize_lanczos_banded
      ritz_analysis_general,xmpi_sum

SOURCE

100 subroutine generate_degeneracy_table_for_poles(debug)
101 !----------------------------------------------------------------------
102 ! This subroutine groups, once and for all, the indices of 
103 ! degenerate eigenstates. This will be useful to compute the poles
104 !
105 !----------------------------------------------------------------------
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'generate_degeneracy_table_for_poles'
111 !End of the abilint section
112 
113 implicit none
114 
115 
116 logical, intent(in) :: debug
117 
118 real(dp) :: degeneracy_tolerance
119 
120 integer  :: nbands
121 integer  :: n, i
122 integer  :: i_set, j_deg
123 integer  :: n_degeneracy
124 
125 real(dp) :: energy
126 
127 integer        :: io_unit
128 character(128) :: filename
129 logical        :: file_exists
130 
131 ! *************************************************************************
132 
133 
134 degeneracy_tolerance = 1.0D-8
135 !--------------------------------------------------------------------------------
136 !
137 ! First, find the largest degeneracy in the eigenvalue spectrum
138 !
139 !--------------------------------------------------------------------------------
140 
141 nbands = size(eig)
142 
143 ! initialize
144 largest_degeneracy      = 1 
145 number_of_denerate_sets = 1
146 
147 n_degeneracy = 1
148 energy       = eig(1)
149 
150 !============================================================
151 ! Notes for the code block below: 
152 !
153 ! The electronic eigenvalues can be grouped in degenerate
154 ! sets. The code below counts how many such sets there are,
155 ! and how many eigenvalues belong to each set.
156 !
157 ! It is important to have a robust algorithm to do this 
158 ! properly. In particular, the edge case where the LAST SET
159 ! is degenerate must be treated with care .
160 
161 ! The algorithm.I'm looking at at the time of this writing 
162 ! has a bug in it and cannot handle a degenerate last set...
163 ! Let's fix that!
164 !============================================================
165 
166 do n = 2, nbands
167   if (abs(eig(n) - energy) < degeneracy_tolerance) then
168     ! A degenerate state! add one
169     n_degeneracy  = n_degeneracy + 1
170   else
171     ! We are no longer degenerate. Update
172     if ( n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy
173 
174     n_degeneracy = 1
175     energy       = eig(n)
176     number_of_denerate_sets = number_of_denerate_sets + 1
177   end if
178 
179   ! If this is the last index, update the largest_degeneracy if necessary
180   if ( n == nbands .and. n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy
181 
182 end do
183 
184 !--------------------------------------------------------------------------------
185 !
186 ! Allocate the array which will contain the indices of the degenerate
187 ! states, and populate it.
188 !--------------------------------------------------------------------------------
189 ABI_ALLOCATE(degeneracy_table,            (number_of_denerate_sets,largest_degeneracy))
190 ABI_ALLOCATE(number_of_degenerate_states, (number_of_denerate_sets))
191 
192 degeneracy_table(:,:) = 0
193 
194 i_set = 1
195 j_deg = 1
196 
197 ! initialize
198 energy = eig(1)
199 
200 degeneracy_table(i_set,j_deg)       = 1
201 number_of_degenerate_states(i_set)  = 1
202 
203 do n = 2, nbands
204 
205   if (abs(eig(n) - energy) < degeneracy_tolerance) then
206     ! A degenerate state! add one
207     j_deg = j_deg + 1
208 
209   else
210     ! We are no longer degenerate. Update
211     j_deg = 1
212     i_set = i_set+1
213     energy = eig(n)
214 
215   end if
216 
217 
218   number_of_degenerate_states(i_set) = j_deg
219   degeneracy_table(i_set,j_deg)      = n
220 
221 end do
222 
223 if (debug .and. mpi_enreg%me == 0) then
224   io_unit  = get_unit()
225   filename = "degeneracy_table.log"
226 
227   i = 0
228   inquire(file=filename,exist=file_exists)
229   do while (file_exists)
230   i = i+1
231   write (filename,'(A,I0,A)') "degeneracy_table_",i,".log"
232   inquire(file=filename,exist=file_exists)
233   end do
234 
235   io_unit = get_unit()
236 
237   open(io_unit,file=filename,status=files_status_new)
238 
239   write(io_unit,10) " "
240   write(io_unit,10) "#==============================================================================================="
241   write(io_unit,10) "#                     Degeneracy table : tabulate the degenerate states                         "
242   write(io_unit,10) "#                     -------------------------------------------------------                   "
243   write(io_unit,10) "#                                                                                               "
244   write(io_unit,14) "#     number_of_denerate_sets = ", number_of_denerate_sets
245   write(io_unit,10) "#                                                                                               "
246   write(io_unit,14) "#      largest_degeneracy     = ", largest_degeneracy
247   write(io_unit,10) "#                                                                                               "
248   write(io_unit,10) "# Eigenvalues (Ha)                                                                              "
249   write(io_unit,10) "#==============================================================================================="
250   write(io_unit,16) eig(:)
251 
252 
253 
254   write(io_unit,10) "#==============================================================================================="
255   write(io_unit,10) "#  i_set    number of states           States                                                   "
256   write(io_unit,10) "#==============================================================================================="
257   flush(io_unit)
258 
259   do i_set = 1, number_of_denerate_sets
260   write(io_unit,12) i_set, number_of_degenerate_states(i_set), degeneracy_table(i_set,:)
261   end do 
262 
263   flush(io_unit)
264 
265   close(io_unit)
266 end if 
267 
268 10 format(A)
269 12 format(I5,10X,I5,15X,1000I5)
270 14 format(A,I5)
271 16 format(1000F12.8,2X)
272 
273 end subroutine generate_degeneracy_table_for_poles