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