Line data Source code
1 : !!****m* ABINIT/m_gwls_ComputePoles
2 : !! NAME
3 : !! m_gwls_ComputePoles
4 : !!
5 : !! FUNCTION
6 : !! .
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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
61 : !!***
62 :
63 : !!****f* m_gwls_ComputePoles/generate_degeneracy_table_for_poles
64 : !! NAME
65 : !! generate_degeneracy_table_for_poles
66 : !!
67 : !! FUNCTION
68 : !! .
69 : !!
70 : !! INPUTS
71 : !!
72 : !! OUTPUT
73 : !!
74 : !! SOURCE
75 :
76 7 : 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 7 : degeneracy_tolerance = 1.0D-8
102 : !--------------------------------------------------------------------------------
103 : !
104 : ! First, find the largest degeneracy in the eigenvalue spectrum
105 : !
106 : !--------------------------------------------------------------------------------
107 :
108 7 : nbands = size(eig)
109 :
110 : ! initialize
111 7 : largest_degeneracy = 1
112 7 : number_of_denerate_sets = 1
113 :
114 7 : n_degeneracy = 1
115 7 : 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 70 : do n = 2, nbands
134 63 : if (abs(eig(n) - energy) < degeneracy_tolerance) then
135 : ! A degenerate state! add one
136 35 : n_degeneracy = n_degeneracy + 1
137 : else
138 : ! We are no longer degenerate. Update
139 28 : if ( n_degeneracy > largest_degeneracy ) largest_degeneracy = n_degeneracy
140 :
141 28 : n_degeneracy = 1
142 28 : energy = eig(n)
143 28 : 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 70 : 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 28 : ABI_MALLOC(degeneracy_table, (number_of_denerate_sets,largest_degeneracy))
157 21 : ABI_MALLOC(number_of_degenerate_states, (number_of_denerate_sets))
158 :
159 133 : degeneracy_table(:,:) = 0
160 :
161 7 : i_set = 1
162 7 : j_deg = 1
163 :
164 : ! initialize
165 7 : energy = eig(1)
166 :
167 7 : degeneracy_table(i_set,j_deg) = 1
168 7 : number_of_degenerate_states(i_set) = 1
169 :
170 70 : do n = 2, nbands
171 :
172 63 : if (abs(eig(n) - energy) < degeneracy_tolerance) then
173 : ! A degenerate state! add one
174 35 : j_deg = j_deg + 1
175 :
176 : else
177 : ! We are no longer degenerate. Update
178 28 : j_deg = 1
179 28 : i_set = i_set+1
180 28 : energy = eig(n)
181 :
182 : end if
183 :
184 :
185 63 : number_of_degenerate_states(i_set) = j_deg
186 70 : degeneracy_table(i_set,j_deg) = n
187 :
188 : end do
189 :
190 7 : if (debug .and. mpi_enreg%me == 0) then
191 4 : io_unit = get_unit()
192 4 : filename = "degeneracy_table.log"
193 :
194 4 : i = 0
195 4 : inquire(file=filename,exist=file_exists)
196 7 : do while (file_exists)
197 3 : i = i+1
198 3 : write (filename,'(A,I0,A)') "degeneracy_table_",i,".log"
199 7 : inquire(file=filename,exist=file_exists)
200 : end do
201 :
202 4 : io_unit = get_unit()
203 :
204 4 : open(io_unit,file=filename,status=files_status_new)
205 :
206 4 : write(io_unit,10) " "
207 4 : write(io_unit,10) "#==============================================================================================="
208 4 : write(io_unit,10) "# Degeneracy table : tabulate the degenerate states "
209 4 : write(io_unit,10) "# ------------------------------------------------------- "
210 4 : write(io_unit,10) "# "
211 4 : write(io_unit,14) "# number_of_denerate_sets = ", number_of_denerate_sets
212 4 : write(io_unit,10) "# "
213 4 : write(io_unit,14) "# largest_degeneracy = ", largest_degeneracy
214 4 : write(io_unit,10) "# "
215 4 : write(io_unit,10) "# Eigenvalues (Ha) "
216 4 : write(io_unit,10) "#==============================================================================================="
217 4 : write(io_unit,16) eig(:)
218 :
219 :
220 :
221 4 : write(io_unit,10) "#==============================================================================================="
222 4 : write(io_unit,10) "# i_set number of states States "
223 4 : write(io_unit,10) "#==============================================================================================="
224 4 : flush(io_unit)
225 :
226 24 : do i_set = 1, number_of_denerate_sets
227 24 : write(io_unit,12) i_set, number_of_degenerate_states(i_set), degeneracy_table(i_set,:)
228 : end do
229 :
230 4 : flush(io_unit)
231 :
232 4 : 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 7 : end subroutine generate_degeneracy_table_for_poles
241 : !!***
242 :
243 : !!****f* m_gwls_ComputePoles/clean_degeneracy_table_for_poles
244 : !! NAME
245 : !! clean_degeneracy_table_for_poles
246 : !!
247 : !! FUNCTION
248 : !! .
249 : !!
250 : !! INPUTS
251 : !!
252 : !! OUTPUT
253 : !!
254 : !! SOURCE
255 :
256 7 : subroutine clean_degeneracy_table_for_poles()
257 :
258 : ! *************************************************************************
259 :
260 7 : if(allocated(degeneracy_table)) then
261 7 : ABI_FREE(degeneracy_table)
262 : end if
263 7 : if(allocated(number_of_degenerate_states)) then
264 7 : ABI_FREE(number_of_degenerate_states)
265 : end if
266 :
267 7 : end subroutine clean_degeneracy_table_for_poles
268 : !!***
269 :
270 : !!****f* m_gwls_ComputePoles/compute_Poles
271 : !! NAME
272 : !! compute_Poles
273 : !!
274 : !! FUNCTION
275 : !! .
276 : !!
277 : !! INPUTS
278 : !!
279 : !! OUTPUT
280 : !!
281 : !!
282 : !!
283 : !! SOURCE
284 :
285 7 : 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 7 : complex(dp), allocatable :: seeds(:,:)
325 :
326 : ! *************************************************************************
327 :
328 7 : compute_Poles = zero
329 :
330 7 : energy_tolerance = 1.0D-8
331 :
332 :
333 7 : if (debug .and. mpi_enreg%me == 0) then
334 0 : io_unit = get_unit()
335 :
336 0 : i = 0
337 :
338 0 : file_exists = .true.
339 0 : do while (file_exists)
340 0 : i = i+1
341 0 : write (filename,'(A,I0.4,A)') "ComputePoles_",i,".log"
342 0 : inquire(file=filename,exist=file_exists)
343 : end do
344 :
345 :
346 0 : open(io_unit,file=filename,status=files_status_new)
347 :
348 0 : write(io_unit,10) " "
349 0 : write(io_unit,10) "#==============================================================================================="
350 0 : write(io_unit,10) "# ComputePoles: debug information for the pole computation "
351 0 : write(io_unit,10) "# -------------------------------------------------------- "
352 0 : write(io_unit,10) "# "
353 0 : write(io_unit,10) "# This file contains data describing the computation of the pole contribution to the "
354 0 : write(io_unit,10) "# correlation self energy. "
355 0 : write(io_unit,10) "# "
356 0 : write(io_unit,10) "#==============================================================================================="
357 0 : write(io_unit,10) " "
358 0 : write(io_unit,10) "#==============================================================================================="
359 0 : write(io_unit,10) "# "
360 0 : write(io_unit,10) "# parameters: "
361 0 : write(io_unit,10) "# "
362 0 : write(io_unit,11) "# external_omega = ",external_omega," Ha "
363 0 : write(io_unit,12) "# kmax_poles = ",kmax_poles
364 0 : write(io_unit,12) "# nbandv = ",nbandv
365 0 : write(io_unit,10) "# "
366 0 : write(io_unit,10) "#==============================================================================================="
367 0 : write(io_unit,10) "# "
368 0 : write(io_unit,10) "# DFT Eigenvalues (Ha) "
369 0 : write(io_unit,10) "#==============================================================================================="
370 0 : write(io_unit,13) eig(:)
371 0 : 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 7 : compute_Poles = zero
385 :
386 7 : pole_is_conduction = .false.
387 7 : pole_is_valence = .false.
388 7 : pole_is_in_gap = .false.
389 :
390 :
391 :
392 : ! Careful here! there may be only nbandv states in memory; nbandv+1 causes segfaults!
393 7 : if ( external_omega <= eig(nbandv)) then
394 7 : pole_is_valence = .true.
395 0 : else if ( external_omega > eig(nbandv)) then
396 0 : pole_is_conduction = .true.
397 : else
398 0 : pole_is_in_gap = .true.
399 : end if
400 :
401 :
402 :
403 :
404 7 : if (debug .and. mpi_enreg%me == 0 ) then
405 0 : write(io_unit,10) "#===================================================================================================="
406 0 : write(io_unit,10) "# "
407 0 : write(io_unit,10) "# Determine where the external energy is: "
408 0 : write(io_unit,10) "# "
409 0 : write(io_unit,14) "# pole_is_valence = ",pole_is_valence
410 0 : write(io_unit,14) "# pole_is_conduction = ",pole_is_conduction
411 0 : write(io_unit,14) "# pole_is_in_gap = ",pole_is_in_gap
412 0 : write(io_unit,10) "# "
413 0 : write(io_unit,10) "#===================================================================================================="
414 0 : flush(io_unit)
415 :
416 : end if
417 :
418 7 : if ( pole_is_in_gap) return
419 :
420 : !--------------------------------------------------------------------------------
421 : !
422 : ! Loop on all degenerate sets
423 : !
424 : !--------------------------------------------------------------------------------
425 :
426 7 : if (debug .and. mpi_enreg%me == 0 ) then
427 0 : write(io_unit,10) "#===================================================================================================="
428 0 : write(io_unit,10) "# "
429 0 : write(io_unit,10) "# Iterating over all degenerate sets of eigenvalues: "
430 0 : write(io_unit,10) "# "
431 0 : write(io_unit,10) "#===================================================================================================="
432 0 : flush(io_unit)
433 :
434 : end if
435 :
436 :
437 21 : do i_set =1, number_of_denerate_sets
438 :
439 21 : n = degeneracy_table(i_set,1)
440 :
441 21 : En_m_omega = eig(n)-external_omega
442 :
443 21 : if (debug .and. mpi_enreg%me == 0) then
444 0 : write(io_unit,12) "# i_set = ", i_set
445 0 : write(io_unit,12) "# n = ",n
446 0 : write(io_unit,16) "# eig(n)-omega = ",En_m_omega," Ha"
447 0 : flush(io_unit)
448 : end if
449 :
450 : !------------------------------------------
451 : ! Test if we need to exit the loop
452 : !------------------------------------------
453 21 : if (pole_is_valence ) then
454 :
455 : ! If the pole is valence, get out when we enter conduction states
456 21 : if ( n > nbandv ) then
457 7 : if (debug.and. mpi_enreg%me == 0) then
458 0 : write(io_unit,10) "#"
459 0 : write(io_unit,10) "# n > nbandv : exit loop!"
460 0 : 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 14 : if (En_m_omega < zero .and. abs(En_m_omega) > energy_tolerance) then
470 : ! careful close to zero!
471 0 : if (debug .and. mpi_enreg%me == 0) then
472 0 : write(io_unit,10) "# "
473 0 : write(io_unit,10) "# eig(n) < omega : cycle!"
474 0 : 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 14 : prefactor = -one
482 :
483 0 : 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 0 : if ( En_m_omega > energy_tolerance ) then
488 0 : if (debug .and. mpi_enreg%me == 0) then
489 0 : write(io_unit,10) "#"
490 0 : write(io_unit,10) "# eig(n) > omega : exit!"
491 0 : 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 0 : if ( n <= nbandv ) then
499 0 : if (debug .and. mpi_enreg%me == 0) then
500 0 : write(io_unit,10) "#"
501 0 : write(io_unit,10) "# n <= nbandv : cycle!"
502 0 : 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 0 : prefactor = one
510 : end if
511 :
512 :
513 : !-------------------------------------------------
514 : ! If we made it this far, we have a contribution!
515 : !-------------------------------------------------
516 :
517 14 : if (abs(En_m_omega) < energy_tolerance ) then
518 :
519 0 : if (debug .and. mpi_enreg%me == 0) then
520 0 : write(io_unit,10) "# "
521 0 : write(io_unit,10) "# En - omega ~ 0: pole at the origin, multiply by 1/2!"
522 0 : 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 0 : prefactor = 0.5_dp*prefactor
528 : end if
529 :
530 :
531 :
532 14 : number_of_seeds = number_of_degenerate_states(i_set)
533 :
534 56 : ABI_MALLOC(seeds, (npw_k,number_of_seeds))
535 :
536 14 : call get_seeds(n, number_of_seeds, seeds) !Missing wrappers
537 :
538 42 : call set_dielectric_function_frequency([En_m_omega,zero])
539 14 : if (debug .and. mpi_enreg%me == 0) then
540 0 : write(io_unit,10) "# Compute pole contribution:"
541 0 : write(io_unit,12) "# number of seeds = ",number_of_seeds
542 0 : write(io_unit,16) "# || seeds || = ",sqrt(sum(abs(seeds(:,:))**2)) !Missing xmpi_sum
543 0 : write(io_unit,16) "# eig(n)-omega = ",En_m_omega, " Ha"
544 0 : write(io_unit,17) "# prefactor = ",prefactor
545 0 : flush(io_unit)
546 : end if
547 14 : if(dtset%zcut > tol12) activate_inf_shift_poles = .true.
548 14 : 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 14 : seeds,debug)
553 14 : if(dtset%zcut > tol12) activate_inf_shift_poles = .false.
554 14 : if (debug .and. mpi_enreg%me == 0) then
555 0 : write(io_unit,16) "# pole contribution = ",prefactor*pole_contribution, " Ha"
556 0 : flush(io_unit)
557 : end if
558 :
559 :
560 14 : compute_Poles = compute_Poles + prefactor*pole_contribution
561 14 : ABI_FREE(seeds)
562 : end do
563 :
564 7 : if (debug .and. mpi_enreg%me == 0) then
565 0 : 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
578 : !!***
579 :
580 : !!****f* m_gwls_ComputePoles/compute_pole_contribution
581 : !! NAME
582 : !! compute_pole_contribution
583 : !!
584 : !! FUNCTION
585 : !! .
586 : !!
587 : !! INPUTS
588 : !!
589 : !! OUTPUT
590 : !!
591 : !!
592 : !!
593 : !! SOURCE
594 :
595 14 : 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(dp), intent(in) :: seeds(npw_k,nseeds)
616 : logical, intent(in) :: debug
617 :
618 : ! local variables
619 :
620 : integer :: mpi_communicator
621 :
622 14 : complex(dp),allocatable :: local_seeds(:,:)
623 14 : complex(dp),allocatable :: Lbasis(:,:) ! array containing the Lanczos basis
624 14 : complex(dp),allocatable :: alpha(:,:,:)
625 14 : complex(dp),allocatable :: beta (:,:,:)
626 14 : real(dp), allocatable :: epsilon_eigenvalues(:)
627 :
628 : real(dp):: matrix_elements
629 :
630 : complex(dp) :: cmplx_value
631 : integer :: l, s
632 : integer :: ierr
633 :
634 : ! *************************************************************************
635 :
636 : ! compute the Lanczos basis
637 70 : ABI_MALLOC(alpha,(nseeds,nseeds,kmax))
638 56 : ABI_MALLOC(beta ,(nseeds,nseeds,kmax))
639 56 : ABI_MALLOC(Lbasis,(npw_k,nseeds*kmax))
640 56 : ABI_MALLOC(local_seeds,(npw_k,nseeds))
641 42 : ABI_MALLOC(epsilon_eigenvalues, (nseeds*kmax))
642 :
643 :
644 14 : mpi_communicator = mpi_enreg%comm_bandfft !Missing maybe something for easy access of LA and FFT comms?
645 :
646 4154 : local_seeds(:,:) = seeds(:,:)
647 :
648 : call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, &
649 14 : & local_seeds,alpha,beta,Lbasis)
650 :
651 14 : write(std_out,*) "alpha:"
652 70 : do l=1,kmax
653 168 : do s=1,nseeds
654 168 : write(std_out,*) alpha(:,s,l)
655 : end do
656 70 : write(std_out,*) " "
657 : end do
658 :
659 14 : write(std_out,*) "beta:"
660 70 : do l=1,kmax
661 168 : do s=1,nseeds
662 168 : write(std_out,*) beta(:,s,l)
663 : end do
664 70 : write(std_out,*) " "
665 : end do
666 :
667 14 : ABI_FREE(local_seeds)
668 :
669 : ! Diagonalize the epsilon matrix, which is banded
670 14 : call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug)
671 :
672 14 : if (debug) then
673 0 : call ritz_analysis_general(mpi_communicator ,epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues)
674 : end if
675 :
676 : compute_pole_contribution = zero
677 :
678 126 : do l = 1, nseeds*kmax
679 :
680 : matrix_elements = zero
681 :
682 392 : do s = 1, nseeds
683 :
684 280 : cmplx_value = complex_vector_product(seeds(:,s),Lbasis(:,l),npw_k)
685 :
686 280 : call xmpi_sum(cmplx_value,mpi_communicator,ierr) ! sum on all processors working on FFT!
687 :
688 392 : matrix_elements = matrix_elements + abs(cmplx_value)**2
689 :
690 :
691 : end do
692 :
693 :
694 : compute_pole_contribution = compute_pole_contribution + &
695 126 : matrix_elements *(one/epsilon_eigenvalues(l)-one)
696 : end do
697 :
698 :
699 :
700 14 : ABI_FREE(alpha)
701 14 : ABI_FREE(beta)
702 14 : ABI_FREE(Lbasis)
703 14 : ABI_FREE(epsilon_eigenvalues)
704 :
705 14 : end function compute_pole_contribution
706 :
707 : end module m_gwls_ComputePoles
708 : !!***
|
