TABLE OF CONTENTS


ABINIT/m_gwls_Projected_AT [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_Projected_AT

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_Projected_AT
24 !----------------------------------------------------------------------------------------------------
25 ! This module contains routines to compute the matrix elements of the so-called A operator,
26 ! which accounts for the Static correlation energy term.
27 !----------------------------------------------------------------------------------------------------
28 ! local modules
29 use m_gwls_utility
30 use m_gwls_wf
31 use m_gwls_TimingLog
32 use m_gwls_hamiltonian
33 use m_gwls_lineqsolver
34 use m_gwls_GWlanczos
35 use m_gwls_LanczosBasis 
36 use m_gwls_LanczosResolvents
37 
38 use m_gwls_GWanalyticPart, only : get_projection_band_indices
39 ! abinit modules
40 use defs_basis
41 use m_abicore
42 use m_xmpi
43 use m_io_tools,  only : get_unit
44 
45 
46 
47 implicit none
48 save
49 private

m_hamiltonian/compute_AT_shift_Lanczos [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_AT_shift_Lanczos

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 73 subroutine compute_AT_shift_Lanczos(nfreq,list_external_omega,model_parameter,lmax, modified_Lbasis,kmax_analytic,list_AT_Lanczos)
 74 !----------------------------------------------------------------------------------------------------
 75 ! This function returns the diagonal matrix elements of the so-called A^T operator, which is pertinent to the
 76 ! computation of the analytic energy term.
 77 !
 78 ! The operator is given by
 79 !
 80 !                 Am(W) = w0/2 Um^dagger  . [ PW/(H-W-w0)+QW/(H-W+w0)] . Um
 81 !
 82 !  where W is the external frequency of the self energy, and w0 is the lorentzian parameter.
 83 !  The operator PW projects on states of energy lower than W, and QW on states or energy higher than W.
 84 !
 85 !
 86 !  For every l, We seek to compute  AT_l = < l |  A^T | l > = < l^* | A | l^* > 
 87 !
 88 !  It will be assumed that the array modified_Lbasis already contains the basis vectors U_m | l^* >.
 89 !
 90 !  This function does not use SQMR, but rather shift Lanczos to extract the values of the matrix elements for 
 91 !  all external frequencies.
 92 !----------------------------------------------------------------------------------------------------
 93 implicit none
 94 
 95 integer,      intent(in) :: nfreq
 96 real(dp),     intent(in) :: list_external_omega(nfreq)
 97 real(dp),     intent(in) :: model_parameter
 98 integer,      intent(in) :: lmax
 99 integer,      intent(in) :: kmax_analytic
100 complex(dpc), intent(in) :: modified_Lbasis(npw_k,lmax)
101 complex(dpc), intent(out):: list_AT_Lanczos(nfreq,lmax)
102 
103 
104 real(dp),     allocatable :: psik_wrk(:,:)
105 real(dp),     allocatable :: psikb_wrk(:,:)
106 real(dp),     allocatable :: psikg_wrk(:,:)
107 
108 real(dp),     allocatable :: psikg(:,:)
109 
110 complex(dpc), allocatable :: seed_vector(:)
111 complex(dpc), allocatable :: list_left_vectors(:,:)
112 complex(dpc), allocatable :: matrix_elements_resolvent(:,:)
113 
114 complex(dpc), allocatable :: list_z_P(:), list_z_Q(:)
115 integer,      allocatable :: frequency_indices_array(:,:)
116 
117 
118 real(dp):: external_omega 
119 logical :: prec
120 integer :: nvec
121 integer :: band_index_below, band_index_above, bib0, bia0, bib, bia
122 
123 integer :: l, lloc, iw_ext, iw_ext_min, iw_ext_max
124 integer :: ierr
125 
126 integer :: number_of_frequency_blocks, ifreq_block
127 
128 integer :: iblk, nbdblock_lanczos
129 integer :: mb 
130 
131 integer :: nz
132 
133 integer :: io_unit
134 integer :: mpi_band_rank
135 
136 ! *************************************************************************
137 
138 mpi_band_rank    = mpi_enreg%me_band
139 
140 !=================================================
141 !
142 ! Find the number of frequency blocks, where
143 ! the projection operators are constant within a
144 ! block.
145 !=================================================
146 
147 if (mpi_enreg%me == 0 ) then
148   io_unit  = get_unit()
149   open(io_unit,file='Frequency_blocks_AT.log',position='append')
150 
151   write(io_unit,10) "#================================================================================"
152   write(io_unit,10) "#                                                                                "
153   write(io_unit,10) "#  This log file documents how the algorithm in compute_AT_shift_Lanczos         "
154   write(io_unit,10) "#  separates the external frequencies in blocks. It is quite easy to put         "
155   write(io_unit,10) "#  bugs in this algorithm, so monitoring is critical.                            "
156   write(io_unit,10) "#                                                                                "
157   write(io_unit,10) "#================================================================================"
158 
159 
160   write(io_unit,10) "                                                                                 "
161   write(io_unit,10) "#================================================================================"
162   write(io_unit,10) "#                                                                                "
163   write(io_unit,10) "#  input parameters:                                                             "
164   write(io_unit,10) "#                                                                                "
165   write(io_unit,15) "#                      nfreq    : ",nfreq                                                     
166   write(io_unit,17) "#          list_external_omega  : ",list_external_omega  
167   write(io_unit,17) "#              model_parameter  : ",model_parameter
168   write(io_unit,15) "#                       lmax    : ",lmax
169   write(io_unit,15) "#                kmax_analytic  : ",kmax_analytic
170   write(io_unit,10) "#                                                                                "
171   write(io_unit,10) "#================================================================================"
172   write(io_unit,10) "#                                                                                "
173   write(io_unit,10) "#  Building the blocks                                                           "
174   write(io_unit,10) "#                                                                                "
175   write(io_unit,10) "#         bib : band index below                                                 "
176   write(io_unit,10) "#         bia : band index above                                                 "
177   write(io_unit,10) "#         nfb : number_of_frequency_blocks                                       "
178   write(io_unit,10) "#                                                                                "
179   write(io_unit,10) "#                                                                                "
180   write(io_unit,10) "#  iw_ext      external_omega (Ha)     bib  bia  nfb                             "
181   write(io_unit,10) "#================================================================================"
182 end if
183 external_omega = list_external_omega(1)
184 
185 
186 ! bib == band_index_below
187 ! bia == band_index_above
188 call get_projection_band_indices(external_omega,bib0, bia0)
189 
190 number_of_frequency_blocks = 1
191 
192 ! loop on all external frequencies
193 do iw_ext = 1 , nfreq
194 
195 external_omega = list_external_omega(iw_ext)
196 ! Define the energy of the state to be corrected
197 
198 ! Find the indices for the projections PW and QW
199 call get_projection_band_indices(external_omega,bib, bia)
200 
201 
202 if (mpi_enreg%me == 0 ) write(io_unit,20) iw_ext, external_omega, bib, bia, number_of_frequency_blocks 
203 
204 if (bib /= bib0 .or. bia /= bia0) then
205 
206   if (mpi_enreg%me == 0 )   write(io_unit,10) "*************    new block!    **********************"
207   bib0 = bib
208   bia0 = bia
209 
210   number_of_frequency_blocks = number_of_frequency_blocks + 1
211 
212 end if
213 end do
214 
215 !=================================================
216 !
217 ! fill the frequency indices array
218 !
219 !=================================================
220 
221 if (mpi_enreg%me == 0 ) then
222   write(io_unit,10) "                                                                                 "
223   write(io_unit,10) "#================================================================================"
224   write(io_unit,10) "#                                                                                "
225   write(io_unit,10) "#  Filling the frequency indices array                                           "
226   write(io_unit,10) "#                                                                                "
227   write(io_unit,10) "#  iw_ext      iw_ext_min     iw_ext_max  ifreq_block                            "
228   write(io_unit,10) "#================================================================================"
229 end if
230 
231 
232 ABI_MALLOC(frequency_indices_array, (2,number_of_frequency_blocks))
233 frequency_indices_array = 0
234 
235 iw_ext_min = 1
236 iw_ext_max = 1
237 
238 ! loop on all external frequencies
239 external_omega = list_external_omega(1)
240 call get_projection_band_indices(external_omega,bib0, bia0)
241 ifreq_block = 1
242 
243 do iw_ext = 1 , nfreq
244 
245 if (mpi_enreg%me == 0 ) write(io_unit,30) iw_ext, iw_ext_min, iw_ext, ifreq_block
246 
247 external_omega = list_external_omega(iw_ext)
248 ! Define the energy of the state to be corrected
249 
250 ! Find the indices for the projections PW and QW
251 call get_projection_band_indices(external_omega,bib, bia)
252 
253 if (bib /= bib0 .or. bia /= bia0) then
254 
255   if (mpi_enreg%me == 0 )   write(io_unit,10) "*************    new block!    **********************"
256   bib0 = bib
257   bia0 = bia
258 
259   ! write previous block
260   frequency_indices_array(1,ifreq_block) = iw_ext_min
261   frequency_indices_array(2,ifreq_block) = iw_ext-1 ! we went one too far
262 
263   ifreq_block  = ifreq_block  + 1
264   iw_ext_min = iw_ext
265 end if 
266 
267 end do
268 
269 ! write last block!
270 frequency_indices_array(1,ifreq_block) = iw_ext_min
271 frequency_indices_array(2,ifreq_block) = nfreq
272 
273 if (mpi_enreg%me == 0 ) then
274   write(io_unit,10) "                                                                                 "
275   write(io_unit,10) "#================================================================================"
276   write(io_unit,10) "#                                                                                "
277   write(io_unit,10) "#                                                                                "
278   write(io_unit,10) "#  Final frequency_indices_array :                                               "
279   write(io_unit,10) "#                                                                                "
280   write(io_unit,10) "#  ifreq_block iw_ext_min     iw_ext_max                                         "
281   write(io_unit,10) "#================================================================================"
282 
283   do ifreq_block = 1 , number_of_frequency_blocks 
284 
285   write(io_unit,40) ifreq_block, frequency_indices_array(1,ifreq_block), frequency_indices_array(2,ifreq_block)
286 
287   end do
288 end if 
289 
290 
291 
292 
293 !=================================================
294 !
295 ! Set up the LanczosResolvents algorithm
296 !
297 !=================================================
298 
299 prec = .false.  ! let's not precondition for now
300 call setup_LanczosResolvents(kmax_analytic, prec)
301 
302 !=================================================
303 !
304 ! iterate on frequency blocks!
305 !
306 !=================================================
307 nvec = 1
308 
309 ABI_MALLOC(list_left_vectors, (npw_g,nvec))
310 ABI_MALLOC(seed_vector, (npw_g))
311 
312 ABI_MALLOC(psik_wrk,  (2,npw_k))
313 ABI_MALLOC(psikb_wrk, (2,npw_kb))
314 
315 ABI_MALLOC(psikg_wrk, (2,npw_g))
316 ABI_MALLOC(psikg,     (2,npw_g))
317 
318 list_AT_Lanczos(:,:) = cmplx_0
319 
320 ! Number of blocks of lanczos vectors
321 nbdblock_lanczos = lmax/blocksize
322 if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
323 
324 
325 do ifreq_block = 1, number_of_frequency_blocks
326 
327 !-------------------------------
328 ! Build the frequency arrays
329 !-------------------------------
330 iw_ext_min = frequency_indices_array(1,ifreq_block) 
331 iw_ext_max = frequency_indices_array(2,ifreq_block)
332 
333 
334 nz = iw_ext_max -iw_ext_min + 1
335 
336 ABI_MALLOC( list_z_P, (nz))
337 ABI_MALLOC( list_z_Q, (nz))
338 ABI_MALLOC(matrix_elements_resolvent, (nz,nvec))
339 
340 
341 external_omega = list_external_omega(iw_ext_min)
342 call get_projection_band_indices(external_omega,band_index_below, band_index_above)
343 
344 
345 do iw_ext = iw_ext_min, iw_ext_max
346 
347 list_z_P(iw_ext-iw_ext_min+1) = list_external_omega(iw_ext)+ model_parameter
348 list_z_Q(iw_ext-iw_ext_min+1) = list_external_omega(iw_ext)- model_parameter
349 
350 end do
351 
352 !-----------------------------------------
353 ! Compute with shift Lanczos, using blocks
354 !-----------------------------------------
355 do iblk = 1, nbdblock_lanczos
356 ! which l is on this row of FFT processors?
357 l = (iblk-1)*blocksize + mpi_band_rank + 1 
358 
359 ! Change the configuration of the data
360 do mb =1, blocksize
361 lloc = (iblk-1)*blocksize+mb
362 if (lloc <= lmax) then 
363   psik_wrk(1,:)  = dble ( modified_Lbasis(:,lloc) )
364   psik_wrk(2,:)  = dimag( modified_Lbasis(:,lloc) )
365 else
366   psik_wrk(:,:)  = zero
367 end if
368 
369 psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:) 
370 
371 end do ! mb 
372 
373 ! change configuration of the data, from LA to FFT
374 call wf_block_distribute(psikb_wrk,  psikg_wrk, 1) ! LA -> FFT
375 
376 
377 if (band_index_below == 0 .and. l <= lmax) then
378   ! There are no DFT states with energy below W!
379   ! PW is thus zero, and QW = I.
380 
381   seed_vector(:)         = cmplx_1*psikg_wrk(1,:) + cmplx_i*psikg_wrk(2,:)
382   list_left_vectors(:,1) = seed_vector(:)
383 
384   call compute_resolvent_column_shift_lanczos(nz, list_z_Q, nvec, list_left_vectors, seed_vector, &
385   &                                                                   matrix_elements_resolvent)
386 
387   list_AT_Lanczos(iw_ext_min:iw_ext_max, l) = 0.5_dp*model_parameter*matrix_elements_resolvent(:,1)
388 
389 
390 else if ( l <= lmax ) then
391 
392   !----------------------------------------------------------------------------------
393   !
394   ! CAREFUL HERE! PW + QW != I. If the external frequency is exactly equal to a 
395   !               DFT eigenvalue, then the projections PW and QW both project OUT
396   !               of the subspace corresponding to this energy! 
397   !
398   !----------------------------------------------------------------------------------
399 
400   !-------------------
401   ! Treat first PW
402   !-------------------
403   psikg(:,:) = psikg_wrk(:,:)
404   call pc_k_valence_kernel(psikg,band_index_below)
405 
406   seed_vector = cmplx_1*psikg(1,:)+cmplx_i*psikg(2,:)
407 
408   list_left_vectors(:,1) = seed_vector(:)
409 
410   call compute_resolvent_column_shift_lanczos(nz, list_z_P, nvec, list_left_vectors, seed_vector, &
411   &                                                                   matrix_elements_resolvent)
412 
413   list_AT_Lanczos(iw_ext_min:iw_ext_max, l) = 0.5_dp*model_parameter*matrix_elements_resolvent(:,1)
414 
415 
416   !-------------------
417   ! Treat second QW
418   !-------------------
419 
420   psikg(:,:) = psikg_wrk(:,:)
421 
422   call pc_k_valence_kernel(psikg,band_index_above)
423 
424   psikg(:,:) = psikg_wrk(:,:)-psikg(:,:)
425 
426   seed_vector = cmplx_1*psikg(1,:)+cmplx_i*psikg(2,:)
427 
428   list_left_vectors(:,1) = seed_vector(:)
429 
430   call compute_resolvent_column_shift_lanczos(nz, list_z_Q, nvec, list_left_vectors, seed_vector, &
431   &                                                                   matrix_elements_resolvent)
432 
433   list_AT_Lanczos(iw_ext_min:iw_ext_max, l) = list_AT_Lanczos(iw_ext_min:iw_ext_max, l) +   &
434   0.5_dp*model_parameter*matrix_elements_resolvent(:,1)
435 
436 end if
437 
438 end do
439 
440 ABI_FREE( list_z_P)
441 ABI_FREE( list_z_Q)
442 ABI_FREE(matrix_elements_resolvent)
443 
444 end do
445 
446 
447 ! sum all results
448 call xmpi_sum(list_AT_Lanczos, mpi_enreg%comm_band,ierr) ! sum on all processors working on bands
449 
450 
451 if (mpi_enreg%me == 0 ) then
452   close(io_unit)
453 end if 
454 
455 ABI_FREE(list_left_vectors)
456 ABI_FREE(seed_vector)
457 ABI_FREE(frequency_indices_array)
458 
459 
460 ABI_FREE(psik_wrk)
461 ABI_FREE(psikb_wrk)
462 
463 ABI_FREE(psikg_wrk)
464 ABI_FREE(psikg)
465 
466 
467 
468 
469 call cleanup_LanczosResolvents
470 
471 10 format(A)
472 15 format(A,I5)
473 17 format(A,1000F8.4)
474 20 format(I5,7X,ES24.16,3I5)
475 30 format(4(I5,10X))
476 40 format(3(I5,10x))
477 
478 end subroutine compute_AT_shift_Lanczos