1! 2! Copyright (C) 2001-2018 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!----------------------------------------------------------------------- 10SUBROUTINE elphon() 11 !----------------------------------------------------------------------- 12 ! 13 ! Electron-phonon calculation from data saved in fildvscf 14 ! 15 USE kinds, ONLY : DP 16 USE constants, ONLY : amu_ry, RY_TO_THZ, RY_TO_CMM1 17 USE cell_base, ONLY : celldm, omega, ibrav, at, bg 18 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, amass 19 USE gvecs, ONLY: doublegrid 20 USE fft_base, ONLY : dfftp, dffts 21 USE fft_interfaces, ONLY : fft_interpolate 22 USE noncollin_module, ONLY : nspin_mag, noncolin, m_loc 23 USE lsda_mod, ONLY : nspin 24 USE uspp, ONLY: okvan 25 USE paw_variables, ONLY : okpaw 26 USE el_phon, ONLY : done_elph 27 USE dynmat, ONLY : dyn, w2 28 USE modes, ONLY : npert, nirr, u 29 USE uspp_param, ONLY : nhm 30 USE control_ph, ONLY : trans, xmldyn 31 USE output, ONLY : fildyn,fildvscf 32 USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, & 33 read_dyn_mat, read_dyn_mat_tail 34 USE units_ph, ONLY : iudyn, lrdrho, iudvscf, iuint3paw, lint3paw 35 USE dfile_star, ONLY : dvscf_star 36 USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, root_bgrp 37 USE mp, ONLY : mp_bcast 38 USE io_global, ONLY: stdout 39 USE lrus, ONLY : int3, int3_nc, int3_paw 40 USE qpoint, ONLY : xq 41 USE dvscf_interpolate, ONLY : ldvscf_interpolate, dvscf_r2q 42 USE ahc, ONLY : elph_ahc 43 ! 44 IMPLICIT NONE 45 ! 46 INTEGER :: irr, imode0, ipert, is, npe 47 ! counter on the representations 48 ! counter on the modes 49 ! the change of Vscf due to perturbations 50 INTEGER :: i,j 51 COMPLEX(DP), ALLOCATABLE :: dvscfin_all(:, :, :) 52 !! dvscfin for all modes. Used when doing dvscf_r2q interpolation. 53 COMPLEX(DP), POINTER :: dvscfin(:,:,:), dvscfins (:,:,:) 54 COMPLEX(DP), allocatable :: phip (:, :, :, :) 55 56 INTEGER :: ntyp_, nat_, ibrav_, nspin_mag_, mu, nu, na, nb, nta, ntb, nqs_ 57 REAL(DP) :: celldm_(6), w1 58 CHARACTER(LEN=3) :: atm(ntyp) 59 60 CALL start_clock ('elphon') 61 62 if(dvscf_star%basis.eq.'cartesian') then 63 write(stdout,*) 'Setting patterns to identity' 64 u=(0.d0,0.d0) 65 do irr=1,3*nat 66 u(irr,irr)=(1.d0,0.d0) 67 enddo 68 endif 69 ! 70 ! If ldvscf_interpolate, use Fourier interpolation instead of reading dVscf 71 ! 72 IF (ldvscf_interpolate) THEN 73 ! 74 WRITE (6, '(5x,a)') "Fourier interpolating dVscf" 75 ALLOCATE(dvscfin_all(dfftp%nnr, nspin_mag, 3 * nat)) 76 CALL dvscf_r2q(xq, u, dvscfin_all) 77 ! 78 ELSE 79 WRITE (6, '(5x,a)') "Reading dVscf from file "//trim(fildvscf) 80 ENDIF 81 ! 82 ! read Delta Vscf and calculate electron-phonon coefficients 83 ! 84 imode0 = 0 85 DO irr = 1, nirr 86 npe=npert(irr) 87 ALLOCATE (dvscfin (dfftp%nnr, nspin_mag , npe) ) 88 IF (okvan) THEN 89 ALLOCATE (int3 ( nhm, nhm, nat, nspin_mag, npe)) 90 IF (okpaw) ALLOCATE (int3_paw (nhm, nhm, nat, nspin_mag, npe)) 91 IF (noncolin) ALLOCATE(int3_nc( nhm, nhm, nat, nspin, npe)) 92 ENDIF 93 DO ipert = 1, npe 94 IF (ldvscf_interpolate) THEN 95 dvscfin(:, :, ipert) = dvscfin_all(:, :, imode0 + ipert) 96 ELSE 97 CALL davcio_drho ( dvscfin(1,1,ipert), lrdrho, iudvscf, & 98 imode0 + ipert, -1 ) 99 ENDIF 100 IF (okpaw .AND. me_bgrp==0) & 101 CALL davcio( int3_paw(:,:,:,:,ipert), lint3paw, & 102 iuint3paw, imode0 + ipert, - 1 ) 103 END DO 104 IF (okpaw) CALL mp_bcast(int3_paw, root_bgrp, intra_bgrp_comm) 105 IF (doublegrid) THEN 106 ALLOCATE (dvscfins (dffts%nnr, nspin_mag , npert(irr)) ) 107 DO is = 1, nspin_mag 108 DO ipert = 1, npe 109 CALL fft_interpolate (dfftp, dvscfin(:,is,ipert), dffts, dvscfins(:,is,ipert)) 110 ENDDO 111 ENDDO 112 ELSE 113 dvscfins => dvscfin 114 ENDIF 115 CALL newdq (dvscfin, npert(irr)) 116 CALL elphel (irr, npert (irr), imode0, dvscfins) 117 ! 118 imode0 = imode0 + npe 119 IF (doublegrid) DEALLOCATE (dvscfins) 120 DEALLOCATE (dvscfin) 121 IF (okvan) THEN 122 DEALLOCATE (int3) 123 IF (okpaw) DEALLOCATE (int3_paw) 124 IF (noncolin) DEALLOCATE(int3_nc) 125 ENDIF 126 ENDDO 127 ! 128 IF (ldvscf_interpolate) DEALLOCATE(dvscfin_all) 129 ! 130 ! In AHC calculation, we do not need the dynamical matrix. So return here. 131 IF (elph_ahc) THEN 132 CALL stop_clock('elphon') 133 RETURN 134 ENDIF 135 ! 136 ! now read the eigenvalues and eigenvectors of the dynamical matrix 137 ! calculated in a previous run 138 ! 139 IF (.NOT.trans) THEN 140 IF (.NOT. xmldyn) THEN 141 WRITE (6, '(5x,a)') "Reading dynamics matrix from file "//trim(fildyn) 142 CALL readmat (iudyn, ibrav, celldm, nat, ntyp, & 143 ityp, omega, amass, tau, xq, w2, dyn) 144 ELSE 145 allocate( phip(3,3,nat,nat) ) 146 CALL read_dyn_mat_param(fildyn, ntyp_, nat_) 147 IF ( ntyp_ /= ntyp .OR. nat_ /= nat ) & 148 CALL errore('elphon','uncorrect nat or ntyp',1) 149 150 CALL read_dyn_mat_header(ntyp, nat, ibrav_, nspin_mag_, & 151 celldm_, at, bg, omega, atm, amass, tau, ityp, & 152 m_loc, nqs_) 153 154 IF (ibrav_.NE.ibrav .OR. ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5 & 155 .OR. (nspin_mag_ /= nspin_mag ) ) CALL errore ('elphon', & 156 'inconsistent data', 1) 157 158 CALL read_dyn_mat(nat,1,xq,phip) 159 ! 160 ! Diagonalize the dynamical matrix 161 ! 162 163 164 DO i=1,3 165 do na=1,nat 166 nta = ityp (na) 167 mu=3*(na-1)+i 168 do j=1,3 169 do nb=1,nat 170 nu=3*(nb-1)+j 171 ntb = ityp (nb) 172 dyn (mu, nu) = phip (i, j, na, nb) / & 173 sqrt( amass(nta)*amass(ntb))/amu_ry 174 enddo 175 enddo 176 enddo 177 enddo 178 179 ! 180 CALL cdiagh (3 * nat, dyn, 3 * nat, w2, dyn) 181 ! 182 ! divide by sqrt(mass) to get displacements 183 ! 184 DO nu = 1, 3 * nat 185 DO mu = 1, 3 * nat 186 na = (mu - 1) / 3 + 1 187 dyn (mu, nu) = dyn (mu, nu) / SQRT ( amu_ry * amass (ityp (na) ) ) 188 ENDDO 189 ENDDO 190 191 CALL read_dyn_mat_tail(nat) 192 193 deallocate( phip ) 194 ENDIF 195 ! 196 ! Write phonon frequency to stdout 197 ! 198 WRITE( stdout, 8000) (xq (i), i = 1, 3) 199 ! 200 DO nu = 1, 3 * nat 201 w1 = SQRT( ABS( w2(nu) ) ) 202 if (w2(nu) < 0.d0) w1 = - w1 203 WRITE( stdout, 8010) nu, w1 * RY_TO_THZ, w1 * RY_TO_CMM1 204 ENDDO 205 ! 206 WRITE( stdout, '(1x,74("*"))') 207 ! 208 ENDIF ! .NOT. trans 209 ! 210 CALL stop_clock ('elphon') 211 ! 2128000 FORMAT(/,5x,'Diagonalizing the dynamical matrix', & 213 & //,5x,'q = ( ',3f14.9,' ) ',//,1x,74('*')) 2148010 FORMAT (5x,'freq (',i5,') =',f15.6,' [THz] =',f15.6,' [cm-1]') 215 ! 216 RETURN 217END SUBROUTINE elphon 218! 219!----------------------------------------------------------------------- 220SUBROUTINE readmat (iudyn, ibrav, celldm, nat, ntyp, ityp, omega, & 221 amass, tau, q, w2, dyn) 222 !----------------------------------------------------------------------- 223 ! 224 USE kinds, ONLY : DP 225 USE constants, ONLY : amu_ry 226 IMPLICIT NONE 227 ! Input 228 INTEGER :: iudyn, ibrav, nat, ntyp, ityp (nat) 229 REAL(DP) :: celldm (6), amass (ntyp), tau (3, nat), q (3), & 230 omega 231 ! output 232 REAL(DP) :: w2 (3 * nat) 233 COMPLEX(DP) :: dyn (3 * nat, 3 * nat) 234 ! local (control variables) 235 INTEGER :: ntyp_, nat_, ibrav_, ityp_ 236 REAL(DP) :: celldm_ (6), amass_, tau_ (3), q_ (3) 237 ! local 238 REAL(DP) :: dynr (2, 3, nat, 3, nat) 239 CHARACTER(len=80) :: line 240 CHARACTER(len=3) :: atm 241 INTEGER :: nt, na, nb, naa, nbb, nu, mu, i, j 242 ! 243 ! 244 REWIND (iudyn) 245 READ (iudyn, '(a)') line 246 READ (iudyn, '(a)') line 247 READ (iudyn, * ) ntyp_, nat_, ibrav_, celldm_ 248 IF ( ntyp.NE.ntyp_ .OR. nat.NE.nat_ .OR.ibrav_.NE.ibrav .OR. & 249 ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5) & 250 CALL errore ('readmat', 'inconsistent data', 1) 251 IF ( ibrav_ == 0 ) THEN 252 READ (iudyn, '(a)') line 253 READ (iudyn, '(a)') line 254 READ (iudyn, '(a)') line 255 READ (iudyn, '(a)') line 256 END IF 257 DO nt = 1, ntyp 258 READ (iudyn, * ) i, atm, amass_ 259 IF ( nt.NE.i .OR. ABS (amass_ - amu_ry*amass (nt) ) > 1.0d-5) & 260 CALL errore ( 'readmat', 'inconsistent data', 1 + nt) 261 ENDDO 262 DO na = 1, nat 263 READ (iudyn, * ) i, ityp_, tau_ 264 IF (na.NE.i.OR.ityp_.NE.ityp (na) ) CALL errore ('readmat', & 265 'inconsistent data', 10 + na) 266 ENDDO 267 READ (iudyn, '(a)') line 268 READ (iudyn, '(a)') line 269 READ (iudyn, '(a)') line 270 READ (iudyn, '(a)') line 271 READ (line (11:80), * ) (q_ (i), i = 1, 3) 272 READ (iudyn, '(a)') line 273 DO na = 1, nat 274 DO nb = 1, nat 275 READ (iudyn, * ) naa, nbb 276 IF (na.NE.naa.OR.nb.NE.nbb) CALL errore ('readmat', 'error reading & 277 &file', nb) 278 READ (iudyn, * ) ( (dynr (1, i, na, j, nb), dynr (2, i, na, j, nb) & 279 , j = 1, 3), i = 1, 3) 280 ENDDO 281 ENDDO 282 ! 283 ! divide the dynamical matrix by the (input) masses (in amu) 284 ! 285 DO nb = 1, nat 286 DO j = 1, 3 287 DO na = 1, nat 288 DO i = 1, 3 289 dynr (1, i, na, j, nb) = dynr (1, i, na, j, nb) / SQRT (amass ( & 290 ityp (na) ) * amass (ityp (nb) ) ) / amu_ry 291 dynr (2, i, na, j, nb) = dynr (2, i, na, j, nb) / SQRT (amass ( & 292 ityp (na) ) * amass (ityp (nb) ) ) / amu_ry 293 ENDDO 294 ENDDO 295 ENDDO 296 ENDDO 297 ! 298 ! solve the eigenvalue problem. 299 ! NOTA BENE: eigenvectors are overwritten on dyn 300 ! 301 CALL cdiagh (3 * nat, dynr, 3 * nat, w2, dyn) 302 ! 303 ! divide by sqrt(mass) to get displacements 304 ! 305 DO nu = 1, 3 * nat 306 DO mu = 1, 3 * nat 307 na = (mu - 1) / 3 + 1 308 dyn (mu, nu) = dyn (mu, nu) / SQRT ( amu_ry * amass (ityp (na) ) ) 309 ENDDO 310 ENDDO 311 ! 312 ! 313 RETURN 314END SUBROUTINE readmat 315! 316!----------------------------------------------------------------------- 317SUBROUTINE elphel (irr, npe, imode0, dvscfins) 318 !----------------------------------------------------------------------- 319 ! 320 ! Calculation of the electron-phonon matrix elements el_ph_mat 321 ! <\psi(k+q)|dV_{SCF}/du^q_{i a}|\psi(k)> 322 ! Original routine written by Francesco Mauri 323 ! Modified by A. Floris and I. Timrov to include Hubbard U (01.10.2018) 324 ! 325 USE kinds, ONLY : DP 326 USE fft_base, ONLY : dffts 327 USE ions_base, ONLY : nat, ityp 328 USE control_flags, ONLY : iverbosity 329 USE wavefunctions, ONLY : evc 330 USE buffers, ONLY : get_buffer, save_buffer 331 USE klist, ONLY : xk, ngk, igk_k 332 USE lsda_mod, ONLY : lsda, current_spin, isk, nspin 333 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 334 USE wvfct, ONLY : nbnd, npwx 335 USE uspp, ONLY : vkb 336 USE el_phon, ONLY : el_ph_mat, el_ph_mat_rec, el_ph_mat_rec_col, & 337 comp_elph, done_elph, elph_nbnd_min, elph_nbnd_max 338 USE modes, ONLY : u, nmodes 339 USE units_ph, ONLY : iubar, lrbar, iundnsscf, iudvpsi, lrdvpsi 340 USE units_lr, ONLY : iuwfc, lrwfc 341 USE control_ph, ONLY : trans, current_iq 342 USE ph_restart, ONLY : ph_writefile 343 USE spin_orb, ONLY : domag 344 USE mp_bands, ONLY : intra_bgrp_comm, ntask_groups 345 USE mp_pools, ONLY : npool 346 USE mp, ONLY : mp_sum, mp_bcast 347 USE mp_world, ONLY : world_comm 348 USE elph_tetra_mod, ONLY : elph_tetra 349 USE eqv, ONLY : dvpsi, evq 350 USE qpoint, ONLY : nksq, ikks, ikqs, nksqtot 351 USE control_lr, ONLY : lgamma 352 USE fft_helper_subroutines 353 USE ldaU, ONLY : lda_plus_u, Hubbard_lmax 354 USE ldaU_ph, ONLY : dnsscf_all_modes, dnsscf 355 USE io_global, ONLY : ionode, ionode_id 356 USE io_files, ONLY : seqopn 357 USE lrus, ONLY : becp1 358 USE phus, ONLY : alphap 359 USE ahc, ONLY : elph_ahc, ib_ahc_gauge_min, ib_ahc_gauge_max 360 361 IMPLICIT NONE 362 ! 363 INTEGER, INTENT(IN) :: irr, npe, imode0 364 COMPLEX(DP), INTENT(IN) :: dvscfins (dffts%nnr, nspin_mag, npe) 365 ! LOCAL variables 366 INTEGER :: npw, npwq, nrec, ik, ikk, ikq, ipert, mode, ibnd, jbnd, ir, ig, & 367 ipol, ios, ierr, nrec_ahc 368 COMPLEX(DP) , ALLOCATABLE :: aux1 (:,:), elphmat (:,:,:), tg_dv(:,:), & 369 tg_psic(:,:), aux2(:,:) 370 INTEGER :: v_siz, incr 371 LOGICAL :: exst 372 COMPLEX(DP), EXTERNAL :: zdotc 373 integer :: ibnd_fst, ibnd_lst 374 ! 375 CALL start_clock('elphel') 376 ! 377 IF (elph_ahc) THEN 378 ibnd_fst = ib_ahc_gauge_min 379 ibnd_lst = ib_ahc_gauge_max 380 elseif(elph_tetra == 0) then 381 ibnd_fst = 1 382 ibnd_lst = nbnd 383 else 384 ibnd_fst = elph_nbnd_min 385 ibnd_lst = elph_nbnd_max 386 end if 387 ! 388 IF (.NOT. comp_elph(irr) .OR. done_elph(irr)) RETURN 389 390 ALLOCATE (aux1 (dffts%nnr, npol)) 391 ALLOCATE (elphmat ( nbnd , nbnd , npe)) 392 ALLOCATE( el_ph_mat_rec (nbnd,nbnd,nksq,npe) ) 393 el_ph_mat_rec=(0.0_DP,0.0_DP) 394 ALLOCATE (aux2(npwx*npol, nbnd)) 395 incr=1 396 IF ( dffts%has_task_groups ) THEN 397 ! 398 v_siz = dffts%nnr_tg 399 ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) 400 ALLOCATE( tg_psic( v_siz, npol ) ) 401 incr = fftx_ntgrp(dffts) 402 ! 403 ENDIF 404 ! 405 ! DFPT+U case 406 ! 407 IF (lda_plus_u .AND. .NOT.trans) THEN 408 ! 409 ! Allocate and read dnsscf_all_modes from file 410 ! 411 ALLOCATE (dnsscf_all_modes(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, nmodes)) 412 dnsscf_all_modes = (0.d0, 0.d0) 413 ! 414 IF (ionode) READ(iundnsscf,*) dnsscf_all_modes 415 CALL mp_bcast(dnsscf_all_modes, ionode_id, world_comm) 416 REWIND(iundnsscf) 417 ! 418 ! Check whether the re-read is correct 419 ! 420 IF (iverbosity==1) CALL elphel_read_dnsscf_check() 421 ! 422 ! Allocate dnsscf 423 ! 424 ALLOCATE (dnsscf(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, npe)) 425 dnsscf = (0.d0, 0.d0) 426 ! 427 ENDIF 428 ! 429 ! Start the loops over the k-points 430 ! 431 DO ik = 1, nksq 432 ! 433 ! ik = counter of k-points with vector k 434 ! ikk= index of k-point with vector k 435 ! ikq= index of k-point with vector k+q 436 ! k and k+q are alternated if q!=0, are the same if q=0 437 ! 438 ikk = ikks(ik) 439 ikq = ikqs(ik) 440 IF (lsda) current_spin = isk (ikk) 441 npw = ngk(ikk) 442 npwq= ngk(ikq) 443 ! 444 CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) 445 ! 446 ! read unperturbed wavefuctions psi(k) and psi(k+q) 447 ! 448 IF (nksq.GT.1) THEN 449 IF (lgamma) THEN 450 CALL get_buffer(evc, lrwfc, iuwfc, ikk) 451 ELSE 452 CALL get_buffer (evc, lrwfc, iuwfc, ikk) 453 CALL get_buffer (evq, lrwfc, iuwfc, ikq) 454 ENDIF 455 ENDIF 456 ! 457 DO ipert = 1, npe 458 nrec = (ipert - 1) * nksq + ik 459 ! 460 ! dvbare_q*psi_kpoint is read from file (if available) or recalculated 461 ! 462 IF (trans) THEN 463 CALL get_buffer (dvpsi, lrbar, iubar, nrec) 464 ELSE 465 mode = imode0 + ipert 466 ! FIXME: .false. or .true. ??? 467 CALL dvqpsi_us (ik, u (1, mode), .FALSE., becp1, alphap) 468 ! 469 ! DFPT+U: calculate the bare derivative of the Hubbard potential in el-ph 470 ! 471 IF (lda_plus_u) CALL dvqhub_barepsi_us (ik, u(1,mode)) 472 ! 473 ENDIF 474 ! 475 ! calculate dvscf_q*psi_k 476 ! 477 IF ( dffts%has_task_groups ) THEN 478 IF (noncolin) THEN 479 CALL tg_cgather( dffts, dvscfins(:,1,ipert), tg_dv(:,1)) 480 IF (domag) THEN 481 DO ipol=2,4 482 CALL tg_cgather( dffts, dvscfins(:,ipol,ipert), tg_dv(:,ipol)) 483 ENDDO 484 ENDIF 485 ELSE 486 CALL tg_cgather( dffts, dvscfins(:,current_spin,ipert), tg_dv(:,1)) 487 ENDIF 488 ENDIF 489 aux2=(0.0_DP,0.0_DP) 490 DO ibnd = ibnd_fst, ibnd_lst, incr 491 IF ( dffts%has_task_groups ) THEN 492 CALL cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd ) 493 CALL apply_dpot(v_siz, tg_psic, tg_dv, 1) 494 CALL cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd) 495 ELSE 496 CALL cft_wave (ik, evc(1, ibnd), aux1, +1) 497 CALL apply_dpot(dffts%nnr, aux1, dvscfins(1,1,ipert), current_spin) 498 CALL cft_wave (ik, aux2(1, ibnd), aux1, -1) 499 ENDIF 500 ENDDO 501 dvpsi=dvpsi+aux2 502 ! 503 CALL adddvscf (ipert, ik) 504 ! 505 ! DFPT+U: add to dvpsi the scf part of the perturbed Hubbard potential 506 ! 507 IF (lda_plus_u) THEN 508 IF (.NOT.trans) dnsscf(:,:,:,:,ipert) = dnsscf_all_modes(:,:,:,:,mode) 509 CALL adddvhubscf (ipert, ik) 510 ENDIF 511 ! 512 ! If doing Allen-Heine-Cardona (AHC) calculation, we need dvpsi 513 ! later. So, write to buffer. 514 ! 515 IF (elph_ahc) THEN 516 nrec_ahc = (ik - 1) * nmodes + ipert + imode0 517 CALL save_buffer(dvpsi(1, ibnd_fst), lrdvpsi, iudvpsi, nrec_ahc) 518 ! 519 ! If elph_ahc, the matrix elements are computed in ahc.f90 520 CYCLE 521 ! 522 ENDIF 523 ! 524 ! calculate elphmat(j,i)=<psi_{k+q,j}|dvscf_q*psi_{k,i}> for this pertur 525 ! 526 DO ibnd = ibnd_fst, ibnd_lst 527 DO jbnd = ibnd_fst, ibnd_lst 528 elphmat (jbnd, ibnd, ipert) = zdotc (npwq, evq (1, jbnd), 1, & 529 dvpsi (1, ibnd), 1) 530 IF (noncolin) & 531 elphmat (jbnd, ibnd, ipert) = elphmat (jbnd, ibnd, ipert)+ & 532 zdotc (npwq, evq(npwx+1,jbnd),1,dvpsi(npwx+1,ibnd), 1) 533 ENDDO 534 ENDDO 535 ENDDO 536 ! 537 ! If elph_ahc, the matrix elements are computed in ahc.f90 538 IF (elph_ahc) CYCLE 539 ! 540 CALL mp_sum (elphmat, intra_bgrp_comm) 541 ! 542 ! save all e-ph matrix elements into el_ph_mat 543 ! 544 DO ipert = 1, npe 545 DO jbnd = ibnd_fst, ibnd_lst 546 DO ibnd = ibnd_fst, ibnd_lst 547 el_ph_mat (ibnd, jbnd, ik, ipert + imode0) = elphmat (ibnd, jbnd, ipert) 548 el_ph_mat_rec (ibnd, jbnd, ik, ipert ) = elphmat (ibnd, jbnd, ipert) 549 ENDDO 550 ENDDO 551 ENDDO 552 ENDDO 553 ! 554 done_elph(irr)=.TRUE. 555 if(elph_tetra == 0 .AND. .NOT. elph_ahc) then 556 IF (npool>1) THEN 557 ALLOCATE(el_ph_mat_rec_col(nbnd,nbnd,nksqtot,npe)) 558 CALL el_ph_collect(npe,el_ph_mat_rec,el_ph_mat_rec_col,nksqtot,nksq) 559 ELSE 560 el_ph_mat_rec_col => el_ph_mat_rec 561 ENDIF 562 CALL ph_writefile('el_phon',current_iq,irr,ierr) 563 IF (npool > 1) DEALLOCATE(el_ph_mat_rec_col) 564 end if 565 DEALLOCATE(el_ph_mat_rec) 566 ! 567 DEALLOCATE (elphmat) 568 DEALLOCATE (aux1) 569 DEALLOCATE (aux2) 570 IF ( dffts%has_task_groups ) THEN 571 DEALLOCATE( tg_dv ) 572 DEALLOCATE( tg_psic ) 573 ENDIF 574 ! 575 IF (lda_plus_u .AND. .NOT.trans) THEN 576 DEALLOCATE (dnsscf_all_modes) 577 DEALLOCATE (dnsscf) 578 ENDIF 579 ! 580 CALL stop_clock('elphel') 581 ! 582 RETURN 583 ! 584END SUBROUTINE elphel 585! 586!------------------------------------------------------------------------ 587SUBROUTINE elphel_read_dnsscf_check() 588 ! 589 ! DFPT+U: This subroutine checks whether dnsscf_all_modes was 590 ! read correctly from file. 591 ! 592 USE kinds, ONLY : DP 593 USE ions_base, ONLY : nat, ityp 594 USE modes, ONLY : u, nmodes 595 USE lsda_mod, ONLY : nspin 596 USE ldaU, ONLY : Hubbard_l, is_hubbard, Hubbard_lmax 597 USE ldaU_ph, ONLY : dnsscf_all_modes 598 USE io_global, ONLY : stdout 599 ! 600 IMPLICIT NONE 601 ! 602 COMPLEX(DP), ALLOCATABLE :: dnsscf_all_modes_cart(:,:,:,:,:) 603 INTEGER :: na_icart, nah, is, m1, m2, na, icart, nt, na_icar, imode 604 ! 605 ALLOCATE(dnsscf_all_modes_cart (2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, nmodes)) 606 dnsscf_all_modes_cart = (0.d0, 0.d0) 607 ! 608 ! Transform dnsscf_all_modes from pattern to cartesian coordinates 609 ! 610 DO na_icart = 1, 3*nat 611 DO imode = 1, nmodes 612 DO nah = 1, nat 613 nt = ityp(nah) 614 IF (is_hubbard(nt)) THEN 615 DO is = 1, nspin 616 DO m1 = 1, 2*Hubbard_l(nt) + 1 617 DO m2 = 1, 2*Hubbard_l(nt) + 1 618 ! 619 dnsscf_all_modes_cart (m1, m2, is, nah, na_icart) = & 620 dnsscf_all_modes_cart (m1, m2, is, nah, na_icart) + & 621 dnsscf_all_modes (m1, m2, is, nah, imode) * & 622 CONJG(u(na_icart,imode)) 623 ! 624 ENDDO 625 ENDDO 626 ENDDO 627 ENDIF 628 ENDDO 629 ENDDO 630 ENDDO 631 ! 632 ! Write dnsscf in cartesian coordinates 633 ! 634 WRITE(stdout,*) 635 WRITE(stdout,*) 'DNS_SCF SYMMETRIZED IN CARTESIAN COORDINATES' 636 ! 637 DO na = 1, nat 638 DO icart = 1, 3 639 WRITE( stdout,'(a,1x,i2,2x,a,1x,i2)') 'displaced atom L =', na, 'ipol=', icart 640 na_icart = 3*(na-1) + icart 641 DO nah = 1, nat 642 nt = ityp(nah) 643 IF (is_hubbard(nt)) THEN 644 DO is = 1, nspin 645 WRITE(stdout,'(a,1x,i2,2x,a,1x,i2)') ' Hubbard atom', nah, 'spin', is 646 DO m1 = 1, 2*Hubbard_l(nt) + 1 647 WRITE(stdout,'(14(f15.10,1x))') dnsscf_all_modes_cart (m1,:,is,nah,na_icart) 648 ENDDO 649 ENDDO 650 ENDIF 651 ENDDO 652 ENDDO 653 ENDDO 654 WRITE(stdout,*) 655 ! 656 DEALLOCATE(dnsscf_all_modes_cart) 657 ! 658 RETURN 659 ! 660END SUBROUTINE elphel_read_dnsscf_check 661!------------------------------------------------------------------------ 662 663!------------------------------------------------------------------------ 664SUBROUTINE elphsum ( ) 665 !----------------------------------------------------------------------- 666 ! 667 ! Sum over BZ of the electron-phonon matrix elements el_ph_mat 668 ! Original routine written by Francesco Mauri, modified by PG 669 ! New version by Malgorzata Wierzbowska 670 ! 671 USE kinds, ONLY : DP 672 USE constants, ONLY : pi, rytoev, ry_to_cmm1, ry_to_ghz, degspin 673 USE ions_base, ONLY : nat, ityp, tau 674 USE cell_base, ONLY : at, bg 675 USE lsda_mod, ONLY: isk, nspin 676 USE klist, ONLY: nks, nkstot, xk, wk, nelec 677 USE start_k, ONLY: nk1, nk2, nk3 678 USE symm_base, ONLY: s, irt, nsym, invs 679 USE noncollin_module, ONLY: nspin_lsda, nspin_mag 680 USE wvfct, ONLY: nbnd, et 681 USE parameters, ONLY : npk 682 USE el_phon, ONLY : el_ph_mat, done_elph, el_ph_nsigma, el_ph_ngauss, & 683 el_ph_sigma 684 USE modes, ONLY : u, nirr 685 USE dynmat, ONLY : dyn, w2 686 USE io_global, ONLY : stdout, ionode, ionode_id 687 USE parallel_include 688 USE mp_pools, ONLY : my_pool_id, npool, kunit 689 USE mp_images, ONLY : intra_image_comm, me_image, nproc_image 690 USE mp, ONLY : mp_bcast 691 USE control_ph, ONLY : tmp_dir_phq, xmldyn, current_iq 692 USE save_ph, ONLY : tmp_dir_save 693 USE io_files, ONLY : prefix, tmp_dir, seqopn, create_directory 694 ! 695 USE lr_symm_base, ONLY : minus_q, nsymq, rtau 696 USE qpoint, ONLY : xq, nksq 697 USE control_lr, ONLY : lgamma 698 ! 699 IMPLICIT NONE 700 ! epsw = 20 cm^-1, in Ry 701 REAL(DP), PARAMETER :: epsw = 20.d0 / ry_to_cmm1 702 REAL(DP), PARAMETER :: eps = 1.0d-6 703 ! 704 INTEGER :: iuna2Fsave = 40 705 ! 706 REAL(DP), allocatable :: gam(:,:), lamb(:,:) 707 ! 708 ! Quantities ending with "fit" are relative to the "dense" grid 709 ! 710 REAL(DP), allocatable :: xkfit(:,:) 711 REAL(DP), allocatable, target :: etfit(:,:), wkfit(:) 712 INTEGER :: nksfit, nk1fit, nk2fit, nk3fit, nkfit, nksfit_real 713 INTEGER, allocatable :: eqkfit(:), eqqfit(:), sfit(:) 714 ! 715 integer :: nq, isq (48), imq 716 ! nq : degeneracy of the star of q 717 ! isq: index of q in the star of a given sym.op. 718 ! imq: index of -q in the star of q (0 if not present) 719 real(DP) :: sxq (3, 48) 720 ! list of vectors in the star of q 721 ! 722 ! workspace used for symmetrisation 723 ! 724 COMPLEX(DP), allocatable :: g1(:,:,:), g2(:,:,:), g0(:,:), gf(:,:,:) 725 COMPLEX(DP), allocatable :: point(:), noint(:) 726 COMPLEX(DP) :: dyn22(3*nat,3*nat), ctemp 727 ! 728 INTEGER :: ik, ikk, ikq, isig, ibnd, jbnd, ipert, jpert, nu, mu, & 729 vu, ngauss1, iuelph, ios, i,k,j, ii, jj 730 INTEGER :: nkBZ, nti, ntj, ntk, nkr, itemp1, itemp2, nn, & 731 qx,qy,qz,iq,jq,kq 732 INTEGER, ALLOCATABLE :: eqBZ(:), sBZ(:) 733 REAL(DP) :: weight, wqa, w0g1, w0g2, degauss1, effit1, dosef, & 734 ef1, lambda, gamma 735 REAL(DP), ALLOCATABLE :: deg(:), effit(:), dosfit(:) 736 REAL(DP) :: etk, etq 737 REAL(DP), EXTERNAL :: dos_ef, efermig, w0gauss 738 character(len=80) :: name 739 LOGICAL :: exst, xmldyn_save 740 ! 741 COMPLEX(DP) :: el_ph_sum (3*nat,3*nat) 742 743 COMPLEX(DP), POINTER :: el_ph_mat_collect(:,:,:,:) 744 REAL(DP), ALLOCATABLE :: xk_collect(:,:) 745 REAL(DP), POINTER :: wkfit_dist(:), etfit_dist(:,:) 746 INTEGER :: nksfit_dist, rest, kunit_save 747 INTEGER :: nks_real, ispin, nksqtot, irr, ierr 748 CHARACTER(LEN=256) :: elph_dir 749 CHARACTER(LEN=6) :: int_to_char 750 ! 751 ! 752 ! 753 ! If the electron phonon matrix elements have not been calculated for 754 ! all representations this routine exit 755 ! 756 DO irr=1,nirr 757 IF (.NOT.done_elph(irr)) RETURN 758 ENDDO 759 760 CALL start_clock('elphsum') 761 762 elph_dir='elph_dir/' 763 IF (ionode) INQUIRE(file=TRIM(elph_dir), EXIST=exst) 764 CALL mp_bcast(exst, ionode_id, intra_image_comm) 765 IF (.NOT.exst) CALL create_directory( elph_dir ) 766 WRITE (6, '(5x,"electron-phonon interaction ..."/)') 767 ngauss1 = 0 768 769 ALLOCATE(xk_collect(3,nkstot)) 770 771 ALLOCATE(deg(el_ph_nsigma)) 772 ALLOCATE(effit(el_ph_nsigma)) 773 ALLOCATE(dosfit(el_ph_nsigma)) 774 775 IF (npool==1) THEN 776! 777! no pool, just copy old variables on the new ones 778! 779 nksqtot=nksq 780 xk_collect(:,1:nks) = xk(:,1:nks) 781 el_ph_mat_collect => el_ph_mat 782 ELSE 783! 784! pools, allocate new variables and collect the results. All the rest 785! remain unchanged. 786! 787 IF (lgamma) THEN 788 nksqtot=nkstot 789 ELSE 790 nksqtot=nkstot/2 791 ENDIF 792 CALL poolcollect(3, nks, xk, nkstot, xk_collect) 793 ALLOCATE(el_ph_mat_collect(nbnd,nbnd,nksqtot,3*nat)) 794 ! FIXME: this routine should be replaced by a generic routine 795 CALL el_ph_collect(3*nat,el_ph_mat,el_ph_mat_collect,nksqtot,nksq) 796 ENDIF 797 ! 798 ! read eigenvalues for the dense grid 799 ! FIXME: this should be done from the xml file, not from a specialized file 800 ! parallel case: only first node reads 801 ! 802 IF ( ionode ) THEN 803 tmp_dir=tmp_dir_save 804 CALL seqopn( iuna2Fsave, 'a2Fsave', 'FORMATTED', exst ) 805 tmp_dir=tmp_dir_phq 806 READ(iuna2Fsave,*) ibnd, nksfit 807 END IF 808 ! 809 CALL mp_bcast (ibnd, ionode_id, intra_image_comm) 810 CALL mp_bcast (nksfit, ionode_id, intra_image_comm) 811 if ( ibnd /= nbnd ) call errore('elphsum','wrong file read',iuna2Fsave) 812 allocate (etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit)) 813 ! 814 IF ( ionode ) THEN 815 READ(iuna2Fsave,*) etfit 816 READ(iuna2Fsave,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit) 817 READ(iuna2Fsave,*) wkfit 818 READ(iuna2Fsave,*) nk1fit, nk2fit, nk3fit 819 CLOSE( UNIT = iuna2Fsave, STATUS = 'KEEP' ) 820 END IF 821 ! 822 ! broadcast all variables read 823 ! 824 CALL mp_bcast (etfit, ionode_id, intra_image_comm) 825 CALL mp_bcast (xkfit, ionode_id, intra_image_comm) 826 CALL mp_bcast (wkfit, ionode_id, intra_image_comm) 827 CALL mp_bcast (nk1fit, ionode_id, intra_image_comm) 828 CALL mp_bcast (nk2fit, ionode_id, intra_image_comm) 829 CALL mp_bcast (nk3fit, ionode_id, intra_image_comm) 830 ! 831 nkfit=nk1fit*nk2fit*nk3fit 832 ! 833 ! efermig and dos_ef require scattered points and eigenvalues 834 ! isk is neither read nor used. phonon with two Fermi energies is 835 ! not yet implemented. 836 ! 837 nksfit_dist = ( nksfit / npool ) 838 rest = ( nksfit - nksfit_dist * npool ) 839 IF ( ( my_pool_id + 1 ) <= rest ) nksfit_dist = nksfit_dist + 1 840 kunit_save=kunit 841 kunit=1 842 843#if defined(__MPI) 844 ALLOCATE(etfit_dist(nbnd,nksfit_dist)) 845 ALLOCATE(wkfit_dist(nksfit_dist)) 846 CALL poolscatter( 1, nksfit, wkfit, nksfit_dist, wkfit_dist ) 847 CALL poolscatter( nbnd, nksfit, etfit, nksfit_dist, etfit_dist ) 848#else 849 wkfit_dist => wkfit 850 etfit_dist => etfit 851#endif 852 ! 853 do isig=1,el_ph_nsigma 854 ! 855 ! recalculate Ef = effit and DOS at Ef N(Ef) = dosfit using dense grid 856 ! for value "deg" of gaussian broadening 857 ! 858 deg(isig) = isig * el_ph_sigma 859 ! 860 effit(isig) = efermig & 861 ( etfit_dist, nbnd, nksfit_dist, nelec, wkfit_dist, & 862 deg(isig), ngauss1, 0, isk) 863 dosfit(isig) = dos_ef ( ngauss1, deg(isig), effit(isig), etfit_dist, & 864 wkfit_dist, nksfit_dist, nbnd) / 2.0d0 865 enddo 866#if defined(__MPI) 867 DEALLOCATE(etfit_dist) 868 DEALLOCATE(wkfit_dist) 869#endif 870 kunit=kunit_save 871 allocate (eqkfit(nkfit), eqqfit(nkfit), sfit(nkfit)) 872 ! 873 ! map k-points in the IBZ to k-points in the complete uniform grid 874 ! 875 nksfit_real=nksfit/nspin_lsda 876 call lint ( nsym, s, .true., at, bg, npk, 0,0,0, & 877 nk1fit,nk2fit,nk3fit, nksfit_real, xkfit, 1, nkfit, eqkfit, sfit) 878 deallocate (sfit, xkfit, wkfit) 879 ! 880 ! find epsilon(k+q) in the dense grid 881 ! 882 call cryst_to_cart (1, xq, at, -1) 883 qx = nint(nk1fit*xq(1)) 884 if (abs(qx-nk1fit*xq(1)) > eps) & 885 call errore('elphsum','q is not a vector in the dense grid',1) 886 if (qx < 0) qx = qx + nk1fit 887 if (qx > nk1fit) qx = qx - nk1fit 888 qy = nint(nk2fit*xq(2)) 889 if (abs(qy-nk2fit*xq(2)) > eps) & 890 call errore('elphsum','q is not a vector in the dense grid',2) 891 if (qy < 0) qy = qy + nk2fit 892 if (qy > nk2fit) qy = qy - nk2fit 893 qz = nint(nk3fit*xq(3)) 894 if (abs(qz-nk3fit*xq(3)) > eps) & 895 call errore('elphsum','q is not a vector in the dense grid',3) 896 if (qz < 0) qz = qz + nk3fit 897 if (qz > nk3fit) qz = qz - nk3fit 898 call cryst_to_cart (1, xq, bg, 1) 899 ! 900 eqqfit(:) = 0 901 do i=1,nk1fit 902 do j=1,nk2fit 903 do k=1,nk3fit 904 ik = k-1 + (j-1)*nk3fit + (i-1)*nk2fit*nk3fit + 1 905 iq = i+qx 906 if (iq > nk1fit) iq = iq - nk1fit 907 jq = j+qy 908 if (jq > nk2fit) jq = jq - nk2fit 909 kq = k+qz 910 if (kq > nk3fit) kq = kq - nk3fit 911 nn = (kq-1)+(jq-1)*nk3fit+(iq-1)*nk2fit*nk3fit + 1 912 eqqfit(ik) = eqkfit(nn) 913 enddo 914 enddo 915 enddo 916 ! 917 ! calculate the electron-phonon coefficient using the dense grid 918 ! 919 nti = nk1fit/nk1 920 ntj = nk2fit/nk2 921 ntk = nk3fit/nk3 922 nkBZ = nk1*nk2*nk3 923 allocate (eqBZ(nkBZ), sBZ(nkBZ)) 924 ! 925 nks_real=nkstot/nspin_lsda 926 IF ( lgamma ) THEN 927 call lint ( nsymq, s, minus_q, at, bg, npk, 0,0,0, & 928 nk1,nk2,nk3, nks_real, xk_collect, 1, nkBZ, eqBZ, sBZ) 929 ELSE 930 call lint ( nsymq, s, minus_q, at, bg, npk, 0,0,0, & 931 nk1,nk2,nk3, nks_real, xk_collect, 2, nkBZ, eqBZ, sBZ) 932 END IF 933 ! 934 allocate (gf(3*nat,3*nat,el_ph_nsigma)) 935 gf = (0.0d0,0.0d0) 936 ! 937 wqa = 1.0d0/nkfit 938 IF (nspin==1) wqa=degspin*wqa 939 ! 940#if defined(__MPI) 941 do ibnd = me_image+1, nbnd, nproc_image 942#else 943 do ibnd = 1, nbnd 944#endif 945 do jbnd = 1, nbnd 946 allocate (g2(nkBZ*nspin_lsda,3*nat,3*nat)) 947 allocate (g1(nksqtot,3*nat,3*nat)) 948 do ik = 1, nksqtot 949 do ii = 1, 3*nat 950 do jj = 1, 3*nat 951 g1(ik,ii,jj)=CONJG(el_ph_mat_collect(jbnd,ibnd,ik,ii))* & 952 el_ph_mat_collect(jbnd,ibnd,ik,jj) 953 enddo ! ipert 954 enddo !jpert 955 enddo ! ik 956 ! 957 allocate (g0(3*nat,3*nat)) 958 do i=1,nk1 959 do j=1,nk2 960 do k=1,nk3 961 do ispin=1,nspin_lsda 962 nn = k-1 + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 963 itemp1 = eqBZ(nn) 964 if (ispin==2) itemp1=itemp1+nksqtot/2 965 g0(:,:) = g1(itemp1,:,:) 966 itemp2 = sBZ(nn) 967 call symm ( g0, u, xq, s, itemp2, rtau, irt, & 968 at, bg, nat) 969 if (ispin==2) nn=nn+nkBZ 970 g2(nn,:,:) = g0(:,:) 971 enddo 972 enddo ! k 973 enddo !j 974 enddo !i 975 deallocate (g0) 976 deallocate (g1) 977 ! 978 allocate ( point(nkBZ), noint(nkfit) ) 979 do jpert = 1, 3 * nat 980 do ipert = 1, 3 * nat 981 do ispin=1,nspin_lsda 982 ! 983 point(1:nkBZ) = & 984 g2(1+nkBZ*(ispin-1):nkBZ+nkBZ*(ispin-1),ipert,jpert) 985 ! 986 CALL clinear(nk1,nk2,nk3,nti,ntj,ntk,point,noint) 987 ! 988 do isig = 1,el_ph_nsigma 989 degauss1 = deg(isig) 990 effit1 = effit(isig) 991 ctemp = 0 992 do ik=1,nkfit 993 etk = etfit(ibnd,eqkfit(ik)+nksfit*(ispin-1)/2) 994 etq = etfit(jbnd,eqqfit(ik)+nksfit*(ispin-1)/2) 995 ctemp = ctemp & 996 + exp(-((effit1-etk)**2 + (effit1-etq)**2)/degauss1**2)*noint(ik) 997 enddo 998 gf(ipert,jpert,isig) = gf(ipert,jpert,isig) + & 999 ctemp * wqa / (degauss1**2) / pi 1000 enddo ! isig 1001 enddo ! ispin 1002 enddo ! ipert 1003 enddo !jpert 1004 deallocate (point, noint) 1005 deallocate (g2) 1006 ! 1007 enddo ! ibnd 1008 enddo ! jbnd 1009 1010#if defined(__MPI) 1011 CALL MPI_ALLREDUCE(MPI_IN_PLACE,gf,3*nat*3*nat*el_ph_nsigma, & 1012 MPI_DOUBLE_COMPLEX,MPI_SUM,intra_image_comm,ierr) 1013#endif 1014 1015 deallocate (eqqfit, eqkfit) 1016 deallocate (etfit) 1017 deallocate (eqBZ, sBZ) 1018! 1019 allocate (gam(3*nat,el_ph_nsigma), lamb(3*nat,el_ph_nsigma)) 1020 lamb(:,:) = 0.0d0 1021 gam (:,:) = 0.0d0 1022 do isig= 1, el_ph_nsigma 1023 do nu = 1,3*nat 1024 gam(nu,isig) = 0.0d0 1025 do mu = 1, 3 * nat 1026 do vu = 1, 3 * nat 1027 gam(nu,isig) = gam(nu,isig) + DBLE(conjg(dyn(mu,nu)) * & 1028 gf(mu,vu,isig) * dyn(vu,nu)) 1029 enddo 1030 enddo 1031 gam(nu,isig) = gam(nu,isig) * pi/2.0d0 1032 ! 1033 ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that appears 1034 ! in the definition of the electron-phonon matrix element g 1035 ! The sqrt(1/M) factor is actually hidden into the normal modes 1036 ! 1037 ! gamma = \pi \sum_k\sum_{i,j} \delta(e_{k,i}-Ef) \delta(e_{k+q,j}-Ef) 1038 ! | \sum_mu z(mu,nu) <psi_{k+q,j}|dvscf_q(mu)*psi_{k,i}> |^2 1039 ! where z(mu,nu) is the mu component of normal mode nu (z = dyn) 1040 ! gamma(nu) is the phonon linewidth of mode nu 1041 ! 1042 ! The factor N(Ef)^2 that appears in most formulations of el-ph interact 1043 ! is absent because we sum, not average, over the Fermi surface. 1044 ! The factor 2 is provided by the sum over spins 1045 ! 1046 if (sqrt(abs(w2(nu))) > epsw) then 1047 ! lambda is the adimensional el-ph coupling for mode nu: 1048 ! lambda(nu)= gamma(nu)/(pi N(Ef) \omega_{q,nu}^2) 1049 lamb(nu,isig) = gam(nu,isig)/pi/w2(nu)/dosfit(isig) 1050 else 1051 lamb(nu,isig) = 0.0d0 1052 endif 1053 gam(nu,isig) = gam(nu,isig)*ry_to_ghz 1054 enddo !nu 1055 enddo ! isig 1056 ! 1057 do isig= 1,el_ph_nsigma 1058 WRITE (6, 9000) deg(isig), ngauss1 1059 WRITE (6, 9005) dosfit(isig), effit(isig) * rytoev 1060 do nu=1,3*nat 1061 WRITE (6, 9010) nu, lamb(nu,isig), gam(nu,isig) 1062 enddo 1063 enddo 1064 ! Isaev: save files in suitable format for processing by lambda.x 1065 name=TRIM(elph_dir)// 'elph.inp_lambda.' //TRIM(int_to_char(current_iq)) 1066 1067 IF (ionode) THEN 1068 open(unit=12, file=TRIM(name), form='formatted', status='unknown', & 1069 iostat=ios) 1070 1071 write(12, "(5x,3f14.6,2i6)") xq(1),xq(2),xq(3), el_ph_nsigma, 3*nat 1072 write(12, "(6e14.6)") (w2(i), i=1,3*nat) 1073 1074 do isig= 1, el_ph_nsigma 1075 WRITE (12, 9000) deg(isig), ngauss1 1076 WRITE (12, 9005) dosfit(isig), effit(isig) * rytoev 1077 do nu=1,3*nat 1078 WRITE (12, 9010) nu, lamb(nu,isig), gam(nu,isig) 1079 enddo 1080 enddo 1081 close (unit=12,status='keep') 1082 ENDIF 1083 ! Isaev end 1084 CALL mp_bcast(ios, ionode_id, intra_image_comm) 1085 IF (ios /= 0) CALL errore('elphsum','problem opening file'//TRIM(name),1) 1086 deallocate (gam) 1087 deallocate (lamb) 1088 write(stdout,*) 1089 ! 1090 ! Prepare interface to q2r and matdyn 1091 ! 1092 call star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, .TRUE. ) 1093 ! 1094 do isig=1,el_ph_nsigma 1095 name=TRIM(elph_dir)//'a2Fq2r.'// TRIM(int_to_char(50 + isig)) & 1096 //'.'//TRIM(int_to_char(current_iq)) 1097 if (ionode) then 1098 iuelph = 4 1099 open(iuelph, file=name, STATUS = 'unknown', FORM = 'formatted', & 1100 iostat=ios) 1101 else 1102 ! 1103 ! this node doesn't write: unit 6 is redirected to /dev/null 1104 ! 1105 iuelph =6 1106 end if 1107 call mp_bcast(ios, ionode_id, intra_image_comm) 1108 IF (ios /= 0) call errore('elphsum','opening output file '// TRIM(name),1) 1109 dyn22(:,:) = gf(:,:,isig) 1110 write(iuelph,*) deg(isig), effit(isig), dosfit(isig) 1111 IF ( imq == 0 ) THEN 1112 write(iuelph,*) 2*nq 1113 ELSE 1114 write(iuelph,*) nq 1115 ENDIF 1116 xmldyn_save=xmldyn 1117 xmldyn=.FALSE. 1118 call q2qstar_ph (dyn22, at, bg, nat, nsym, s, invs, & 1119 irt, rtau, nq, sxq, isq, imq, iuelph) 1120 xmldyn=xmldyn_save 1121 if (ionode) CLOSE( UNIT = iuelph, STATUS = 'KEEP' ) 1122 enddo 1123 deallocate (gf) 1124 DEALLOCATE( deg ) 1125 DEALLOCATE( effit ) 1126 DEALLOCATE( dosfit ) 1127 DEALLOCATE( xk_collect ) 1128 IF (npool /= 1) DEALLOCATE(el_ph_mat_collect) 1129 1130 call stop_clock('elphsum') 1131 1132 ! 11339000 FORMAT(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) 11349005 FORMAT(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=', & 1135 & f10.6,' eV') 11369006 FORMAT(5x,'double delta at Ef =',f10.6) 11379010 FORMAT(5x,'lambda(',i5,')=',f8.4,' gamma=',f8.2,' GHz') 1138 ! 1139 RETURN 1140END SUBROUTINE elphsum 1141 1142!----------------------------------------------------------------------- 1143SUBROUTINE elphsum_simple 1144 !----------------------------------------------------------------------- 1145 ! 1146 ! Sum over BZ of the electron-phonon matrix elements el_ph_mat 1147 ! Original routine written by Francesco Mauri 1148 ! Rewritten by Matteo Calandra 1149 !----------------------------------------------------------------------- 1150 USE kinds, ONLY : DP 1151 USE constants, ONLY : pi, ry_to_cmm1, ry_to_ghz, rytoev 1152 USE ions_base, ONLY : nat 1153 USE cell_base, ONLY : at, bg 1154 USE symm_base, ONLY : s, irt, nsym, invs 1155 USE klist, ONLY : xk, nelec, nks, wk 1156 USE wvfct, ONLY : nbnd, et 1157 USE el_phon, ONLY : el_ph_mat, el_ph_nsigma, el_ph_ngauss, el_ph_sigma 1158 USE mp_pools, ONLY : inter_pool_comm, npool 1159 USE mp_images, ONLY : intra_image_comm 1160 USE output, ONLY : fildyn 1161 USE dynmat, ONLY : dyn, w2 1162 USE modes, ONLY : u, nirr 1163 USE control_ph, only : current_iq, qplot 1164 USE lsda_mod, only : isk 1165 USE el_phon, ONLY : done_elph, gamma_disp 1166 USE io_global, ONLY : stdout, ionode, ionode_id 1167 USE mp, ONLY: mp_sum, mp_bcast 1168 1169 USE lr_symm_base, ONLY : rtau, nsymq, irotmq, minus_q 1170 USE qpoint, ONLY : xq, nksq, ikks, ikqs 1171 ! 1172 IMPLICIT NONE 1173 REAL(DP), PARAMETER :: eps = 20_dp/ry_to_cmm1 ! eps = 20 cm^-1, in Ry 1174 ! 1175 INTEGER :: ik, ikk, ikq, isig, ibnd, jbnd, ipert, jpert, nu, mu, & 1176 vu, ngauss1, iuelph, ios, irr 1177 INTEGER :: nmodes 1178 REAL(DP) :: weight, w0g1, w0g2, w0gauss, wgauss, degauss1, dosef, & 1179 ef1, phase_space, lambda, gamma 1180 REAL(DP), EXTERNAL :: dos_ef, efermig 1181 character(len=80) :: filelph 1182 CHARACTER(len=256) :: file_elphmat 1183 ! 1184 COMPLEX(DP) :: el_ph_sum (3*nat,3*nat), dyn_corr(3*nat,3*nat) 1185 1186 INTEGER, EXTERNAL :: find_free_unit 1187 CHARACTER(LEN=6) :: int_to_char 1188 1189 1190 DO irr=1,nirr 1191 IF (.NOT.done_elph(irr)) RETURN 1192 ENDDO 1193 1194 nmodes=3*nat 1195 1196 filelph=TRIM(fildyn)//'.elph.'//TRIM(int_to_char(current_iq)) 1197 1198 ! parallel case: only first node writes 1199 IF ( ionode ) THEN 1200 ! 1201 iuelph = find_free_unit() 1202 OPEN (unit = iuelph, file = TRIM(filelph), status = 'unknown', err = & 1203 100, iostat = ios) 1204 REWIND (iuelph) 1205 ELSE 1206 iuelph = 0 1207 ! 1208 END IF 1209100 CONTINUE 1210 CALL mp_bcast(ios,ionode_id,intra_image_comm) 1211 CALL errore ('elphsum_simple', 'opening file '//filelph, ABS (ios) ) 1212 1213 IF (ionode) THEN 1214 WRITE (iuelph, '(3f15.8,2i8)') xq, el_ph_nsigma, 3 * nat 1215 WRITE (iuelph, '(6e14.6)') (w2 (nu) , nu = 1, nmodes) 1216 ENDIF 1217 1218 1219 ngauss1=0 1220 DO isig = 1, el_ph_nsigma 1221 ! degauss1 = 0.01 * isig 1222 degauss1 = el_ph_sigma * isig 1223 el_ph_sum(:,:) = (0.d0, 0.d0) 1224 phase_space = 0.d0 1225 ! 1226 ! Recalculate the Fermi energy Ef=ef1 and the DOS at Ef, dosef = N(Ef) 1227 ! for this gaussian broadening 1228 ! 1229 ! Note that the weights of k+q points must be set to zero for the 1230 ! following call to yield correct results 1231 ! 1232 1233 ef1 = efermig (et, nbnd, nks, nelec, wk, degauss1, el_ph_ngauss, 0, isk) 1234 dosef = dos_ef (el_ph_ngauss, degauss1, ef1, et, wk, nks, nbnd) 1235 ! N(Ef) is the DOS per spin, not summed over spin 1236 dosef = dosef / 2.d0 1237 ! 1238 ! Sum over bands with gaussian weights 1239 ! 1240 1241 DO ik = 1, nksq 1242 1243 ! 1244 ! see subroutine elphel for the logic of indices 1245 ! 1246 ikk = ikks(ik) 1247 ikq = ikqs(ik) 1248 DO ibnd = 1, nbnd 1249 w0g1 = w0gauss ( (ef1 - et (ibnd, ikk) ) / degauss1, ngauss1) & 1250 / degauss1 1251 DO jbnd = 1, nbnd 1252 w0g2 = w0gauss ( (ef1 - et (jbnd, ikq) ) / degauss1, ngauss1) & 1253 / degauss1 1254 ! note that wk(ikq)=wk(ikk) 1255 weight = wk (ikk) * w0g1 * w0g2 1256 DO jpert = 1, 3 * nat 1257 DO ipert = 1, 3 * nat 1258 el_ph_sum (ipert, jpert) = el_ph_sum (ipert, jpert) + weight * & 1259 CONJG (el_ph_mat (jbnd, ibnd, ik, ipert) ) * & 1260 el_ph_mat (jbnd, ibnd, ik, jpert) 1261 ENDDO 1262 ENDDO 1263 phase_space = phase_space+weight 1264 ENDDO 1265 ENDDO 1266 1267 ENDDO 1268 1269 ! el_ph_sum(mu,nu)=\sum_k\sum_{i,j}[ <psi_{k+q,j}|dvscf_q(mu)*psi_{k,i}> 1270 ! x <psi_{k+q,j}|dvscf_q(nu)*psi_{k,i}> 1271 ! x \delta(e_{k,i}-Ef) \delta(e_{k+q,j} 1272 ! 1273 ! collect contributions from all pools (sum over k-points) 1274 ! 1275 1276 call mp_sum ( el_ph_sum , inter_pool_comm ) 1277 call mp_sum ( phase_space , inter_pool_comm ) 1278 1279 ! 1280 ! symmetrize el_ph_sum(mu,nu) : it transforms as the dynamical matrix 1281 ! 1282 1283 CALL symdyn_munu_new (el_ph_sum, u, xq, s, invs, rtau, irt, at, & 1284 bg, nsymq, nat, irotmq, minus_q) 1285 ! 1286 WRITE (stdout, *) 1287 WRITE (stdout, 9000) degauss1, ngauss1 1288 WRITE (stdout, 9005) dosef, ef1 * rytoev 1289 WRITE (stdout, 9006) phase_space 1290 IF (ionode) THEN 1291 WRITE (iuelph, 9000) degauss1, ngauss1 1292 WRITE (iuelph, 9005) dosef, ef1 * rytoev 1293 ENDIF 1294 1295 DO nu = 1, nmodes 1296 gamma = 0.d0 1297 DO mu = 1, 3 * nat 1298 DO vu = 1, 3 * nat 1299 gamma = gamma + DBLE (CONJG (dyn (mu, nu) ) * el_ph_sum (mu, vu)& 1300 * dyn (vu, nu) ) 1301 ENDDO 1302 ENDDO 1303 gamma = pi * gamma / 2.d0 1304 ! 1305 ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that appears 1306 ! in the definition of the electron-phonon matrix element g 1307 ! The sqrt(1/M) factor is actually hidden into the normal modes 1308 ! 1309 ! gamma = \pi \sum_k\sum_{i,j} \delta(e_{k,i}-Ef) \delta(e_{k+q,j}-Ef) 1310 ! | \sum_mu z(mu,nu) <psi_{k+q,j}|dvscf_q(mu)*psi_{k,i}> |^2 1311 ! where z(mu,nu) is the mu component of normal mode nu (z = dyn) 1312 ! gamma(nu) is the phonon linewidth of mode nu 1313 ! 1314 ! The factor N(Ef)^2 that appears in most formulations of el-ph interact 1315 ! is absent because we sum, not average, over the Fermi surface. 1316 ! The factor 2 is provided by the sum over spins 1317 ! 1318 IF (SQRT (ABS (w2 (nu) ) ) > eps) THEN 1319 ! lambda is the adimensional el-ph coupling for mode nu: 1320 ! lambda(nu)= gamma(nu)/(pi N(Ef) \omega_{q,nu}^2) 1321 lambda = gamma / pi / w2 (nu) / dosef 1322 ELSE 1323 lambda = 0.d0 1324 ENDIF 1325 WRITE (stdout, 9010) nu, lambda, gamma * ry_to_gHz 1326 IF (ionode) WRITE (iuelph, 9010) nu, lambda, gamma * ry_to_gHz 1327 IF (qplot) gamma_disp(nu,isig,current_iq) = gamma * ry_to_gHz 1328 ENDDO 1329 ENDDO 1330 1331 13329000 FORMAT(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) 13339005 FORMAT(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=', & 1334 & f10.6,' eV') 13359006 FORMAT(5x,'double delta at Ef =',f10.6) 13369010 FORMAT(5x,'lambda(',i5,')=',f8.4,' gamma=',f8.2,' GHz') 1337 ! 1338 ! 1339 IF (ionode) CLOSE (unit = iuelph) 1340 RETURN 1341 1342 1343 1344 ! call star_q(x_q(1,iq), at, bg, nsym , s , invs , nq, sxq, & 1345 ! isq, imq, .FALSE. ) 1346 1347 1348END SUBROUTINE elphsum_simple 1349 1350!----------------------------------------------------------------------- 1351SUBROUTINE elphfil_epa(iq) 1352 !----------------------------------------------------------------------- 1353 ! 1354 ! Writes electron-phonon matrix elements to a file 1355 ! which is subsequently processed by the epa code 1356 ! Original routine written by Georgy Samsonidze 1357 ! 1358 !----------------------------------------------------------------------- 1359 USE cell_base, ONLY : ibrav, alat, omega, tpiba, at, bg 1360 USE disp, ONLY : nq1, nq2, nq3, nqs, x_q, wq, lgamma_iq 1361 USE dynmat, ONLY : dyn, w2 1362 USE el_phon, ONLY : el_ph_mat, done_elph 1363 USE fft_base, ONLY : dfftp, dffts, dfftb 1364 USE gvect, ONLY : ngm_g, ecutrho 1365 USE io_global, ONLY : ionode, ionode_id 1366 USE ions_base, ONLY : nat, nsp, atm, ityp, tau 1367 USE kinds, ONLY : DP 1368 USE klist, ONLY : xk, wk, nelec, nks, nkstot, ngk 1369 USE lsda_mod, ONLY : nspin, isk 1370 USE modes, ONLY : nirr, nmodes, npert, npertx, u, t, tmq, & 1371 name_rap_mode, num_rap_mode 1372 USE lr_symm_base, ONLY : irgq, nsymq, irotmq, rtau, gi, gimq, & 1373 minus_q, invsymq 1374 USE mp, ONLY : mp_bcast, mp_sum 1375 USE mp_images, ONLY : intra_image_comm 1376 USE mp_pools, ONLY : npool, intra_pool_comm 1377 USE qpoint, ONLY : nksq, nksqtot, ikks, ikqs, eigqts 1378 USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3 1379 USE symm_base, ONLY : s, invs, ft, nrot, nsym, nsym_ns, & 1380 nsym_na, ft, sr, sname, t_rev, irt, time_reversal, & 1381 invsym, nofrac, allfrac, nosym, nosym_evc, no_t_rev 1382 USE wvfct, ONLY : nbnd, et, wg 1383 USE gvecw, ONLY : ecutwfc 1384 USE io_files, ONLY : prefix 1385 1386 IMPLICIT NONE 1387 1388 INTEGER, INTENT(IN) :: iq 1389 1390 INTEGER :: iuelph, ios, irr, ii, jj, kk, ll 1391 character :: cdate*9, ctime*9, sdate*32, stime*32, & 1392 stitle*32, myaccess*10, mystatus*7 1393 CHARACTER(LEN=80) :: filelph 1394 1395 REAL(DP), ALLOCATABLE :: xk_collect(:,:), wk_collect(:) 1396 REAL(DP), ALLOCATABLE :: et_collect(:,:), wg_collect(:,:) 1397 INTEGER, ALLOCATABLE :: ngk_collect(:) 1398 INTEGER, ALLOCATABLE :: ikks_collect(:), ikqs_collect(:) 1399 COMPLEX(DP), ALLOCATABLE :: el_ph_mat_collect(:,:,:,:) 1400 INTEGER :: ftau(3,48) 1401 INTEGER, EXTERNAL :: find_free_unit, atomic_number 1402 1403 filelph = TRIM(prefix) // '.epa.k' 1404 1405 DO irr = 1, nirr 1406 IF (.NOT. done_elph(irr)) RETURN 1407 ENDDO 1408 1409 IF (iq .EQ. 1) THEN 1410 myaccess = 'sequential' 1411 mystatus = 'replace' 1412 ELSE 1413 myaccess = 'append' 1414 mystatus = 'old' 1415 ENDIF 1416 IF (ionode) THEN 1417 iuelph = find_free_unit() 1418 OPEN(unit = iuelph, file = TRIM(filelph), form = 'unformatted', & 1419 access = myaccess, status = mystatus, iostat = ios) 1420 ELSE 1421 iuelph = 0 1422 ENDIF 1423 CALL mp_bcast(ios, ionode_id, intra_image_comm) 1424 CALL errore('elphfil_epa', 'opening file ' // filelph, ABS(ios)) 1425 1426 IF (iq .EQ. 1) THEN 1427 CALL date_and_tim(cdate, ctime) 1428 WRITE(sdate, '(A2,"-",A3,"-",A4,21X)') cdate(1:2), cdate(3:5), cdate(6:9) 1429 WRITE(stime, '(A8,24X)') ctime(1:8) 1430 WRITE(stitle, '("EPA-Complex",21X)') 1431 CALL cryst_to_cart(nqs, x_q, at, -1) 1432 ! write header 1433 IF (ionode) THEN 1434 WRITE(iuelph) stitle, sdate, stime 1435 WRITE(iuelph) ibrav, nat, nsp, nrot, nsym, nsym_ns, nsym_na, & 1436 ngm_g, nspin, nbnd, nmodes, nqs 1437 WRITE(iuelph) nq1, nq2, nq3, nk1, nk2, nk3, k1, k2, k3 1438 WRITE(iuelph) time_reversal, invsym, nofrac, allfrac, nosym, & 1439 nosym_evc, no_t_rev 1440 WRITE(iuelph) alat, omega, tpiba, nelec, ecutrho, ecutwfc 1441 WRITE(iuelph) dfftp%nr1, dfftp%nr2, dfftp%nr3 1442 WRITE(iuelph) dffts%nr1, dffts%nr2, dffts%nr3 1443 WRITE(iuelph) dfftb%nr1, dfftb%nr2, dfftb%nr3 1444 WRITE(iuelph) ((at(ii, jj), ii = 1, 3), jj = 1, 3) 1445 WRITE(iuelph) ((bg(ii, jj), ii = 1, 3), jj = 1, 3) 1446 WRITE(iuelph) (atomic_number(atm(ii)), ii = 1, nsp) 1447 WRITE(iuelph) (ityp(ii), ii = 1, nat) 1448 WRITE(iuelph) ((tau(ii, jj), ii = 1, 3), jj = 1, nat) 1449 WRITE(iuelph) ((x_q(ii, jj), ii = 1, 3), jj = 1, nqs) 1450 WRITE(iuelph) (wq(ii), ii = 1, nqs) 1451 WRITE(iuelph) (lgamma_iq(ii), ii = 1, nqs) 1452 ENDIF 1453 CALL cryst_to_cart(nqs, x_q, bg, 1) 1454 ENDIF 1455 1456 ! collect data for current q-point 1457 ALLOCATE(xk_collect(3, nkstot)) 1458 ALLOCATE(wk_collect(nkstot)) 1459 ALLOCATE(et_collect(nbnd, nkstot)) 1460 ALLOCATE(wg_collect(nbnd, nkstot)) 1461 ALLOCATE(ngk_collect(nkstot)) 1462 ALLOCATE(ikks_collect(nksqtot)) 1463 ALLOCATE(ikqs_collect(nksqtot)) 1464 ALLOCATE(el_ph_mat_collect(nbnd, nbnd, nksqtot, nmodes)) 1465 IF (npool > 1) THEN 1466 CALL poolcollect(3, nks, xk, nkstot, xk_collect) 1467 CALL poolcollect(1, nks, wk, nkstot, wk_collect) 1468 CALL poolcollect(nbnd, nks, et, nkstot, et_collect) 1469 CALL poolcollect(nbnd, nks, wg, nkstot, wg_collect) 1470 CALL ipoolcollect(1, nks, ngk, nkstot, ngk_collect) 1471 CALL jpoolcollect(1, nksq, ikks, nksqtot, ikks_collect) 1472 CALL jpoolcollect(1, nksq, ikqs, nksqtot, ikqs_collect) 1473 CALL el_ph_collect(nmodes, el_ph_mat, el_ph_mat_collect, nksqtot, nksq) 1474 ELSE 1475 xk_collect(1:3, 1:nks) = xk(1:3, 1:nks) 1476 wk_collect(1:nks) = wk(1:nks) 1477 et_collect(1:nbnd, 1:nks) = et(1:nbnd, 1:nks) 1478 wg_collect(1:nbnd, 1:nks) = wg(1:nbnd, 1:nks) 1479 ngk_collect(1:nks) = ngk(1:nks) 1480 ikks_collect(1:nksq) = ikks(1:nksq) 1481 ikqs_collect(1:nksq) = ikqs(1:nksq) 1482 el_ph_mat_collect(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes) = & 1483 el_ph_mat(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes) 1484 ENDIF 1485 CALL cryst_to_cart(nkstot, xk_collect, at, -1) 1486 ! write data for current q-point 1487 IF (ionode) THEN 1488 WRITE(iuelph) nsymq, irotmq, nirr, npertx, nkstot, nksqtot 1489 WRITE(iuelph) minus_q, invsymq 1490 WRITE(iuelph) (irgq(ii), ii = 1, 48) 1491 WRITE(iuelph) (npert(ii), ii = 1, nmodes) 1492 WRITE(iuelph) (((rtau(ii, jj, kk), ii = 1, 3), jj = 1, 48), & 1493 kk = 1, nat) 1494 WRITE(iuelph) ((gi(ii, jj), ii = 1, 3), jj = 1, 48) 1495 WRITE(iuelph) (gimq(ii), ii = 1, 3) 1496 WRITE(iuelph) ((u(ii, jj), ii = 1, nmodes), jj = 1, nmodes) 1497 WRITE(iuelph) ((((t(ii, jj, kk, ll), ii = 1, npertx), & 1498 jj = 1, npertx), kk = 1, 48), ll = 1, nmodes) 1499 WRITE(iuelph) (((tmq(ii, jj, kk), ii = 1, npertx), & 1500 jj = 1, npertx), kk = 1, nmodes) 1501 WRITE(iuelph) (name_rap_mode(ii), ii = 1, nmodes) 1502 WRITE(iuelph) (num_rap_mode(ii), ii = 1, nmodes) 1503 WRITE(iuelph) (((s(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48) 1504 WRITE(iuelph) (invs(ii), ii = 1, 48) 1505 ! FIXME: should disappear 1506 ftau(1,1:48) = NINT(ft(1,1:48)*dfftp%nr1) 1507 ftau(2,1:48) = NINT(ft(2,1:48)*dfftp%nr2) 1508 ftau(3,1:48) = NINT(ft(3,1:48)*dfftp%nr3) 1509 WRITE(iuelph) ((ftau(ii, jj), ii = 1, 3), jj = 1, 48) 1510 ! end FIXME 1511 WRITE(iuelph) ((ft(ii, jj), ii = 1, 3), jj = 1, 48) 1512 WRITE(iuelph) (((sr(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48) 1513 WRITE(iuelph) (sname(ii), ii = 1, 48) 1514 WRITE(iuelph) (t_rev(ii), ii = 1, 48) 1515 WRITE(iuelph) ((irt(ii, jj), ii = 1, 48), jj = 1, nat) 1516 WRITE(iuelph) ((xk_collect(ii, jj), ii = 1, 3), jj = 1, nkstot) 1517 WRITE(iuelph) (wk_collect(ii), ii = 1, nkstot) 1518 WRITE(iuelph) ((et_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot) 1519 WRITE(iuelph) ((wg_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot) 1520 WRITE(iuelph) (isk(ii), ii = 1, nkstot) 1521 WRITE(iuelph) (ngk_collect(ii), ii = 1, nkstot) 1522 WRITE(iuelph) (ikks_collect(ii), ii = 1, nksqtot) 1523 WRITE(iuelph) (ikqs_collect(ii), ii = 1, nksqtot) 1524 WRITE(iuelph) (eigqts(ii), ii = 1, nat) 1525 WRITE(iuelph) (w2(ii), ii = 1, nmodes) 1526 WRITE(iuelph) ((dyn(ii, jj), ii = 1, nmodes), jj = 1, nmodes) 1527 WRITE(iuelph) ((((el_ph_mat_collect(ii, jj, kk, ll), ii = 1, nbnd), & 1528 jj = 1, nbnd), kk = 1, nksqtot), ll = 1, nmodes) 1529 CLOSE (unit = iuelph, status = 'keep') 1530 ENDIF 1531 CALL cryst_to_cart(nkstot, xk_collect, bg, 1) 1532 DEALLOCATE(xk_collect) 1533 DEALLOCATE(wk_collect) 1534 DEALLOCATE(et_collect) 1535 DEALLOCATE(wg_collect) 1536 DEALLOCATE(ngk_collect) 1537 DEALLOCATE(ikks_collect) 1538 DEALLOCATE(ikqs_collect) 1539 DEALLOCATE(el_ph_mat_collect) 1540 1541 RETURN 1542 1543END SUBROUTINE elphfil_epa 1544 1545!---------------------------------------------------------------------------- 1546SUBROUTINE ipoolcollect( length, nks, f_in, nkstot, f_out ) 1547 !---------------------------------------------------------------------------- 1548 ! 1549 ! ... as poolcollect, for an integer vector 1550 ! 1551 USE mp_pools, ONLY : my_pool_id, npool, kunit, & 1552 inter_pool_comm, intra_pool_comm 1553 USE mp, ONLY : mp_sum 1554 ! 1555 IMPLICIT NONE 1556 ! 1557 INTEGER, INTENT(IN) :: length, nks, nkstot 1558 ! first dimension of arrays 1559 ! number of k-points per pool 1560 ! total number of k-points 1561 INTEGER, INTENT(IN) :: f_in (length,nks) 1562 ! pool-distributed function 1563 INTEGER, INTENT(OUT) :: f_out(length,nkstot) 1564 ! pool-collected function 1565 ! 1566 INTEGER :: nbase, rest, nks1 1567 ! 1568 nks1 = kunit * ( nkstot / kunit / npool ) 1569 ! 1570 rest = ( nkstot - nks1 * npool ) / kunit 1571 ! 1572 IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + kunit 1573 ! 1574 IF (nks1.ne.nks) & 1575 call errore('ipoolcollect','inconsistent number of k-points',1) 1576 ! 1577 ! ... calculates nbase = the position in the list of the first point that 1578 ! ... belong to this npool - 1 1579 ! 1580 nbase = nks * my_pool_id 1581 ! 1582 IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest * kunit 1583 ! 1584 ! copy the original points in the correct position of the list 1585 ! 1586 f_out=0 1587 f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks) 1588 ! 1589 CALL mp_sum( f_out, inter_pool_comm ) 1590 ! 1591 RETURN 1592 ! 1593END SUBROUTINE ipoolcollect 1594 1595!---------------------------------------------------------------------------- 1596SUBROUTINE jpoolcollect( length, nks, f_in, nkstot, f_out ) 1597 !---------------------------------------------------------------------------- 1598 ! 1599 ! ... as ipoolcollect, without kunit and with an index shift 1600 ! 1601 USE mp_pools, ONLY : my_pool_id, npool, kunit, & 1602 inter_pool_comm, intra_pool_comm 1603 USE mp, ONLY : mp_sum 1604 ! 1605 IMPLICIT NONE 1606 ! 1607 INTEGER, INTENT(IN) :: length, nks, nkstot 1608 ! first dimension of arrays 1609 ! number of k-points per pool 1610 ! total number of k-points 1611 INTEGER, INTENT(IN) :: f_in (length,nks) 1612 ! pool-distributed function 1613 INTEGER, INTENT(OUT) :: f_out(length,nkstot) 1614 ! pool-collected function 1615 ! 1616 INTEGER :: nbase, rest, nks1 1617 ! 1618 nks1 = ( nkstot / npool ) 1619 ! 1620 rest = ( nkstot - nks1 * npool ) 1621 ! 1622 IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + 1 1623 ! 1624 IF (nks1.ne.nks) & 1625 call errore('jpoolcollect','inconsistent number of k-points',1) 1626 ! 1627 ! ... calculates nbase = the position in the list of the first point that 1628 ! ... belong to this npool - 1 1629 ! 1630 nbase = nks * my_pool_id 1631 ! 1632 IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest 1633 ! 1634 ! copy the original points in the correct position of the list 1635 ! 1636 f_out=0 1637 f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks) + nbase * kunit 1638 ! 1639 CALL mp_sum( f_out, inter_pool_comm ) 1640 ! 1641 RETURN 1642 ! 1643END SUBROUTINE jpoolcollect 1644 1645!----------------------------------------------------------------------- 1646FUNCTION dos_ef (ngauss, degauss, ef, et, wk, nks, nbnd) 1647 !----------------------------------------------------------------------- 1648 ! 1649 USE kinds, ONLY : DP 1650 USE mp_pools, ONLY : inter_pool_comm 1651 USE mp, ONLY : mp_sum 1652 IMPLICIT NONE 1653 REAL(DP) :: dos_ef 1654 INTEGER :: ngauss, nbnd, nks 1655 REAL(DP) :: et (nbnd, nks), wk (nks), ef, degauss 1656 ! 1657 INTEGER :: ik, ibnd 1658 REAL(DP), EXTERNAL :: w0gauss 1659 ! 1660 ! Compute DOS at E_F (states per Ry per unit cell) 1661 ! 1662 dos_ef = 0.0d0 1663 DO ik = 1, nks 1664 DO ibnd = 1, nbnd 1665 dos_ef = dos_ef + wk (ik) * w0gauss ( (et (ibnd, ik) - ef) & 1666 / degauss, ngauss) / degauss 1667 ENDDO 1668 ENDDO 1669 ! 1670 ! Collects partial sums on k-points from all pools 1671 ! 1672 CALL mp_sum ( dos_ef, inter_pool_comm ) 1673 ! 1674 RETURN 1675END FUNCTION dos_ef 1676 1677!a2F 1678subroutine lint ( nsym, s, minus_q, at, bg, npk, k1,k2,k3, & 1679 nk1,nk2,nk3, nks, xk, kunit, nkBZ, eqBZ, sBZ) 1680 !----------------------------------------------------------------------- 1681 ! 1682 ! Find which k-points of a uniform grid are in the IBZ 1683 ! 1684 use kinds, only : DP 1685 implicit none 1686 integer, intent (IN) :: nks, nsym, s(3,3,48), npk, k1, k2, k3, & 1687 nk1, nk2, nk3, kunit, nkBZ 1688 logical, intent (IN) :: minus_q 1689 real(kind=DP), intent(IN):: at(3,3), bg(3,3), xk(3,npk) 1690 integer, INTENT(OUT) :: eqBZ(nkBZ), sBZ(nkBZ) 1691 ! 1692 real(kind=DP) :: xkr(3), deltap(3), deltam(3) 1693 real(kind=DP), parameter:: eps=1.0d-5 1694 real(kind=DP), allocatable :: xkg(:,:), xp(:,:) 1695 integer :: i,j,k, ns, n, nk 1696 integer :: nkh 1697 ! 1698 ! Re-generate a uniform grid of k-points xkg 1699 ! 1700 allocate (xkg( 3,nkBZ)) 1701 ! 1702 if(kunit < 1 .or. kunit > 2) call errore('lint','bad kunit value',kunit) 1703 ! 1704 ! kunit=2: get only "true" k points, not k+q points, from the list 1705 ! 1706 nkh = nks/kunit 1707 allocate (xp(3,nkh)) 1708 if (kunit == 1) then 1709 xp(:,1:nkh) = xk(:,1:nkh) 1710 else 1711 do j=1,nkh 1712 xp(:,j) = xk(:,2*j-1) 1713 enddo 1714 end if 1715 do i=1,nk1 1716 do j=1,nk2 1717 do k=1,nk3 1718 n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 1719 xkg(1,n) = dble(i-1)/nk1 + dble(k1)/2/nk1 1720 xkg(2,n) = dble(j-1)/nk2 + dble(k2)/2/nk2 1721 xkg(3,n) = dble(k-1)/nk3 + dble(k3)/2/nk3 1722 end do 1723 end do 1724 end do 1725 1726 call cryst_to_cart (nkh,xp,at,-1) 1727 1728 do nk=1,nkBZ 1729 do n=1,nkh 1730 do ns=1,nsym 1731 do i=1,3 1732 xkr(i) = s(i,1,ns) * xp(1,n) + & 1733 s(i,2,ns) * xp(2,n) + & 1734 s(i,3,ns) * xp(3,n) 1735 end do 1736 do i=1,3 1737 deltap(i) = xkr(i)-xkg(i,nk) - nint (xkr(i)-xkg(i,nk) ) 1738 deltam(i) = xkr(i)+xkg(i,nk) - nint (xkr(i)+xkg(i,nk) ) 1739 end do 1740 if ( sqrt ( deltap(1)**2 + & 1741 deltap(2)**2 + & 1742 deltap(3)**2 ) < eps .or. ( minus_q .and. & 1743 sqrt ( deltam(1)**2 + & 1744 deltam(2)**2 + & 1745 deltam(3)**2 ) < eps ) ) then 1746 eqBZ(nk) = n 1747 sBZ(nk) = ns 1748 go to 15 1749 end if 1750 end do 1751 end do 1752 call errore('lint','cannot locate k point xk',nk) 175315 continue 1754 end do 1755 1756 do n=1,nkh 1757 do nk=1,nkBZ 1758 if (eqBZ(nk) == n) go to 20 1759 end do 1760 ! this failure of the algorithm may indicate that the displaced grid 1761 ! (with k1,k2,k3.ne.0) does not have the full symmetry of the lattice 1762 call errore('lint','cannot remap grid on k-point list',n) 176320 continue 1764 end do 1765 1766 deallocate(xkg) 1767 deallocate(xp) 1768 1769 return 1770end subroutine lint 1771