TABLE OF CONTENTS
ABINIT/read_gkk [ Functions ]
NAME
read_gkk
FUNCTION
This routine reads in elphon matrix elements and completes them using the appropriate symmetries
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer, MG) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
elph_ds = datastructure containing elphon matrix elements Cryst<crystal_t>=Info on the crystal unit cell. Ifc<ifc_type>=Object containing the interatomic force constants. FSfullpqtofull = mapping of k+q to k n1wf = number of 1WF files to be read and analyzed nband = number of bands per kpoint unitgkk = unit of GKK file for reading
OUTPUT
elph_ds = modified gkq
gkk_qpt = el-ph matrix elements for irreducible qpoints and
kpoints (as a function of the reduced symmetry for the qpoint)
gkk_flag = flag array:
-1 -> element is missing
0 -> element is from symmetric qpt (Now done in complete_gkk)
1 -> element is from symmetric pert
2 -> element is kptsym of gkk file
3 -> element was read from gkk file
PARENTS
get_all_gkq
CHILDREN
completeperts,get_rank_1kpt,hdr_bcast,hdr_fort_read,hdr_free,ifc_fourq littlegroup_pert,littlegroup_q,mati3inv,normsq_gkq,phdispl_cart2red prt_gkk_yambo,wrap2_pmhalf,wrtout,xmpi_bcast
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine read_gkk(elph_ds,Cryst,ifc,Bst,FSfullpqtofull,gkk_flag,n1wf,nband,ep_prt_yambo,unitgkk) 55 56 use defs_basis 57 use defs_datatypes 58 use defs_abitypes 59 use defs_elphon 60 use m_errors 61 use m_profiling_abi 62 use m_xmpi 63 use m_kptrank 64 use m_hdr 65 66 use m_numeric_tools, only : wrap2_pmhalf 67 use m_geometry, only : phdispl_cart2red 68 use m_crystal, only : crystal_t 69 use m_ifc, only : ifc_type, ifc_fourq 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'read_gkk' 75 use interfaces_14_hidewrite 76 use interfaces_32_util 77 use interfaces_41_geometry 78 use interfaces_77_ddb, except_this_one => read_gkk 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: n1wf,nband,unitgkk,ep_prt_yambo 86 type(crystal_t),intent(in) :: Cryst 87 type(ifc_type),intent(in) :: ifc 88 type(ebands_t),intent(in) :: Bst 89 type(elph_type),intent(inout) :: elph_ds 90 !arrays 91 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 92 integer,intent(out) :: gkk_flag(elph_ds%nbranch,elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol,elph_ds%nqpt_full) 93 94 !Local variables------------------------------- 95 !scalars 96 integer :: nsppol,nbranch,nFSband,minFSband,comm,use_sym 97 integer :: fform,i1wf,ikpt_phon,iatom1,iatom2 98 integer :: ib,ib1,ib2,ibb,ibranch,idir,idir1,idir2,ierr,ii,ikpt1 99 integer :: ipert,ipert1,ipert2,iqptirred,iqptfull,isppol,isym1 100 integer :: itim1,jkpt_phon,new 101 integer :: nsym1,qtimrev,syuse 102 integer :: tdonecompl,test_flag,verify 103 integer :: nqptirred_local 104 integer :: master, me 105 integer :: symrankkpt, ikpt1_phon, ik_this_proc 106 real(dp) :: res,ss,timsign 107 character(len=500) :: msg 108 type(hdr_type) :: hdr1 109 !arrays 110 integer :: FSirrtok(3,elph_ds%k_phon%nkpt) 111 integer :: symaf1(Cryst%nsym),symq(4,2,Cryst%nsym) 112 integer :: symrc1(3,3,Cryst%nsym),symrl1(3,3,Cryst%nsym) 113 integer :: tmpflg(3,Cryst%natom+2,3,Cryst%natom+2) 114 real(dp) :: displ_cart(2,3*Cryst%natom,3*Cryst%natom) 115 real(dp) :: displ_red(2,3*Cryst%natom,3*Cryst%natom) 116 real(dp) :: eigvec(2,3*Cryst%natom,3*Cryst%natom),kpt(3),phfrq_tmp(3*Cryst%natom),redkpt(3) 117 real(dp) :: qptirred_local(3,n1wf) 118 real(dp) :: tnons1(3,Cryst%nsym) 119 real(dp),allocatable :: eigen1(:,:,:),gkk_qpt_tmp(:,:,:,:) 120 real(dp),allocatable :: h1_mat_el(:,:,:,:,:),h1_mat_el_sq(:,:,:,:,:) 121 real(dp),allocatable :: qdata(:,:,:),qdata_tmp(:,:,:,:) 122 123 ! ************************************************************************* 124 125 ABI_UNUSED(Bst%bantot) 126 127 use_sym = 1 128 nsppol = elph_ds%nsppol 129 nbranch = elph_ds%nbranch 130 if (ep_prt_yambo==1) then 131 nFSband = nband 132 minFSband = 1 133 else 134 nFSband = elph_ds%nFSband 135 minFSband = elph_ds%minFSband 136 end if 137 138 !init values for parallelization 139 comm = xmpi_world 140 me = xmpi_comm_rank(comm) 141 master = 0 142 143 ABI_STAT_ALLOCATE(h1_mat_el,(2, nFSband**2, nbranch, elph_ds%k_phon%my_nkpt, nsppol), ierr) 144 ABI_CHECK(ierr==0, 'trying to allocate array h1_mat_el') 145 h1_mat_el= zero 146 147 ABI_STAT_ALLOCATE(h1_mat_el_sq,(2, nFSband**2, nbranch**2,elph_ds%k_phon%my_nkpt, nsppol), ierr) 148 ABI_CHECK(ierr==0, 'trying to allocate array h1_mat_el_sq') 149 h1_mat_el_sq = zero 150 151 ABI_ALLOCATE(elph_ds%qirredtofull,(elph_ds%nqptirred)) 152 153 !MG array to store the e-ph quantities calculated over the input Q-grid 154 ABI_ALLOCATE(qdata_tmp,(elph_ds%nqptirred,nbranch,nsppol,3)) 155 qdata_tmp=zero 156 157 nqptirred_local=0 !zero number of irred q-points found 158 qptirred_local(:,:)=zero 159 160 gkk_flag = -1 161 162 if (elph_ds%gkqwrite ==0) then 163 elph_ds%gkk_qpt = zero 164 165 else if (elph_ds%gkqwrite == 1) then 166 ABI_STAT_ALLOCATE(gkk_qpt_tmp,(2,elph_ds%ngkkband**2,nbranch**2,nsppol), ierr) 167 ABI_CHECK(ierr==0, 'trying to allocate array gkk_qpt_tmp') 168 gkk_qpt_tmp = zero 169 do iqptirred=1,elph_ds%nqptirred*elph_ds%k_phon%nkpt 170 write (elph_ds%unitgkq,REC=iqptirred) gkk_qpt_tmp 171 end do 172 ABI_DEALLOCATE(gkk_qpt_tmp) 173 174 else 175 write (msg,'(a,i0)')' Wrong values for gkqwrite = ',elph_ds%gkqwrite 176 MSG_BUG(msg) 177 end if !gkqwrite 178 179 !=========================================================== 180 !Loop over all files we have 181 !read in header for perturbation 182 !should check that all files are complete, have same header 183 !(taking into account the symmetries for the qpoint), 184 !represent the correct qpoints ... 185 !MG: this task should be performed in mrggkk 186 !=========================================================== 187 188 ABI_ALLOCATE(eigen1,(2,nband,nband)) 189 do i1wf=1,n1wf 190 191 if (master == me) then 192 write (msg,'(2a,i4,a,i4)')ch10,' read_gkk : reading 1WF header # ',i1wf,' /',n1wf 193 call wrtout(std_out,msg,'COLL') 194 195 ! Could check for compatibility of natom, kpt grids, ecut, qpt with DDB grid... 196 ! MG: Also this task should be done in mrggkk 197 198 call hdr_fort_read(hdr1, unitgkk, fform) 199 if (fform == 0) then 200 write (msg,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0' 201 MSG_ERROR(msg) 202 end if 203 204 write(msg,'(a,i4)')' read_gkk : have read 1WF header #',i1wf 205 call wrtout(std_out,msg,'COLL') 206 write (msg,'(2a,i4,a)')ch10,' read_gkk : # of kpt for this perturbation: ',hdr1%nkpt,ch10 207 call wrtout(std_out,msg,'COLL') 208 209 end if 210 211 ! broadcast data to all nodes: 212 call hdr_bcast(hdr1, master, me, comm) 213 214 ! Find qpoint in full grid 215 new=1 216 do iqptfull=1,elph_ds%nqpt_full 217 kpt(:) = hdr1%qptn(:) - elph_ds%qpt_full(:,iqptfull) 218 call wrap2_pmhalf(kpt(1),redkpt(1),res) 219 call wrap2_pmhalf(kpt(2),redkpt(2),res) 220 call wrap2_pmhalf(kpt(3),redkpt(3),res) 221 ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2 222 if(ss < tol6) then 223 new = 0 224 exit !exit with iqptfull 225 end if 226 end do !iqptfull 227 228 if (new == 1) then 229 ! Test should be at the end: dont care if there are additional 230 ! qpts in gkk file which are not on the main grid. Ignore them. 231 write (msg,'(4a,3es16.6,2a)')ch10,& 232 & ' read_gkk : WARNING- ',ch10,& 233 & ' qpoint = ',hdr1%qptn(:),ch10,& 234 & ' not found in the input q-grid. Ignoring this point ' 235 call wrtout(ab_out,msg,'COLL') 236 call wrtout(std_out,msg,'COLL') 237 if (me == master) then 238 do isppol=1,hdr1%nsppol 239 do ikpt1=1,hdr1%nkpt 240 read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband) 241 end do 242 end do 243 end if 244 245 cycle !cycle the loop on i1wf 246 end if !end if (new ==1) 247 248 249 ! Check whether other pieces of the DDB have used this qpt already 250 new=1 251 do iqptirred=1,nqptirred_local 252 kpt(:) = qptirred_local(:,iqptirred) - hdr1%qptn(:) 253 call wrap2_pmhalf(kpt(1),redkpt(1),res) 254 call wrap2_pmhalf(kpt(2),redkpt(2),res) 255 call wrap2_pmhalf(kpt(3),redkpt(3),res) 256 ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2 257 if(ss < tol6) then 258 new=0 259 exit !MG We can use this information to avoid recalculating the dynamical matrix 260 end if !but we need to use a fixed format in GKK! 261 end do !iqptirred 262 263 if (new==1) then !we have a new valid irreducible qpoint, add it! 264 nqptirred_local = nqptirred_local+1 265 if (nqptirred_local > elph_ds%nqptirred) then 266 write (msg, '(a,a,a,i6,i6)') & 267 & 'found too many qpoints in GKK file wrt anaddb input ', ch10, & 268 & 'nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local 269 MSG_ERROR(msg) 270 end if 271 qptirred_local(:,nqptirred_local) = hdr1%qptn(:) 272 iqptirred = nqptirred_local 273 tdonecompl = 0 274 h1_mat_el = zero 275 end if 276 277 ! now iqptirred is the index of the present qpoint in the array qptirred_local 278 ! and iqptfull is the index in the full qpt_full array for future reference 279 elph_ds%qirredtofull(iqptirred) = iqptfull 280 281 write (msg,'(a,i5,a,3es16.8)')& 282 & ' read_gkk : full zone qpt number ',iqptfull,' is ',elph_ds%qpt_full(:,iqptfull) 283 call wrtout(std_out,msg,'COLL') 284 285 ! if this perturbation has already been filled (overcomplete gkk) 286 ! check only 1st kpoint and spinpol, then check others 287 verify = 0 288 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) /= -1) then 289 ! 290 do isppol=1,nsppol 291 do ik_this_proc=1,elph_ds%k_phon%my_nkpt 292 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 293 write (std_out,*)" hdr1%pertcase,ik_this_proc,iqptirred",hdr1%pertcase,ik_this_proc,iqptirred 294 MSG_ERROR('Partially filled perturbation ') 295 end if 296 end do ! ikpt_phon 297 end do ! isppol 298 ! 299 MSG_WARNING(' gkk perturbation is already filled') 300 write(std_out,*)' hdr1%pertcase,iqptirred,iqptfull = ',hdr1%pertcase,iqptirred,iqptfull,& 301 & gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) 302 verify = 1 303 write (125,*) '# matrix elements for symmetric perturbation' 304 ! Instead of reading eigen1 into void, verify == 1 checks them later on wrt values in memory 305 end if !gkk_flag 306 307 ! Examine the symmetries of the q wavevector 308 ! these will be used to complete the perturbations for other atoms and idir 309 if (ep_prt_yambo==1) then 310 ! If one wants to print GKKs along phonon modes, it mean mixing of 311 ! perturbations with differnt jauge. Symmetries must then be disable. 312 call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0,use_sym=0) 313 else 314 call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0) 315 end if 316 317 ! Determine dynamical matrix, phonon frequencies and displacement vector for qpoint 318 !call wrtout(std_out,' read_gkk: calling inpphon to calculate the dynamical matrix','COLL') 319 320 call ifc_fourq(ifc,cryst,qptirred_local(:,iqptirred),phfrq_tmp,displ_cart,out_eigvec=eigvec) 321 322 ! Get displacement vectors for all branches in reduced coordinates 323 ! used in scalar product with H(1)_atom,idir matrix elements 324 ! Calculate $displ_red = displ_cart \cdot gprimd$ for each phonon branch 325 326 call phdispl_cart2red(Cryst%natom,Cryst%gprimd,displ_cart,displ_red) 327 328 ! prefactors for gk+q,n\prime;k,n matrix element 329 ! COMMENT : in decaft there is a weird term in the mass factor, of M-zval(species) 330 ! dont know why. Not needed to reproduce decaft results, though... 331 ! weight is squared in evaluation of 332 ! gamma_{k,q,j} = 2 \pi omega_{q,j} sum_{nu,nu\prime} |g^{q,j}_{k+q,nu\prime; k,nu}|^2 333 ! normally cancels with the 2 \pi omega_{q,j} factor in front of the sum... 334 335 ! hdr1%pertcase = idir + (ipert-1)*3 where ipert=iatom in the interesting cases 336 idir = mod (hdr1%pertcase-1,3)+1 337 ipert = int(dble(hdr1%pertcase-idir)/three)+1 338 339 write (msg,'(4a,i3,a,i3,a,i4,a)')ch10,& 340 & ' read_gkk : calling littlegroup_pert to examine the symmetries of the full perturbation ',ch10,& 341 & ' idir = ',idir,' ipert = ',ipert,' and Q point = ',iqptirred,ch10 342 call wrtout(std_out,msg,'COLL') 343 344 ! Examine the symmetries of the full perturbation these will be used to complete the kpoints 345 ! DOESNT USE TIME REVERSAL IN littlegroup_pert except for gamma 346 347 syuse=0 348 349 call littlegroup_pert(Cryst%gprimd,idir,Cryst%indsym,ab_out,ipert,Cryst%natom,Cryst%nsym,nsym1,2,Cryst%symafm,symaf1,& 350 & symq,Cryst%symrec,Cryst%symrel,symrl1,syuse,Cryst%tnons,tnons1) 351 352 do isym1=1,nsym1 353 call mati3inv(symrl1(:,:,isym1),symrc1(:,:,isym1)) 354 end do 355 FSirrtok = 0 356 357 ! ======================================================== 358 ! Loop over irred kpts in file, and fill the default gkk 359 ! ======================================================== 360 361 ! MG NOTE : in the present implementation, if nsppol /=1 the code stops in rchkGSheader! 362 do isppol=1,hdr1%nsppol !Loop over spins is trivial? Not tested. 363 write (std_out,*) ' read_gkk : isppol = ', isppol 364 365 do ikpt1=1,hdr1%nkpt !Loop over irred kpoints, WARNING nkpt depends on qpoint and symmetry! 366 ! 367 ! this is the main read of the gkk matrix elements from the file (eigen1 arrays) 368 ! it has to be done exactly nsppol*nkpt times, and the kpt_phon are completed 369 ! where appropriate in the loop below (normally succeeding only once for each kpt) 370 ! 371 if (master == me) then 372 read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband) 373 end if 374 375 ! MPI broadcast data to all nodes: 376 call xmpi_bcast(eigen1, master, comm, ierr) 377 378 ! find place of irred k in k_phon 379 ! the kpoints in the file (kptns) could be ordered arbitrarily 380 call get_rank_1kpt (hdr1%kptns(:,ikpt1)-qptirred_local(:,iqptirred), & 381 & symrankkpt, elph_ds%k_phon%kptrank_t) 382 ikpt1_phon = elph_ds%k_phon%kptrank_t%invrank(symrankkpt) 383 if (ikpt1_phon < 0) then 384 write (msg,'(a,3es16.6,a)')& 385 & ' irred k ',hdr1%kptns(:,ikpt1),' was not found in full grid' 386 MSG_ERROR(msg) 387 end if 388 ! find correspondence between this kpt_phon and the others 389 ! symrc1 conserves perturbation as well as qpoint 390 ! add to FSirrtok list 391 do isym1=1,nsym1 392 do itim1=0,qtimrev 393 timsign=one-two*itim1 394 kpt(:) = timsign*matmul(symrc1(:,:,isym1), elph_ds%k_phon%kpt(:,ikpt1_phon)) 395 396 call get_rank_1kpt (kpt,symrankkpt,elph_ds%k_phon%kptrank_t) 397 jkpt_phon = elph_ds%k_phon%kptrank_t%invrank(symrankkpt) 398 399 if (jkpt_phon > 0) then 400 FSirrtok(1,jkpt_phon) = ikpt1_phon 401 FSirrtok(2,jkpt_phon) = isym1 402 FSirrtok(3,jkpt_phon) = itim1 403 else 404 write (msg,'(a,3es16.6,a,i5,a,i4,a)')& 405 & ' sym equivalent of kpt ',hdr1%kptns(:,ikpt1),' by sym ',& 406 & isym1,' and itime ',itim1,' was not found' 407 MSG_ERROR(msg) 408 end if 409 end do !itim1 410 end do !isim1 411 412 413 ! 414 ! Here check if the symmetry-copied gkk at new k point is equal to the one found in the file for non-irreducible point 415 ! NB This is DEBUG code 416 ! 417 if (verify == 1 .and. elph_ds%k_phon%my_kpt(ikpt1_phon) == me) then 418 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 419 if (elph_ds%k_phon%my_ikpt(ik_this_proc) == ikpt1_phon) exit 420 end do 421 do ib1=1,nFSband 422 do ib2=1,nFSband 423 ibb = (ib1-1)*nFSband+ib2 424 write (125,'(2(2E16.6,2x))') h1_mat_el(:,ibb,hdr1%pertcase,ik_this_proc,isppol),& 425 & eigen1(:,minFSband-1+ib2,minFSband-1+ib1) 426 end do 427 end do 428 end if !verify end DEBUG code 429 430 431 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 432 ! should I be dealing with this k-point? 433 jkpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 434 435 ! does present ikpt1 contribute to this k-point? 436 if (FSirrtok(1,jkpt_phon) /= ikpt1_phon) cycle 437 438 ! if this kpoint has already been filled (overcomplete gkk) 439 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) /= -1) then 440 MSG_WARNING("gkk element is already filled") 441 write(std_out,*)' hdr1%pertcase,ik_this_proc,isppol,iqptirred = ',& 442 & hdr1%pertcase,ik_this_proc,isppol,iqptirred,& 443 & gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) 444 ! exit 445 end if !gkk_flag 446 447 ! =============================================================== 448 ! TODO: if there is a phase factor in swapping k-points, insert it here in copy to h1_mat_el 449 ! as a function of symops in FSirrtok 450 ! complete gkk for symmetric ikpt_phon with sym1 which conserve 451 ! the full perturbation+qpoint 452 ! Not tested explicitly, but the results for Pb using reduced kpts look good 453 ! should do same RF calculation with nsym=1 and check 454 ! =============================================================== 455 456 ! save this kpoint 457 do ib1=1,nFSband 458 do ib2=1,nFSband 459 ibb = (ib1-1)*nFSband+ib2 460 461 ! real 462 res=eigen1(1,minFSband-1+ib2,minFSband-1+ib1) 463 h1_mat_el(1,ibb,hdr1%pertcase,ik_this_proc,isppol) = res 464 465 ! imag 466 res=eigen1(2,minFSband-1+ib2,minFSband-1+ib1) 467 h1_mat_el(2,ibb,hdr1%pertcase,ik_this_proc,isppol) = res 468 end do !ib2 469 end do !ib1 470 ! if jkpt is equal to ikpt1_phon (if clause above) flag == 3 471 if (FSirrtok(2,jkpt_phon) == 1) then 472 gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 3 473 ! if jkpt_phon comes from ikpt1_phon flag == 2 with some symop 474 else 475 gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 2 476 end if 477 478 end do !jkpt_phon 479 480 ! =============================================================== 481 ! we now have contribution to g(k+q,k; \kappa,\alpha) from one 482 ! kpoint,and one perturbation, 483 ! NB: each perturbation will contribute to all the modes later! 484 ! 485 ! SHOULD ONLY DO THIS FOR THE SYMS NEEDED 486 ! TO COMPLETE THE PERTURBATIONS!!! 487 ! ================================================================ 488 489 end do !ikpt1 490 end do !isppol 491 492 ! 14 Jan 2014 removed test on verify - in new scheme full BZ is read in and should be used to avoid phase errors 493 ! if (verify == 1) cycle 494 495 ! Checks on irred grid provided and on gkk_flag accumulated up to now 496 if (elph_ds%tuniformgrid == 1) then ! check if irred kpoints found reconstitute the FS kpts 497 do ikpt_phon=1,elph_ds%k_phon%nkpt 498 if (FSirrtok(1,ikpt_phon) == 0) then 499 write(msg,'(a,3es16.6,2a)')& 500 & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,& 501 & ' is not the symmetric of one of those found in the GKK file' 502 MSG_ERROR(msg) 503 end if 504 end do !ikpt_phon 505 506 ! normally at this point we have used all the gkk for all kpoints on the FS 507 ! for the given irred perturbation: check 508 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 509 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 510 511 if (gkk_flag(hdr1%pertcase, hdr1%pertcase, ik_this_proc, 1, elph_ds%qirredtofull(iqptirred)) == -1) then 512 write (msg,'(a,i3,a,3es18.6,2a,i3,a,i3,a,3es18.6,a,a,i4,a,a)')& 513 & ' For irreducible qpt ', iqptirred,' = ',qptirred_local(:,iqptirred),ch10, & 514 & ' the gkk element : pertcase = ',hdr1%pertcase,' ik_this_proc = ',ik_this_proc, & 515 & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,& 516 & ' and isppol ',1,ch10,& 517 & ' was not found by symmetry operations on the irreducible kpoints given' 518 MSG_ERROR(msg) 519 end if 520 end do !ikpt_phon 521 end if ! end elph_ds%tuniformgrid == 1 checks 522 523 write(msg,'(a,i0)')' read_gkk : Done completing the kpoints for pertcase ',hdr1%pertcase 524 call wrtout(std_out,msg,'COLL') 525 526 tmpflg(:,:,:,:) = 0 527 528 do idir1=1,3 529 do iatom1=1,Cryst%natom 530 ipert1 = (iatom1-1)*3+idir1 531 do idir2=1,3 532 do iatom2=1,Cryst%natom 533 ipert2 = (iatom2-1)*3+idir2 534 if (gkk_flag(ipert1,ipert1,1,1,elph_ds%qirredtofull(iqptirred)) >= 0 .and. & 535 & gkk_flag(ipert2,ipert2,1,1,elph_ds%qirredtofull(iqptirred)) >= 0) then 536 tmpflg(idir1,iatom1,idir2,iatom2) = 1 537 end if 538 end do 539 end do 540 end do 541 end do 542 543 544 ! =============================================== 545 ! Full test: need all perturbations explicitly 546 ! =============================================== 547 548 test_flag = 0 549 if (sum(tmpflg(:,1:Cryst%natom,:,1:Cryst%natom)) == (3*Cryst%natom)**2 .and. tdonecompl == 0) test_flag = 1 550 551 write(std_out,*)'read_gkk: tdonecompl = ', tdonecompl 552 553 ! de-activate completion of perts by symmetry for now. 554 ! Must be called when all irreducible perturbations are in memory!!!! 555 if (test_flag == 1 .and. tdonecompl == 0) then 556 557 ! write(std_out,*) ' read_gkk : enter fxgkkphase before completeperts' 558 ! call fxgkkphase(elph_ds,gkk_flag,h1_mat_el,iqptirred) 559 560 if (ep_prt_yambo==1) then 561 if (elph_ds%k_phon%my_nkpt /= elph_ds%k_phon%nkpt) then 562 write (msg, '(a)') 'prt_gkk_yambo can not handle parallel anaddb yet' 563 MSG_ERROR(msg) 564 end if 565 call prt_gkk_yambo(displ_cart,displ_red,elph_ds%k_phon%kpt,h1_mat_el,iqptirred,& 566 & Cryst%natom,nFSband,elph_ds%k_phon%my_nkpt,phfrq_tmp,hdr1%qptn) 567 end if 568 569 ! ======================================================================== 570 ! Now use more general symops to complete the other equivalent 571 ! perturbations: the kpoints are also shuffled by these symops 572 ! afterwards h1_mat_el_sq contains gamma_\tau\alpha,\tau'\alpha' in reduced coordinates 573 ! 574 ! \gamma_{\tau'\alpha',\tau\alpha} = 575 ! <psi_{k+q,ib2}| H(1)_{\tau'\alpha'}| psi_{k,ib1}>* \cdot 576 ! <psi_{k+q,ib2}| H(1)_{\tau \alpha }| psi_{k,ib1}> 577 ! 578 ! ======================================================================== 579 580 call completeperts(Cryst,nbranch,nFSband,elph_ds%k_phon%my_nkpt,nsppol,& 581 & gkk_flag(:,:,:,:,elph_ds%qirredtofull(iqptirred)),h1_mat_el,h1_mat_el_sq,qptirred_local(:,iqptirred),symq,qtimrev) 582 583 tdonecompl = 1 584 end if 585 586 ! ============================================================== 587 ! if we have all the perturbations for this qpoint, proceed 588 ! with scalar product, norm squared, and add weight factors 589 ! 590 ! SHOULD HAVE A TEST SO h1_mat_el IS NOT OVERWRITTEN 591 ! BEFORE PREVIOUS QPOINT IS FINISHED!!!!! 592 ! ============================================================== 593 594 test_flag = 1 595 do isppol=1,nsppol 596 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 597 do ibranch=1,nbranch 598 if (gkk_flag (ibranch,ibranch,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 599 test_flag = 0 600 exit 601 end if 602 end do 603 end do 604 end do 605 606 if (test_flag /= 0) then 607 call wrtout(std_out,' read_gkk : enter normsq_gkq',"COLL") 608 609 ! MG temporary array to save ph-linewidths before Fourier interpolation 610 ABI_ALLOCATE(qdata,(nbranch,nsppol,3)) 611 qdata(:,:,:)=zero 612 613 call normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,& 614 & h1_mat_el_sq,iqptirred,phfrq_tmp,qptirred_local,qdata) 615 616 ! save gkk_qpt, eventually to disk, for bands up to ngkkband, 617 ! NB: if the sum over bands has been performed ngkkband is 1 instead of nFSband 618 if (elph_ds%gkqwrite == 0) then 619 elph_ds%gkk_qpt(:,:,:,:,:,iqptirred) = h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,:,:) 620 else 621 ! write all kpoints to disk 622 write (std_out,*) 'size of record to be written: ', 8 * 2*elph_ds%ngkkband*elph_ds%ngkkband*& 623 & elph_ds%nbranch*elph_ds%nbranch*elph_ds%k_phon%my_nkpt*elph_ds%nsppol 624 inquire(unit=elph_ds%unitgkq, recl=isppol) 625 write (std_out,*) 'recl =', isppol 626 write (std_out,*) 'iqptirred ', iqptirred 627 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 628 write (elph_ds%unitgkq,REC=((iqptirred-1)*elph_ds%k_phon%my_nkpt+ik_this_proc)) & 629 & h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,ik_this_proc,:) 630 end do 631 end if 632 633 qdata_tmp(iqptirred,:,:,:)=qdata(:,:,:) 634 ABI_DEALLOCATE(qdata) 635 end if 636 637 call hdr_free(hdr1) 638 639 end do !of i1wf 640 641 !got all the gkk perturbations 642 643 ABI_DEALLOCATE(eigen1) 644 ABI_DEALLOCATE(h1_mat_el) 645 ABI_DEALLOCATE(h1_mat_el_sq) 646 647 if (nqptirred_local /= elph_ds%nqptirred) then 648 write (msg, '(3a,i0,i0)') & 649 & ' Found wrong number of qpoints in GKK file wrt anaddb input ', ch10, & 650 & ' nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local 651 MSG_ERROR(msg) 652 end if 653 654 !normally at this point we have the gkk for all kpoints on the FS 655 !for all the perturbations. Otherwise a 1WF file is missing. 656 !NOTE: still havent checked the qpoint grid completeness 657 do iqptirred=1,elph_ds%nqptirred 658 do isppol=1,nsppol 659 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 660 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 661 do ipert=1,nbranch 662 if (gkk_flag(ipert,ipert,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 663 write (msg,'(a,i5,1x,i5,1x,i5,1x,i5,a,a)')& 664 & ' gkk element',ipert,ikpt_phon,isppol,iqptirred,' was not found by symmetry operations ',& 665 & ' on the irreducible perturbations and qpoints given' 666 MSG_ERROR(msg) 667 end if 668 end do !ipert 669 end do !ik_this_proc 670 end do !isppol 671 end do !iqptirred 672 673 call wrtout(std_out,'read_gkk : done completing the perturbations (and checked!)','COLL') 674 675 !MG save phonon frequencies, ph-linewidths and lambda(q,n) values before Fourier interpolation 676 ABI_ALLOCATE(elph_ds%qgrid_data,(elph_ds%nqptirred,nbranch,nsppol,3)) 677 678 do iqptirred=1,elph_ds%nqptirred 679 elph_ds%qgrid_data(iqptirred,:,:,:)=qdata_tmp(iqptirred,:,:,:) 680 end do 681 682 ABI_DEALLOCATE(qdata_tmp) 683 684 end subroutine read_gkk