TABLE OF CONTENTS
ABINIT/gwls_Projected_AT [ Modules ]
NAME
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 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 gwls_wf 36 use gwls_TimingLog 37 use gwls_hamiltonian 38 use gwls_lineqsolver 39 use gwls_GWlanczos 40 use gwls_LanczosBasis 41 use gwls_LanczosResolvents 42 43 use gwls_GWanalyticPart, only : get_projection_band_indices 44 ! abinit modules 45 use defs_basis 46 use m_profiling_abi 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