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