1! 2! Copyright (C) 2020-2020 Quantum ESPRESSO group 3! 4! This file is distributed under the terms of the GNU General Public 5! License. See the file `LICENSE' in the root directory of the 6! present distribution, or http://www.gnu.org/copyleft.gpl.txt . 7! 8!------------------------------------------------------------------------------ 9PROGRAM postahc 10!------------------------------------------------------------------------------ 11!! 12!! This program 13!! - Reads the electron-phonon quantities calculated by ph.x with the 14!! electron_phonon='ahc' option. 15!! - Calculate the phonon-induced electron self-energy in the full matrix 16!! form at a given temperature. 17!! 18!! Input data (namelist "input") is described in Doc/INPUT_POSTAHC. 19!! 20!------------------------------------------------------------------------------ 21 USE kinds, ONLY : DP 22 USE constants, ONLY : ry_to_kelvin, AMU_RY, RY_TO_CMM1 23 USE mp, ONLY : mp_bcast, mp_sum 24 USE mp_world, ONLY : world_comm 25 USE mp_global, ONLY : mp_startup, mp_global_end 26 USE io_global, ONLY : ionode_id, ionode, stdout 27 USE environment, ONLY : environment_start, environment_end 28 ! 29 IMPLICIT NONE 30 ! 31 INTEGER, PARAMETER :: nat_max = 1000 32 !! Max number of atoms. 33 ! 34 ! Input variables 35 ! 36 LOGICAL :: skip_upperfan 37 !! Skip calculation of upper Fan self-energy 38 LOGICAL :: skip_dw 39 !! Skip calculation of Debye-Waller self-energy 40 INTEGER :: nk 41 !! Number of k points 42 INTEGER :: nbnd 43 !! Number of bands in the NSCF calculation 44 INTEGER :: ahc_nbnd 45 !! Number of bands for which the electron self-energy is to be computed. 46 INTEGER :: ahc_nbndskip 47 !! Number of bands to exclude when computing the self-energy. The 48 !! self-energy is computed for ibnd from ahc_nbndskip + 1 49 !! to ahc_nbndskip + ahc_nbnd. 50 INTEGER :: nat 51 !! Number of atoms 52 INTEGER :: nq 53 !! Number of q points 54 CHARACTER(LEN=256) :: filename 55 !! Name of binary files 56 CHARACTER(LEN=256) :: ahc_dir 57 !! Directory where the output binary files for AHC e-ph coupling are located 58 CHARACTER(LEN=256) :: flvec 59 !! output file for normalized phonon displacements generated by matdyn.x 60 !! The normalized phonon displacements are the eigenvectors divided by the 61 !! square root of the mass then normalized. As such they are not orthogonal. 62 REAL(DP) :: eta 63 !! Infinitesimal to add in the denominator for self-energy, in Ry 64 REAL(DP) :: temp_kelvin 65 !! Temperature at which the electron self-energy is calculated, in Kelvins. 66 REAL(DP) :: efermi 67 !! Fermi level energy, in Ry 68 REAL(DP) :: amass_amu(nat_max) 69 !! Mass of atoms one for each atom (not type), in atomic mass unit. 70 ! 71 ! Local variables 72 ! 73 LOGICAL :: lgamma 74 !! .true. if q == Gamma 75 INTEGER :: ios 76 !! io status 77 INTEGER :: iat 78 !! Counter for atoms 79 INTEGER :: ibnd 80 !! Counter for bands 81 INTEGER :: jbnd 82 !! Counter for bands 83 INTEGER :: ik 84 !! Counter for k points 85 INTEGER :: iq 86 !! Counter for q points 87 INTEGER :: imode 88 !! Counter for modes 89 INTEGER :: iun 90 !! Unit for reading file 91 INTEGER :: recl 92 !! Record length for reading file 93 INTEGER :: count 94 !! Counter for degeneracy 95 INTEGER :: nmodes 96 !! Number of modes. 3 * nat 97 REAL(DP) :: temperature 98 !! temp_kelvin transformed from Kelvin to Ry 99 REAL(DP) :: unorm 100 !! Norm of u multiplied with amass 101 REAL(DP) :: rval 102 !! Temporary real variables 103 REAL(DP) :: omega_zero_cutoff = 1.d-4 104 !! Cutoff of phonon frequency. Modes with omega smaller than 105 !! omega_zero_cutoff is neglected. 106 REAL(DP) :: e_degen_cutoff = 2.d-5 107 !! degeneracy cutoff. Ignore couping between degenerate states with energy 108 !! difference less than e_degen_cutoff. 109 COMPLEX(DP) :: selfen_avg_temp(5) 110 !! Diagonal self-energy averaged over degenerate states 111 REAL(DP), ALLOCATABLE :: inv_omega(:) 112 !! (nmodes) 1 / omega 113 REAL(DP), ALLOCATABLE :: occph(:) 114 !! (nmodes) Bose-Einstein occupation of phonon 115 REAL(DP), ALLOCATABLE :: wtq(:) 116 !! (nq) Weight of q points. Set to 1/nq. 117 REAL(DP), ALLOCATABLE :: amass(:) 118 !! Mass of atoms in Ry 119 REAL(DP), ALLOCATABLE :: xq(:, :) 120 !! (3, nq) q point vectors in Cartesian basis 121 REAL(DP), ALLOCATABLE :: omega(:, :) 122 !! (nmodes, nq) Phonon frequency 123 REAL(DP), ALLOCATABLE :: etk(:, :) 124 !! (nbnd, nk) Energy at k 125 REAL(DP), ALLOCATABLE :: etk_all(:, :) 126 !! (ahc_nbnd, nk) Energy at k 127 COMPLEX(DP), ALLOCATABLE :: u(:, :, :) 128 !! (nmodes, nmodes, nq) Phonon modes 129 COMPLEX(DP), ALLOCATABLE :: ahc_dw(:, :, :, :, :) 130 !! Debye-Waller matrix element 131 COMPLEX(DP), ALLOCATABLE :: selfen_dw(:, :, :) 132 !! Debye-Waller self-energy 133 COMPLEX(DP), ALLOCATABLE :: selfen_upfan(:, :, :) 134 !! Upper Fan self-energy 135 COMPLEX(DP), ALLOCATABLE :: selfen_lofan(:, :, :) 136 !! Lower Fan self-energy 137 COMPLEX(DP), ALLOCATABLE :: selfen_fan(:, :, :) 138 !! Fan self-energy (lower + upper) 139 COMPLEX(DP), ALLOCATABLE :: selfen_tot(:, :, :) 140 !! Total self-energy 141 COMPLEX(DP), ALLOCATABLE :: selfen_diag(:, :) 142 !! Diagonal self-energy 143 COMPLEX(DP), ALLOCATABLE :: selfen_diag_avg(:, :) 144 !! Diagonal self-energy averaged over degenerate states 145 ! 146 INTEGER, EXTERNAL :: find_free_unit 147 CHARACTER(LEN=256), EXTERNAL :: trimcheck 148 CHARACTER(len=6), EXTERNAL :: int_to_char 149 REAL(DP), EXTERNAL :: wgauss 150 ! 151 ! --------------------------------------------------------------------------- 152 ! 153 NAMELIST / input / skip_upperfan, skip_dw, nk, nbnd, nat, nq, ahc_nbnd, & 154 ahc_nbndskip, ahc_dir, flvec, eta, temp_kelvin, efermi, amass_amu 155 ! 156 CALL mp_startup() 157 CALL environment_start('POSTAHC') 158 ! 159 ! --------------------------------------------------------------------------- 160 ! 161 ! Reading input 162 ! 163 IF (ionode) CALL input_from_file() 164 ! 165 ! Default values of input arguments 166 ! 167 skip_upperfan = .FALSE. 168 skip_dw = .FALSE. 169 nk = -1 170 nbnd = -1 171 nat = -1 172 nq = -1 173 ahc_nbnd = -1 174 ahc_nbndskip = 0 175 ahc_dir = ' ' 176 flvec = ' ' 177 eta = -9999.d0 178 temp_kelvin = -9999.d0 179 efermi = -9999.d0 180 amass_amu(:) = -9999.d0 181 ! 182 ! Read input file 183 ! 184 IF (ionode) READ (5, input, IOSTAT=ios) 185 CALL mp_bcast(ios, ionode_id, world_comm) 186 CALL errore('postahc','error reading input namelist', ABS(ios)) 187 ! 188 ! Broadcast input arguments 189 ! 190 CALL mp_bcast(skip_upperfan, ionode_id, world_comm) 191 CALL mp_bcast(skip_dw, ionode_id, world_comm) 192 CALL mp_bcast(nk, ionode_id, world_comm) 193 CALL mp_bcast(nbnd, ionode_id, world_comm) 194 CALL mp_bcast(nat, ionode_id, world_comm) 195 CALL mp_bcast(nq, ionode_id, world_comm) 196 CALL mp_bcast(ahc_nbnd, ionode_id, world_comm) 197 CALL mp_bcast(ahc_nbndskip, ionode_id, world_comm) 198 CALL mp_bcast(ahc_dir, ionode_id, world_comm) 199 CALL mp_bcast(flvec, ionode_id, world_comm) 200 CALL mp_bcast(eta, ionode_id, world_comm) 201 CALL mp_bcast(temp_kelvin, ionode_id, world_comm) 202 CALL mp_bcast(efermi, ionode_id, world_comm) 203 CALL mp_bcast(amass_amu, ionode_id, world_comm) 204 ! 205 ! Check input argument validity 206 ! 207 IF (nk < 0) CALL errore('postahc', 'nk must be specified', 1) 208 IF (nbnd < 0) CALL errore('postahc', 'nbnd must be specified', 1) 209 IF (nat < 0) CALL errore('postahc', 'nat must be specified', 1) 210 IF (nq < 0) CALL errore('postahc', 'nq must be specified', 1) 211 IF (ahc_nbnd < 0) CALL errore('postahc', 'ahc_nbnd must be specified', 1) 212 IF (ahc_dir == '') CALL errore('postahc', 'ahc_dir must be specified', 1) 213 IF (flvec == '') CALL errore('postahc', 'flvec must be specified', 1) 214 IF (eta < -9990.d0) CALL errore('postahc', 'eta must be specified', 1) 215 IF (efermi < -9990.d0) CALL errore('postahc', 'efermi must be specified', 1) 216 IF (temp_kelvin < -9990.d0) CALL errore('postahc', & 217 'temp_kelvin must be specified', 1) 218 DO iat = 1, nat 219 IF (amass_amu(iat) < -9990.d0) CALL errore('postahc', & 220 'amass_amu(iat) must be specified for iat = 1 to nat', 1) 221 ENDDO 222 ! 223 ! Parallelization not implemented. Other processors goto the end of the code. 224 ! 225 IF (.NOT. ionode) GOTO 1001 226 ! 227 ! Setup local variables, unit converstion 228 ! 229 ahc_dir = trimcheck(ahc_dir) 230 nmodes = 3 * nat 231 temperature = temp_kelvin / ry_to_kelvin 232 ! 233 ALLOCATE(amass(nmodes)) 234 ALLOCATE(wtq(nq)) 235 ALLOCATE(xq(3, nq)) 236 ALLOCATE(omega(nmodes, nq)) 237 ALLOCATE(u(nmodes, nmodes, nq)) 238 ALLOCATE(occph(nmodes)) 239 ALLOCATE(inv_omega(nmodes)) 240 ALLOCATE(etk_all(nbnd, nk)) 241 ALLOCATE(etk(ahc_nbnd, nk)) 242 ALLOCATE(ahc_dw(ahc_nbnd, ahc_nbnd, nmodes, 3, nk)) 243 ALLOCATE(selfen_dw(ahc_nbnd, ahc_nbnd, nk)) 244 ALLOCATE(selfen_upfan(ahc_nbnd, ahc_nbnd, nk)) 245 ALLOCATE(selfen_lofan(ahc_nbnd, ahc_nbnd, nk)) 246 ALLOCATE(selfen_fan(ahc_nbnd, ahc_nbnd, nk)) 247 ALLOCATE(selfen_tot(ahc_nbnd, ahc_nbnd, nk)) 248 ALLOCATE(selfen_diag(ahc_nbnd, 5)) 249 ALLOCATE(selfen_diag_avg(ahc_nbnd, 5)) 250 ! 251 etk_all = 0.d0 252 etk = 0.d0 253 selfen_dw = (0.d0, 0.d0) 254 selfen_upfan = (0.d0, 0.d0) 255 selfen_lofan = (0.d0, 0.d0) 256 selfen_fan = (0.d0, 0.d0) 257 selfen_tot = (0.d0, 0.d0) 258 ! 259 DO iat = 1, nat 260 amass((iat-1) * 3 + 1:(iat-1) * 3 + 3) = amass_amu(iat) * AMU_RY 261 ENDDO 262 wtq = 1.d0 / REAL(nq, KIND=DP) 263 ! 264 ! --------------------------------------------------------------------------- 265 ! 266 ! Read flvec file 267 ! 268 CALL read_flvec(flvec, nq, nmodes, xq, omega, u) 269 ! 270 omega = omega / RY_TO_CMM1 271 ! 272 DO iq = 1, nq 273 DO imode = 1, nmodes 274 unorm = SUM(CONJG(u(:, imode, iq)) * u(:, imode, iq) * amass(:)) 275 u(:, imode, iq) = u(:, imode, iq) / SQRT(unorm) 276 ENDDO 277 ENDDO 278 ! 279 ! --------------------------------------------------------------------------- 280 ! 281 ! Read ahc_dw which does not depend on iq. 282 ! 283 iun = find_free_unit() 284 ! 285 filename = TRIM(ahc_dir) // 'ahc_dw.bin' 286 ! 287 INQUIRE(IOLENGTH=recl) ahc_dw 288 OPEN(iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', & 289 ACCESS='DIRECT', RECL=recl, IOSTAT=ios) 290 IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1) 291 READ(iun, REC=1) ahc_dw 292 CLOSE(iun) 293 ! 294 ! Read ahc_etk_iq 295 ! 296 DO iq = 1, nq 297 filename = TRIM(ahc_dir) // 'ahc_etk_iq' // TRIM(int_to_char(iq)) // '.bin' 298 ! 299 INQUIRE(IOLENGTH=recl) etk_all 300 OPEN(iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', & 301 ACCESS='DIRECT', RECL=recl, IOSTAT=ios) 302 IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1) 303 READ(iun, REC=1) etk_all 304 CLOSE(iun) 305 ENDDO 306 ! 307 ! Skip ahc_nbndskip bands in etk 308 ! 309 etk = etk_all(ahc_nbndskip+1:ahc_nbndskip+ahc_nbnd, :) 310 ! 311 ! --------------------------------------------------------------------------- 312 ! 313 ! Loop over q points and calculate self-energy 314 ! 315 WRITE(stdout, *) 316 WRITE(stdout,'(5x,a)') 'Calculating electron self-energy. Loop over q points' 317 ! 318 DO iq = 1, nq 319 ! 320 WRITE(stdout, '(i8)', ADVANCE='no') iq 321 IF(MOD(iq, 10) == 0 ) WRITE(stdout,*) 322 FLUSH(stdout) 323 ! 324 lgamma = .FALSE. 325 IF ( ALL( ABS(xq(:, iq)) < 1.d-4 ) ) lgamma = .TRUE. 326 ! 327 ! Set Bose-Einstein occupation of phonons and set inv_omega 328 ! 329 DO imode = 1, nmodes 330 IF (omega(imode, iq) < omega_zero_cutoff) THEN 331 inv_omega(imode) = 0.d0 332 occph(imode) = 0.d0 333 ELSE 334 inv_omega(imode) = 1.d0 / omega(imode, iq) 335 IF (temperature < 1.d-4) THEN 336 occph = 0.d0 337 ELSE 338 rval = wgauss( omega(imode, iq) / temperature, -99 ) 339 occph(imode) = 1.d0 / (4.d0 * rval - 2.d0) - 0.5d0 340 ENDIF 341 ENDIF 342 ENDDO 343 ! 344 IF (.NOT. skip_dw) CALL calc_debye_waller(iq, selfen_dw) 345 ! 346 IF (.NOT. skip_upperfan) CALL calc_upper_fan(iq, selfen_upfan) 347 ! 348 CALL calc_lower_fan(iq, selfen_lofan) 349 ! 350 ENDDO ! iq 351 ! 352 ! --------------------------------------------------------------------------- 353 ! 354 selfen_fan = selfen_lofan + selfen_upfan 355 selfen_tot = selfen_fan + selfen_dw 356 ! 357 ! Write self-energy to stdout 358 ! 359 WRITE(stdout, *) 360 WRITE(stdout, *) 361 IF (skip_dw) THEN 362 WRITE(stdout, '(5x,a)') 'Skip Debye-Waller: Debye-Waller self-energy & 363 &is set to zero' 364 ENDIF 365 IF (skip_upperfan) THEN 366 WRITE(stdout, '(5x,a)') 'Skip upper Fan: upper Fan self-energy is set & 367 &to zero' 368 ENDIF 369 ! 370 WRITE(stdout, *) 371 WRITE(stdout, '(5x,a)') 'Real part of diagonal electron self-energy in Ry' 372 WRITE(stdout, '(5x,a)') 'Self-energy of degenerate states are averaged.' 373 WRITE(stdout, '(5x,a)') 'Total_Fan = Upper_Fan + Lower_Fan' 374 WRITE(stdout, '(5x,a)') 'Total = Total_Fan + DW' 375 WRITE(stdout, *) 376 WRITE(stdout, '(5x,a)') 'Begin postahc output' 377 WRITE(stdout, '(5x,a)') ' ik ibnd Total DW Total_Fan& 378 & Upper_Fan Lower_Fan' 379 DO ik = 1, nk 380 DO ibnd = 1, ahc_nbnd 381 selfen_diag(ibnd, 1) = selfen_tot(ibnd, ibnd, ik) 382 selfen_diag(ibnd, 2) = selfen_dw(ibnd, ibnd, ik) 383 selfen_diag(ibnd, 3) = selfen_fan(ibnd, ibnd, ik) 384 selfen_diag(ibnd, 4) = selfen_upfan(ibnd, ibnd, ik) 385 selfen_diag(ibnd, 5) = selfen_lofan(ibnd, ibnd, ik) 386 ENDDO 387 ! 388 ! Average over degenerate states 389 ! 390 DO ibnd = 1, ahc_nbnd 391 ! 392 selfen_avg_temp = (0.d0, 0.d0) 393 count = 0 394 ! 395 DO jbnd = 1, ahc_nbnd 396 ! 397 IF (ABS(etk(ibnd, ik) - etk(jbnd, ik)) < e_degen_cutoff) THEN 398 count = count + 1 399 selfen_avg_temp = selfen_avg_temp + selfen_diag(jbnd, :) 400 ENDIF 401 ! 402 ENDDO 403 ! 404 selfen_diag_avg(ibnd, :) = selfen_avg_temp / REAL(count, DP) 405 ! 406 ENDDO ! ibnd 407 ! 408 ! Write averaged self-energy 409 ! 410 DO ibnd = 1, ahc_nbnd 411 WRITE(stdout, '(5x, 2I6, 5F12.7)') ik, ibnd, REAL(selfen_diag_avg(ibnd, :)) 412 ENDDO 413 ! 414 ENDDO ! ik 415 ! 416 WRITE(stdout, '(5x,a)') 'End postahc output' 417 WRITE(stdout, *) 418 WRITE(stdout, '(5x,a)') 'Full off-diagonal complex self-energy & 419 &matrix is written in files' 420 WRITE(stdout, '(5x,a)') 'selfen_real.dat and selfen_imag.dat. & 421 &These data can differ from' 422 WRITE(stdout, '(5x,a)') 'the output above because the self-energy & 423 &of degenerate states are' 424 WRITE(stdout, '(5x,a)') 'NOT averaged in the selfen_*.dat output.' 425 ! 426 ! Write self-energy to selfen.dat 427 ! 428 iun = find_free_unit() 429 ! 430 OPEN(iun, FILE='selfen_real.dat', FORM='FORMATTED') 431 WRITE(iun, '(a)') '# Real part of diagonal electron self-energy & 432 &Re[sigma(ibnd, jbnd, ik)] in Ry' 433 WRITE(iun, '(a)') '# ik ibnd jbnd Total DW Total_Fan& 434 & Upper_Fan Lower_Fan' 435 DO ik = 1, nk 436 DO jbnd = 1, ahc_nbnd 437 DO ibnd = 1, ahc_nbnd 438 WRITE(iun, '(3I6, 5F12.7)') ik, ibnd, jbnd, & 439 REAL(selfen_tot(ibnd, jbnd, ik)), REAL(selfen_dw(ibnd, jbnd, ik)), & 440 REAL(selfen_fan(ibnd, jbnd, ik)), REAL(selfen_upfan(ibnd, jbnd, ik)), & 441 REAL(selfen_lofan(ibnd, jbnd, ik)) 442 ENDDO 443 ENDDO 444 ENDDO 445 CLOSE(iun) 446 ! 447 OPEN(iun, FILE='selfen_imag.dat', FORM='FORMATTED') 448 WRITE(iun, '(a)') '# Imaginary part of diagonal electron self-energy & 449 &Im[sigma(ibnd, jbnd, ik)] in Ry' 450 WRITE(iun, '(a)') '# ik ibnd jbnd Total DW Total_Fan& 451 & Upper_Fan Lower_Fan' 452 DO ik = 1, nk 453 DO jbnd = 1, ahc_nbnd 454 DO ibnd = 1, ahc_nbnd 455 WRITE(iun, '(3I6, 5F12.7)') ik, ibnd, jbnd, & 456 AIMAG(selfen_tot(ibnd, jbnd, ik)), AIMAG(selfen_dw(ibnd, jbnd, ik)), & 457 AIMAG(selfen_fan(ibnd, jbnd, ik)), AIMAG(selfen_upfan(ibnd, jbnd, ik)), & 458 AIMAG(selfen_lofan(ibnd, jbnd, ik)) 459 ENDDO 460 ENDDO 461 ENDDO 462 CLOSE(iun) 463 ! 464 ! --------------------------------------------------------------------------- 465 ! 466 DEALLOCATE(amass) 467 DEALLOCATE(wtq) 468 DEALLOCATE(xq) 469 DEALLOCATE(omega) 470 DEALLOCATE(u) 471 DEALLOCATE(occph) 472 DEALLOCATE(inv_omega) 473 DEALLOCATE(etk_all) 474 DEALLOCATE(etk) 475 DEALLOCATE(ahc_dw) 476 DEALLOCATE(selfen_dw) 477 DEALLOCATE(selfen_upfan) 478 DEALLOCATE(selfen_lofan) 479 DEALLOCATE(selfen_fan) 480 DEALLOCATE(selfen_tot) 481 DEALLOCATE(selfen_diag) 482 DEALLOCATE(selfen_diag_avg) 483 ! 484 ! --------------------------------------------------------------------------- 485 ! 486 WRITE(stdout, *) 487 WRITE(stdout, *) 488 CALL print_clock('debye_waller') 489 CALL print_clock('lower_fan') 490 CALL print_clock('upper_fan') 491 ! 4921001 CONTINUE 493 ! 494 CALL environment_end('POSTAHC') 495 CALL mp_global_end() 496 ! 497CONTAINS 498!------------------------------------------------------------------------------ 499SUBROUTINE calc_debye_waller(iq, selfen_dw) 500 !---------------------------------------------------------------------------- 501 !! 502 !! Compute Debye-Waller self-energy at iq 503 !! 504 !! Implements Eq.(8) of PHonon/Doc/dfpt_self_energy.pdf 505 !! 506 !! Here, the "operator-generalized acoustic sum rule" is used to represent 507 !! Debye-Waller self-energy as a simple matrix element. 508 !! See Eq.(13) of the following reference: 509 !! Jae-Mo Lihm and Cheol-Hwan Park, Phys. Rev. B, 101, 121102(R) (2020). 510 !! 511 !---------------------------------------------------------------------------- 512 USE kinds, ONLY : DP 513 ! 514 INTEGER, INTENT(IN) :: iq 515 COMPLEX(DP), INTENT(INOUT) :: selfen_dw(ahc_nbnd, ahc_nbnd, nk) 516 ! 517 INTEGER :: imode, jmode, kmode 518 !! Counter for modes 519 INTEGER :: idir, jdir 520 !! Counter for directions 521 COMPLEX(DP), ALLOCATABLE :: selfen_dw_iq(:, :, :) 522 !! Debye-Waller self-energy at iq 523 COMPLEX(DP), ALLOCATABLE :: coeff_dw(:, :, :) 524 !! Coefficients for Debye-Waller 525 ! 526 CALL start_clock('debye_waller') 527 ! 528 ALLOCATE(selfen_dw_iq(ahc_nbnd, ahc_nbnd, nk)) 529 ALLOCATE(coeff_dw(nmodes, 3, nmodes)) 530 coeff_dw = (0.d0, 0.d0) 531 selfen_dw_iq = (0.d0, 0.d0) 532 ! 533 ! coeff_dw(3*iat + idir, jdir, kmode) 534 ! = 0.5 * ( CONJG(u(iat*3+idir, kmode)) * u(iat*3+jdir, kmode) ).real 535 ! * (occph(kmode) + 0.5) * inv_omega(kmode) 536 ! 537 DO kmode = 1, nmodes 538 DO jdir = 1, 3 539 DO idir = 1, 3 540 DO iat = 1, nat 541 imode = 3 * (iat - 1) + idir 542 jmode = 3 * (iat - 1) + jdir 543 coeff_dw(imode, jdir, kmode) = coeff_dw(imode, jdir, kmode) & 544 + 0.5d0 * CONJG(u(imode, kmode, iq)) * u(jmode, kmode, iq) & 545 * (occph(kmode) + 0.5d0) * inv_omega(kmode) 546 ENDDO 547 ENDDO 548 ENDDO 549 ENDDO 550 coeff_dw = REAL(coeff_dw, KIND=DP) 551 ! 552 DO kmode = 1, nmodes 553 DO jdir = 1, 3 554 DO imode = 1, nmodes 555 selfen_dw_iq = selfen_dw_iq + ahc_dw(:, :, imode, jdir, :) & 556 * coeff_dw(imode, jdir, kmode) 557 ENDDO 558 ENDDO 559 ENDDO 560 ! 561 selfen_dw = selfen_dw + selfen_dw_iq * wtq(iq) 562 ! 563 DEALLOCATE(coeff_dw) 564 DEALLOCATE(selfen_dw_iq) 565 ! 566 CALL stop_clock('debye_waller') 567 ! 568!------------------------------------------------------------------------------ 569END SUBROUTINE calc_debye_waller 570!------------------------------------------------------------------------------ 571! 572!------------------------------------------------------------------------------ 573SUBROUTINE calc_upper_fan(iq, selfen_upfan) 574 !---------------------------------------------------------------------------- 575 !! 576 !! Compute upper Fan (high-energy band contribution from Sternheimer 577 !! equation) self-energy at iq 578 !! 579 !! Implements Eq.(11) of PHonon/Doc/dfpt_self_energy.pdf 580 !! 581 !---------------------------------------------------------------------------- 582 USE kinds, ONLY : DP 583 ! 584 INTEGER, INTENT(IN) :: iq 585 COMPLEX(DP), INTENT(INOUT) :: selfen_upfan(ahc_nbnd, ahc_nbnd, nk) 586 ! 587 CHARACTER(LEN=256) :: filename 588 !! Name of ahc_upfan_iq*.bin file 589 INTEGER :: iun 590 !! Unit for reading file 591 INTEGER :: recl 592 !! Record length for reading file 593 INTEGER :: ibnd, jbnd 594 !! Counter for bands 595 INTEGER :: imode, jmode, kmode 596 !! Counter for modes 597 REAL(DP) :: coeff 598 !! real coefficient 599 COMPLEX(DP), ALLOCATABLE :: selfen_upfan_iq(:, :, :) 600 !! Upper Fan self-energy at iq 601 COMPLEX(DP), ALLOCATABLE :: ahc_upfan(:, :, :, :, :) 602 !! Upper Fan matrix element 603 COMPLEX(DP), ALLOCATABLE :: ahc_upfan_mode(:, :, :, :) 604 !! Upper Fan matrix element in mode basis 605 INTEGER, EXTERNAL :: find_free_unit 606 ! 607 CALL start_clock('upper_fan') 608 ! 609 ALLOCATE(selfen_upfan_iq(ahc_nbnd, ahc_nbnd, nk)) 610 ALLOCATE(ahc_upfan(ahc_nbnd, ahc_nbnd, nmodes, nmodes, nk)) 611 ALLOCATE(ahc_upfan_mode(ahc_nbnd, ahc_nbnd, nmodes, nk)) 612 selfen_upfan_iq = (0.d0, 0.d0) 613 ahc_upfan_mode = (0.d0, 0.d0) 614 ! 615 ! Read ahc_upfan_iq*.bin 616 ! 617 filename = TRIM(ahc_dir) // 'ahc_upfan_iq' // TRIM(int_to_char(iq)) // '.bin' 618 ! 619 iun = find_free_unit() 620 ! 621 INQUIRE(IOLENGTH=recl) ahc_upfan 622 OPEN(iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', & 623 ACCESS='DIRECT', RECL=recl, IOSTAT=ios) 624 IF (ios /= 0) CALL errore('postahc', 'Error reading ' // TRIM(filename), 1) 625 ! 626 READ(iun, REC=1) ahc_upfan 627 CLOSE(iun) 628 ! 629 ! rotate ahc_upfan from Cartesian to eigenmode basis 630 ! 631 DO ik = 1, nk 632 DO imode = 1, nmodes 633 DO kmode = 1, nmodes 634 DO jmode = 1, nmodes 635 ahc_upfan_mode(:, :, imode, ik) = ahc_upfan_mode(:, :, imode, ik) & 636 + ahc_upfan(:, :, jmode, kmode, ik) & 637 * CONJG(u(jmode, imode, iq)) * u(kmode, imode, iq) 638 ENDDO 639 ENDDO 640 ENDDO 641 ENDDO 642 ! 643 ! Compute selfen_upfan_iq 644 ! 645 DO ik = 1, nk 646 DO imode = 1, nmodes 647 ! 648 coeff = 0.5d0 * inv_omega(imode) * (occph(imode) + 0.5d0) 649 ! 650 DO jbnd = 1, ahc_nbnd 651 DO ibnd = 1, ahc_nbnd 652 selfen_upfan_iq(ibnd, jbnd, ik) = selfen_upfan_iq(ibnd, jbnd, ik) & 653 + ( ahc_upfan_mode(ibnd, jbnd, imode, ik) & 654 + CONJG(ahc_upfan_mode(jbnd, ibnd, imode, ik)) ) & 655 * coeff 656 ENDDO 657 ENDDO 658 ! 659 ENDDO 660 ENDDO 661 ! 662 selfen_upfan = selfen_upfan + selfen_upfan_iq * wtq(iq) 663 ! 664 DEALLOCATE(ahc_upfan) 665 DEALLOCATE(selfen_upfan_iq) 666 ! 667 CALL stop_clock('upper_fan') 668 ! 669!------------------------------------------------------------------------------ 670END SUBROUTINE calc_upper_fan 671!------------------------------------------------------------------------------ 672! 673!------------------------------------------------------------------------------ 674SUBROUTINE calc_lower_fan(iq, selfen_lofan) 675 !---------------------------------------------------------------------------- 676 !! 677 !! Compute lower Fan (low-energy band contribution) self-energy at iq 678 !! 679 !! Implements Eq.(10) of PHonon/Doc/dfpt_self_energy.pdf 680 !! 681 !---------------------------------------------------------------------------- 682 USE kinds, ONLY : DP 683 ! 684 INTEGER, INTENT(IN) :: iq 685 COMPLEX(DP), INTENT(INOUT) :: selfen_lofan(ahc_nbnd, ahc_nbnd, nk) 686 ! 687 CHARACTER(LEN=256) :: filename 688 !! Name of ahc_gkk_iq*.bin file 689 INTEGER :: iun 690 !! Unit for reading file 691 INTEGER :: recl 692 !! Record length for reading file 693 INTEGER :: ibq, ibk, jbk 694 !! Counter for bands 695 INTEGER :: imode, jmode 696 !! Counter for modes 697 REAL(DP) :: sign 698 !! coefficients 699 COMPLEX(DP) :: de1, de2 700 !! coefficients 701 REAL(DP), ALLOCATABLE :: etq(:, :) 702 !! (nbnd, nk) Energy at k+q 703 REAL(DP), ALLOCATABLE :: occq(:, :) 704 !! (nbnd, nk) Fermi-Dirac occupation of energy at k+q 705 COMPLEX(DP), ALLOCATABLE :: selfen_lofan_iq(:, :, :) 706 !! Lower Fan self-energy at iq 707 COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :) 708 !! electron-phonon matrix element 709 COMPLEX(DP), ALLOCATABLE :: ahc_gkk_mode(:, :, :, :) 710 !! electron-phonon matrix element in mode basis 711 COMPLEX(DP), ALLOCATABLE :: coeff(:, :, :, :) 712 !! coefficient for lower Fan self-energy 713 COMPLEX(DP), ALLOCATABLE :: coeff1(:, :, :, :) 714 !! coefficient for lower Fan self-energy 715 COMPLEX(DP), ALLOCATABLE :: coeff2(:, :, :, :) 716 !! coefficient for lower Fan self-energy 717 INTEGER, EXTERNAL :: find_free_unit 718 ! 719 CALL start_clock('lower_fan') 720 ! 721 ALLOCATE(selfen_lofan_iq(ahc_nbnd, ahc_nbnd, nk)) 722 ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nk)) 723 ALLOCATE(ahc_gkk_mode(nbnd, ahc_nbnd, nmodes, nk)) 724 ALLOCATE(coeff(nbnd, ahc_nbnd, nmodes, nk)) 725 ALLOCATE(coeff1(nbnd, ahc_nbnd, nmodes, nk)) 726 ALLOCATE(coeff2(nbnd, ahc_nbnd, nmodes, nk)) 727 ALLOCATE(etq(nbnd, nk)) 728 ALLOCATE(occq(nbnd, nk)) 729 ! 730 selfen_lofan_iq = (0.d0, 0.d0) 731 ahc_gkk = (0.d0, 0.d0) 732 ahc_gkk_mode = (0.d0, 0.d0) 733 coeff = (0.d0, 0.d0) 734 coeff1 = (0.d0, 0.d0) 735 coeff2 = (0.d0, 0.d0) 736 etq = 0.d0 737 occq = 0.d0 738 ! 739 ! Read files: ahc_etq, ahc_gkk 740 ! 741 iun = find_free_unit() 742 ! 743 filename = TRIM(ahc_dir) // 'ahc_gkk_iq' // TRIM(int_to_char(iq)) // '.bin' 744 ! 745 INQUIRE(IOLENGTH=recl) ahc_gkk 746 OPEN(iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', & 747 ACCESS='DIRECT', RECL=recl, IOSTAT=ios) 748 IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1) 749 READ(iun, REC=1) ahc_gkk 750 CLOSE(iun) 751 ! 752 filename = TRIM(ahc_dir) // 'ahc_etq_iq' // TRIM(int_to_char(iq)) // '.bin' 753 ! 754 INQUIRE(IOLENGTH=recl) etq 755 OPEN(iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', & 756 ACCESS='DIRECT', RECL=recl, IOSTAT=ios) 757 IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1) 758 READ(iun, REC=1) etq 759 CLOSE(iun) 760 ! 761 ! Fermi-Dirac occupation at k+q 762 ! 763 DO ik = 1, nk 764 DO ibnd = 1, nbnd 765 IF (temperature < 1.d-4) THEN 766 IF (etq(ibnd, ik) < efermi) THEN 767 occq(ibnd, ik) = 1.0d0 768 ELSEIF (etq(ibnd, ik) > efermi) THEN 769 occq(ibnd, ik) = 0.0d0 770 ELSE 771 occq(ibnd, ik) = 0.5d0 772 ENDIF 773 ELSE 774 occq(ibnd, ik) = wgauss( (efermi - etq(ibnd, ik)) / temperature, -99 ) 775 ENDIF 776 ENDDO 777 ENDDO 778 ! 779 ! rotate ahc_gkk from Cartesian to eigenmode basis 780 ! 781 DO ik = 1, nk 782 DO imode = 1, nmodes 783 DO jmode = 1, nmodes 784 ahc_gkk_mode(:, :, imode, ik) = ahc_gkk_mode(:, :, imode, ik) & 785 + ahc_gkk(:, :, jmode, ik) * u(jmode, imode, iq) 786 ENDDO 787 ENDDO 788 ENDDO 789 ! 790 ! Compute coefficients 791 ! 792 ! sign = +1 if etk(ibk, ik) > efermi 793 ! = -1 otherwise 794 ! 795 ! coeff1(ibq, ibk, imode, ik) 796 ! = ( 1 - occq(ibq, ik) + occph(imode) ) 797 ! / ( etk(ibk, ik) - etq(ibq, ik) - omega(imode) + 1j * eta * sign ) 798 ! 799 ! coeff2(ibq, ibk, imode, ik) 800 ! = ( occq(ibq, ik) + occph(imode) ) 801 ! / ( etk(ibk, ik) - etq(ibq, ik) + omega(imode) + 1j * eta * sign ) 802 ! 803 ! coeff = (coeff1 + coeff2) * 0.5 * inv_omega(imode) 804 ! 805 DO ik = 1, nk 806 DO imode = 1, nmodes 807 DO ibk = 1, ahc_nbnd 808 ! 809 sign = 1.d0 810 IF (etk(ibk, ik) < efermi) sign = -1.d0 811 ! 812 DO ibq = 1, nbnd 813 de1 = CMPLX(etk(ibk, ik) - etq(ibq, ik) - omega(imode, iq),& 814 eta * sign, KIND=DP) 815 coeff1(ibq, ibk, imode, ik) = (1.d0 - occq(ibq, ik) + occph(imode)) & 816 / de1 817 ! 818 de2 = CMPLX(etk(ibk, ik) - etq(ibq, ik) + omega(imode, iq),& 819 eta * sign, KIND=DP) 820 coeff2(ibq, ibk, imode, ik) = ( occq(ibq, ik) + occph(imode) ) & 821 / de2 822 ENDDO 823 ENDDO 824 ENDDO 825 ENDDO 826 ! 827 coeff = coeff1 + coeff2 828 ! 829 DO ik = 1, nk 830 DO imode = 1, nmodes 831 coeff(:, :, imode, ik) = coeff(:, :, imode, ik) * 0.5d0 * inv_omega(imode) 832 ENDDO 833 ENDDO 834 ! 835 ! Remove coupling with oneself 836 ! 837 IF (lgamma) THEN 838 DO ibk = 1, ahc_nbnd 839 coeff(ibk + ahc_nbndskip, ibk, :, :) = (0.d0, 0.d0) 840 ENDDO 841 ENDIF 842 ! 843 ! Remove coupling between degenerate states 844 ! 845 IF (lgamma) THEN 846 DO ik = 1, nk 847 DO ibk = 1, ahc_nbnd 848 DO ibq = 1, nbnd 849 IF ( ABS(etk(ibk, ik) - etq(ibq, ik)) < e_degen_cutoff ) THEN 850 coeff(ibq, ibk, :, ik) = (0.d0, 0.d0) 851 ENDIF 852 ENDDO 853 ENDDO 854 ENDDO 855 ENDIF 856 ! 857 ! Compute selfen_lofan_iq 858 ! 859 DO ik = 1, nk 860 DO imode = 1, nmodes 861 ! 862 DO jbk = 1, ahc_nbnd 863 DO ibk = 1, ahc_nbnd 864 DO ibq = 1, nbnd 865 ! 866 selfen_lofan_iq(ibk, jbk, ik) = selfen_lofan_iq(ibk, jbk, ik) & 867 + CONJG(ahc_gkk_mode(ibq, ibk, imode, ik)) & 868 * ahc_gkk_mode(ibq, jbk, imode, ik) & 869 * ( coeff(ibq, ibk, imode, ik) & 870 + coeff(ibq, jbk, imode, ik) ) * 0.5d0 871 ! 872 ENDDO 873 ENDDO 874 ENDDO 875 ! 876 ENDDO 877 ENDDO 878 ! 879 selfen_lofan = selfen_lofan + selfen_lofan_iq * wtq(iq) 880 ! 881 DEALLOCATE(coeff) 882 DEALLOCATE(coeff1) 883 DEALLOCATE(coeff2) 884 DEALLOCATE(ahc_gkk) 885 DEALLOCATE(ahc_gkk_mode) 886 DEALLOCATE(selfen_lofan_iq) 887 DEALLOCATE(etq) 888 DEALLOCATE(occq) 889 ! 890 CALL stop_clock('lower_fan') 891 ! 892!------------------------------------------------------------------------------ 893END SUBROUTINE calc_lower_fan 894!------------------------------------------------------------------------------ 895! 896!------------------------------------------------------------------------------ 897SUBROUTINE read_flvec(flvec, nq, nmodes, xq, omega, u) 898 !---------------------------------------------------------------------------- 899 !! 900 !! Read flvec (.modes) file generated by matdyn.x 901 !! 902 !! Output 903 !! - xq : list of q point vectors in Cartesian coordinate 904 !! - omega : phonon frequency in Ry 905 !! - u : phonon modes, normalized to satisfy u^dagger * M * u = identity 906 !! (Eq.(1) of PHonon/Doc/dfpt_self_energy.pdf) 907 !! 908 !---------------------------------------------------------------------------- 909 USE kinds, ONLY : DP 910 ! 911 CHARACTER(LEN=256), INTENT(IN) :: flvec 912 !! Name of the modes file 913 INTEGER, INTENT(IN) :: nq 914 !! Number of q points 915 INTEGER, INTENT(IN) :: nmodes 916 !! Number of modes 917 REAL(DP), INTENT(INOUT) :: xq(3, nq) 918 !! q point vectors 919 REAL(DP), INTENT(INOUT) :: omega(nmodes, nq) 920 !! Phonon frequency 921 COMPLEX(DP), INTENT(INOUT) :: u(nmodes, nmodes, nq) 922 !! Phonon modes 923 ! 924 INTEGER :: i 925 !! dummy variable for reading flvec 926 REAL(DP) :: omega_ 927 !! dummy variable for reading flvec 928 INTEGER :: iq 929 !! Counter for q points 930 INTEGER :: imode 931 !! Counter for modes 932 INTEGER :: iat 933 !! Counter for atoms 934 INTEGER :: nat 935 !! number of atoms. nmodes / 3 936 INTEGER :: iun 937 !! Unit for reading flvec 938 INTEGER :: ios 939 !! io status 940 ! 941 INTEGER, EXTERNAL :: find_free_unit 942 ! 943 nat = nmodes / 3 944 ! 945 iun = find_free_unit() 946 OPEN(UNIT=iun, FILE=TRIM(flvec), STATUS='OLD', FORM='FORMATTED', IOSTAT=ios) 947 IF (ios /= 0) CALL errore('postahc', & 948 'problem reading flvec file ' // TRIM(flvec), 1) 949 ! 950 DO iq = 1, nq 951 READ(iun, '(a)') 952 READ(iun, '(a)') 953 READ(iun, '(1x,6x,3f12.4)') (xq(i, iq), i=1, 3) 954 READ(iun, '(a)') 955 DO imode = 1, nmodes 956 READ(iun, 9010) i, omega_, omega(imode, iq) 957 DO iat = 1, nat 958 READ(iun, 9020) (u(i, imode, iq), i=(iat-1)*3+1, (iat-1)*3+3) 959 ENDDO 960 ENDDO 961 READ(iun, '(a)') 962 ENDDO 9639010 format(5x, 6x, i5, 3x, f15.6, 8x, f15.6, 7x) 9649020 format(1x, 1x, 3(f10.6, 1x, f10.6, 3x), 1x) 965 ! 966 CLOSE(iun, IOSTAT=ios) 967 IF (ios /= 0) CALL errore('postahc', & 968 'problem closing flvec file ' // TRIM(flvec), 1) 969 ! 970!------------------------------------------------------------------------------ 971END SUBROUTINE read_flvec 972!------------------------------------------------------------------------------ 973! 974!------------------------------------------------------------------------------ 975END PROGRAM postahc 976!------------------------------------------------------------------------------ 977