TABLE OF CONTENTS
ABINIT/m_gwls_Projected_AT [ Modules ]
NAME
m_gwls_Projected_AT
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2024 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