TABLE OF CONTENTS


ABINIT/m_gwls_Projected_AT [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_Projected_AT

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

m_hamiltonian/compute_AT_shift_Lanczos [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_AT_shift_Lanczos

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      cleanup_lanczosresolvents,compute_resolvent_column_shift_lanczos
      get_projection_band_indices,pc_k_valence_kernel,setup_lanczosresolvents
      wf_block_distribute,xmpi_sum

SOURCE

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