TABLE OF CONTENTS


ABINIT/m_gwls_GenerateEpsilon [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_GenerateEpsilon

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

m_gwls_GenerateEpsilon/driver_generate_dielectric_matrix [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  driver_generate_dielectric_matrix

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 80 subroutine driver_generate_dielectric_matrix(epsilon_matrix_function,nseeds,kmax,&
 81 epsilon_eigenvalues,Lbasis,debug)
 82 !----------------------------------------------------------------------
 83 ! This routine computes the Lanczos approximate representation of the
 84 ! implicit dielectic operator and then diagonalizes the banded
 85 ! Lanczos matrix.
 86 !----------------------------------------------------------------------
 87 interface
 88   subroutine epsilon_matrix_function(v_out,v_in,l)
 89   use defs_basis
 90 
 91   integer,     intent(in)  :: l
 92   complex(dp), intent(out) :: v_out(l)
 93   complex(dp), intent(in)  :: v_in(l)
 94 
 95   end subroutine epsilon_matrix_function
 96 end interface
 97 
 98 integer,       intent(in) :: nseeds, kmax
 99 logical,       intent(in) :: debug
100 
101 real   (dp),  intent(out) :: epsilon_eigenvalues(nseeds*kmax)
102 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax)  ! array containing the Lanczos basis
103 
104 
105 ! local variables
106 
107 complex(dpc), allocatable :: seeds(:,:)
108 
109 complex(dpc),allocatable :: alpha(:,:,:)
110 complex(dpc),allocatable :: beta (:,:,:)
111 
112 integer :: mpi_communicator
113 
114 ! *************************************************************************
115 
116 
117 ! The epsilon operator will act in LA mode.
118 mpi_communicator = mpi_enreg%comm_bandfft
119 
120 
121 !Create seeds
122 ABI_MALLOC(seeds,(npw_k,nseeds))
123 call get_seeds(first_seed, nseeds, seeds)
124 
125 ! compute the Lanczos basis
126 ABI_MALLOC(alpha,(nseeds,nseeds,kmax))
127 ABI_MALLOC(beta ,(nseeds,nseeds,kmax))
128 
129 call block_lanczos_algorithm(mpi_communicator,epsilon_matrix_function,kmax,nseeds,npw_k,        &
130 seeds,alpha,beta,Lbasis)
131 
132 ! Diagonalize the epsilon matrix, which is banded
133 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug)
134 
135 if (debug) then
136   call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues)
137 end if
138 
139 ABI_FREE(seeds)
140 ABI_FREE(alpha)
141 ABI_FREE(beta)
142 
143 end subroutine driver_generate_dielectric_matrix

m_gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  Driver_GeneratePrintDielectricEigenvalues

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

371 subroutine Driver_GeneratePrintDielectricEigenvalues(dtset)
372 !----------------------------------------------------------------------
373 ! Compute the eigenvalues of the various dielectric operators
374 !----------------------------------------------------------------------
375 type(dataset_type),intent(in) :: dtset
376 
377 integer  ::  kmax_exact, kmax_model, kmax
378 real(dp) :: second_model_parameter
379 
380 
381 integer  ::  lm, k, lmax, l1, l2
382 integer  ::  io_unit
383 integer  ::  io_unit2
384 
385 real(dp) :: time1, time2, time
386 ! local variables
387 character(128)  :: output_filename
388 character(256)  :: timing_string
389 
390 
391 complex(dpc), allocatable :: Lbasis_exact(:,:)
392 complex(dpc), allocatable :: Lbasis_model(:,:)
393 
394 complex(dpc), allocatable :: sub_Lbasis_exact(:,:)
395 complex(dpc), allocatable :: sub_Lbasis_model(:,:)
396 
397 complex(dpc), allocatable :: dummy(:,:)
398 complex(dpc), allocatable :: dummy2(:,:)
399 complex(dpc), allocatable :: dummy3(:,:)
400 
401 complex(dpc), allocatable :: alpha_exact(:,:,:)
402 complex(dpc), allocatable :: beta_exact (:,:,:)
403 
404 complex(dpc), allocatable :: alpha_model(:,:,:)
405 complex(dpc), allocatable :: beta_model (:,:,:)
406 
407 real(dp), allocatable :: eig_exact(:)
408 real(dp), allocatable :: eig_model(:)
409 
410 complex(dpc), allocatable :: model_epsilon_matrix(:,:)
411 complex(dpc), allocatable :: vector(:)
412 
413 real(dpc) :: tr_eps_1, tr_eps_2, tr_eps_3
414 
415 
416 integer   ::  lwork, lrwork, liwork, info
417 complex(dpc), allocatable :: work(:)
418 real(dp)    , allocatable :: rwork(:)
419 integer     , allocatable :: iwork(:)
420 
421 integer        :: debug_unit
422 character(50)  :: debug_filename
423 
424 ! *************************************************************************
425 
426 
427 
428 
429 kmax_exact   = dtset%gwls_stern_kmax
430 kmax_model   = dtset%gwls_kmax_complement
431 
432 !second_model_parameter  = dtset%gwls_second_model_parameter
433 second_model_parameter  = zero
434 
435 
436 ! global stuff
437 nseeds       = dtset%gwls_nseeds
438 first_seed   = dtset%gwls_first_seed
439 e            = dtset%gwls_band_index
440 
441 
442 
443 ABI_MALLOC(Lbasis_exact,(npw_k,kmax_exact*nseeds))
444 ABI_MALLOC(Lbasis_model,(npw_k,kmax_model*nseeds))
445 
446 ABI_MALLOC(alpha_exact, (nseeds,nseeds,kmax_exact))
447 ABI_MALLOC(beta_exact , (nseeds,nseeds,kmax_exact))
448 ABI_MALLOC(alpha_model, (nseeds,nseeds,kmax_model))
449 ABI_MALLOC(beta_model , (nseeds,nseeds,kmax_model))
450 
451 
452 ! set omega=0 for exact dielectric operator
453 call set_dielectric_function_frequency([zero,zero])
454 
455 
456 
457 call cpu_time(time1)
458 output_filename = 'EIGENVALUES_EXACT.dat'
459 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_k, nseeds, kmax_exact, &
460 output_filename, Lbasis_exact, alpha_exact, beta_exact)
461 
462 
463 
464 call cpu_time(time2)
465 time = time2-time1
466 write(timing_string,'(A)')  "Time to compute the EXACT Static Dielectric Matrix  :   "
467 call write_timing_log(timing_string,time)
468 
469 
470 
471 call cpu_time(time1)
472 call setup_Pk_model(zero,second_model_parameter)
473 output_filename = 'EIGENVALUES_MODEL.dat'
474 call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_model_operator, nseeds, kmax_model, &
475 &output_filename, Lbasis_model, alpha_model, beta_model)
476 call cpu_time(time2)
477 time = time2-time1
478 write(timing_string,'(A)')  "Time to compute the MODEL Static Dielectric Matrix  :   "
479 call write_timing_log(timing_string,time)
480 
481 
482 call cpu_time(time1)
483 if (kmax_exact <= kmax_model) then
484   kmax  = kmax_exact
485 else
486   kmax  = kmax_model
487 end if
488 lmax = nseeds*kmax
489 
490 ! Build model operator matrix elements in the exact basis
491 ABI_MALLOC(model_epsilon_matrix, (lmax,lmax))
492 ABI_MALLOC(vector, (npw_k))
493 
494 model_epsilon_matrix = cmplx_0
495 
496 do l2 =1 , lmax
497 call matrix_function_epsilon_model_operator(vector ,Lbasis_exact(:,l2),npw_k)
498 
499 
500 do l1 =1, lmax
501 
502 model_epsilon_matrix(l1, l2) = complex_vector_product(Lbasis_exact(:,l1),vector,npw_k)
503 
504 end do
505 
506 
507 end do
508 
509 ABI_FREE(vector)
510 
511 call cpu_time(time2)
512 time = time2-time1
513 write(timing_string,'(A)')  "Compute MODEL matrix elements in EXACT basis        :   "
514 call write_timing_log(timing_string,time)
515 
516 
517 
518 
519 ! Compare the traces
520 
521 
522 io_unit = get_unit()
523 open(file='DIELECTRIC_TRACE.dat',status=files_status_new,unit=io_unit)
524 write(io_unit,10) '#----------------------------------------------------------------------------'
525 write(io_unit,10) '#                                                                            '
526 write(io_unit,10) '#                       Partial traces                                       '
527 write(io_unit,10) '#           ==========================================                       '
528 write(io_unit,10) '#                                                                            '
529 write(io_unit,10) '#  Tabulate the trace of various operators, as function of the number        '
530 write(io_unit,10) '#  of lanczos steps performed.                                               '
531 write(io_unit,10) '#                                                                            '
532 write(io_unit,10) '#                                                                            '
533 write(io_unit,10) '#  NOTES:                                                                    '
534 write(io_unit,10) '#         Tr[1-eps^{-1}] is evaluated in the Lanczos basis of eps            '
535 write(io_unit,10) '#         Tr[1-eps_m^{-1}] is evaluated in the Lanczos basis of eps_m        '
536 write(io_unit,10) '#         Tr[eps_m^{-1}-eps^{-1}] is evaluated in the Lanczos basis of eps   '
537 write(io_unit,10) '#                                                                            '
538 write(io_unit,10) '#                                                                            '
539 write(io_unit,10) '#  k       Tr[1-eps^{-1}]        Tr[1-eps_m^{-1}]     Tr[eps_m^{-1}-eps^{-1}]'
540 write(io_unit,10) '#----------------------------------------------------------------------------'
541 flush(io_unit)
542 
543 
544 io_unit2 = get_unit()
545 open(file='RPA_ENERGY.dat',status=files_status_new,unit=io_unit2)
546 write(io_unit2,10) '#----------------------------------------------------------------------------'
547 write(io_unit2,10) '#                                                                            '
548 write(io_unit2,10) '#                       RPA TOTAL ENERGY                                     '
549 write(io_unit2,10) '#           ==========================================                       '
550 write(io_unit2,10) '#                                                                            '
551 write(io_unit2,10) '#  It can be shown that the correlation energy, within the RPA, is given     '
552 write(io_unit2,10) '#  by:                                                                       '
553 write(io_unit2,10) '#       E_c = int_0^{infty} dw/(2pi) Tr[ ln(eps{iw)}+1-eps(iw)]              '
554 write(io_unit2,10) '#                                                                            '
555 write(io_unit2,10) '#  As a gauge of what can be expected as far as convergence is concerned,    '
556 write(io_unit2,10) '#  the following will be printed.                                            '
557 write(io_unit2,10) '#                                                                            '
558 write(io_unit2,10) '#  I_1 = Tr[  ln(eps)  + 1 - eps   ]                                         '
559 write(io_unit2,10) '#  I_2 = Tr[ ln(eps_m) + 1 - eps_m ]                                         '
560 write(io_unit2,10) '#  I_3 = Tr[ ln(eps^{-1}_m . eps ) + eps_m - eps ]                           '
561 write(io_unit2,10) '#                                                                            '
562 write(io_unit2,10) '#                                                                            '
563 write(io_unit2,10) '#                                                                            '
564 write(io_unit2,10) '#  k       I_1                   I_2                  I_3                    '
565 write(io_unit2,10) '#----------------------------------------------------------------------------'
566 flush(io_unit2)
567 
568 
569 ! Iterate every 10 values of k max, or else the linear algebra gets too expensive...
570 do k = 4, kmax, 4
571 
572 ABI_MALLOC(sub_Lbasis_exact,(npw_k,k*nseeds))
573 ABI_MALLOC(sub_Lbasis_model,(npw_k,k*nseeds))
574 
575 
576 ABI_MALLOC(eig_exact,(k*nseeds))
577 ABI_MALLOC(eig_model,(k*nseeds))
578 ABI_MALLOC(dummy,(k*nseeds,k*nseeds))
579 ABI_MALLOC(dummy2,(k*nseeds,k*nseeds))
580 
581 sub_Lbasis_exact(:,:) = Lbasis_exact(:,1:k*nseeds)
582 sub_Lbasis_model(:,:) = Lbasis_model(:,1:k*nseeds)
583 
584 ! Diagonalize the epsilon matrix, which is banded
585 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_exact(:,:,1:k),beta_exact(:,:,1:k),sub_Lbasis_exact,eig_exact,.false.)
586 call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_model(:,:,1:k),beta_model(:,:,1:k),sub_Lbasis_model,eig_model,.false.)
587 
588 tr_eps_1 = sum(one-one/eig_exact(:))
589 tr_eps_2 = sum(one-one/eig_model(:))
590 
591 tr_eps_3 = -sum(one/eig_exact(:))
592 
593 dummy(:,:) = model_epsilon_matrix(1:k*nseeds, 1:k*nseeds)
594 
595 call driver_invert_positive_definite_hermitian_matrix(dummy,k*nseeds)
596 
597 do lm = 1, k*nseeds
598 
599 tr_eps_3 = tr_eps_3 +  dble(dummy(lm,lm))
600 end do
601 
602 
603 write(io_unit,20) k, tr_eps_1, tr_eps_2, tr_eps_3
604 flush(io_unit)
605 
606 
607 
608 tr_eps_1 = sum(log(eig_exact(:))+one-eig_exact(:))
609 tr_eps_2 = sum(log(eig_model(:))+one-eig_model(:))
610 
611 dummy2(:,:) = zero
612 tr_eps_3    = zero
613 do lm = 1, k*nseeds
614 dummy2(lm,lm) = eig_exact(lm)
615 tr_eps_3 = tr_eps_3 + dble(model_epsilon_matrix(lm,lm)) - dble(eig_exact(lm))
616 end do
617 
618 ABI_MALLOC(dummy3,(k*nseeds,k*nseeds))
619 call ZGEMM(      'N',   & ! Hermitian conjugate the first array
620 'N',   & ! Leave second array as is
621 k*nseeds,   & ! the number of rows of the  matrix op( A )
622 k*nseeds,   & ! the number of columns of the  matrix op( B )
623 k*nseeds,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
624 cmplx_1,   & ! alpha constant
625 dummy2,   & ! matrix A
626 k*nseeds,   & ! LDA
627 dummy,   & ! matrix B
628 k*nseeds,   & ! LDB
629 cmplx_0,   & ! beta constant
630 dummy3,   & ! matrix C
631 k*nseeds)     ! LDC
632 
633 dummy2(:,:) = dummy3(:,:)
634 ABI_FREE(dummy3)
635 
636 ! find eigenvalues
637 !call heevd(dummy2, eig_exact)
638 
639 lwork  = k*nseeds+1
640 lrwork = k*nseeds
641 liwork = 1
642 
643 ABI_MALLOC(work,(lwork))
644 ABI_MALLOC(rwork,(lrwork))
645 ABI_MALLOC(iwork,(liwork))
646 
647 call zheevd('N', 'U',k*nseeds, dummy2, k*nseeds, eig_exact, work, lwork, rwork, lrwork, iwork, liwork, info)
648 if ( info /= 0) then
649   debug_unit = get_unit()
650   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
651 
652   open(debug_unit,file=trim(debug_filename),status='unknown')
653 
654   write(debug_unit,'(A)')      '*************************************************************************************'
655   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZHEEVD(1), gwls_GenerateEpsilon'
656   write(debug_unit,'(A)')      '*************************************************************************************'
657 
658   close(debug_unit)
659 
660 end if
661 
662 
663 
664 
665 
666 tr_eps_3 =  tr_eps_3 + sum(log(eig_exact))
667 
668 write(io_unit2,20) k, tr_eps_1, tr_eps_2, tr_eps_3
669 flush(io_unit2)
670 
671 
672 
673 ABI_FREE(work)
674 ABI_FREE(rwork)
675 ABI_FREE(iwork)
676 
677 
678 
679 ABI_FREE(sub_Lbasis_exact)
680 ABI_FREE(sub_Lbasis_model)
681 
682 ABI_FREE(eig_exact)
683 ABI_FREE(eig_model)
684 ABI_FREE(dummy)
685 ABI_FREE(dummy2)
686 end do
687 
688 close(io_unit)
689 close(io_unit2)
690 
691 call cpu_time(time2)
692 time = time2-time1
693 write(timing_string,'(A)')  "Time to compute the TRACES of the Dielectric Matrices:   "
694 call write_timing_log(timing_string,time)
695 
696 ABI_FREE(Lbasis_exact)
697 ABI_FREE(Lbasis_model)
698 ABI_FREE(alpha_exact)
699 ABI_FREE(beta_exact )
700 ABI_FREE(alpha_model)
701 ABI_FREE(beta_model )
702 ABI_FREE(model_epsilon_matrix)
703 
704 
705 
706 10 format(A)
707 20 format(I5,3ES24.16)
708 
709 end subroutine Driver_GeneratePrintDielectricEigenvalues

m_gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues [ Functions ]

[ Top ] [ m_gwls_GenerateEpsilon ] [ Functions ]

NAME

  GeneratePrintDielectricEigenvalues

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

159 subroutine GeneratePrintDielectricEigenvalues(epsilon_matrix_function,nseeds,kmax,output_filename,Lbasis,alpha,beta)
160 !----------------------------------------------------------------------
161 ! This routine computes the Lanczos approximate representation of the
162 ! implicit dielectic operator and then diagonalizes the banded
163 ! Lanczos matrix.
164 !----------------------------------------------------------------------
165 interface
166   subroutine epsilon_matrix_function(v_out,v_in,l)
167   use defs_basis
168 
169   integer,     intent(in)  :: l
170   complex(dp), intent(out) :: v_out(l)
171   complex(dp), intent(in)  :: v_in(l)
172 
173   end subroutine epsilon_matrix_function
174 end interface
175 
176 integer,       intent(in) :: nseeds, kmax
177 
178 character(*),  intent(in) :: output_filename
179 
180 
181 complex(dpc), intent(out) :: Lbasis(npw_k,nseeds*kmax)
182 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax)
183 complex(dpc), intent(out) :: beta (nseeds,nseeds,kmax)
184 
185 
186 ! local variables
187 
188 
189 complex(dpc),allocatable :: seeds(:,:)
190 
191 complex(dpc),allocatable :: Lbasis_diag(:,:)
192 
193 
194 real(dp),    allocatable :: psik(:,:)
195 real(dp),    allocatable :: psir(:,:,:,:)
196 
197 real(dp),    allocatable :: epsilon_eigenvalues(:)
198 
199 
200 integer :: mpi_communicator
201 integer :: io_unit
202 integer :: lmax
203 integer :: l
204 integer :: ir1, ir2, ir3
205 integer :: n1, n2, n3
206 
207 real(dp) :: R, G
208 real(dp) :: sigma_R, sigma_G
209 real(dp) :: x, y, z
210 
211 real(dp),allocatable :: G_array(:)
212 real(dp),allocatable :: R_array(:,:,:)
213 
214 logical :: debug
215 
216 ! *************************************************************************
217 
218 
219 debug = .false.
220 lmax = kmax*nseeds
221 mpi_communicator = mpi_enreg%comm_bandfft
222 !Create seeds
223 ABI_MALLOC(seeds,(npw_k,nseeds))
224 call get_seeds(first_seed, nseeds, seeds)
225 
226 ! compute the Lanczos basis
227 ABI_MALLOC(Lbasis_diag,(npw_k,lmax))
228 ABI_MALLOC(epsilon_eigenvalues,(lmax))
229 
230 ABI_MALLOC(psik,(2,npw_k))
231 ABI_MALLOC(psir,(2,n4,n5,n6))
232 ABI_MALLOC(G_array,(npw_k))
233 ABI_MALLOC(R_array,(n4,n5,n6))
234 
235 psir = zero
236 R_array = zero
237 
238 n1 = n4-1
239 n2 = n5-1
240 n3 = n6
241 
242 ! Generate the Lanczos basis and banded eigenvalue representation
243 call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k,        seeds,alpha,beta,Lbasis)
244 
245 Lbasis_diag = Lbasis
246 
247 ! Diagonalize the epsilon matrix, which is banded
248 call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis_diag,epsilon_eigenvalues,debug)
249 
250 call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,lmax,npw_k,Lbasis_diag,epsilon_eigenvalues)
251 
252 io_unit = get_unit()
253 open(file=output_filename,status=files_status_new,unit=io_unit)
254 write(io_unit,10) '#----------------------------------------------------------------------------'
255 write(io_unit,10) '#                                                                            '
256 write(io_unit,10) '#                       Partial eigenvalues                                  '
257 write(io_unit,10) '#           ==========================================                       '
258 write(io_unit,10) '#                                                                            '
259 write(io_unit,10) '#  Tabulate the eigenvalues of the dielectic matrix, as well as some         '
260 write(io_unit,10) '#  information regarding the eigenstates.                                    '
261 write(io_unit,10) '#                                                                            '
262 write(io_unit,10) '#                                                                            '
263 write(io_unit,10) '#  definitions:                                                              '
264 write(io_unit,10) '#                l      index of the eigenvalue                              '
265 write(io_unit,10) '#                eig    eigenvalue                                           '
266 write(io_unit,10) '#                                                                            '
267 write(io_unit,10) '#                                                                            '
268 write(io_unit,10) '#                (the following vectors are expressed in crystal units)      '
269 write(io_unit,10) '#                                                                            '
270 write(io_unit,10) '#                R       =   < l | |r| | l >                                 '
271 write(io_unit,10) '#                sigma_R = sqrt{< l | (|r|-R)^2 | l > }                      '
272 write(io_unit,10) '#                                                                            '
273 write(io_unit,10) '#                G       =   < l | |G| | l >                                 '
274 write(io_unit,10) '#                sigma_G = sqrt{< l | (|G|-G)^2 | l > }                      '
275 write(io_unit,10) '#                                                                            '
276 write(io_unit,10) '#                                                                            '
277 write(io_unit,10) '#  l            eig                r         sigma_r       G         sigma_G '
278 write(io_unit,10) '#----------------------------------------------------------------------------'
279 flush(io_unit)
280 
281 
282 G_array(:) = kg_k(1,:)**2+ kg_k(2,:)**2+ kg_k(3,:)**2
283 
284 G_array(:) = sqrt(G_array(:))
285 
286 R_array  = zero
287 
288 do ir3=1,n3
289 
290 if (ir3 <= n3/2 ) then
291   z = (one*ir3)/(one*n3)
292 else
293   z = (one*ir3)/(one*n3)-one
294 end if
295 
296 do ir2=1,n2
297 
298 if (ir2 <= n2/2 ) then
299   y = (one*ir2)/(one*n2)
300 else
301   y = (one*ir2)/(one*n2)-one
302 end if
303 
304 do ir1=1,n1
305 
306 if (ir1 <= n1/2 ) then
307   x = (one*ir1)/(one*n1)
308 else
309   x = (one*ir1)/(one*n1)-one
310 end if
311 
312 
313 R_array(ir1,ir2,ir3)  = sqrt(x**2+y**2+z**2)
314 end do
315 end do
316 end do
317 
318 
319 
320 do l=1, lmax
321 
322 psik(1,:) = dble (Lbasis_diag(:,l))
323 psik(2,:) = dimag(Lbasis_diag(:,l))
324 
325 call g_to_r(psir ,psik)
326 
327 
328 G       = sum(G_array(:)*(psik(1,:)**2+psik(2,:)**2))
329 R       = sum(R_array(:,:,:)*(psir(1,:,:,:)**2+psir(2,:,:,:)**2) )*ucvol/nfft
330 
331 sigma_G = sqrt(sum((G_array(:)    -G)**2*(psik(1,:)**2    +psik(2,:)**2)))
332 sigma_R = sqrt(sum((R_array(:,:,:)-R)**2*(psir(1,:,:,:)**2+psir(2,:,:,:)**2))*ucvol/nfft)
333 
334 
335 
336 write(io_unit,20) l, epsilon_eigenvalues(l), R,sigma_R, G,sigma_G
337 
338 end do
339 
340 close(io_unit)
341 
342 ABI_FREE(seeds)
343 ABI_FREE(Lbasis_diag)
344 ABI_FREE(psik)
345 ABI_FREE(psir)
346 ABI_FREE(G_array)
347 ABI_FREE(R_array)
348 
349 ABI_FREE(epsilon_eigenvalues)
350 
351 
352 10 format(A)
353 20 format(I5,ES24.16,4F12.8)
354 
355 end subroutine GeneratePrintDielectricEigenvalues