TABLE OF CONTENTS


ABINIT/m_gwls_GenerateEpsilon [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_GenerateEpsilon

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 module m_gwls_GenerateEpsilon
29 !----------------------------------------------------------------------------------------------------
30 ! This module contains routines to compute and store the dielectric matrix, which plays a
31 ! central role in the self energy computations. In particular, global arrays are used to store
32 ! the static dielectric matrix.
33 !----------------------------------------------------------------------------------------------------
34 ! local modules
35 use m_gwls_utility
36 use m_gwls_wf
37 use m_gwls_hamiltonian
38 use m_gwls_lineqsolver
39 use m_gwls_TimingLog
40 use m_gwls_polarisability
41 use m_gwls_model_polarisability
42 use m_gwls_GWlanczos
43 ! Abinit modules
44 use m_abicore
45 use defs_basis
46 use defs_datatypes
47 use defs_abitypes
48 
49 use m_io_tools,    only : get_unit
50 
51 implicit none
52 save
53 private

m_gwls_GenerateEpsilon/driver_generate_dielectric_matrix [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  driver_generate_dielectric_matrix

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      cpu_time,diagonalize_lanczos_banded
      driver_invert_positive_definite_hermitian_matrix
      generateprintdielectriceigenvalues
      matrix_function_epsilon_model_operator
      set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm
      zheevd

SOURCE

 97 subroutine driver_generate_dielectric_matrix(epsilon_matrix_function,nseeds,kmax,&
 98 epsilon_eigenvalues,Lbasis,debug)
 99 !----------------------------------------------------------------------
100 ! This routine computes the Lanczos approximate representation of the
101 ! implicit dielectic operator and then diagonalizes the banded
102 ! Lanczos matrix.
103 !----------------------------------------------------------------------
104 
105 !This section has been created automatically by the script Abilint (TD).
106 !Do not modify the following lines by hand.
107 #undef ABI_FUNC
108 #define ABI_FUNC 'driver_generate_dielectric_matrix'
109 !End of the abilint section
110 
111 implicit none
112 interface
113   subroutine epsilon_matrix_function(v_out,v_in,l)
114   use defs_basis
115 
116   integer,     intent(in)  :: l
117   complex(dp), intent(out) :: v_out(l)
118   complex(dp), intent(in)  :: v_in(l)
119 
120   end subroutine epsilon_matrix_function
121 end interface
122 
123 integer,       intent(in) :: nseeds, kmax
124 logical,       intent(in) :: debug
125 
126 real   (dp),  intent(out) :: epsilon_eigenvalues(nseeds*kmax)
127 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax)  ! array containing the Lanczos basis
128 
129 
130 ! local variables
131 
132 complex(dpc), allocatable :: seeds(:,:)
133 
134 complex(dpc),allocatable :: alpha(:,:,:)
135 complex(dpc),allocatable :: beta (:,:,:)
136 
137 integer :: mpi_communicator
138 
139 ! *************************************************************************
140 
141 
142 ! The epsilon operator will act in LA mode.
143 mpi_communicator = mpi_enreg%comm_bandfft
144 
145 
146 !Create seeds
147 ABI_ALLOCATE(seeds,(npw_k,nseeds))
148 call get_seeds(first_seed, nseeds, seeds)
149 
150 ! compute the Lanczos basis
151 ABI_ALLOCATE(alpha,(nseeds,nseeds,kmax))
152 ABI_ALLOCATE(beta ,(nseeds,nseeds,kmax))
153 
154 call block_lanczos_algorithm(mpi_communicator,epsilon_matrix_function,kmax,nseeds,npw_k,        &
155 seeds,alpha,beta,Lbasis)
156 
157 ! Diagonalize the epsilon matrix, which is banded
158 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug)
159 
160 if (debug) then
161   call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues)
162 end if
163 
164 ABI_DEALLOCATE(seeds)
165 ABI_DEALLOCATE(alpha)
166 ABI_DEALLOCATE(beta)
167 
168 end subroutine driver_generate_dielectric_matrix

m_gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  Driver_GeneratePrintDielectricEigenvalues

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_sternheimer

CHILDREN

      cpu_time,diagonalize_lanczos_banded
      driver_invert_positive_definite_hermitian_matrix
      generateprintdielectriceigenvalues
      matrix_function_epsilon_model_operator
      set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm
      zheevd

SOURCE

426 subroutine Driver_GeneratePrintDielectricEigenvalues(dtset)
427 !----------------------------------------------------------------------
428 ! Compute the eigenvalues of the various dielectric operators
429 !----------------------------------------------------------------------
430 
431 !This section has been created automatically by the script Abilint (TD).
432 !Do not modify the following lines by hand.
433 #undef ABI_FUNC
434 #define ABI_FUNC 'Driver_GeneratePrintDielectricEigenvalues'
435 !End of the abilint section
436 
437 type(dataset_type),intent(in) :: dtset
438 
439 integer  ::  kmax_exact, kmax_model, kmax
440 real(dp) :: second_model_parameter
441 
442 
443 integer  ::  lm, k, lmax, l1, l2
444 integer  ::  io_unit
445 integer  ::  io_unit2
446 
447 real(dp) :: time1, time2, time
448 ! local variables
449 character(128)  :: output_filename
450 character(256)  :: timing_string
451 
452 
453 complex(dpc), allocatable :: Lbasis_exact(:,:)
454 complex(dpc), allocatable :: Lbasis_model(:,:)
455 
456 complex(dpc), allocatable :: sub_Lbasis_exact(:,:)
457 complex(dpc), allocatable :: sub_Lbasis_model(:,:)
458 
459 complex(dpc), allocatable :: dummy(:,:)
460 complex(dpc), allocatable :: dummy2(:,:)
461 complex(dpc), allocatable :: dummy3(:,:)
462 
463 complex(dpc), allocatable :: alpha_exact(:,:,:)
464 complex(dpc), allocatable :: beta_exact (:,:,:)
465 
466 complex(dpc), allocatable :: alpha_model(:,:,:)
467 complex(dpc), allocatable :: beta_model (:,:,:)
468 
469 real(dp), allocatable :: eig_exact(:)
470 real(dp), allocatable :: eig_model(:)
471 
472 complex(dpc), allocatable :: model_epsilon_matrix(:,:)
473 complex(dpc), allocatable :: vector(:)
474 
475 real(dpc) :: tr_eps_1, tr_eps_2, tr_eps_3
476 
477 
478 integer   ::  lwork, lrwork, liwork, info
479 complex(dpc), allocatable :: work(:)
480 real(dp)    , allocatable :: rwork(:)
481 integer     , allocatable :: iwork(:)
482 
483 integer        :: debug_unit
484 character(50)  :: debug_filename
485 
486 ! *************************************************************************
487 
488 
489 
490 
491 kmax_exact   = dtset%gwls_stern_kmax
492 kmax_model   = dtset%gwls_kmax_complement
493 
494 !second_model_parameter  = dtset%gwls_second_model_parameter
495 second_model_parameter  = zero
496 
497 
498 ! global stuff
499 nseeds       = dtset%gwls_nseeds
500 first_seed   = dtset%gwls_first_seed
501 e            = dtset%gwls_band_index
502 
503 
504 
505 ABI_ALLOCATE(Lbasis_exact,(npw_k,kmax_exact*nseeds))
506 ABI_ALLOCATE(Lbasis_model,(npw_k,kmax_model*nseeds))
507 
508 ABI_ALLOCATE(alpha_exact, (nseeds,nseeds,kmax_exact))
509 ABI_ALLOCATE(beta_exact , (nseeds,nseeds,kmax_exact))
510 ABI_ALLOCATE(alpha_model, (nseeds,nseeds,kmax_model))
511 ABI_ALLOCATE(beta_model , (nseeds,nseeds,kmax_model))
512 
513 
514 ! set omega=0 for exact dielectric operator
515 call set_dielectric_function_frequency([zero,zero])
516 
517 
518 
519 call cpu_time(time1)
520 output_filename = 'EIGENVALUES_EXACT.dat'
521 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_k, kmax_exact, &
522 output_filename, Lbasis_exact, alpha_exact, beta_exact)
523 
524 
525 
526 call cpu_time(time2)
527 time = time2-time1
528 write(timing_string,'(A)')  "Time to compute the EXACT Static Dielectric Matrix  :   "
529 call write_timing_log(timing_string,time)
530 
531 
532 
533 call cpu_time(time1)
534 call setup_Pk_model(zero,second_model_parameter)
535 output_filename = 'EIGENVALUES_MODEL.dat'
536 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_model_operator, kmax_model, output_filename,&
537 Lbasis_model, alpha_model, beta_model)
538 call cpu_time(time2)
539 time = time2-time1
540 write(timing_string,'(A)')  "Time to compute the MODEL Static Dielectric Matrix  :   "
541 call write_timing_log(timing_string,time)
542 
543 
544 call cpu_time(time1)
545 if (kmax_exact <= kmax_model) then
546   kmax  = kmax_exact
547 else
548   kmax  = kmax_model
549 end if
550 lmax = nseeds*kmax
551 
552 ! Build model operator matrix elements in the exact basis
553 ABI_ALLOCATE(model_epsilon_matrix, (lmax,lmax))
554 ABI_ALLOCATE(vector, (npw_k))
555 
556 model_epsilon_matrix = cmplx_0
557 
558 do l2 =1 , lmax
559 call matrix_function_epsilon_model_operator(vector ,Lbasis_exact(:,l2),npw_k)
560 
561 
562 do l1 =1, lmax
563 
564 model_epsilon_matrix(l1, l2) = complex_vector_product(Lbasis_exact(:,l1),vector,npw_k)
565 
566 end do
567 
568 
569 end do
570 
571 ABI_DEALLOCATE(vector)
572 
573 call cpu_time(time2)
574 time = time2-time1
575 write(timing_string,'(A)')  "Compute MODEL matrix elements in EXACT basis        :   "
576 call write_timing_log(timing_string,time)
577 
578 
579 
580 
581 ! Compare the traces 
582 
583 
584 io_unit = get_unit()
585 open(file='DIELECTRIC_TRACE.dat',status=files_status_new,unit=io_unit)
586 write(io_unit,10) '#----------------------------------------------------------------------------'
587 write(io_unit,10) '#                                                                            '
588 write(io_unit,10) '#                       Partial traces                                       '
589 write(io_unit,10) '#           ==========================================                       '
590 write(io_unit,10) '#                                                                            '
591 write(io_unit,10) '#  Tabulate the trace of various operators, as function of the number        '
592 write(io_unit,10) '#  of lanczos steps performed.                                               '
593 write(io_unit,10) '#                                                                            '
594 write(io_unit,10) '#                                                                            '
595 write(io_unit,10) '#  NOTES:                                                                    '
596 write(io_unit,10) '#         Tr[1-eps^{-1}] is evaluated in the Lanczos basis of eps            '
597 write(io_unit,10) '#         Tr[1-eps_m^{-1}] is evaluated in the Lanczos basis of eps_m        '
598 write(io_unit,10) '#         Tr[eps_m^{-1}-eps^{-1}] is evaluated in the Lanczos basis of eps   '
599 write(io_unit,10) '#                                                                            '
600 write(io_unit,10) '#                                                                            '
601 write(io_unit,10) '#  k       Tr[1-eps^{-1}]        Tr[1-eps_m^{-1}]     Tr[eps_m^{-1}-eps^{-1}]'
602 write(io_unit,10) '#----------------------------------------------------------------------------'
603 flush(io_unit)
604 
605 
606 io_unit2 = get_unit()
607 open(file='RPA_ENERGY.dat',status=files_status_new,unit=io_unit2)
608 write(io_unit2,10) '#----------------------------------------------------------------------------'
609 write(io_unit2,10) '#                                                                            '
610 write(io_unit2,10) '#                       RPA TOTAL ENERGY                                     '
611 write(io_unit2,10) '#           ==========================================                       '
612 write(io_unit2,10) '#                                                                            '
613 write(io_unit2,10) '#  It can be shown that the correlation energy, within the RPA, is given     '
614 write(io_unit2,10) '#  by:                                                                       '
615 write(io_unit2,10) '#       E_c = int_0^{infty} dw/(2pi) Tr[ ln(eps{iw)}+1-eps(iw)]              '
616 write(io_unit2,10) '#                                                                            '
617 write(io_unit2,10) '#  As a gauge of what can be expected as far as convergence is concerned,    '
618 write(io_unit2,10) '#  the following will be printed.                                            '
619 write(io_unit2,10) '#                                                                            '
620 write(io_unit2,10) '#  I_1 = Tr[  ln(eps)  + 1 - eps   ]                                         '
621 write(io_unit2,10) '#  I_2 = Tr[ ln(eps_m) + 1 - eps_m ]                                         '
622 write(io_unit2,10) '#  I_3 = Tr[ ln(eps^{-1}_m . eps ) + eps_m - eps ]                           '
623 write(io_unit2,10) '#                                                                            '
624 write(io_unit2,10) '#                                                                            '
625 write(io_unit2,10) '#                                                                            '
626 write(io_unit2,10) '#  k       I_1                   I_2                  I_3                    '
627 write(io_unit2,10) '#----------------------------------------------------------------------------'
628 flush(io_unit2)
629 
630 
631 ! Iterate every 10 values of k max, or else the linear algebra gets too expensive...
632 do k = 4, kmax, 4
633 
634 ABI_ALLOCATE(sub_Lbasis_exact,(npw_k,k*nseeds))
635 ABI_ALLOCATE(sub_Lbasis_model,(npw_k,k*nseeds))
636 
637 
638 ABI_ALLOCATE(eig_exact,(k*nseeds))
639 ABI_ALLOCATE(eig_model,(k*nseeds))
640 ABI_ALLOCATE(dummy,(k*nseeds,k*nseeds))
641 ABI_ALLOCATE(dummy2,(k*nseeds,k*nseeds))
642 
643 sub_Lbasis_exact(:,:) = Lbasis_exact(:,1:k*nseeds)
644 sub_Lbasis_model(:,:) = Lbasis_model(:,1:k*nseeds)
645 
646 ! Diagonalize the epsilon matrix, which is banded
647 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_exact(:,:,1:k),beta_exact(:,:,1:k),sub_Lbasis_exact,eig_exact,.false.)
648 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_model(:,:,1:k),beta_model(:,:,1:k),sub_Lbasis_model,eig_model,.false.)
649 
650 tr_eps_1 = sum(one-one/eig_exact(:))
651 tr_eps_2 = sum(one-one/eig_model(:))
652 
653 tr_eps_3 = -sum(one/eig_exact(:))
654 
655 dummy(:,:) = model_epsilon_matrix(1:k*nseeds, 1:k*nseeds)
656 
657 call driver_invert_positive_definite_hermitian_matrix(dummy,k*nseeds)
658 
659 do lm = 1, k*nseeds
660 
661 tr_eps_3 = tr_eps_3 +  dble(dummy(lm,lm))
662 end do
663 
664 
665 write(io_unit,20) k, tr_eps_1, tr_eps_2, tr_eps_3
666 flush(io_unit)
667 
668 
669 
670 tr_eps_1 = sum(log(eig_exact(:))+one-eig_exact(:))
671 tr_eps_2 = sum(log(eig_model(:))+one-eig_model(:))
672 
673 dummy2(:,:) = zero
674 tr_eps_3    = zero
675 do lm = 1, k*nseeds
676 dummy2(lm,lm) = eig_exact(lm)
677 tr_eps_3 = tr_eps_3 + dble(model_epsilon_matrix(lm,lm)) - dble(eig_exact(lm))
678 end do
679 
680 ABI_ALLOCATE(dummy3,(k*nseeds,k*nseeds))
681 call ZGEMM(      'N',   & ! Hermitian conjugate the first array
682 'N',   & ! Leave second array as is
683 k*nseeds,   & ! the number of rows of the  matrix op( A )
684 k*nseeds,   & ! the number of columns of the  matrix op( B )
685 k*nseeds,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
686 cmplx_1,   & ! alpha constant
687 dummy2,   & ! matrix A
688 k*nseeds,   & ! LDA
689 dummy,   & ! matrix B
690 k*nseeds,   & ! LDB
691 cmplx_0,   & ! beta constant
692 dummy3,   & ! matrix C
693 k*nseeds)     ! LDC
694 
695 dummy2(:,:) = dummy3(:,:)
696 ABI_DEALLOCATE(dummy3)
697 
698 ! find eigenvalues
699 !call heevd(dummy2, eig_exact)
700 
701 lwork  = k*nseeds+1
702 lrwork = k*nseeds
703 liwork = 1
704 
705 ABI_ALLOCATE(work,(lwork))
706 ABI_ALLOCATE(rwork,(lrwork))
707 ABI_ALLOCATE(iwork,(liwork))
708 
709 call zheevd('N', 'U',k*nseeds, dummy2, k*nseeds, eig_exact, work, lwork, rwork, lrwork, iwork, liwork, info)
710 if ( info /= 0) then        
711   debug_unit = get_unit()
712   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
713 
714   open(debug_unit,file=trim(debug_filename),status='unknown')
715 
716   write(debug_unit,'(A)')      '*************************************************************************************'
717   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZHEEVD(1), gwls_GenerateEpsilon'
718   write(debug_unit,'(A)')      '*************************************************************************************'
719 
720   close(debug_unit)
721 
722 end if
723 
724 
725 
726 
727 
728 tr_eps_3 =  tr_eps_3 + sum(log(eig_exact))
729 
730 write(io_unit2,20) k, tr_eps_1, tr_eps_2, tr_eps_3
731 flush(io_unit2)
732 
733 
734 
735 ABI_DEALLOCATE(work)
736 ABI_DEALLOCATE(rwork)
737 ABI_DEALLOCATE(iwork)
738 
739 
740 
741 ABI_DEALLOCATE(sub_Lbasis_exact)
742 ABI_DEALLOCATE(sub_Lbasis_model)
743 
744 ABI_DEALLOCATE(eig_exact)
745 ABI_DEALLOCATE(eig_model)
746 ABI_DEALLOCATE(dummy)
747 ABI_DEALLOCATE(dummy2)
748 end do
749 
750 close(io_unit)
751 close(io_unit2)
752 
753 call cpu_time(time2)
754 time = time2-time1
755 write(timing_string,'(A)')  "Time to compute the TRACES of the Dielectric Matrices:   "
756 call write_timing_log(timing_string,time)
757 
758 ABI_DEALLOCATE(Lbasis_exact)
759 ABI_DEALLOCATE(Lbasis_model)
760 ABI_DEALLOCATE(alpha_exact)
761 ABI_DEALLOCATE(beta_exact )
762 ABI_DEALLOCATE(alpha_model)
763 ABI_DEALLOCATE(beta_model )
764 ABI_DEALLOCATE(model_epsilon_matrix)
765 
766 
767 
768 10 format(A)
769 20 format(I5,3ES24.16)
770 
771 end subroutine Driver_GeneratePrintDielectricEigenvalues

m_gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  GeneratePrintDielectricEigenvalues

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gwls_GenerateEpsilon

CHILDREN

      cpu_time,diagonalize_lanczos_banded
      driver_invert_positive_definite_hermitian_matrix
      generateprintdielectriceigenvalues
      matrix_function_epsilon_model_operator
      set_dielectric_function_frequency,setup_pk_model,write_timing_log,zgemm
      zheevd

SOURCE

195 subroutine GeneratePrintDielectricEigenvalues(epsilon_matrix_function,kmax,output_filename,Lbasis,alpha,beta)
196 !----------------------------------------------------------------------
197 ! This routine computes the Lanczos approximate representation of the
198 ! implicit dielectic operator and then diagonalizes the banded
199 ! Lanczos matrix.
200 !----------------------------------------------------------------------
201 
202 !This section has been created automatically by the script Abilint (TD).
203 !Do not modify the following lines by hand.
204 #undef ABI_FUNC
205 #define ABI_FUNC 'GeneratePrintDielectricEigenvalues'
206 !End of the abilint section
207 
208 implicit none
209 interface
210   subroutine epsilon_matrix_function(v_out,v_in,l)
211   use defs_basis
212 
213   integer,     intent(in)  :: l
214   complex(dp), intent(out) :: v_out(l)
215   complex(dp), intent(in)  :: v_in(l)
216 
217   end subroutine epsilon_matrix_function
218 end interface
219 
220 integer,       intent(in) :: kmax
221 
222 character(*),  intent(in) :: output_filename
223 
224 
225 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax)  
226 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax)
227 complex(dpc), intent(out) :: beta (nseeds,nseeds,kmax)
228 
229 
230 ! local variables
231 
232 
233 complex(dpc),allocatable :: seeds(:,:)
234 
235 complex(dpc),allocatable :: Lbasis_diag(:,:)  
236 
237 
238 real(dp),    allocatable :: psik(:,:)
239 real(dp),    allocatable :: psir(:,:,:,:)
240 
241 real(dp),    allocatable :: epsilon_eigenvalues(:)
242 
243 
244 integer :: mpi_communicator
245 integer :: io_unit
246 integer :: lmax
247 integer :: l
248 integer :: ir1, ir2, ir3
249 integer :: n1, n2, n3 
250 
251 real(dp) :: R, G
252 real(dp) :: sigma_R, sigma_G
253 real(dp) :: x, y, z
254 
255 real(dp),allocatable :: G_array(:)
256 real(dp),allocatable :: R_array(:,:,:)
257 
258 logical :: debug
259 
260 ! *************************************************************************
261 
262 
263 debug = .false.
264 lmax = kmax*nseeds
265 mpi_communicator = mpi_enreg%comm_bandfft
266 !Create seeds
267 ABI_ALLOCATE(seeds,(npw_k,nseeds))
268 call get_seeds(first_seed, nseeds, seeds)
269 
270 ! compute the Lanczos basis
271 ABI_ALLOCATE(Lbasis_diag,(npw_k,lmax))
272 ABI_ALLOCATE(epsilon_eigenvalues,(lmax))
273 
274 ABI_ALLOCATE(psik,(2,npw_k))
275 ABI_ALLOCATE(psir,(2,n4,n5,n6))
276 ABI_ALLOCATE(G_array,(npw_k))
277 ABI_ALLOCATE(R_array,(n4,n5,n6))
278 
279 psir = zero
280 R_array = zero
281 
282 n1 = n4-1
283 n2 = n5-1
284 n3 = n6
285 
286 ! Generate the Lanczos basis and banded eigenvalue representation
287 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k,        seeds,alpha,beta,Lbasis)
288 
289 Lbasis_diag = Lbasis
290 
291 ! Diagonalize the epsilon matrix, which is banded
292 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis_diag,epsilon_eigenvalues,debug)
293 
294 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,lmax,npw_k,Lbasis_diag,epsilon_eigenvalues)
295 
296 io_unit = get_unit()
297 open(file=output_filename,status=files_status_new,unit=io_unit)
298 write(io_unit,10) '#----------------------------------------------------------------------------'
299 write(io_unit,10) '#                                                                            '
300 write(io_unit,10) '#                       Partial eigenvalues                                  '
301 write(io_unit,10) '#           ==========================================                       '
302 write(io_unit,10) '#                                                                            '
303 write(io_unit,10) '#  Tabulate the eigenvalues of the dielectic matrix, as well as some         '
304 write(io_unit,10) '#  information regarding the eigenstates.                                    '
305 write(io_unit,10) '#                                                                            '
306 write(io_unit,10) '#                                                                            '
307 write(io_unit,10) '#  definitions:                                                              '
308 write(io_unit,10) '#                l      index of the eigenvalue                              '
309 write(io_unit,10) '#                eig    eigenvalue                                           '
310 write(io_unit,10) '#                                                                            '
311 write(io_unit,10) '#                                                                            '
312 write(io_unit,10) '#                (the following vectors are expressed in crystal units)      '
313 write(io_unit,10) '#                                                                            '
314 write(io_unit,10) '#                R       =   < l | |r| | l >                                 '
315 write(io_unit,10) '#                sigma_R = sqrt{< l | (|r|-R)^2 | l > }                      '
316 write(io_unit,10) '#                                                                            '
317 write(io_unit,10) '#                G       =   < l | |G| | l >                                 '
318 write(io_unit,10) '#                sigma_G = sqrt{< l | (|G|-G)^2 | l > }                      '
319 write(io_unit,10) '#                                                                            '
320 write(io_unit,10) '#                                                                            '
321 write(io_unit,10) '#  l            eig                r         sigma_r       G         sigma_G '
322 write(io_unit,10) '#----------------------------------------------------------------------------'
323 flush(io_unit)
324 
325 
326 G_array(:) = kg_k(1,:)**2+ kg_k(2,:)**2+ kg_k(3,:)**2
327 
328 G_array(:) = sqrt(G_array(:))
329 
330 R_array  = zero
331 
332 do ir3=1,n3
333 
334 if (ir3 <= n3/2 ) then
335   z = (one*ir3)/(one*n3)
336 else
337   z = (one*ir3)/(one*n3)-one
338 end if
339 
340 do ir2=1,n2
341 
342 if (ir2 <= n2/2 ) then
343   y = (one*ir2)/(one*n2)
344 else
345   y = (one*ir2)/(one*n2)-one
346 end if
347 
348 do ir1=1,n1
349 
350 if (ir1 <= n1/2 ) then
351   x = (one*ir1)/(one*n1)
352 else
353   x = (one*ir1)/(one*n1)-one
354 end if
355 
356 
357 R_array(ir1,ir2,ir3)  = sqrt(x**2+y**2+z**2)
358 end do
359 end do
360 end do
361 
362 
363 
364 do l=1, lmax
365 
366 psik(1,:) = dble (Lbasis_diag(:,l))
367 psik(2,:) = dimag(Lbasis_diag(:,l))
368 
369 call g_to_r(psir ,psik)
370 
371 
372 G       = sum(G_array(:)*(psik(1,:)**2+psik(2,:)**2))
373 R       = sum(R_array(:,:,:)*(psir(1,:,:,:)**2+psir(2,:,:,:)**2) )*ucvol/nfft
374 
375 sigma_G = sqrt(sum((G_array(:)    -G)**2*(psik(1,:)**2    +psik(2,:)**2)))
376 sigma_R = sqrt(sum((R_array(:,:,:)-R)**2*(psir(1,:,:,:)**2+psir(2,:,:,:)**2))*ucvol/nfft)
377 
378 
379 
380 write(io_unit,20) l, epsilon_eigenvalues(l), R,sigma_R, G,sigma_G
381 
382 end do
383 
384 close(io_unit)
385 
386 ABI_DEALLOCATE(seeds)
387 ABI_DEALLOCATE(Lbasis_diag)
388 ABI_DEALLOCATE(psik)
389 ABI_DEALLOCATE(psir)
390 ABI_DEALLOCATE(G_array)
391 ABI_DEALLOCATE(R_array)
392 
393 ABI_DEALLOCATE(epsilon_eigenvalues)
394 
395 
396 10 format(A)
397 20 format(I5,ES24.16,4F12.8)
398 
399 end subroutine GeneratePrintDielectricEigenvalues