TABLE OF CONTENTS
ABINIT/elphon [ Functions ]
NAME
elphon
FUNCTION
This routine extracts the electron phonon coupling matrix elements and calculates related properties - Tc, phonon linewidths...
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer,BXu,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
anaddb_dtset=dataset with input variables
anaddb_dtset%a2fsmear = smearing for alpha2F function
anaddb_dtset%brav = type of Bravais lattice
anaddb_dtset%elphsmear = smearing width for gaussian integration
or buffer in energy for calculations with tetrahedra (telphint=0)
anaddb_dtset%elph_fermie = input value of Fermi energy
0 means use value from wfk file
anaddb_dtset%enunit = governs the units to be used for the output of
the phonon frequencies and e-ph quantities
anaddb_dtset%gkk2write= flag to write out gkk2 matrix elements to disk
anaddb_dtset%gkk_rptwrite= flag to write out real space gkk_rpt matrix elements to disk
anaddb_dtset%gkqwrite= flag to write out gkq matrix elements to disk
anaddb_dtset%ep_b_min= first band taken into account in FS integration (if telphint==2)
anaddb_dtset%ep_b_max= last band taken into account in FS integration (if telphint==2)
anaddb_dtset%prtfsurf = integer flag for the output of the Fermi surface (XCrysden file format)
anaddb_dtset%prtnest = integer flag for the calculation of the nesting function
anaddb_dtset%ifcflag = flag for IFC matrices in anaddb calling routine
the IFCs are presumed to be known!
anaddb_dtset%ifltransport= flag for transport properties (no=0: yes_LOVA=1; yes_nonLOVA=2 )
anaddb_dtset%kptrlatt=kpoint grid generating vectors, as in abinit
anaddb_dtset%kptrlatt_fine=kpoint grid generating vectors, for fine grid used in FS integration
anaddb_dtset%mustar = parameter for Coulombic pseudo-potential in McMillan T_c calculation
anaddb_dtset%ngqpt(3)=integers defining the number of points in the qpt sampling
anaddb_dtset%nqpath=number of vertices in the path in reciprocal space, for band structure
and phonon linewidth output
anaddb_dtset%nqshft= number of shift vectors for defining the sampling of q points
anaddb_dtset%ntemper = number of temperature points to calculate, from tempermin to
tempermin+ntemper*temperinc
anaddb_dtset%qpath=vertices in the path in reciprocal space, for band structure
and phonon linewidth output
anaddb_dtset%q1shft(3,4) =qpoint shifts considered
anaddb_dtset%telphint = flag for integration over the FS with 0=tetrahedra 1=gaussians
anaddb_dtset%tempermin = minimum temperature at which resistivity etc are calculated (in K)
anaddb_dtset%temperinc = interval temperature grid on which resistivity etc are calculated (in K)
anaddb_dtset%ep_keepbands = flag to keep gamma matrix dependence on electronic bands
Cryst<crystal_t>=data type gathering info on the crystalline structure.
Ifc<ifc_type>=Object containing the interatomic force constants.
atmfrc = inter-atomic force constants from anaddb
rpt(3,nprt) =canonical positions of R points in the unit cell
nrpt =number of real space points used to integrate IFC (for interpolation of dynamical matrices)
wghatm(natom,natom,nrpt) =Weight for the pair of atoms and the R vector
filnam(7)=character strings giving file names
comm=MPI communicator.
OUTPUT
NOTES
inspired to a large extent by epcouple.f from the DecAFT package by J. Kay Dewhurst most inputs taken from mkifc.f in anaddb anaddb_dtset%ifcflag must be 1 such that the IFC are calculated in atmfrc prior to calling elphon brav not taken into account propely in all of the code. (MG?) could choose to make a full 3 dimensional kpt array (:,:,:). Easier for many operations
PARENTS
anaddb
CHILDREN
complete_gamma,complete_gamma_tr,copy_kptrank,d2c_weights,ebands_free ebands_update_occ,eliashberg_1d,elph_ds_clean,elph_k_procs elph_tr_ds_clean,ep_fs_weights,ep_setupqpt,ftgam,ftgam_init get_all_gkk2,get_all_gkq,get_all_gkr,get_fs_bands,get_nv_fs_en get_nv_fs_temp,get_rank_1kpt,get_tau_k,get_veloc_tr,hdr_bcast hdr_fort_read,hdr_free,integrate_gamma,integrate_gamma_tr integrate_gamma_tr_lova,mka2f,mka2f_tr,mka2f_tr_lova,mka2fqgrid mkfskgrid,mknesting,mkph_linwid,mkqptequiv,order_fs_kpts,outelph printvtk,rchkgsheader,read_el_veloc,symkpt,timein,wrap2_pmhalf,wrtout xmpi_bcast
NOTES
SOURCE
93 #if defined HAVE_CONFIG_H 94 #include "config.h" 95 #endif 96 97 #include "abi_common.h" 98 99 subroutine elphon(anaddb_dtset,Cryst,Ifc,filnam,comm) 100 101 use defs_basis 102 use defs_datatypes 103 use defs_abitypes 104 use defs_elphon 105 use m_profiling_abi 106 use m_kptrank 107 use m_errors 108 use m_xmpi 109 use m_hdr 110 use m_ebands 111 112 use m_io_tools, only : open_file, is_open 113 use m_numeric_tools, only : wrap2_pmhalf 114 use m_pptools, only : printvtk 115 use m_dynmat, only : ftgam_init, ftgam 116 use m_crystal, only : crystal_t 117 use m_ifc, only : ifc_type 118 use m_nesting, only : mknesting 119 use m_anaddb_dataset, only : anaddb_dataset_type 120 121 !This section has been created automatically by the script Abilint (TD). 122 !Do not modify the following lines by hand. 123 #undef ABI_FUNC 124 #define ABI_FUNC 'elphon' 125 use interfaces_14_hidewrite 126 use interfaces_18_timing 127 use interfaces_41_geometry 128 use interfaces_77_ddb, except_this_one => elphon 129 !End of the abilint section 130 131 implicit none 132 133 !Arguments ------------------------------------ 134 !scalars 135 type(anaddb_dataset_type),intent(inout) :: anaddb_dtset 136 type(crystal_t),intent(in) :: Cryst 137 type(ifc_type),intent(inout) :: Ifc 138 integer,intent(in) :: comm 139 !arrays 140 character(len=fnlen),intent(in) :: filnam(7) 141 142 !Local variables------------------------------- 143 !scalars 144 integer,parameter :: timrev2=2,space_group0=0,master=0 145 integer :: ikpt_fine,ierr,unitgkk, unit_epts,iband,ibandp,ii 146 integer :: ikpt,jkpt,kkpt, ik1,ik2,ik3,nk1, nk2, nk3 147 integer :: iqpt,isppol,n1wf,nband,natom,onegkksize 148 integer :: timrev,unitfskgrid,qtor,idir,iFSkpq,symrankkpt,ikpt_irr 149 integer :: ep_prt_wtk ! eventually to be made into an input variable 150 integer :: fform,ie,ie1,ie2,i_start,i_end 151 integer :: ssp,s1,s2,tmp_nenergy, top_vb,nproc,me 152 integer :: nkpt_tmp 153 real(dp) :: max_occ,realdp_ex,res !,ss 154 real(dp) :: tcpu, twall, tcpui, twalli 155 real(dp) :: e1, e2, btocm3,diff, omega_max 156 real(dp) :: e_vb_max, e_cb_min, etemp_vb 157 logical :: make_gkk2,use_afm,use_tr 158 character(len=500) :: message 159 character(len=fnlen) :: fname,elph_base_name,ddkfilename,gkk_fname 160 character(len=fnlen) :: nestname 161 type(elph_tr_type) :: elph_tr_ds 162 type(elph_type) :: elph_ds 163 type(hdr_type) :: hdr,hdr1 164 type(ebands_t) :: Bst 165 !arrays 166 integer :: s1ofssp(4), s2ofssp(4) 167 integer :: qptrlatt(3,3),kptrlatt_fine(3,3) 168 integer,allocatable :: indkpt1(:) 169 integer,allocatable :: FSfullpqtofull(:,:) 170 integer,allocatable :: qpttoqpt(:,:,:) 171 integer,allocatable :: pair2red(:,:), red2pair(:,:) 172 !real(dp) :: acell_in(3),rprim_in(3,3),rprim(3,3),acell(3), 173 real(dp) :: kpt(3),shiftk(3) 174 real(dp),allocatable :: wtk_fullbz(:),wtk_folded(:) 175 real(dp),allocatable :: a2f_1d(:),dos_phon(:) 176 real(dp),allocatable :: eigenGS(:,:,:),eigenGS_fine(:,:,:) 177 real(dp),allocatable :: v_surf(:,:,:,:,:,:) 178 real(dp),allocatable :: tmp_veloc_sq1(:,:), tmp_veloc_sq2(:,:) 179 real(dp),allocatable :: coskr(:,:), sinkr(:,:) 180 181 ! ************************************************************************* 182 183 write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,& 184 & ' Properties based on electron-phonon coupling ',ch10 185 call wrtout(std_out,message,'COLL') 186 call wrtout(ab_out,message,'COLL') 187 188 call timein(tcpui,twalli) 189 write(message, '(a,f11.3,a,f11.3,a)' )& 190 & '-begin elphon at tcpu',tcpui,' and twall',twalli,' sec' 191 call wrtout(std_out,message,'COLL') 192 193 nproc = xmpi_comm_size(comm); me = xmpi_comm_rank(comm) 194 195 write(message, '(a,i0,a,i0)' )'- running on ', nproc,' cpus me = ', me 196 call wrtout(std_out,message,'PERS') 197 write(std_out,*) message 198 199 !================================== 200 !Initialization of some variables 201 !================================== 202 203 if (master == me) then 204 gkk_fname = filnam(5) 205 if (open_file(gkk_fname,message,newunit=unitgkk,form="unformatted",status="old",action="read") /=0) then 206 MSG_ERROR(message) 207 end if 208 end if 209 210 elph_base_name=trim(filnam(2))//"_ep" 211 ddkfilename=trim(filnam(7)) 212 213 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 214 215 natom = Cryst%natom 216 elph_ds%mustar = anaddb_dtset%mustar ! input mustar 217 elph_ds%nbranch = 3*natom ! number of phonon modes = 3 * natom 218 elph_ds%natom = natom ! 219 elph_ds%ep_keepbands = anaddb_dtset%ep_keepbands ! flag to sum over bands 220 elph_ds%a2fsmear = anaddb_dtset%a2fsmear ! smearing for Eliashberg functions 221 elph_ds%elphsmear = anaddb_dtset%elphsmear ! smearing for Eliashberg functions 222 elph_ds%ep_b_min = anaddb_dtset%ep_b_min 223 elph_ds%ep_b_max = anaddb_dtset%ep_b_max 224 elph_ds%telphint = anaddb_dtset%telphint 225 elph_ds%kptrlatt = anaddb_dtset%kptrlatt 226 elph_ds%kptrlatt_fine= anaddb_dtset%kptrlatt_fine 227 elph_ds%tempermin = anaddb_dtset%tempermin 228 elph_ds%temperinc = anaddb_dtset%temperinc 229 elph_ds%ntemper = anaddb_dtset%ntemper 230 elph_ds%use_k_fine = anaddb_dtset%use_k_fine 231 elph_ds%ep_int_gkk = anaddb_dtset%ep_int_gkk 232 elph_ds%ep_nspline = anaddb_dtset%ep_nspline 233 elph_ds%ep_scalprod = anaddb_dtset%ep_scalprod 234 elph_ds%prtbltztrp = anaddb_dtset%prtbltztrp 235 236 elph_ds%tuniformgrid = 1 237 elph_ds%na2f = 400 ! maximum number of Matsubara frequencies. 238 elph_ds%ep_lova = 0 ! 1 for lova and 0 for general 239 elph_ds%nenergy = 8 240 btocm3 = 1.4818474347690475d-25 241 242 !The nenergy needs to be 1) large enough to converge the integral, 2) greater 243 !than the max phonon energy. 244 !elph_ds%nenergy = INT(8*(anaddb_dtset%tempermin+anaddb_dtset%ntemper*anaddb_dtset%temperinc)/ & 245 !& (anaddb_dtset%tempermin+anaddb_dtset%temperinc)) ! number of energy levels 246 247 write(message,'(a,i6)')' The initial number of energy levels above/below Ef is set to be :',elph_ds%nenergy 248 call wrtout(std_out,message,'COLL') 249 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 250 251 !The precise number used depends on the value of Tc: 252 !they span $w_n = (2n+1) \pi T_c$ where $abs(w_n) < w_{cutoff}$ 253 !ie $|n| < n_{cutoff} = ( \frac{w_{cutoff}}{\pi T_c} ) / 2$ 254 255 !save gkk data for full kpoints to file on disk 256 257 elph_ds%gkqwrite = anaddb_dtset%gkqwrite 258 elph_ds%gkk_rptwrite = anaddb_dtset%gkk_rptwrite 259 elph_ds%gkk2write = anaddb_dtset%gkk2write 260 261 !This should never be turned off: symmetrization of elphon matrix elements in complete_gkk. See get_all_gkq 262 elph_ds%symgkq=anaddb_dtset%symgkq 263 264 elph_ds%elph_base_name = trim(elph_base_name) 265 266 !MG: @Matthieu: Why this? Now we should always use the value of rprim and acell reported in IFC 267 !rprim_in = Ifc%rprim 268 !acell_in = Ifc%acell 269 270 !normalize input rprim and acell. 271 !do ii=1,3 272 ! ss = sqrt(rprim_in(1,ii)**2+rprim_in(2,ii)**2+rprim_in(3,ii)**2) 273 ! rprim(:,ii) = rprim_in(:,ii)/ss 274 ! acell(ii) = acell_in(ii) * ss 275 !end do 276 277 !make dimension-ful rprimd and gprimd for transformation of derivatives to cartesian coordinates. 278 !call mkrdim(acell,rprim,rprimd) 279 !call matr3inv(rprimd,gprimd) 280 281 !rprimd = cryst%rprimd 282 !gprimd = cryst%gprimd 283 284 !=================== 285 !Check some inputs 286 !=================== 287 if (Cryst%nsym==1) then 288 write (message,'(7a)')ch10,& 289 & ' elphon: COMMENT- ',ch10,& 290 & ' Symmetries are not used! ',ch10,& 291 & ' Full matrix elements must be supplied for all perturbations and qpoints!',ch10 292 call wrtout(std_out,message,'COLL') 293 call wrtout(ab_out,message,'COLL') 294 if ( ANY( ABS(Cryst%tnons(:,1)) > tol10) ) then 295 MSG_ERROR('nsym==1 but the symmetry is not the identity') 296 end if 297 end if 298 299 if (anaddb_dtset%ifcflag/=1) then 300 write(message,'(a,i0)')& 301 & ' ifcflag should be set to 1 since the IFC matrices are supposed to exist but ifcflag= ',anaddb_dtset%ifcflag 302 MSG_ERROR(message) 303 end if 304 305 call timein(tcpu,twall) 306 write(message, '(a,f11.3,a,f11.3,a)' )& 307 & '-elphon begin setup after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 308 call wrtout(std_out,message,'COLL') 309 tcpui = tcpu 310 twalli = twall 311 312 !================================= 313 !Set up the full grid of qpoints 314 !================================= 315 !use time reversal symmetry always when possible for kpoint reduction, 316 !and suppose it has been used in WF generation 317 !not used for the moment: values are always taken from input files. 318 timrev = 1 319 call ep_setupqpt(elph_ds,cryst,anaddb_dtset,qptrlatt,timrev) 320 321 !==================================== 322 !Read the GS header of the GKK file 323 !this will give the phon grid of k 324 !and the Fermi surface integration weights 325 !==================================== 326 call wrtout (std_out,' elphon: reading and checking the GS header of the GKK file','COLL') 327 328 if (master == me) then 329 call rchkGSheader(hdr,natom,nband,unitgkk) 330 end if 331 332 !the following is for the non master nodes 333 call hdr_bcast(hdr,master,me,comm) 334 call xmpi_bcast(nband, master,comm,ierr) 335 elph_ds%nband = nband 336 337 elph_ds%nsppol =hdr%nsppol 338 elph_ds%nspinor=hdr%nspinor 339 340 !in spinor or spin polarized case, orbitals have occupation <= 1 instead of 2 341 max_occ = one 342 if (hdr%nspinor == 2) max_occ = half ! this accounts for the doubling of the num of bands, even though spin channels are not well defined 343 if (elph_ds%nsppol > 1) max_occ = one 344 write (std_out,*) ' max_occ factor ', max_occ 345 346 elph_ds%occ_factor = one 347 if (hdr%nspinor == 1 .and. hdr%nsppol == 1) then 348 elph_ds%occ_factor = one 349 else if (hdr%nspinor == 2) then 350 elph_ds%occ_factor = two 351 else if (hdr%nsppol == 2) then 352 elph_ds%occ_factor = one 353 end if 354 355 !================================================== 356 !Read GS eigenvalues for each irreducible kpt and 357 !number of 1WF files contributing to the GKK file 358 !================================================== 359 360 ABI_ALLOCATE(eigenGS,(nband,hdr%nkpt,elph_ds%nsppol)) 361 362 if (master == me) then 363 do isppol=1,elph_ds%nsppol 364 do ikpt=1,hdr%nkpt 365 read(unitgkk) eigenGS(:,ikpt,isppol) 366 end do 367 end do 368 369 ! read number of 1WF files contributing to the GKK file 370 read(unitgkk) n1wf 371 write(message,'(a,i0)')' elphon : number of perturbations in the gkk file = ',n1wf 372 call wrtout(std_out,message,'COLL') 373 end if 374 call xmpi_bcast(n1wf, master, comm, ierr) 375 call xmpi_bcast(eigenGS, master, comm, ierr) 376 377 !================================================== 378 !Set elph_ds%fermie: either comes from anaddb input file or from wfk file 379 !================================================== 380 elph_ds%fermie = hdr%fermie 381 !elph_ds%nelect = hdr_get_nelect_byocc(Hdr) 382 elph_ds%nelect = Hdr%nelect 383 if (abs(anaddb_dtset%elph_fermie) > tol10) then 384 elph_ds%fermie = anaddb_dtset%elph_fermie 385 write(message,'(a,E20.12)')' Fermi level set by the user at :',elph_ds%fermie 386 call wrtout(std_out,message,'COLL') 387 Bst = ebands_from_hdr(Hdr,nband,eigenGS) 388 else if (abs(anaddb_dtset%ep_extrael) > tol10) then 389 if (abs(anaddb_dtset%ep_extrael) > 1.0d2) then 390 write(message,'(a,E20.12)')' Doping set by the user is (negative for el doping) :',& 391 & anaddb_dtset%ep_extrael 392 call wrtout(std_out,message,'COLL') 393 anaddb_dtset%ep_extrael = anaddb_dtset%ep_extrael*cryst%ucvol*btocm3*(-1.0d0) 394 end if 395 write(message,'(a,E20.12)')' Additional electrons per unit cell set by the user at :',& 396 & anaddb_dtset%ep_extrael 397 call wrtout(std_out,message,'COLL') 398 elph_ds%nelect = elph_ds%nelect + anaddb_dtset%ep_extrael 399 bst = ebands_from_hdr(Hdr,nband,eigenGS,nelect=elph_ds%nelect) 400 401 ! set Bst to use FD occupations: 402 Bst%occopt = 3 403 ! Bst%tsmear = 0.00001_dp ! is this small etol9 Bst%tsmeatol90001_dp ! last used 404 Bst%tsmear = tol9 ! is this small etol9 Bst%tsmeatol90001_dp ! last used 405 ! Calculate occupation numbers. 406 call ebands_update_occ(Bst,-99.99_dp) 407 write(message,'(a,E20.12)')' Fermi level is now calculated to be :',Bst%fermie 408 call wrtout(std_out,message,'COLL') 409 elph_ds%fermie = BSt%fermie 410 else 411 bst = ebands_from_hdr(Hdr,nband,eigenGS) 412 end if 413 call wrtout(std_out,message,'COLL') 414 415 !==================================================================== 416 !Setup of the phon k-grid : 417 !1) get bands near Ef 418 !==================================================================== 419 call get_fs_bands(eigenGS,hdr,elph_ds%fermie,anaddb_dtset%ep_b_min, anaddb_dtset%ep_b_max,& 420 & elph_ds%minFSband,elph_ds%maxFSband,elph_ds%k_phon%nkptirr) 421 422 elph_ds%nFSband = elph_ds%maxFSband - elph_ds%minFSband + 1 423 424 !Modify the band gap by sissor shift of the CB 425 if (abs(anaddb_dtset%band_gap) < 10.0d0) then 426 anaddb_dtset%band_gap = anaddb_dtset%band_gap*0.036749309 ! eV2Ha 427 do isppol=1,elph_ds%nsppol 428 429 !First find where the gap is 430 etemp_vb = 999.0d0 431 top_vb = elph_ds%minFSband 432 do iband = elph_ds%minFSband, elph_ds%maxFSband 433 e_vb_max = maxval(eigenGS(iband,:,isppol)) 434 if (dabs(e_vb_max-elph_ds%fermie) < etemp_vb) then 435 etemp_vb = dabs(e_vb_max-elph_ds%fermie) 436 top_vb = iband 437 end if 438 end do 439 do iband = top_vb, elph_ds%maxFSband 440 e_vb_max = maxval(eigenGS(iband,:,isppol)) 441 if (dabs(e_vb_max-maxval(eigenGS(top_vb,:,isppol))) < tol6) then 442 etemp_vb = dabs(e_vb_max-elph_ds%fermie) 443 top_vb = iband 444 end if 445 end do 446 e_vb_max = maxval(eigenGS(top_vb,:,isppol)) 447 e_cb_min = minval(eigenGS(top_vb+1,:,isppol)) 448 write(message,'(a,E20.12,2x,E20.12)')' elphon : original fermi energy = ', elph_ds%fermie 449 call wrtout(std_out,message,'COLL') 450 write(message,'(a,E20.12,2x,E20.12)')' elphon : top of VB, bottom of CB = ',e_vb_max, e_cb_min 451 call wrtout(std_out,message,'COLL') 452 453 do iband = top_vb+1, elph_ds%maxFSband 454 eigenGS(iband,:,isppol) = eigenGS(iband,:,isppol) + (anaddb_dtset%band_gap-(e_cb_min-e_vb_max)) 455 end do 456 end do !nsppol 457 458 !! recalculate Fermi level 459 !elph_ds%nelect = hdr_get_nelect_byocc(Hdr) 460 elph_ds%nelect = Hdr%nelect 461 if (abs(anaddb_dtset%elph_fermie) > tol10) then 462 elph_ds%fermie = anaddb_dtset%elph_fermie 463 write(message,'(a,E20.12)')' Fermi level set by the user at :',elph_ds%fermie 464 call wrtout(std_out,message,'COLL') 465 bst = ebands_from_hdr(Hdr,nband,eigenGS) 466 else if (abs(anaddb_dtset%ep_extrael) > tol10) then 467 write(message,'(a,E20.12)')' Additional electrons per unit cell set by the user at :',anaddb_dtset%ep_extrael 468 call wrtout(std_out,message,'COLL') 469 elph_ds%nelect = elph_ds%nelect + anaddb_dtset%ep_extrael 470 bst = ebands_from_hdr(Hdr,nband,eigenGS,nelect=elph_ds%nelect) 471 472 ! set Bst to use FD occupations: 473 Bst%occopt = 3 474 ! Bst%tsmear = 0.00001_dp ! is this small etol9 Bst%tsmeatol90001_dp ! last used 475 Bst%tsmear = tol9 ! is this small etol9 Bst%tsmeatol90001_dp ! last used 476 ! Calculate occupation numbers. 477 call ebands_update_occ(Bst,-99.99_dp) 478 write(message,'(a,E20.12)')' Fermi level is now calculated to be :',Bst%fermie 479 call wrtout(std_out,message,'COLL') 480 elph_ds%fermie = BSt%fermie 481 else 482 bst = ebands_from_hdr(Hdr,nband,eigenGS) 483 end if 484 call wrtout(std_out,message,'COLL') 485 end if !modify band_gap 486 487 if (elph_ds%ep_keepbands == 0) then !we are summing over bands 488 elph_ds%ngkkband = 1 489 else if (elph_ds%ep_keepbands == 1) then 490 ! keep the band dependency btw elph_ds%minFSband and elph_ds%maxFSband 491 elph_ds%ngkkband = elph_ds%nFSband 492 else 493 write(message,'(a,i0)')' ep_keepbands must be 0 or 1 while it is: ',elph_ds%ep_keepbands 494 MSG_BUG(message) 495 end if 496 497 write(message,'(a,i0,2x,i0)')' elphon : minFSband, maxFSband = ',elph_ds%minFSband,elph_ds%maxFSband 498 call wrtout(std_out,message,'COLL') 499 500 501 ABI_ALLOCATE(elph_ds%k_phon%kptirr,(3,elph_ds%k_phon%nkptirr)) 502 ABI_ALLOCATE(elph_ds%k_phon%irredtoGS,(elph_ds%k_phon%nkptirr)) 503 504 !==================================================================== 505 !2) order irred k-points 506 !==================================================================== 507 if (master == me) then 508 call order_fs_kpts(hdr%kptns, hdr%nkpt, elph_ds%k_phon%kptirr,elph_ds%k_phon%nkptirr,elph_ds%k_phon%irredtoGS) 509 end if 510 call xmpi_bcast(elph_ds%k_phon%nkptirr, master, comm, ierr) 511 call xmpi_bcast(elph_ds%k_phon%kptirr, master, comm, ierr) 512 call xmpi_bcast(elph_ds%k_phon%irredtoGS, master, comm, ierr) 513 514 !========================================== 515 !3) reconstruct full kgrid from irred kpoints, 516 !========================================== 517 call mkFSkgrid (elph_ds%k_phon, Cryst%nsym, Cryst%symrec, timrev) 518 519 ! check that kptrlatt is coherent with kpt found here 520 nkpt_tmp = elph_ds%kptrlatt(1,1)*elph_ds%kptrlatt(2,2)*elph_ds%kptrlatt(3,3) 521 if (sum(abs(elph_ds%kptrlatt(:,:))) /= nkpt_tmp) then 522 MSG_WARNING(' the input kptrlatt is not diagonal... ') 523 end if 524 if (anaddb_dtset%ifltransport > 1 .and. nkpt_tmp /= elph_ds%k_phon%nkpt) then 525 write(message,'(a,i0,a,i0)')& 526 & ' the input kptrlatt is inconsistent ', nkpt_tmp, " /= ", elph_ds%k_phon%nkpt 527 MSG_ERROR(message) 528 end if 529 530 if (anaddb_dtset%ifltransport==3 ) then 531 !==================================================================== 532 ! The real irred kpt, now only used by get_tau_k 533 !==================================================================== 534 535 ABI_ALLOCATE(indkpt1,(elph_ds%k_phon%nkpt)) 536 ABI_ALLOCATE(wtk_fullbz,(elph_ds%k_phon%nkpt)) 537 ABI_ALLOCATE(wtk_folded,(elph_ds%k_phon%nkpt)) 538 539 wtk_fullbz(:) = one/dble(elph_ds%k_phon%nkpt) !weights normalized to unity 540 call symkpt(0,cryst%gmet,indkpt1,0,elph_ds%k_phon%kpt,elph_ds%k_phon%nkpt,elph_ds%k_phon%new_nkptirr,& 541 & Cryst%nsym,Cryst%symrec,timrev,wtk_fullbz,wtk_folded) 542 543 write (message,'(2a,i0)')ch10,' Number of irreducible k-points = ',elph_ds%k_phon%new_nkptirr 544 call wrtout(std_out,message,'COLL') 545 546 ABI_ALLOCATE(elph_ds%k_phon%new_kptirr,(3,elph_ds%k_phon%new_nkptirr)) 547 ABI_ALLOCATE(elph_ds%k_phon%new_wtkirr,(elph_ds%k_phon%new_nkptirr)) 548 ABI_ALLOCATE(elph_ds%k_phon%new_irredtoGS,(elph_ds%k_phon%new_nkptirr)) 549 550 ikpt_irr = 0 551 do ikpt=1,elph_ds%k_phon%nkpt 552 if (wtk_folded(ikpt) /= zero) then 553 ikpt_irr = ikpt_irr + 1 554 elph_ds%k_phon%new_kptirr(:,ikpt_irr) = elph_ds%k_phon%kpt(:,ikpt) 555 elph_ds%k_phon%new_wtkirr(ikpt_irr) = wtk_folded(ikpt) 556 elph_ds%k_phon%new_irredtoGS(ikpt_irr) = ikpt 557 end if 558 end do 559 if (ikpt_irr .ne. elph_ds%k_phon%new_nkptirr) then 560 write (message,'(a)')' The number of irred nkpt does not match! ' 561 MSG_ERROR(message) 562 end if 563 564 ABI_DEALLOCATE(indkpt1) 565 ABI_DEALLOCATE(wtk_fullbz) 566 ABI_DEALLOCATE(wtk_folded) 567 end if 568 569 !==================================================================== 570 !4) setup weights for integration (gaussian or tetrahedron method) 571 !==================================================================== 572 elph_ds%k_phon%nband = elph_ds%nFSband 573 elph_ds%k_phon%nsppol = elph_ds%nsppol 574 elph_ds%k_phon%nsym = Cryst%nsym 575 ABI_ALLOCATE(elph_ds%k_phon%wtk,(elph_ds%nFSband,elph_ds%k_phon%nkpt,elph_ds%k_phon%nsppol)) 576 577 call ep_fs_weights(anaddb_dtset%ep_b_min, anaddb_dtset%ep_b_max, eigenGS, anaddb_dtset%elphsmear, & 578 & elph_ds%fermie, cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, elph_ds%minFSband, nband, elph_ds%nFSband, & 579 & elph_ds%nsppol, anaddb_dtset%telphint, elph_ds%k_phon) 580 581 !distribute k-points among processors, if any 582 call elph_k_procs(nproc, elph_ds%k_phon) 583 584 !===================================================== 585 !get kpt info from the fine grid part 586 !===================================================== 587 if (anaddb_dtset%use_k_fine == 1) then 588 589 if (abs(anaddb_dtset%band_gap) < 10.0d0) then 590 write (message,'(a)')' Not coded yet when use_k_fine and band_gap are both used' 591 MSG_ERROR(message) 592 end if 593 594 if (master == me) then 595 if (open_file("densergrid_GKK",message,newunit=unitfskgrid,form="unformatted",status="old") /=0) then 596 MSG_ERROR(message) 597 end if 598 !read the header of file 599 call hdr_fort_read(hdr1, unitfskgrid, fform) 600 ABI_CHECK(fform/=0,'denser grid GKK header was mis-read. fform == 0') 601 end if 602 call hdr_bcast(hdr1,master,me,comm) 603 604 ABI_ALLOCATE(eigenGS_fine,(nband,hdr1%nkpt,elph_ds%nsppol)) 605 606 if (master == me) then 607 do isppol=1,elph_ds%nsppol 608 do ikpt=1,hdr1%nkpt 609 read(unitfskgrid) eigenGS_fine(:,ikpt,isppol) 610 end do 611 end do 612 close(unitfskgrid) 613 end if 614 call xmpi_bcast(eigenGS_fine, master, comm, ierr) 615 616 ! Reinit the structure storing the eigevalues. 617 ! Be careful. This part has not been tested. 618 call ebands_free(Bst) 619 bst = ebands_from_hdr(hdr1,nband,eigenGS_fine) 620 621 elph_ds%k_fine%nkptirr = hdr1%nkpt 622 ABI_ALLOCATE(elph_ds%k_fine%kptirr,(3,elph_ds%k_fine%nkptirr)) 623 ABI_ALLOCATE(elph_ds%k_fine%irredtoGS,(elph_ds%k_fine%nkptirr)) 624 625 call order_fs_kpts(hdr1%kptns, hdr1%nkpt, elph_ds%k_fine%kptirr,& 626 & elph_ds%k_fine%nkptirr,elph_ds%k_fine%irredtoGS) 627 628 call hdr_free(hdr1) 629 630 call mkFSkgrid (elph_ds%k_fine, Cryst%nsym, Cryst%symrec, timrev) 631 632 elph_ds%k_fine%nband = elph_ds%nFSband 633 elph_ds%k_fine%nsppol = elph_ds%nsppol 634 elph_ds%k_fine%nsym = Cryst%nsym 635 636 ABI_ALLOCATE(elph_ds%k_fine%wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol)) 637 638 kptrlatt_fine = elph_ds%kptrlatt_fine 639 640 call ep_fs_weights(anaddb_dtset%ep_b_min, anaddb_dtset%ep_b_max, & 641 & eigenGS_fine, anaddb_dtset%elphsmear, & 642 & elph_ds%fermie, cryst%gprimd, elph_ds%k_fine%irredtoGS, kptrlatt_fine, & 643 & max_occ, elph_ds%minFSband, nband, elph_ds%nFSband, & 644 & elph_ds%nsppol, anaddb_dtset%telphint, elph_ds%k_fine) 645 646 else ! not using k_fine 647 elph_ds%k_fine%nband = elph_ds%k_phon%nband 648 elph_ds%k_fine%nsppol = elph_ds%k_phon%nsppol 649 elph_ds%k_fine%nsym = elph_ds%k_phon%nsym 650 651 elph_ds%k_fine%nkpt = elph_ds%k_phon%nkpt 652 elph_ds%k_fine%nkptirr = elph_ds%k_phon%nkptirr 653 654 elph_ds%k_fine%my_nkpt = elph_ds%k_phon%my_nkpt 655 656 ABI_ALLOCATE(elph_ds%k_fine%my_kpt,(elph_ds%k_fine%nkpt)) 657 elph_ds%k_fine%my_kpt = elph_ds%k_phon%my_kpt 658 659 ABI_ALLOCATE(elph_ds%k_fine%my_ikpt,(elph_ds%k_fine%my_nkpt)) 660 elph_ds%k_fine%my_ikpt = elph_ds%k_phon%my_ikpt 661 662 ABI_ALLOCATE(elph_ds%k_fine%kptirr,(3,elph_ds%k_fine%nkptirr)) 663 elph_ds%k_fine%kptirr = elph_ds%k_phon%kptirr 664 ABI_ALLOCATE(elph_ds%k_fine%wtkirr,(elph_ds%k_fine%nkptirr)) 665 elph_ds%k_fine%wtkirr = elph_ds%k_phon%wtkirr 666 667 ABI_ALLOCATE(elph_ds%k_fine%wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%k_fine%nsppol)) 668 elph_ds%k_fine%wtk = elph_ds%k_phon%wtk 669 ABI_ALLOCATE(elph_ds%k_fine%kpt,(3,elph_ds%k_fine%nkpt)) 670 elph_ds%k_fine%kpt = elph_ds%k_phon%kpt 671 672 call copy_kptrank(elph_ds%k_phon%kptrank_t, elph_ds%k_fine%kptrank_t) 673 674 ABI_ALLOCATE(elph_ds%k_fine%irr2full,(elph_ds%k_fine%nkptirr)) 675 elph_ds%k_fine%irr2full = elph_ds%k_phon%irr2full 676 ABI_ALLOCATE(elph_ds%k_fine%full2irr,(3,elph_ds%k_fine%nkpt)) 677 elph_ds%k_fine%full2irr = elph_ds%k_phon%full2irr 678 ABI_ALLOCATE(elph_ds%k_fine%full2full,(2,elph_ds%k_fine%nsym,elph_ds%k_fine%nkpt)) 679 elph_ds%k_fine%full2full = elph_ds%k_phon%full2full 680 681 ABI_ALLOCATE(elph_ds%k_fine%irredtoGS,(elph_ds%k_fine%nkptirr)) 682 elph_ds%k_fine%irredtoGS = elph_ds%k_phon%irredtoGS 683 684 ! call elph_k_copy(elph_ds%k_phon, elph_ds%k_fine) 685 686 kptrlatt_fine = elph_ds%kptrlatt 687 688 ABI_ALLOCATE(eigenGS_fine,(nband,elph_ds%k_fine%nkptirr,elph_ds%nsppol)) 689 690 eigenGS_fine = eigenGS 691 end if ! k_fine or not 692 693 if (elph_ds%kptrlatt_fine(1,1) == 0) then ! when there is not input for kptrlatt_fine 694 elph_ds%kptrlatt_fine = kptrlatt_fine 695 end if 696 697 call timein(tcpu,twall) 698 write(message, '(a,f11.3,a,f11.3,a)' )& 699 & '-elphon k and q grids have been setup after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 700 call wrtout(std_out,message,'COLL') 701 tcpui = tcpu 702 twalli = twall 703 704 !==================================================================== 705 !5) calculate DOS at Ef 706 !==================================================================== 707 ABI_ALLOCATE(elph_ds%n0,(elph_ds%nsppol)) 708 709 !SPPOL sum over spin channels to get total DOS 710 !channels decoupled => use separate values for DOS_up(Ef) resp down 711 do isppol=1,elph_ds%nsppol 712 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 713 end do 714 715 if (elph_ds%nsppol == 1) then 716 write (std_out,*) ' elphon : the estimated DOS(E_Fermi) = ', elph_ds%n0(1), ' states/Ha/spin ' 717 write (std_out,*) ' elphon : the total FS weight and # of kpoints = ',sum(elph_ds%k_fine%wtk),elph_ds%k_fine%nkpt 718 else if (elph_ds%nsppol == 2) then 719 write (std_out,*) ' elphon : the spin up DOS(E_Fermi) = ', elph_ds%n0(1), ' states/Ha/spin ' 720 write (std_out,*) ' elphon : the spin down DOS(E_Fermi) = ', elph_ds%n0(2), ' states/Ha/spin ' 721 write (std_out,*) ' elphon : total DOS(E_Fermi) = ', elph_ds%n0(1)+elph_ds%n0(2), ' states/Ha ' 722 write (std_out,*) ' elphon : the spin up FS weight and # of kpoints = ',& 723 & sum(elph_ds%k_fine%wtk(:,:,1)),elph_ds%k_fine%nkpt 724 write (std_out,*) ' elphon : the spin down FS weight and # of kpoints = ',& 725 & sum(elph_ds%k_fine%wtk(:,:,2)),elph_ds%k_fine%nkpt 726 else 727 write (message,'(a,i0)') 'bad value for nsppol ', elph_ds%nsppol 728 MSG_ERROR(message) 729 end if 730 731 ABI_ALLOCATE(elph_ds%gkk_intweight,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,elph_ds%nsppol)) 732 733 if (elph_ds%ep_keepbands == 0) then 734 ! use trivial integration weights for single band, 735 ! since average over bands is done in normsq_gkk 736 elph_ds%gkk_intweight(1,:,:) = one 737 738 else if (elph_ds%ep_keepbands == 1) then 739 ! use elph_ds%k_fine%wtk since average over bands is not done in normsq_gkk 740 if (elph_ds%use_k_fine == 1) then 741 call d2c_weights(elph_ds) 742 end if 743 elph_ds%gkk_intweight(:,:,:) = elph_ds%k_phon%wtk(:,:,:) 744 else 745 write(message,'(a,i0)')' ep_keepbands must be 0 or 1 while it is : ',elph_ds%ep_keepbands 746 MSG_ERROR(message) 747 end if 748 749 ep_prt_wtk = 0 750 if (ep_prt_wtk == 1) then 751 do iband=1, elph_ds%ngkkband 752 do ikpt_fine=1, elph_ds%k_fine%nkpt 753 write (300,*) ikpt_fine, elph_ds%gkk_intweight(iband,ikpt_fine,1) 754 end do 755 end do 756 end if 757 758 759 call timein(tcpu,twall) 760 write(message, '(a,f11.3,a,f11.3,a)' )& 761 & '-elphon weights and DOS setup after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 762 call wrtout(std_out,message,'COLL') 763 tcpui = tcpu 764 twalli = twall 765 766 !Output of the Fermi Surface 767 if (anaddb_dtset%prtfsurf == 1 .and. master == me) then 768 fname=trim(elph_ds%elph_base_name) // '_BXSF' 769 if (ebands_write_bxsf(Bst, Cryst, fname) /= 0) then 770 MSG_WARNING("Cannot produce file for Fermi surface, check log file for more info") 771 end if 772 end if 773 774 !========================================================= 775 !Get equivalence between a kpt_phon pair and a qpt in qpt_full 776 !only works if the qpt grid is complete (identical to 777 !the kpt one, with a basic shift of (0,0,0) 778 !========================================================= 779 780 !mapping of k + q onto k' for k and k' in full BZ 781 ABI_ALLOCATE(FSfullpqtofull,(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)) 782 783 !qpttoqpt(itim,isym,iqpt) = qpoint index which transforms to iqpt under isym and with time reversal itim. 784 ABI_ALLOCATE(qpttoqpt,(2,Cryst%nsym,elph_ds%nqpt_full)) 785 786 call wrtout(std_out,'elphon: calling mkqptequiv to set up the FS qpoint set',"COLL") 787 788 call mkqptequiv (FSfullpqtofull,Cryst,elph_ds%k_phon%kpt,elph_ds%k_phon%nkpt,& 789 & elph_ds%nqpt_full,qpttoqpt,elph_ds%qpt_full) 790 791 !========================================== 792 !Set up dataset for phonon interpolations 793 !========================================== 794 795 !transfer ifltransport flag to structure 796 elph_tr_ds%ifltransport=anaddb_dtset%ifltransport 797 !transfer name of files file for ddk 798 elph_tr_ds%ddkfilename=ddkfilename 799 800 !reduce qpt_full to correct zone 801 do iqpt=1,elph_ds%nqpt_full 802 call wrap2_pmhalf(elph_ds%qpt_full(1,iqpt),kpt(1),res) 803 call wrap2_pmhalf(elph_ds%qpt_full(2,iqpt),kpt(2),res) 804 call wrap2_pmhalf(elph_ds%qpt_full(3,iqpt),kpt(3),res) 805 elph_ds%qpt_full(:,iqpt)=kpt 806 end do 807 808 !test density of k+q grid: the following should be close to n0 squared 809 !FIXME: generalize for sppol 810 res = zero 811 do ikpt_fine = 1, elph_ds%k_phon%nkpt 812 do iqpt = 1, elph_ds%nqpt_full 813 kpt = elph_ds%k_phon%kpt(:,ikpt_fine) + elph_ds%qpt_full(:,iqpt) 814 call get_rank_1kpt (kpt,symrankkpt,elph_ds%k_phon%kptrank_t) 815 iFSkpq = elph_ds%k_phon%kptrank_t%invrank(symrankkpt) 816 do iband = 1, elph_ds%ngkkband 817 do ibandp = 1, elph_ds%ngkkband 818 res = res + elph_ds%gkk_intweight(iband,ikpt_fine,1)*elph_ds%gkk_intweight(ibandp,iFSkpq,1) 819 end do 820 end do 821 end do 822 end do 823 res = res / elph_ds%k_phon%nkpt/elph_ds%k_phon%nkpt 824 write (std_out,*) 'elphon: integrated value of intweight for given k and q grid : ', res, res / elph_ds%n0(1)**2 825 826 res = zero 827 do ikpt_fine = 1, elph_ds%k_phon%nkpt 828 do iqpt = 1, elph_ds%k_phon%nkpt 829 kpt = elph_ds%k_phon%kpt(:,ikpt_fine) + elph_ds%k_phon%kpt(:,iqpt) 830 call get_rank_1kpt (kpt,symrankkpt,elph_ds%k_phon%kptrank_t) 831 iFSkpq = elph_ds%k_phon%kptrank_t%invrank(symrankkpt) 832 do iband = 1, elph_ds%ngkkband 833 do ibandp = 1, elph_ds%ngkkband 834 res = res + elph_ds%gkk_intweight(iband,ikpt_fine,1)*elph_ds%gkk_intweight(ibandp,iFSkpq,1) 835 end do 836 end do 837 end do 838 end do 839 res = res / elph_ds%k_phon%nkpt/elph_ds%k_phon%nkpt 840 write (std_out,*) 'elphon: integrated value of intweight for double k grid : ', res, res / elph_ds%n0(1)**2 841 842 !=================================================== 843 !Allocate all important arrays for FS integrations 844 !=================================================== 845 846 !Record sizes for matrices on disk: complex and real versions (for real and recip space resp!) 847 onegkksize = 2*elph_ds%nbranch*elph_ds%nbranch*& 848 & elph_ds%ngkkband*elph_ds%ngkkband*& 849 & elph_ds%nsppol*kind(realdp_ex) 850 851 elph_tr_ds%onegkksize=onegkksize 852 853 write (message,'(4a)')& 854 & ' elphon : preliminary setup completed ',ch10,& 855 & ' calling get_all_gkq to read in all the e-ph matrix elements',ch10 856 call wrtout(std_out,message,'COLL') 857 858 !flag to do scalar product in gkq before interpolation: 859 !should also used in interpolate_gkk and mkph_linwid 860 if (elph_ds%ep_scalprod==0) then 861 write (std_out,*) ' elphon: will NOT perform scalar product with phonon' 862 write (std_out,*) ' displacement vectors in read_gkk. ep_scalprod==0' 863 else if (elph_ds%ep_scalprod==1) then 864 write (std_out,*) ' elphon: will perform scalar product with phonon' 865 write (std_out,*) ' displacement vectors in read_gkk. ep_scalprod==1' 866 else 867 MSG_ERROR('illegal value for ep_scalprod') 868 end if 869 870 call timein(tcpu,twall) 871 write(message, '(a,f11.3,a,f11.3,a)' )& 872 & '-elphon begin gkq construction after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 873 call wrtout(std_out,message,'COLL') 874 tcpui = tcpu 875 twalli = twall 876 877 call get_all_gkq (elph_ds,Cryst,ifc,Bst,FSfullpqtofull,nband,n1wf,onegkksize,& 878 & qpttoqpt,anaddb_dtset%ep_prt_yambo,unitgkk,elph_tr_ds%ifltransport) 879 880 if (master == me) then 881 close (unitgkk) 882 end if 883 884 call timein(tcpu,twall) 885 write(message, '(a,f11.3,a,f11.3,a)' )& 886 & '-elphon end gkq construction after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 887 call wrtout(std_out,message,'COLL') 888 tcpui = tcpu 889 twalli = twall 890 891 if (elph_tr_ds%ifltransport==1 .or. elph_tr_ds%ifltransport==2 .or. elph_tr_ds%ifltransport==3)then 892 893 ! check inputs 894 ! TODO: should be done at earlier stage of initialization and checking 895 if (elph_ds%ngkkband /= elph_ds%nFSband) then 896 write (message,'(a)') 'need to keep electron band dependency in memory for transport calculations' 897 MSG_ERROR(message) 898 end if 899 900 ! bxu, moved the allocation from get_veloc_tr to elphon 901 if (anaddb_dtset%use_k_fine == 1) then 902 ABI_ALLOCATE(elph_tr_ds%el_veloc,(elph_ds%k_fine%nkpt,nband,3,elph_ds%nsppol)) 903 else 904 ABI_ALLOCATE(elph_tr_ds%el_veloc,(elph_ds%k_phon%nkpt,nband,3,elph_ds%nsppol)) 905 end if 906 ABI_ALLOCATE(elph_tr_ds%FSelecveloc_sq,(3,elph_ds%nsppol)) 907 908 ! this only needs to be read in once - the fermi level average is later done many times with get_veloc_tr 909 if (me == master) then 910 if (anaddb_dtset%use_k_fine == 1) then 911 call read_el_veloc(nband,elph_ds%k_fine%nkpt,elph_ds%k_fine%kpt,elph_ds%nsppol,elph_tr_ds) 912 else 913 call read_el_veloc(nband,elph_ds%k_phon%nkpt,elph_ds%k_phon%kpt,elph_ds%nsppol,elph_tr_ds) 914 end if 915 end if 916 call xmpi_bcast (elph_tr_ds%el_veloc, master, comm, ierr) 917 918 call get_veloc_tr(elph_ds,elph_tr_ds) 919 end if 920 921 !Output of the Fermi velocities 922 !to be used for Mayavi visualization 923 if (anaddb_dtset%prtfsurf == 1 .and. master == me) then 924 fname = trim(elph_ds%elph_base_name) // '_VTK' 925 926 ! FIXME 927 ! shiftk is defined neither in the anaddb nor in the hdr data type 928 ! an incorrect FS will be produced in case of a shifted k-grid used during the GS calculation 929 ! check if we are using a unshifthed kgrid, obviously doesnt work in case 930 ! of multiple shifts containg a zero translation but in this case prtbxsf should work 931 shiftk=one 932 do ii=1,hdr%nkpt 933 if (all(hdr%kptns(:,ii) == zero)) shiftk=zero 934 end do 935 936 use_afm=(hdr%nsppol==1.and.hdr%nspden==2) 937 ! MG FIXME warning time reversal is always assumed to be present. 938 ! the header should report this information. 939 940 use_tr=(timrev==1) 941 942 nk1 = elph_ds%kptrlatt_fine(1,1) 943 nk2 = elph_ds%kptrlatt_fine(2,2) 944 nk3 = elph_ds%kptrlatt_fine(3,3) 945 946 ABI_ALLOCATE(v_surf,(nband,nk1+1,nk2+1,nk3+1,3,elph_ds%nsppol)) 947 v_surf = zero 948 do isppol=1,elph_ds%nsppol 949 do iband=1,nband 950 do ikpt = 1, nk1+1 951 do jkpt = 1, nk2+1 952 do kkpt = 1, nk3+1 953 ik1 = ikpt 954 ik2 = jkpt 955 ik3 = kkpt 956 if (ikpt > nk1) ik1 = ikpt - nk1 957 if (jkpt > nk2) ik2 = jkpt - nk2 958 if (kkpt > nk3) ik3 = kkpt - nk3 959 ikpt_fine = (ik1-1)*nk2*nk3 + (ik2-1)*nk3 + ik3 960 ! v_surf(iband,ikpt,jkpt,kkpt,:,isppol)=elph_tr_ds%el_veloc(ikpt_fine,iband,:,isppol)*elph_ds%k_fine%wtk(iband,ikpt_fine,isppol) 961 v_surf(iband,ikpt,jkpt,kkpt,:,isppol)=elph_tr_ds%el_veloc(ikpt_fine,iband,:,isppol) 962 end do 963 end do 964 end do 965 end do 966 end do 967 968 call printvtk(eigenGS,v_surf,zero,elph_ds%fermie,Cryst%gprimd,& 969 & elph_ds%kptrlatt_fine,nband,hdr%nkpt,hdr%kptns,& 970 & Cryst%nsym,use_afm,Cryst%symrec,Cryst%symafm,use_tr,elph_ds%nsppol,shiftk,1,fname,ierr) 971 972 ABI_DEALLOCATE(v_surf) 973 974 end if !anaddb_dtset%prtfsurf 975 976 !============================================================================ 977 !Evaluate lambda and omega_log using the weighted sum over the irred q-points 978 !found in the GKK file. All the data we need are stored in elph_ds%qgrid_data 979 !============================================================================ 980 981 if (master == me) then 982 fname=trim(elph_ds%elph_base_name) // '_QPTS' 983 call outelph(elph_ds,anaddb_dtset%enunit,fname) 984 end if 985 986 !======================================================== 987 !Get FS averaged gamma matrices and Fourier transform to real space 988 !======================================================== 989 990 ABI_ALLOCATE(coskr, (elph_ds%nqpt_full,Ifc%nrpt)) 991 ABI_ALLOCATE(sinkr, (elph_ds%nqpt_full,Ifc%nrpt)) 992 call ftgam_init(ifc%gprim, elph_ds%nqpt_full,Ifc%nrpt, elph_ds%qpt_full, Ifc%rpt, coskr, sinkr) 993 994 call timein(tcpu,twall) 995 write(message, '(a,f11.3,a,f11.3,a)' )& 996 & '-elphon begin integration of gkq after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 997 call wrtout(std_out,message,'COLL') 998 tcpui = tcpu 999 twalli = twall 1000 1001 call integrate_gamma(elph_ds,FSfullpqtofull) 1002 1003 if (elph_ds%symgkq ==1) then 1004 ! complete the gamma_qpt here instead of the gkk previously 1005 call complete_gamma(Cryst,elph_ds%nbranch,elph_ds%nsppol,elph_ds%nqptirred,elph_ds%nqpt_full,& 1006 & elph_ds%ep_scalprod,elph_ds%qirredtofull,qpttoqpt,elph_ds%gamma_qpt) 1007 end if 1008 1009 !Now FT to real space too 1010 !NOTE: gprim (not gprimd) is used for all FT interpolations, 1011 !to be consistent with the dimensions of the rpt, which come from anaddb. 1012 ABI_ALLOCATE(elph_ds%gamma_rpt, (2,elph_ds%nbranch**2,elph_ds%nsppol,Ifc%nrpt)) 1013 elph_ds%gamma_rpt = zero 1014 1015 qtor = 1 ! q --> r 1016 do isppol=1,elph_ds%nsppol 1017 call ftgam(Ifc%wghatm,elph_ds%gamma_qpt(:,:,isppol,:),elph_ds%gamma_rpt(:,:,isppol,:),natom,& 1018 & elph_ds%nqpt_full,Ifc%nrpt,qtor, coskr, sinkr) 1019 end do 1020 1021 call timein(tcpu,twall) 1022 write(message, '(a,f11.3,a,f11.3,a)' )& 1023 & '-elphon end integration and completion of gkq after tcpu',tcpu-tcpui,' and twall',twall-twalli,' sec' 1024 call wrtout(std_out,message,'COLL') 1025 tcpui = tcpu 1026 twalli = twall 1027 1028 1029 !========================================================== 1030 !calculate transport matrix elements, integrated over FS 1031 !========================================================== 1032 1033 if (elph_tr_ds%ifltransport == 1)then ! LOVA 1034 1035 call integrate_gamma_tr_lova(elph_ds,FSfullpqtofull,elph_tr_ds) 1036 1037 call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,& 1038 & elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_trout,elph_ds%qirredtofull,qpttoqpt) 1039 1040 call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,& 1041 & elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_trin,elph_ds%qirredtofull,qpttoqpt) 1042 1043 ABI_ALLOCATE(elph_tr_ds%gamma_rpt_trout,(2,9,elph_ds%nbranch**2,elph_ds%nsppol,Ifc%nrpt)) 1044 elph_tr_ds%gamma_rpt_trout = zero 1045 1046 ABI_ALLOCATE(elph_tr_ds%gamma_rpt_trin,(2,9,elph_ds%nbranch**2,elph_ds%nsppol,Ifc%nrpt)) 1047 elph_tr_ds%gamma_rpt_trin = zero 1048 1049 ! Now FT to real space too 1050 qtor = 1 ! q --> r 1051 do isppol=1,elph_ds%nsppol 1052 do idir=1,9 1053 call ftgam(Ifc%wghatm,elph_tr_ds%gamma_qpt_trout(:,idir,:,isppol,:),& 1054 & elph_tr_ds%gamma_rpt_trout(:,idir,:,isppol,:),natom,& 1055 & elph_ds%nqpt_full,Ifc%nrpt,qtor, coskr, sinkr) 1056 1057 call ftgam(Ifc%wghatm,elph_tr_ds%gamma_qpt_trin(:,idir,:,isppol,:),& 1058 & elph_tr_ds%gamma_rpt_trin(:,idir,:,isppol,:),natom,& 1059 & elph_ds%nqpt_full,Ifc%nrpt,qtor, coskr, sinkr) 1060 end do 1061 end do 1062 1063 else if (elph_tr_ds%ifltransport==2) then ! non-LOVA case 1064 1065 ! Get Ef, DOS(Ef), veloc(Ef) for looped temperatures 1066 call get_nv_fs_temp(elph_ds,BSt,eigenGS_fine,cryst%gprimd,max_occ,elph_tr_ds) 1067 1068 ! Get DOS(E), veloc(E) for looped energy levels 1069 call get_nv_fs_en(cryst,ifc,elph_ds,eigenGS_fine,max_occ,elph_tr_ds,omega_max) 1070 1071 ! Save the E, N(E), v^2(E), dE 1072 if (master == me) then 1073 fname = trim(elph_ds%elph_base_name) // '_EPTS' 1074 if (open_file(fname,message,newunit=unit_epts,status="unknown") /=0) then 1075 MSG_ERROR(message) 1076 end if 1077 do isppol = 1, elph_ds%nsppol 1078 write(unit_epts,"(a,i6)") '# E, N(E), v^2(E), dE for spin channel ', isppol 1079 do ie1 = 1, elph_ds%nenergy 1080 write(unit_epts,"(4E20.12)") elph_tr_ds%en_all(isppol,ie1), elph_tr_ds%dos_n(ie1,isppol),& 1081 & elph_tr_ds%veloc_sq(1,isppol,ie1), elph_tr_ds%de_all(isppol,ie1) 1082 end do 1083 end do 1084 close(unit=unit_epts) 1085 end if 1086 1087 ABI_ALLOCATE(tmp_veloc_sq1,(3,elph_ds%nsppol)) 1088 ABI_ALLOCATE(tmp_veloc_sq2,(3,elph_ds%nsppol)) 1089 ABI_ALLOCATE(elph_tr_ds%tmp_gkk_intweight1,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,elph_ds%nsppol)) 1090 ABI_ALLOCATE(elph_tr_ds%tmp_gkk_intweight2,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,elph_ds%nsppol)) 1091 ABI_ALLOCATE(elph_tr_ds%tmp_velocwtk1,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,3,elph_ds%nsppol)) 1092 ABI_ALLOCATE(elph_tr_ds%tmp_velocwtk2,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,3,elph_ds%nsppol)) 1093 ABI_ALLOCATE(elph_tr_ds%tmp_vvelocwtk1,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,3,3,elph_ds%nsppol)) 1094 ABI_ALLOCATE(elph_tr_ds%tmp_vvelocwtk2,(elph_ds%ngkkband,elph_ds%k_phon%nkpt,3,3,elph_ds%nsppol)) 1095 1096 tmp_veloc_sq1 = zero 1097 tmp_veloc_sq2 = zero 1098 elph_tr_ds%tmp_gkk_intweight1 = zero 1099 elph_tr_ds%tmp_gkk_intweight2 = zero 1100 elph_tr_ds%tmp_velocwtk1 = zero 1101 elph_tr_ds%tmp_velocwtk2 = zero 1102 elph_tr_ds%tmp_vvelocwtk1 = zero 1103 elph_tr_ds%tmp_vvelocwtk2 = zero 1104 1105 if (elph_ds%ep_lova .eq. 1) then 1106 tmp_nenergy = 1 1107 else if (elph_ds%ep_lova .eq. 0) then 1108 tmp_nenergy = elph_ds%nenergy 1109 else 1110 write(message,'(a,i0)')' ep_lova must be 0 or 1 while it is : ', elph_ds%ep_lova 1111 MSG_ERROR(message) 1112 end if 1113 1114 ! This only works for ONE temperature!! for test only 1115 elph_ds%n0(:) = elph_tr_ds%dos_n0(1,:) 1116 1117 ! bxu, no need for complete sets of ie1 and ie2 1118 ! Only save those within the range of omega_max from Ef 1119 ABI_ALLOCATE(pair2red,(tmp_nenergy,tmp_nenergy)) 1120 pair2red = zero 1121 1122 elph_ds%n_pair = 0 1123 do ie1=1,tmp_nenergy 1124 e1 = elph_tr_ds%en_all(1,ie1) 1125 e2 = e1 - omega_max 1126 if (e2 .lt. elph_tr_ds%en_all(1,1)) then 1127 i_start = 1 1128 else 1129 i_start = 1 1130 diff = dabs(e2-elph_tr_ds%en_all(1,1)) 1131 do ie2 = 2, tmp_nenergy 1132 if (dabs(e2-elph_tr_ds%en_all(1,ie2)) .lt. diff) then 1133 diff = dabs(e2-elph_tr_ds%en_all(1,ie2)) 1134 i_start = ie2 1135 end if 1136 end do 1137 end if 1138 e2 = e1 + omega_max 1139 if (e2 .gt. elph_tr_ds%en_all(1,tmp_nenergy)) then 1140 i_end = tmp_nenergy 1141 else 1142 i_end = 1 1143 diff = dabs(e2-elph_tr_ds%en_all(1,1)) 1144 do ie2 = 2, tmp_nenergy 1145 if (dabs(e2-elph_tr_ds%en_all(1,ie2)) .lt. diff) then 1146 diff = dabs(e2-elph_tr_ds%en_all(1,ie2)) 1147 i_end = ie2 1148 end if 1149 end do 1150 end if 1151 do ie2 = i_start, i_end 1152 elph_ds%n_pair = elph_ds%n_pair + 1 1153 pair2red(ie1,ie2) = elph_ds%n_pair 1154 end do 1155 end do 1156 1157 ! symmetrize paire2red 1158 elph_ds%n_pair = 0 1159 do ie1 = 1, tmp_nenergy 1160 do ie2 = 1, tmp_nenergy 1161 if (pair2red(ie1,ie2) .ne. 0 .or. pair2red(ie2,ie1) .ne. 0) then 1162 elph_ds%n_pair = elph_ds%n_pair + 1 1163 pair2red(ie1,ie2) = elph_ds%n_pair 1164 end if 1165 end do 1166 end do 1167 1168 write(message,'(a,i3,a)')' There are ', elph_ds%n_pair, ' energy pairs. ' 1169 call wrtout(std_out,message,'COLL') 1170 1171 ABI_ALLOCATE(red2pair,(2,elph_ds%n_pair)) 1172 red2pair = zero 1173 elph_ds%n_pair = 0 1174 do ie1 = 1, tmp_nenergy 1175 do ie2 = 1, tmp_nenergy 1176 if (pair2red(ie1,ie2) .ne. 0 .or. pair2red(ie2,ie1) .ne. 0) then 1177 elph_ds%n_pair = elph_ds%n_pair + 1 1178 red2pair(1,elph_ds%n_pair) = ie1 1179 red2pair(2,elph_ds%n_pair) = ie2 1180 end if 1181 end do 1182 end do 1183 1184 ! moved from integrate_gamma_tr to here 1185 ABI_ALLOCATE(elph_tr_ds%gamma_qpt_tr,(2,9,elph_ds%nbranch**2,elph_ds%nsppol,elph_ds%nqpt_full)) 1186 ABI_ALLOCATE(elph_tr_ds%gamma_rpt_tr,(2,9,elph_ds%nbranch**2,elph_ds%nsppol,Ifc%nrpt,4,elph_ds%n_pair)) 1187 elph_tr_ds%gamma_rpt_tr = zero 1188 1189 s1ofssp = (/1,1,-1,-1/) 1190 s2ofssp = (/1,-1,1,-1/) 1191 1192 ! Get gamma 1193 do ie=1,elph_ds%n_pair 1194 ie1 = red2pair(1,ie) 1195 ie2 = red2pair(2,ie) 1196 1197 tmp_veloc_sq1(:,:)=elph_tr_ds%veloc_sq(:,:,ie1) 1198 elph_tr_ds%tmp_gkk_intweight1(:,:,:) = elph_tr_ds%tmp_gkk_intweight(:,:,:,ie1) 1199 elph_tr_ds%tmp_velocwtk1(:,:,:,:) = elph_tr_ds%tmp_velocwtk(:,:,:,:,ie1) 1200 elph_tr_ds%tmp_vvelocwtk1(:,:,:,:,:) = elph_tr_ds%tmp_vvelocwtk(:,:,:,:,:,ie1) 1201 1202 tmp_veloc_sq2(:,:)=elph_tr_ds%veloc_sq(:,:,ie2) 1203 elph_tr_ds%tmp_gkk_intweight2(:,:,:) = elph_tr_ds%tmp_gkk_intweight(:,:,:,ie2) 1204 elph_tr_ds%tmp_velocwtk2(:,:,:,:) = elph_tr_ds%tmp_velocwtk(:,:,:,:,ie2) 1205 elph_tr_ds%tmp_vvelocwtk2(:,:,:,:,:) = elph_tr_ds%tmp_vvelocwtk(:,:,:,:,:,ie2) 1206 1207 do ssp=1,4 ! (s,s'=+/-1, condense the indices) 1208 s1=s1ofssp(ssp) 1209 s2=s2ofssp(ssp) 1210 elph_tr_ds%gamma_qpt_tr = zero 1211 1212 call integrate_gamma_tr(elph_ds,FSfullpqtofull,s1,s2, & 1213 & tmp_veloc_sq1,tmp_veloc_sq2,elph_tr_ds) 1214 1215 call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,& 1216 & elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_tr,elph_ds%qirredtofull,qpttoqpt) 1217 1218 ! Now FT to real space too 1219 qtor = 1 ! q --> r 1220 do isppol=1,elph_ds%nsppol 1221 do idir=1,9 1222 call ftgam(Ifc%wghatm,elph_tr_ds%gamma_qpt_tr(:,idir,:,isppol,:),& 1223 & elph_tr_ds%gamma_rpt_tr(:,idir,:,isppol,:,ssp,ie),natom,& 1224 & elph_ds%nqpt_full,Ifc%nrpt,qtor,coskr, sinkr) 1225 end do 1226 end do 1227 1228 end do !ss 1229 end do !ie 1230 1231 ABI_DEALLOCATE(tmp_veloc_sq1) 1232 ABI_DEALLOCATE(tmp_veloc_sq2) 1233 end if ! ifltransport 1234 1235 ABI_DEALLOCATE(qpttoqpt) 1236 ABI_DEALLOCATE(FSfullpqtofull) 1237 1238 1239 !============================================================== 1240 !Calculate phonon linewidths, interpolating on chosen qpoints 1241 !============================================================== 1242 1243 call mkph_linwid(Cryst,ifc,elph_ds,anaddb_dtset%nqpath,anaddb_dtset%qpath) 1244 1245 !============================================================== 1246 !the nesting factor calculation 1247 !FIXME: this could go higher up, before the call to get_all_gkq 1248 !you only need the kpt and weight info 1249 !============================================================== 1250 if (any(anaddb_dtset%prtnest==[1,2])) then 1251 1252 nestname = trim(elph_ds%elph_base_name) // "_NEST" 1253 call mknesting(elph_ds%k_phon%nkpt,elph_ds%k_phon%kpt,elph_ds%kptrlatt,elph_ds%nFSband,& 1254 & elph_ds%k_phon%wtk,anaddb_dtset%nqpath,anaddb_dtset%qpath,elph_ds%nqpt_full, & 1255 & elph_ds%qpt_full,nestname,cryst%gprimd,cryst%gmet,anaddb_dtset%prtnest,qptrlatt) 1256 end if 1257 1258 !====================================================== 1259 !Calculate alpha^2 F integrating over fine kpt_phon grid 1260 !====================================================== 1261 1262 ABI_ALLOCATE(a2f_1d,(elph_ds%na2f)) 1263 ABI_ALLOCATE(dos_phon,(elph_ds%na2f)) 1264 1265 call mka2f(Cryst,Ifc,a2f_1d,dos_phon,elph_ds,elph_ds%kptrlatt_fine,elph_ds%mustar) 1266 1267 !calculate transport spectral function and coefficients 1268 if (elph_tr_ds%ifltransport==1 )then ! LOVA 1269 1270 call mka2f_tr_lova(cryst,ifc,elph_ds,elph_ds%ntemper,elph_ds%tempermin,elph_ds%temperinc,elph_tr_ds) 1271 1272 else if (elph_tr_ds%ifltransport==2 )then ! non LOVA 1273 1274 call mka2f_tr(cryst,ifc,elph_ds,elph_ds%ntemper,elph_ds%tempermin,elph_ds%temperinc,pair2red,elph_tr_ds) 1275 1276 ABI_DEALLOCATE(pair2red) 1277 ABI_DEALLOCATE(red2pair) 1278 1279 else if (elph_tr_ds%ifltransport==3 )then ! get k-dependent tau 1280 1281 call get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ) 1282 !call trans_rta(elph_ds,elph_tr_ds,cryst%gprimd,eigenGS,max_occ,cryst%ucvol) 1283 end if ! ifltransport 1284 1285 ABI_DEALLOCATE(eigenGS) 1286 ABI_DEALLOCATE(eigenGS_fine) 1287 1288 1289 !evaluate a2F only using the input Q-grid (without using interpolated matrices) 1290 !SCOPE: test the validity of the Fourier interpolation 1291 call wrtout(std_out,' elphon : calling mka2fQgrid',"COLL") 1292 1293 fname=trim(elph_ds%elph_base_name) // '_A2F_QGRID' 1294 call mka2fQgrid(elph_ds,fname) 1295 1296 !============================================= 1297 !Eliashberg equation in 1-D (isotropic case) 1298 !============================================= 1299 1300 call eliashberg_1d(a2f_1d,elph_ds,anaddb_dtset%mustar) 1301 1302 ABI_DEALLOCATE(a2f_1d) 1303 ABI_DEALLOCATE(dos_phon) 1304 1305 !MJV: 20070805 should exit here. None of the rest is tested or used yet to my knowledge 1306 1307 !======================================================================== 1308 !Now gkk contains the matrix elements of dH(1)/dxi i=1,2,3 1309 !for kpoints on the FS but qpoints only in the given grid {Q}. 1310 ! 1311 !1.) Need to complete the gkk elements for q and k\prime=k+q not 1312 !in the set of {k+Q} by Fourier interpolation on the Q. 1313 ! 1314 !2.) Need to complete the dynamical matrices and phonon freqs for 1315 !all q between points on the FS. 1316 ! 1317 !3.) With the eigenvectors e_ph of the dyn mats, do the scalar product 1318 !e_ph . gkk, which implies the gkk are turned to the eigenbasis of 1319 !the phonons. Before the (non eigen-) modes are ordered 1320 !atom1 xred1 atom1 xred2 atom1 xred3 1321 !atom2 xred1 atom2 xred2 atom2 xred3 ... 1322 !======================================================================= 1323 1324 make_gkk2=.false. 1325 1326 if (.not. make_gkk2) then 1327 call wrtout(std_out,' elphon : skipping full g(k,k") interpolation ',"COLL") 1328 else 1329 1330 ! ========================================================== 1331 ! FT of recip space gkk matrices to real space (gkk_rpt) 1332 ! NOTE: could be made into FFT, couldnt it? If shifts are 1333 ! used with a homogeneous grid 1334 ! ========================================================== 1335 write (message,'(2a,i0)')ch10,& 1336 & ' elphon : Fourier transform (q --> r) of the gkk matrices using nrpt = ',Ifc%nrpt 1337 call wrtout(std_out,message,'COLL') 1338 1339 call get_all_gkr(elph_ds,ifc%gprim,natom,Ifc%nrpt,onegkksize,Ifc%rpt,elph_ds%qpt_full,Ifc%wghatm) 1340 1341 ! ========================================================= 1342 ! complete gkk2 for all qpts between points 1343 ! on full kpt grid (interpolation from real space values) 1344 ! ========================================================= 1345 1346 write(message,'(2a)')ch10,& 1347 & ' elphon : Calling get_all_gkk2 to calculate gkk2 for q points over the full k grid' 1348 call wrtout(std_out,message,'COLL') 1349 1350 call get_all_gkk2(cryst,ifc,elph_ds,elph_ds%k_phon%kptirr,elph_ds%k_phon%kpt) 1351 end if 1352 1353 !===================================================== 1354 !Here should be the anisotropic Eliashberg equations. 1355 !===================================================== 1356 1357 !clean and deallocate junk 1358 call ebands_free(Bst) 1359 call elph_ds_clean(elph_ds) 1360 call elph_tr_ds_clean(elph_tr_ds) 1361 call hdr_free(hdr) 1362 1363 ABI_DEALLOCATE(coskr) 1364 ABI_DEALLOCATE(sinkr) 1365 1366 if (is_open(elph_ds%unitgkq)) close(elph_ds%unitgkq) 1367 1368 end subroutine elphon