Branch data 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 : : !!***
|
