1! 2! Copyright (C) 2001-2008 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 ep_matrix_element_wannier() 11 !----------------------------------------------------------------------- 12 ! 13 ! Electron-phonon calculation from data saved in fildvscf 14 ! 15 USE kinds, ONLY : DP 16 USE cell_base, ONLY : celldm, omega, ibrav 17 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, amass 18 USE gvecs, ONLY: doublegrid 19 USE fft_base, ONLY : dfftp, dffts 20 USE fft_interfaces, ONLY : fft_interpolate 21 USE noncollin_module, ONLY : nspin_mag, noncolin 22 USE dynmat, ONLY : dyn, w2 23 USE modes, ONLY : npert, nirr, u 24 USE control_ph, ONLY : trans 25 USE units_ph, ONLY : iudyn, lrdrho, iudvscf 26 USE io_global, ONLY : stdout 27 USE mp_pools, ONLY : me_pool, root_pool 28 USE klist, ONLY : xk 29 USE el_phon, ONLY: elph_mat, kpq, g_kpq, igqg, xk_gamma 30 USE uspp, ONLY: okvan 31 USE paw_variables, ONLY : okpaw 32 USE uspp_param, ONLY : nhm 33 USE lsda_mod, ONLY : nspin 34 USE lrus, ONLY : int3, int3_nc, int3_paw 35 USE qpoint, ONLY : xq, nksq, ikks 36 ! 37 IMPLICIT NONE 38 ! 39 LOGICAL :: read_dvscf_cart, ascii_dvscf 40 INTEGER :: irr, imode0, ipert, is, ik 41 ! counter on the representations 42 ! counter on the modes 43 ! the change of Vscf due to perturbations 44 COMPLEX(DP), POINTER :: dvscfin(:,:,:), dvscfins (:,:,:) 45 46 47 CALL start_clock ('elphon') 48 49 50 ascii_dvscf=.false. 51 if(elph_mat) read_dvscf_cart=.true. 52 if(read_dvscf_cart) then 53 write(stdout,*) 54 write(stdout,*) 'Reading dvscf in cartesian coordinates !' 55 write(stdout,*) 56 57 u=(0.d0,0.d0) 58 do irr=1,3*nat 59 u(irr,irr)=(1.d0,0.d0) 60 enddo 61 62 63! if(ascii_dvscf) then 64! ALLOCATE (dvrot ( nrxx , nspin , 3*nat) ) 65! fildvscf_asc=trim(tmp_dir)//trim(prefix)//"."//trim(fildvscf)//'1' 66! open(unit=7899,file=fildvscf_asc,status='unknown') 67! dvrot=(0.0,0.0) 68! do na=1,nat 69! do ipol=1,3 70! irr=(na-1)*3+ipol 71! do k = 1, dfftp%nr3 72! do j = 1, dfftp%nr2 73! do i = 1, dfftp%nr1 74! read(7899,*) n, rep,imp 75! dvrot(n,1,irr)=CMPLX(rep,imp,kind=dp) 76! enddo 77! enddo 78! enddo 79! enddo 80! enddo 81! close(7899) 82! endif 83 84 endif 85 86 87 88 ! 89 ! read Delta Vscf and calculate electron-phonon coefficients 90 ! 91 imode0 = 0 92 DO irr = 1, nirr 93 ALLOCATE (dvscfin (dfftp%nnr, nspin_mag , npert(irr)) ) 94 IF (okvan) THEN 95 ALLOCATE (int3 ( nhm, nhm, nat, nspin_mag, npert(irr))) 96 IF (okpaw) ALLOCATE (int3_paw (nhm, nhm, nat, nspin_mag, npert(irr))) 97 IF (noncolin) ALLOCATE(int3_nc( nhm, nhm, nat, nspin, npert(irr))) 98 ENDIF 99 100! if(ascii_dvscf) then 101! DO ipert = 1, npert(irr) 102! dvscfin(1:dfftp%nnr,:,ipert)=dvrot(1:dfftp%nnr,:,imode0+ipert) 103! enddo 104! else 105 DO ipert = 1, npert (irr) 106 CALL davcio_drho ( dvscfin(1,1,ipert), lrdrho, iudvscf, & 107 imode0 + ipert, -1 ) 108 END DO 109! endif 110 IF (doublegrid) THEN 111 ALLOCATE (dvscfins (dffts%nnr, nspin_mag , npert(irr)) ) 112 DO is = 1, nspin_mag 113 DO ipert = 1, npert(irr) 114 CALL fft_interpolate (dfftp, dvscfin(:,is,ipert), dffts, dvscfins(:,is,ipert)) 115 ENDDO 116 ENDDO 117 ELSE 118 dvscfins => dvscfin 119 ENDIF 120 CALL newdq (dvscfin, npert(irr)) 121 CALL elphel_refolded (npert (irr), imode0, dvscfins) 122 ! 123 imode0 = imode0 + npert (irr) 124 IF (doublegrid) DEALLOCATE (dvscfins) 125 DEALLOCATE (dvscfin) 126 127 IF (okvan) THEN 128 DEALLOCATE (int3) 129 IF (okpaw) DEALLOCATE (int3_paw) 130 IF (noncolin) DEALLOCATE(int3_nc) 131 ENDIF 132 133 ENDDO 134 ! 135 ! now read the eigenvalues and eigenvectors of the dynamical matrix 136 ! calculated in a previous run 137 ! 138! IF (.NOT.trans) CALL readmat (iudyn, ibrav, celldm, nat, ntyp, & 139! ityp, omega, amass, tau, xq, w2, dyn) 140 141 IF (.NOT.trans) CALL readmat_findq (iudyn, ibrav, celldm, nat, ntyp, & 142 ityp, omega, amass, tau, xq, w2, dyn) 143 ! 144 145 146 deallocate(xk_gamma) 147 deallocate(kpq,g_kpq,igqg) 148 149! CALL stop_clock ('elphon') 150 RETURN 151END SUBROUTINE ep_matrix_element_wannier 152 153!----------------------------------------------------------------------- 154SUBROUTINE elphsum_wannier(q_index) 155 !----------------------------------------------------------------------- 156 ! 157 ! Sum over BZ of the electron-phonon matrix elements el_ph_mat 158 ! Original routine written by Francesco Mauri 159 ! Adapted to wannier functions by Matteo Calandra 160 ! Dev. Comment: missing calc_sigma_yet 161 !----------------------------------------------------------------------- 162 USE kinds, ONLY : DP 163 USE ions_base, ONLY : nat, ityp, tau,amass,tau, ntyp => nsp, atm 164 USE cell_base, ONLY : at, bg, ibrav, celldm 165 USE symm_base, ONLY : s, sr, irt, nsym, time_reversal, invs, copy_sym, inverse_s 166 USE klist, ONLY : xk, nelec 167 USE wvfct, ONLY : nbnd, et 168 USE el_phon 169 USE mp_pools, ONLY : me_pool, root_pool, inter_pool_comm, npool 170 USE mp_bands, ONLY : intra_bgrp_comm 171 USE io_global, ONLY : stdout,ionode 172 USE io_files, ONLY : prefix 173 USE dynmat, ONLY : dyn, w2 174 USE modes, ONLY : u 175 USE lsda_mod, only : isk,nspin, current_spin,lsda 176 USE mp, ONLY: mp_sum 177 178 USE lr_symm_base, ONLY : irotmq, irgq, gimq, gi 179 USE qpoint, ONLY : xq, nksq, ikks, ikqs 180 USE control_lr, ONLY : lgamma 181 USE noncollin_module, ONLY : noncolin 182 ! 183 IMPLICIT NONE 184 ! 185 LOGICAL :: lborn 186 INTEGER :: q_index 187 ! 188 ! 189 logical :: minus_qloc,sym (48) 190 integer :: nq, imq, isq(48) 191 INTEGER :: ik, ikk, ikq, ibnd, jbnd, ipert, jpert, nu, mu, & 192 ios, iuelphmat,i,j,nt,k 193 INTEGER :: iu_sym,nmodes,nsymq 194 INTEGER :: io_file_unit 195 ! for star_q 196 REAL(DP) :: rtauloc(3,48,nat) 197 ! end of star_q definitions 198 real(DP) :: sxq (3, 48) 199 REAL(DP) xk_dummy(3) 200 character(len=80) :: filelph 201 CHARACTER(len=256) :: file_elphmat 202 ! 203 COMPLEX(DP) :: el_ph_sum (3*nat,3*nat), dyn_corr(3*nat,3*nat) 204 205 INTEGER, EXTERNAL :: find_free_unit 206 CHARACTER (LEN=6), EXTERNAL :: int_to_char 207 208 nmodes=3*nat 209 210 write(filelph,'(A5,f9.6,A1,f9.6,A1,f9.6)') 'elph.',xq(1),'.',xq(2),'.',xq(3) 211 file_elphmat=trim(adjustl(prefix))//'_elph.mat.q_'// TRIM( int_to_char( q_index ) ) 212 213 lborn=.false. 214 ! parallel case: only first node writes 215 IF ( .not.ionode ) THEN 216 iuelphmat = 0 217 ELSE 218 ! 219 ! First I dump information for the electron-phonon interaction 220 ! 221 222 ! 223 iuelphmat = find_free_unit() 224 OPEN (unit = iuelphmat, file = file_elphmat, status = 'unknown', err = & 225 111, iostat = ios, form='unformatted') 226111 CALL errore ('elphsum_wannier', 'opening file'//file_elphmat, ABS (ios) ) 227 REWIND (iuelphmat) 228 xk_dummy(:)=xq(:) 229 call cryst_to_cart(1,xk_dummy,at,-1) 230 WRITE (iuelphmat) xk_dummy 231 WRITE (iuelphmat) noncolin, nspin, lborn 232 WRITE (iuelphmat) nelec 233 WRITE (iuelphmat) elph_nbnd_min,elph_nbnd_max,nbnd 234 WRITE (iuelphmat) nmodes, nksq, nat, ntyp 235 WRITE (iuelphmat) ibrav,(celldm(j),j=1,6) 236 WRITE(iuelphmat) (atm(j),j=1,ntyp),(amass(j),j=1,ntyp), & 237 (ityp(j),j=1,nat),((tau(j,i),j=1,3),i=1,nat) 238 WRITE (iuelphmat) (w2 (nu) , nu = 1, nmodes) 239 WRITE (iuelphmat) ((u(ipert,jpert),ipert=1,nmodes),jpert=1,nmodes) 240 WRITE (iuelphmat) ((dyn(ipert,jpert),ipert=1,3*nat),jpert=1,3*nat) 241 do ik=1,nksq 242 ikk=ikks(ik) 243 ikq=ikqs(ik) 244 xk_dummy(:)=xk(:,ikk) 245 call cryst_to_cart(1,xk_dummy,at,-1) 246 WRITE (iuelphmat) (xk_dummy(ipert),ipert=1,3) 247 WRITE (iuelphmat) (et(ibnd,ikk),ibnd=elph_nbnd_min,elph_nbnd_max) 248 249 do nu=1,nmodes 250 WRITE (iuelphmat) & 251 ((el_ph_mat (jbnd, ibnd, ik, nu),jbnd=elph_nbnd_min,elph_nbnd_max),& 252 ibnd=elph_nbnd_min,elph_nbnd_max) 253 enddo 254 enddo 255 256 257 258 259 ! 260 ! Then I dump symmetry operations 261 ! 262 minus_qloc = .true. 263 sym = .false. 264 sym(1:nsym) = .true. 265 266 call smallg_q (xq, 0, at, bg, nsym, s, sym, minus_qloc) 267 nsymq = copy_sym(nsym, sym) 268 ! recompute the inverses as the order of sym.ops. has changed 269 CALL inverse_s ( ) 270 271 ! part 2: this redoes most of the above, plus it computes irgq, gi, gimq 272 CALL smallgq (xq, at, bg, s, nsym, irgq, nsymq, irotmq, & 273 minus_qloc, gi, gimq) 274 275 sym(1:nsym)=.true. 276 call sgam_ph (at, bg, nsym, s, irt, tau, rtauloc, nat, sym) 277 call star_q(xq, at, bg, nsym , s , invs , nq, sxq, & 278 isq, imq, .FALSE. ) 279 280 281 do j=1,3 282 write(iuelphmat) (at(i,j),i=1,3) 283 enddo 284 do j=1,3 285 write(iuelphmat) (bg(i,j),i=1,3) 286 enddo 287 write(iuelphmat) nsym,nq,imq 288 do i=1,nsym 289 write(iuelphmat) i,invs(i),isq(i) 290 do j=1,3 291 do k=1,3 292 write(iuelphmat) k,j, s(k,j,i) 293 enddo 294 enddo 295 do j=1,nat 296 write(iuelphmat) j, irt(i,j) 297 enddo 298 do j=1,3 299 do k=1,nat 300 write(iuelphmat) j,i, rtauloc(j,i,k) 301 enddo 302 enddo 303 do j=1,3 304 write(iuelphmat) j, sxq(j,i) 305 enddo 306 enddo 307 308 close(iuelphmat) 309 endif 310 311 ! 312 ! 313 314 315 RETURN 316 317 318END SUBROUTINE ELPHSUM_WANNIER 319 320! 321!----------------------------------------------------------------------- 322SUBROUTINE elphel_refolded (npe, imode0, dvscfins) 323 !----------------------------------------------------------------------- 324 ! 325 ! Calculation of the electron-phonon matrix elements el_ph_mat 326 ! <\psi(k+q)|dV_{SCF}/du^q_{i a}|\psi(k)> 327 ! Original routine written by Francesco Mauri 328 ! 329 USE kinds, ONLY : DP 330 USE fft_base, ONLY : dffts 331 USE wavefunctions, ONLY: evc 332 USE io_files, ONLY: prefix, diropn 333 USE klist, ONLY: xk, ngk, igk_k 334 USE lsda_mod, ONLY: lsda, current_spin, isk 335 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 336 USE buffers, ONLY : get_buffer 337 USE wvfct, ONLY: nbnd, npwx 338 USE uspp, ONLY : vkb 339 USE el_phon, ONLY : el_ph_mat, iunwfcwann, igqg, kpq, g_kpq, & 340 xk_gamma, npwq_refolded, lrwfcr 341 USE modes, ONLY : u 342 USE units_ph, ONLY : iubar, lrbar 343 USE control_ph, ONLY : trans 344 USE mp_bands, ONLY: intra_bgrp_comm 345 USE mp_pools, ONLY: me_pool, root_pool 346 USE mp, ONLY: mp_sum 347 USE ions_base, ONLY : nat 348 USE io_global, ONLY : stdout 349 350 USE eqv, ONLY : dvpsi!, evq 351 USE qpoint, ONLY : nksq, ikks, ikqs 352 USE control_lr, ONLY : lgamma 353 USE lrus, ONLY : becp1 354 USE phus, ONLY : alphap 355 356 IMPLICIT NONE 357 ! 358 INTEGER :: npe, imode0 359 COMPLEX(DP) :: dvscfins (dffts%nnr, nspin_mag, npe) 360 COMPLEX(DP), allocatable :: evq(:,:) 361 362 363 ! LOCAL variables 364 logical :: exst 365 INTEGER :: npw, npwq 366 INTEGER :: nrec, ik, ikk, ikq, ikqg,ipert, mode, ibnd, jbnd, ir, ig, & 367 ios 368 369 COMPLEX(DP) , ALLOCATABLE :: aux1 (:,:), elphmat (:,:,:) 370 COMPLEX(DP), EXTERNAL :: zdotc 371 INTEGER, EXTERNAL :: find_free_unit 372 ! 373 allocate (evq(npol*npwx,nbnd)) 374 ALLOCATE (aux1 (dffts%nnr, npol)) 375 ALLOCATE (elphmat ( nbnd , nbnd , 3*nat)) 376 377 378! iunwfcwann=find_free_unit() 379! CALL diropn (iunwfcwann, 'wfc', lrwfc, exst, dvscf_dir) 380! IF (.NOT.exst) THEN 381! CALL errore ('elphel_refolded', 'file '//trim(prefix)//'.wfc not found in Rotated_DVSCF', 1) 382! END IF 383 ! 384 ! Start the loops over the k-points 385 ! 386 387 DO ik = 1, nksq 388 ! 389 ! ik = counter of k-points with vector k 390 ! ikk= index of k-point with vector k 391 ! ikq= index of k-point with vector k+q 392 ! k and k+q are alternated if q!=0, are the same if q=0 393 ! 394 ikk = ikks(ik) 395 ikq = ikqs(ik) 396 ikqg = kpq(ik) 397 npw = ngk(ikk) 398 npwq= ngk(ikq) 399 IF (lsda) current_spin = isk (ikk) 400 ! 401 CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) 402 ! 403 ! read unperturbed wavefuctions psi(k) and psi(k+q) 404 ! 405 evc=(0.d0,0.d0) 406 407! Warning error in reading wfc, this could explain. 408! We read here the wfc at the Gamma point, that is 409! that saved by Wannier. 410 411! CALL davcio (evc, lrwfc, iunwfcwann, ik, - 1) 412! CALL davcio (evq, lrwfc, iunwfcwann, ikqg, - 1) 413 414 415 416! IF (nksq.GT.1) THEN 417! IF (lgamma) THEN 418! CALL davcio (evc, lrwfc, iunwfcwann, ikk, - 1) 419! ELSE 420! CALL davcio (evc, lrwfc, iunwfcwann, ik, - 1) 421! CALL davcio (evq, lrwfc, iunwfcwann, ikqg, - 1) 422! ENDIF 423! ENDIF 424 ! 425 426 call read_wfc_rspace_and_fwfft( evc , ik , lrwfcr , iunwfcwann , npw , igk_k(1,ikk) ) 427 428 429 call calculate_and_apply_phase(ik, ikqg, igqg, npwq_refolded, g_kpq,xk_gamma, evq, .true.) 430 431 432 DO ipert = 1, npe 433 nrec = (ipert - 1) * nksq + ik 434 ! 435 ! dvbare_q*psi_kpoint is read from file (if available) or recalculated 436 ! 437 IF (trans) THEN 438 CALL get_buffer (dvpsi, lrbar, iubar, nrec) 439 ELSE 440 mode = imode0 + ipert 441 ! FIXME : .false. or .true. ??? 442 CALL dvqpsi_us (ik, u (1, mode), .FALSE., becp1, alphap) 443 ENDIF 444 ! 445 ! calculate dvscf_q*psi_k 446 ! 447 448 DO ibnd = 1, nbnd 449 CALL cft_wave (ik, evc(1, ibnd), aux1, +1) 450 CALL apply_dpot(dffts%nnr, aux1, dvscfins(1,1,ipert), current_spin) 451 CALL cft_wave (ik, dvpsi(1, ibnd), aux1, -1) 452 END DO 453 CALL adddvscf (ipert, ik) 454 ! 455 ! calculate elphmat(j,i)=<psi_{k+q,j}|dvscf_q*psi_{k,i}> for this pertur 456 ! 457 DO ibnd =1, nbnd 458 DO jbnd = 1, nbnd 459 elphmat (jbnd, ibnd, ipert) = zdotc (npwq_refolded, evq (1, jbnd), 1, & 460 dvpsi (1, ibnd), 1) 461 IF (noncolin) & 462 elphmat (jbnd, ibnd, ipert) = elphmat (jbnd, ibnd, ipert)+ & 463 zdotc (npwq_refolded, evq(npwx+1,jbnd),1,dvpsi(npwx+1,ibnd), 1) 464 ENDDO 465 ENDDO 466 ENDDO 467 ! 468 CALL mp_sum (elphmat, intra_bgrp_comm) 469 ! 470 ! save all e-ph matrix elements into el_ph_mat 471 ! 472 DO ipert = 1, npe 473 DO jbnd = 1, nbnd 474 DO ibnd = 1, nbnd 475 el_ph_mat (ibnd, jbnd, ik, ipert + imode0) = elphmat (ibnd, jbnd, ipert) 476 ENDDO 477 ENDDO 478 ENDDO 479 ENDDO 480 481 482! CLOSE( UNIT = iunwfcwann, STATUS = 'KEEP' ) 483 ! 484 DEALLOCATE (elphmat) 485 DEALLOCATE (aux1) 486 DEALLOCATE(evq) 487 ! 488 RETURN 489END SUBROUTINE elphel_refolded 490! 491subroutine get_equivalent_kpq(xk,xq,kpq,g_kpq, igqg) 492 !==================================================================! 493 ! ! 494 ! Set up the k+q shell for electron-phonon coupling ! 495 ! ! 496 ! This routine finds the G vectors such that ! 497 ! k+q+G=k' with k and k' belonging to nksq ! 498 ! for each k, the G vector is stored in g_kpq ! 499 ! k'=kpq(ik) ! 500 ! and finally igqg(ik) is the index that allows to find ! 501 ! the g vector g_kpq in the list of all the G vectors ! 502 ! ! 503 ! Matteo Calandra ! 504 !=================================================================== 505 USE kinds, ONLY : DP 506 USE io_global, ONLY : stdout 507 USE cell_base, ONLY : at, bg 508 USE qpoint, ONLY : nksq, ikks 509 USE gvect, ONLY: g, gg 510 USE qpoint, ONLY : nksq 511 USE mp_bands, ONLY : intra_bgrp_comm 512 USE mp, ONLY : mp_sum 513 ! WARNING g_kpq mesh is an integer 514 implicit none 515 516 ! Variables that are private 517 518 integer :: iqx,iqy,iqz,i,j,k,n,nn,iq,ik, ig 519 integer :: kpq(nksq),g_kpq(3,nksq),igqg(nksq) 520 integer, allocatable :: ig_check(:) 521 real(kind=dp) :: gg_ 522 real(kind=dp) :: xq(3), xk(3,*) 523 real(kind=dp) :: xkpq(3),Gvec(3),xq_crys(3) 524 real(kind=dp), allocatable :: xk_crys(:,:) 525 526 ! 527 ! nksq = number of k point per pool withour k+q 528 ! 529 530 ! 531 ! The xk_point entering here must be the k and not 532 ! the k+q 533 ! 534 535 xq_crys=xq 536 537 call cryst_to_cart (1, xq_crys, at, -1) 538 539 540 allocate(xk_crys(3,nksq)) 541 542 do ik=1,nksq 543 xk_crys(1:3,ik)=xk(1:3,ik) 544 enddo 545 call cryst_to_cart (nksq, xk_crys, at, -1) 546 547 ! 548 ! kpt_latt are the BZ vectors in crystalline coordinates 549 ! xq is the q vector in crystalline coordinates 550 ! 551 552 do iq=1,nksq 553 xkpq(:)=xk_crys(:,iq)+xq_crys(:) 554 do i=1,nksq 555 do iqx=-4,4 556 do iqy=-4,4 557 do iqz=-4,4 558 Gvec(1)=real(iqx,dp)+xkpq(1) 559 Gvec(2)=real(iqy,dp)+xkpq(2) 560 Gvec(3)=real(iqz,dp)+xkpq(3) 561 562 if(dabs(xk_crys(1,i)-Gvec(1)).lt.1.d-6.and. & 563 dabs(xk_crys(2,i)-Gvec(2)).lt.1.d-6.and. & 564 dabs(xk_crys(3,i)-Gvec(3)).lt.1.d-6) then 565 kpq(iq)=i 566 g_kpq(1,iq)=-iqx 567 g_kpq(2,iq)=-iqy 568 g_kpq(3,iq)=-iqz 569 goto 99 570 endif 571 enddo 572 enddo 573 enddo 574 enddo 575 CALL errore ('get_equivalent_kpq', 'cannot find index k+q ', 2 ) 576 stop 57799 continue 578 enddo 579 580 ! 581 ! here between all the g-vectors I find the index of that 582 ! related to the translation in the Brillouin zone. 583 ! Warning: if G does not belong to the processor igqg is zero. 584 ! 585 586 igqg=0 587 do ik=1,nksq 588 Gvec(:) = REAL( g_kpq(:,ik),dp ) 589 call cryst_to_cart (1, Gvec, bg, 1) 590 gg_ = Gvec(1)*Gvec(1) + Gvec(2)*Gvec(2) + Gvec(3)*Gvec(3) 591 igqg(ik)=0 592 ig=1 593 do while (gg(ig) <= gg_ + 1.d-6) 594 if ( (abs(g(1,ig)-Gvec(1)) < 1.d-6) .and. & 595 (abs(g(2,ig)-Gvec(2)) < 1.d-6) .and. & 596 (abs(g(3,ig)-Gvec(3)) < 1.d-6) ) then 597 igqg(ik) = ig 598 599 endif 600 ig= ig +1 601 end do 602 end do 603 604 allocate(ig_check(nksq)) 605 ig_check=igqg 606 CALL mp_sum( ig_check, intra_bgrp_comm ) 607 do ik=1,nksq 608 if(ig_check(ik).eq.0) & 609 CALL errore('get_equivalent_kpq', & 610 ' g_kpq vector is not in the list of Gs', 100*ik ) 611 612 enddo 613 deallocate(xk_crys) 614 615 616end subroutine get_equivalent_kpq 617 618subroutine calculate_and_apply_phase(ik, ikqg, igqg, npwq_refolded, g_kpq, xk_gamma, evq, lread) 619 USE kinds, ONLY : DP 620 USE fft_base, ONLY : dffts 621 USE fft_interfaces, ONLY : fwfft, invfft 622 USE wvfct, ONLY: nbnd, npwx 623 USE gvect, ONLY : ngm, g 624 USE gvecw, ONLY : gcutw 625 USE cell_base, ONLY : bg 626 USE qpoint, ONLY : nksq 627 USE wavefunctions, ONLY : evc 628 USE noncollin_module, ONLY : npol, noncolin 629 USE el_phon, ONLY:iunwfcwann, lrwfcr 630 631 IMPLICIT NONE 632 633 LOGICAL :: lread 634 INTEGER :: ik, ikqg, npwq_refolded 635 INTEGER :: igqg(nksq) 636 INTEGER :: g_kpq(3,nksq) 637 REAL (DP) :: xk_gamma(3,nksq) 638 complex(dp) :: evq(npwx*npol,nbnd) 639! internal 640 INTEGER :: npw_, m,i 641 INTEGER, allocatable :: igk_(:), igkq_(:) 642 REAL(DP) :: xkqg(3), g_(3), g_scra(3,ngm) 643 REAL(DP), ALLOCATABLE :: gk(:) 644 COMPLEX (DP), allocatable :: psi_scratch(:) 645 complex(DP), allocatable :: phase(:) 646 647 allocate(igk_(npwx), igkq_(npwx), gk(npwx) ) 648 allocate (psi_scratch ( dffts%nnr) ) 649 allocate (phase(dffts%nnr)) 650 FLUSH (6) 651 652 g_scra=g 653 654 g_(:)=real( g_kpq(:,ik), dp ) 655 call cryst_to_cart (1, g_, bg, 1) 656 xkqg(:)=xk_gamma(:,ikqg)+g_(:) 657 658 npw_=0 659 npwq_refolded=0 660 igk_=0 661 igkq_=0 662 663 664 call gk_sort (xk_gamma(1,ikqg), ngm, g_scra, gcutw, npw_, igk_, gk) 665 666 if(lread) then 667 call read_wfc_rspace_and_fwfft( evq , ikqg , lrwfcr , iunwfcwann , npw_ , igk_ ) 668 endif 669 670 call gk_sort (xkqg, ngm, g_scra, gcutw, npwq_refolded, igkq_, gk) 671 672 phase(:) = (0.d0,0.d0) 673 674 if ( igqg(ik)>0) then 675 phase( dffts%nl(igqg(ik)) ) = (1.d0,0.d0) 676 endif 677 678 679 CALL invfft ('Wave', phase, dffts) 680 ! call cft3s (phase, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2) 681 phase(:)=conjg(phase(:)) 682 683 684 if(npwq_refolded.ne.npw_) call errore('calculate_and_apply_phase', 'Warning : npwq_refolded \= npw_',-1) 685 686 do m=1,nbnd 687 psi_scratch = (0.d0, 0.d0) 688 psi_scratch(dffts%nl (igk_ (1:npw_) ) ) = evq (1:npw_, m) 689! psi_scratch(dffts%nl (igk_ (1:npw) ) ) = evq (1:npw, m) 690 CALL invfft ('Wave', psi_scratch, dffts) 691 ! call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2) 692 psi_scratch(1:dffts%nnr) = psi_scratch(1:dffts%nnr) * phase(1:dffts%nnr) 693 ! call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2) 694 CALL fwfft ('Wave', psi_scratch, dffts) 695 evq(1:npwq_refolded,m) = psi_scratch(dffts%nl (igkq_(1:npwq_refolded) ) ) 696 enddo 697 698 if(noncolin) then 699 do m=1,nbnd 700 psi_scratch = (0.d0, 0.d0) 701 psi_scratch(dffts%nl (igk_ (1:npw_) ) ) = evq (npwx+1:npwx+npw_, m) 702 ! psi_scratch(dffts%nl (igk_ (1:npw) ) ) = evq (1:npw, m) 703 CALL invfft ('Wave', psi_scratch, dffts) 704 ! call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2) 705 psi_scratch(1:dffts%nnr) = psi_scratch(1:dffts%nnr) * phase(1:dffts%nnr) 706 ! call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2) 707 CALL fwfft ('Wave', psi_scratch, dffts) 708 ! evq(npwx+1:npwx+npwq_refolded,m) = psi_scratch(dffts%nl (igkq_(1:npwq_refolded) ) ) 709 evq((npwx+1):(npwx+npwq_refolded),m) = psi_scratch(dffts%nl (igkq_(1:npwq_refolded) ) ) 710 enddo 711 endif 712 713 deallocate(psi_scratch) 714 DEALLOCATE(phase) 715 deallocate(gk, igk_, igkq_) 716 717 return 718end subroutine calculate_and_apply_phase 719 720 721!----------------------------------------------------------------------- 722SUBROUTINE readmat_findq (iudyn, ibrav, celldm, nat, ntyp, ityp, omega, & 723 amass, tau, q, w2, dyn) 724 !----------------------------------------------------------------------- 725 ! 726 USE kinds, ONLY : DP 727 USE constants, ONLY : amu_ry 728 USE control_ph, ONLY : xmldyn 729 USE output, ONLY : fildyn 730 USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, & 731 read_dyn_mat, read_dyn_mat_tail 732 IMPLICIT NONE 733 ! Input 734 INTEGER :: iudyn, ibrav, nat, ntyp, ityp (nat) 735 REAL(DP) :: celldm (6), amass (ntyp), tau (3, nat), q (3), & 736 omega 737 ! output 738 REAL(DP) :: w2 (3 * nat) 739 COMPLEX(DP) :: dyn (3 * nat, 3 * nat) 740 ! local (control variables) 741 INTEGER :: ntyp_, nat_, ibrav_, ityp_ 742 REAL(DP) :: celldm_ (6), amass_, tau_ (3), q_ (3) 743 ! local 744 INTEGER :: nspin_mag, nqs 745 REAL(DP) :: at(3,3) 746 REAL(DP) :: bg(3,3) 747 REAL(DP) :: m_loc(3,nat) 748 INTEGER :: ityp__ (nat) 749 REAL(DP) :: amass__ (ntyp) 750 INTEGER :: iq 751 REAL(DP) :: xq(3) 752 COMPLEX(DP) :: u(3*nat,3*nat) 753 REAL(DP) :: dynr (2, 3, nat, 3, nat), err_q(3) 754 COMPLEX(DP) :: dynr_c(3,3,nat,nat) 755 756 CHARACTER(len=80) :: line 757 CHARACTER(len=3) :: atm 758 INTEGER :: nt, na, nb, naa, nbb, nu, mu, i, j 759 LOGICAL :: lfound 760 761 ! 762 ! 763 IF(xmldyn) THEN 764 CALL read_dyn_mat_param(fildyn, ntyp_, nat_ ) 765 CALL read_dyn_mat_header(ntyp_, nat_, ibrav_, nspin_mag, & 766 celldm_, at, bg, omega, atm, amass__, tau_, ityp__, m_loc, & 767 nqs ) 768 IF ( ntyp.NE.ntyp_ .OR. nat.NE.nat_ .OR.ibrav_.NE.ibrav .OR. & 769 ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5) & 770 CALL errore ('readmat', 'inconsistent data a', 1) 771 DO nt = 1, ntyp 772 IF ( ABS (amass__ (nt) - amass (nt) ) > 1.0d-5) & 773 CALL errore ( 'readmat', 'inconsistent data b', 1 + nt) 774 ENDDO 775 DO na = 1, nat 776 IF (ityp__ (na).NE.ityp (na) ) CALL errore ('readmat', & 777 'inconsistent data c', na) 778 ENDDO 779 780 ELSE 781 REWIND (iudyn) 782 READ (iudyn, '(a)') line 783 READ (iudyn, '(a)') line 784 READ (iudyn, * ) ntyp_, nat_, ibrav_, celldm_ 785 IF ( ntyp.NE.ntyp_ .OR. nat.NE.nat_ .OR.ibrav_.NE.ibrav .OR. & 786 ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5) & 787 CALL errore ('readmat', 'inconsistent data', 1) 788 DO nt = 1, ntyp 789 READ (iudyn, * ) i, atm, amass_ 790 IF ( nt.NE.i .OR. ABS (amass_ - amu_ry*amass (nt) ) > 1.0d-5) & 791 CALL errore ( 'readmat', 'inconsistent data', 1 + nt) 792 ENDDO 793 DO na = 1, nat 794 READ (iudyn, * ) i, ityp_, tau_ 795 IF (na.NE.i.OR.ityp_.NE.ityp (na) ) CALL errore ('readmat', & 796 'inconsistent data', 10 + na) 797 ENDDO 798 799 ENDIF 800 801 lfound=.false. 802 iq=0 803 804 do while(.not.lfound) 805 806 IF(xmldyn) THEN 807 808 iq = iq+1 809 CALL read_dyn_mat(nat,iq,xq,dynr_c) 810 ! CALL read_dyn_mat_tail(nat,omega,u) 811 err_q(1:3)=dabs(xq(1:3)-q(1:3)) 812 813 if(err_q(1).lt.1.d-7.and.err_q(2).lt.1.d-7.and.err_q(3).lt.1.d-7) lfound=.true. 814 815 DO nb = 1, nat 816 DO j = 1, 3 817 DO na = 1, nat 818 DO i = 1, 3 819 dynr (1, i, na, j, nb) = REAL(dynr_c(i, j, na, nb)) 820 dynr (2, i, na, j, nb) = AIMAG(dynr_c(i, j, na, nb)) 821 ENDDO 822 ENDDO 823 ENDDO 824 ENDDO 825 826 ELSE 827 READ (iudyn, '(a)') line 828 READ (iudyn, '(a)') line 829 READ (iudyn, '(a)') line 830 READ (iudyn, '(a)') line 831 832 READ (line (11:80), * ) (q_ (i), i = 1, 3) 833 834 err_q(1:3)=dabs(q_(1:3)-q(1:3)) 835 836 if(err_q(1).lt.1.d-7.and.err_q(2).lt.1.d-7.and.err_q(3).lt.1.d-7) lfound=.true. 837 838 839 840 READ (iudyn, '(a)') line 841 DO na = 1, nat 842 DO nb = 1, nat 843 READ (iudyn, * ) naa, nbb 844 IF (na.NE.naa.OR.nb.NE.nbb) CALL errore ('readmat', 'error reading & 845 &file', nb) 846 READ (iudyn, * ) ( (dynr (1, i, na, j, nb), dynr (2, i, na, j, nb) & 847 , j = 1, 3), i = 1, 3) 848 ENDDO 849 ENDDO 850 851 ENDIF 852 853 if(lfound) then 854 ! 855 ! divide the dynamical matrix by the (input) masses (in amu) 856 ! 857 DO nb = 1, nat 858 DO j = 1, 3 859 DO na = 1, nat 860 DO i = 1, 3 861 dynr (1, i, na, j, nb) = dynr (1, i, na, j, nb) / SQRT (amass ( & 862 ityp (na) ) * amass (ityp (nb) ) ) / amu_ry 863 dynr (2, i, na, j, nb) = dynr (2, i, na, j, nb) / SQRT (amass ( & 864 ityp (na) ) * amass (ityp (nb) ) ) / amu_ry 865 ENDDO 866 ENDDO 867 ENDDO 868 ENDDO 869 ! 870 ! solve the eigenvalue problem. 871 ! NOTA BENE: eigenvectors are overwritten on dyn 872 ! 873 CALL cdiagh (3 * nat, dynr, 3 * nat, w2, dyn) 874 ! 875 ! divide by sqrt(mass) to get displacements 876 ! 877 DO nu = 1, 3 * nat 878 DO mu = 1, 3 * nat 879 na = (mu - 1) / 3 + 1 880 dyn (mu, nu) = dyn (mu, nu) / SQRT ( amu_ry * amass (ityp (na) ) ) 881 ENDDO 882 ENDDO 883 ! 884 ! 885 endif 886enddo 887 888 889 RETURN 890END SUBROUTINE readmat_findq 891 892