TABLE OF CONTENTS
ABINIT/m_eph_driver [ Modules ]
NAME
m_eph_driver
FUNCTION
Driver for EPH calculations
COPYRIGHT
Copyright (C) 2009-2024 ABINIT group (MG, MVer, GA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_eph_driver 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_xmpi 28 use m_xomp 29 use m_hdr 30 use m_crystal 31 use m_ebands 32 use m_dtset 33 use m_efmas_defs 34 use m_dtfil 35 use m_ddb 36 use m_ddb_hdr 37 use m_dvdb 38 use m_ifc 39 use m_phonons 40 use m_nctk 41 use m_wfk 42 use m_distribfft 43 use netcdf 44 45 use defs_datatypes, only : pseudopotential_type, ebands_t 46 use defs_abitypes, only : MPI_type 47 use m_io_tools, only : file_exists, open_file 48 use m_time, only : cwtime, cwtime_report 49 use m_fstrings, only : strcat, sjoin, ftoa, itoa 50 use m_fftcore, only : print_ngfft 51 use m_rta, only : rta_driver, ibte_driver 52 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 53 use m_pawang, only : pawang_type 54 use m_pawrad, only : pawrad_type 55 use m_pawtab, only : pawtab_type 56 use m_paw_an, only : paw_an_type, paw_an_free !, paw_an_nullify, paw_an_init, 57 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 58 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init 59 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_symrhoij 60 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 61 use m_phgamma, only : eph_phgamma 62 use m_efmas, only : efmasdeg_free_array, efmasval_free_array, efmas_ncread 63 use m_gkk, only : eph_gkk, ncwrite_v1qnu 64 use m_phpi, only : eph_phpi 65 use m_sigmaph, only : sigmaph 66 use m_pspini, only : pspini 67 use m_ephtk, only : ephtk_update_ebands 68 use m_gstore, only : gstore_t 69 use m_migdal_eliashberg, only : migdal_eliashberg_iso !, migdal_eliashberg_aniso 70 use m_berry_curvature, only : berry_curvature 71 use m_cumulant, only : cumulant_driver 72 use m_frohlich, only : frohlich_t, frohlichmodel_zpr, frohlichmodel_polaronmass 73 74 implicit none 75 76 private
m_eph_driver/eph [ Functions ]
[ Top ] [ m_eph_driver ] [ Functions ]
NAME
eph
FUNCTION
Main routine to compute electron phonon coupling matrix elements and calculate related properties - superconducting Tc, phonon linewidths, electronic renormalization due to phonons and temperature effects...
INPUTS
acell(3)=Length scales of primitive translations (bohr) codvsn=Code version dtfil<datafiles_type>=Variables related to files. dtset<dataset_type>=All input variables for this dataset. pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. Before entering the first time in the routine, a significant part of Psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=Dimensionless real space primitive translations. xred(3,natom)=Reduced atomic coordinates.
NOTES
ON THE USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
130 subroutine eph(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred) 131 132 !Arguments ------------------------------------ 133 !scalars 134 character(len=8),intent(in) :: codvsn 135 type(datafiles_type),intent(in) :: dtfil 136 type(dataset_type),intent(inout) :: dtset 137 type(pawang_type),intent(inout) :: pawang 138 type(pseudopotential_type),intent(inout) :: psps 139 !arrays 140 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom) 141 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 142 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 143 144 !Local variables ------------------------------ 145 !scalars 146 integer,parameter :: master = 0, selectz0 = 0, nsphere0 = 0, prtsrlr0 = 0, with_cplex1 = 1, with_cplex2 = 2 147 integer :: ii,comm,nprocs,my_rank,psp_gencond,mgfftf,nfftf 148 integer :: iblock_dielt_zeff, iblock_dielt, iblock_quadrupoles, ddb_nqshift, ierr, npert_miss 149 integer :: omp_ncpus, work_size, nks_per_proc, mtyp, mpert, lwsym !msize, 150 integer :: iatdir, iq2dir, iq1dir, quad_unt, iatom, jj, qptopt 151 integer :: ncid ! ,ncerr 152 real(dp):: eff, mempercpu_mb, max_wfsmem_mb, nonscal_mem 153 real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff 154 real(dp) :: cpu,wall,gflops 155 logical :: use_wfk, use_wfq, use_dvdb, use_sigeph 156 character(len=500) :: msg 157 character(len=fnlen) :: wfk0_path, wfq_path, ddb_filepath, dvdb_filepath, sigeph_filepath, path 158 type(hdr_type) :: wfk0_hdr, wfq_hdr 159 type(crystal_t) :: cryst, cryst_ddb 160 type(ebands_t) :: ebands, ebands_kq 161 type(ddb_type) :: ddb, ddb_lw 162 type(ddb_hdr_type) :: ddb_hdr 163 type(dvdb_t) :: dvdb 164 type(ifc_type) :: ifc 165 type(pawfgr_type) :: pawfgr 166 type(mpi_type) :: mpi_enreg 167 type(phdos_t) :: phdos 168 type(gstore_t) :: gstore 169 type(frohlich_t) :: frohlich 170 !arrays 171 integer :: ngfftc(18), ngfftf(18), count_wminmax(2), units(2) 172 real(dp),parameter :: k0(3)=zero 173 real(dp) :: wminmax(2), dielt(3,3), zeff(3,3,dtset%natom), zeff_raw(3,3,dtset%natom) 174 real(dp) :: qdrp_cart(3,3,3,dtset%natom) 175 real(dp),allocatable :: ddb_qshifts(:,:), kpt_efmas(:,:) 176 type(efmasdeg_type),allocatable :: efmasdeg(:) 177 type(efmasval_type),allocatable :: efmasval(:,:) 178 !type(pawfgrtab_type),allocatable :: pawfgrtab(:) 179 !type(paw_ij_type),allocatable :: paw_ij(:) 180 !type(paw_an_type),allocatable :: paw_an(:) 181 182 !************************************************************************ 183 184 ! This part performs the initialization of the basic objects used to perform e-ph calculations: 185 ! 186 ! 1) Crystal structure `cryst` 187 ! 2) Ground state band energies: `ebands` 188 ! 3) Interatomic force constants: `ifc` 189 ! 4) DVDB database with the dvscf potentials 190 ! 5) Pseudos and PAW basic objects. 191 ! 192 ! Once we have these objects, we can call specialized routines for e-ph calculations. 193 ! Notes: 194 ! 195 ! * Any modification to the basic objects mentioned above should be done here (e.g. change of efermi) 196 ! * This routines shall not allocate big chunks of memory. The CPU-demanding sections should be 197 ! performed in the subdriver that will employ different MPI distribution schemes optimized for that particular task. 198 199 DBG_ENTER('COLL') 200 201 if (psps%usepaw == 1) then 202 ABI_ERROR("PAW not implemented") 203 ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/)) 204 end if 205 206 ! abirules! 207 if (.False.) write(std_out,*)acell,codvsn,rprim,xred 208 209 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 210 units = [std_out, ab_out] 211 212 #ifndef HAVE_MPI_IBCAST 213 do ii=1,5 214 ABI_WARNING("Your MPI library does not provide MPI_IBCAST. Calculations parallelized over perturbations will be slow") 215 end do 216 #endif 217 218 ! Initialize filenames 219 wfk0_path = dtfil%fnamewffk 220 wfq_path = dtfil%fnamewffq 221 ddb_filepath = dtfil%filddbsin 222 223 ! Use the ddb file as prefix if getdvdb or irddvb are not given in the input. 224 dvdb_filepath = dtfil%fildvdbin 225 if (dvdb_filepath == ABI_NOFILE) then 226 dvdb_filepath = dtfil%filddbsin; ii = len_trim(dvdb_filepath); dvdb_filepath(ii-2:ii+1) = "DVDB" 227 end if 228 229 sigeph_filepath = dtfil%filsigephin 230 231 use_wfk = all(dtset%eph_task /= [0, 5, -5, 6, +15, -15, -16, 16]) 232 use_wfq = ((dtset%irdwfq /= 0 .or. dtset%getwfq /= 0 .or. dtset%getwfq_filepath /= ABI_NOFILE) .and. dtset%eph_frohlichm /= 1) 233 234 ! If eph_task is needed and ird/get variables are not provided, assume WFQ == WFK 235 if (any(dtset%eph_task == [2, -2, 3]) .and. .not. use_wfq) then 236 wfq_path = wfk0_path 237 use_wfq = .True. 238 write(msg, "(4a)")& 239 "eph_task requires WFQ but none among (irdwfq, getwfq, getwfk_filepath) is specified in input.", ch10, & 240 "Will read WFQ wavefunctions from WFK file:", trim(wfk0_path) 241 ABI_COMMENT(msg) 242 end if 243 244 use_dvdb = (dtset%eph_task /= 0 .and. dtset%eph_frohlichm /= 1 .and. abs(dtset%eph_task) /= 7) 245 use_sigeph = (dtset%eph_task == 9) 246 247 if (my_rank == master) then 248 ! GA: Let ddb object handle the error at reading time 249 !if (.not. file_exists(ddb_filepath)) ABI_ERROR(sjoin("Cannot find DDB file:", ddb_filepath)) 250 if (use_dvdb .and. .not. file_exists(dvdb_filepath)) ABI_ERROR(sjoin("Cannot find DVDB file:", dvdb_filepath)) 251 if (use_sigeph .and. .not. file_exists(sigeph_filepath)) ABI_ERROR(sjoin("Cannot find SIGEPH file:", sigeph_filepath)) 252 253 ! Accept WFK file in Fortran or netcdf format. 254 if (use_wfk .and. nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then 255 ABI_ERROR(sjoin("Cannot find GS WFK file:", wfk0_path, ". Error:", msg)) 256 end if 257 ! WFQ file 258 if (use_wfq) then 259 if (nctk_try_fort_or_ncfile(wfq_path, msg) /= 0) then 260 ABI_ERROR(sjoin("Cannot find GS WFQ file:", wfq_path, ". Error:", msg)) 261 end if 262 end if 263 end if ! master 264 265 ! Broadcast filenames (needed because they might have been changed if we are using netcdf files) 266 if (use_wfk) then 267 call xmpi_bcast(wfk0_path, master, comm, ierr) 268 call wrtout(units, sjoin("- Reading GS states from WFK file:", wfk0_path)) 269 end if 270 if (use_wfq) then 271 call xmpi_bcast(wfq_path, master, comm, ierr) 272 call wrtout(units, sjoin("- Reading GS states from WFQ file:", wfq_path) ) 273 end if 274 call wrtout(units, sjoin("- Reading DDB from file:", ddb_filepath)) 275 if (use_dvdb) call wrtout(units, sjoin("- Reading DVDB from file:", dvdb_filepath)) 276 if (dtset%eph_frohlichm /= 0) call wrtout(units, sjoin("- Reading EFMAS information from file:", dtfil%fnameabi_efmas)) 277 call wrtout(units, ch10//ch10) 278 279 ! autoparal section 280 ! TODO: This just to activate autoparal in AbiPy. Lot of things should be improved. 281 if (dtset%max_ncpus /= 0) then 282 write(ab_out,'(a)')"--- !Autoparal" 283 write(ab_out,"(a)")"# Autoparal section for EPH runs" 284 write(ab_out,"(a)") "info:" 285 write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal 286 write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 287 write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt 288 write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol 289 write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor 290 write(ab_out,"(a,i0)")" mband: ",dtset%mband 291 write(ab_out,"(a,i0)")" eph_task: ",dtset%eph_task 292 293 work_size = dtset%nkpt * dtset%nsppol 294 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI. 295 nonscal_mem = zero 296 max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp 297 298 ! List of configurations. 299 ! Assuming an OpenMP implementation with perfect speedup! 300 write(ab_out,"(a)")"configurations:" 301 302 do ii=1,dtset%max_ncpus 303 nks_per_proc = work_size / ii 304 nks_per_proc = nks_per_proc + mod(work_size, ii) 305 eff = (one * work_size) / (ii * nks_per_proc) 306 ! Add the non-scalable part and increase by 10% to account for other datastructures. 307 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 308 309 do omp_ncpus=1,1 !xomp_get_max_threads() 310 write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 311 write(ab_out,"(a,i0)")" mpi_ncpus: ",ii 312 write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus 313 write(ab_out,"(a,f12.9)")" efficiency: ",eff 314 write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 315 end do 316 end do 317 write(ab_out,'(a)')"..." 318 ABI_STOP("Stopping now!") 319 end if 320 321 call cwtime(cpu, wall, gflops, "start") 322 323 if (use_wfk) then 324 ! Construct crystal and ebands from the GS WFK file. 325 ebands = wfk_read_ebands(wfk0_path, comm, out_hdr=wfk0_hdr) 326 call wfk0_hdr%vs_dtset(dtset) 327 328 cryst = wfk0_hdr%get_crystal() 329 call cryst%print(header="crystal structure from WFK file") 330 331 ! Here we change the GS bands (Fermi level, scissors operator ...) 332 ! All the modifications to ebands should be done here. 333 call ephtk_update_ebands(dtset, ebands, "Ground state energies") 334 335 ! Need to update the WFK header to reflect the changes in ebands. 336 ! because we may need to write the header to ncfile 337 ! NB: eigenvalues are not stored in the header. 338 339 wfk0_hdr%occopt = ebands%occopt 340 call get_eneocc_vect(ebands, "occ", wfk0_hdr%occ) 341 wfk0_hdr%fermie = ebands%fermie 342 wfk0_hdr%nelect = ebands%nelect 343 end if 344 345 if (use_wfq) then 346 ! Read WFQ and construct ebands on the shifted grid. 347 ebands_kq = wfk_read_ebands(wfq_path, comm, out_hdr=wfq_hdr) 348 ! GKA TODO: Have to construct a header with the proper set of q-shifted k-points then compare against dtset. 349 !call wfq_hdr%vs_dtset(dtset) 350 call wfq_hdr%free() 351 call ephtk_update_ebands(dtset, ebands_kq, "Ground state energies (K+Q)") 352 end if 353 354 call cwtime_report(" eph%init", cpu, wall, gflops) 355 356 ! ======================================= 357 ! Output useful info on electronic bands 358 ! ======================================= 359 call cwtime(cpu, wall, gflops, "start") 360 361 if (my_rank == master) then 362 ! Fermi Surface 363 if (dtset%prtfsurf /= 0) then 364 path = strcat(dtfil%filnam_ds(4), "_BXSF") 365 call wrtout(units, sjoin("- Writing Fermi surface to file:", path)) 366 if (ebands_write_bxsf(ebands, cryst, path) /= 0) then 367 msg = "Cannot produce file for Fermi surface, check log file for more info" 368 ABI_WARNING(msg) 369 call wrtout(ab_out, msg) 370 end if 371 end if 372 373 ! Nesting factor (requires qpath) 374 if (dtset%prtnest /= 0 .and. dtset%ph_nqpath > 0) then 375 path = strcat(dtfil%filnam_ds(4), "_NEST") 376 call wrtout(ab_out, sjoin("- Writing nesting factor to file:", path)) 377 if (ebands_write_nesting(ebands, cryst, path, dtset%prtnest, & 378 dtset%tsmear, dtset%fermie_nest, dtset%ph_qpath(:,1:dtset%ph_nqpath), msg) /= 0) then 379 ABI_WARNING(msg) 380 call wrtout(ab_out,msg) 381 end if 382 end if 383 if (use_wfk) call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4)) 384 end if 385 386 call cwtime_report(" eph%ebands_postprocess:", cpu, wall, gflops) 387 388 ! Read the DDB file. 389 if (use_wfk) then 390 call ddb%from_file(ddb_filepath, ddb_hdr, cryst_ddb, comm, prtvol=dtset%prtvol) 391 392 call ddb%set_brav(dtset%brav) 393 394 ! DDB cryst comes from DFPT --> no time-reversal if q /= 0 395 ! Change the value so that we use the same as the GS part. 396 cryst_ddb%timrev = cryst%timrev 397 if (cryst%compare(cryst_ddb, header=" Comparing WFK crystal with DDB crystal") /= 0) then 398 ABI_ERROR("Crystal structure from WFK and DDB do not agree! Check messages above!") 399 end if 400 call cryst_ddb%free() 401 else 402 ! Get crystal from DDB. 403 ! Warning: We may loose precision in rprimd and xred because DDB in text format does not have enough significant digits. 404 call ddb%from_file(ddb_filepath, ddb_hdr, cryst, comm, prtvol=dtset%prtvol) 405 406 call ddb%set_brav(dtset%brav) 407 408 end if 409 410 ! Set the q-shift for the DDB (well we mainly use gamma-centered q-meshes) 411 ddb_nqshift = 1 412 ABI_CALLOC(ddb_qshifts, (3, ddb_nqshift)) 413 ddb_qshifts(:,1) = dtset%ddb_shiftq(:) 414 415 ! Get Dielectric Tensor 416 iblock_dielt = ddb%get_dielt(dtset%rfmeth, dielt) 417 418 ! Get Dielectric Tensor and Effective Charges 419 ! (initialized to one_3D and zero if the derivatives are not available in the DDB file) 420 iblock_dielt_zeff = ddb%get_dielt_zeff(cryst, dtset%rfmeth, dtset%chneut, selectz0, dielt, zeff, zeff_raw=zeff_raw) 421 if (my_rank == master) then 422 if (iblock_dielt_zeff == 0) then 423 call wrtout(units, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_filepath)) 424 call wrtout(units, " Values initialized with zeros.") 425 else 426 call wrtout(units, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_filepath)) 427 end if 428 end if 429 430 ! Read the quadrupoles 431 !iblock_quadrupoles = ddb%get_quadrupoles(ddb_hdr%ddb_version,1, 3, qdrp_cart) 432 433 iblock_quadrupoles = 0 434 qdrp_cart = zero 435 436 if (my_rank == master) then 437 ! MG: Temporary hack to read the quadrupole tensor from a text file 438 ! Will be removed when the EPH code will be able to read Q* from the DDB. 439 440 if (file_exists("quadrupoles_cart.out")) then 441 call wrtout(std_out, " Reading quadrupoles from quadrupoles_cart.out") 442 quad_unt = 71 443 open(unit=quad_unt,file="quadrupoles_cart.out",action="read") 444 do ii=1,2 445 read(quad_unt,*) msg 446 write(std_out, *)" msg: ", trim(msg) 447 end do 448 449 do ii=1,3 450 do jj=1,3*3*ddb%natom 451 read(quad_unt,'(4(i5,3x),2(1x,f20.10))') iq2dir,iatom,iatdir,iq1dir,qdrp_cart(iq1dir,iq2dir,iatdir,iatom) 452 write(std_out, *) iq2dir,iatom,iatdir,iq1dir,qdrp_cart(iq1dir,iq2dir,iatdir,iatom) 453 end do 454 read(quad_unt,'(a)') msg 455 end do 456 close(quad_unt) 457 iblock_quadrupoles = 1 458 end if 459 end if 460 call xmpi_bcast(iblock_quadrupoles, master, comm, ierr) 461 call xmpi_bcast(qdrp_cart, master, comm, ierr) 462 463 ! Here we get the quadrupoles from the DDB file (this should become the official API). 464 ! Section Copied from Anaddb. 465 466 if (iblock_quadrupoles == 0) then 467 mtyp = ddb_hdr%mblktyp 468 mpert = ddb_hdr%mpert 469 !msize = 3*mpert*3*mpert; if (mtyp==3) msize=msize*3*mpert 470 !call wrtout(std_out, sjoin(" Trying to read Q* from DDB file, mtyp:", itoa(mtyp))) 471 472 ! MR: a new ddb is necessary for the longwave quantities due to incompability of it with authomatic reshapes 473 ! that ddb%val and ddb%flg experience when passed as arguments of some routines 474 475 ! Get Quadrupole tensor 476 qdrp_cart=zero 477 if (mtyp==33) then 478 lwsym=1 479 call ddb_lw_copy(ddb, ddb_lw, mpert, dtset%natom, dtset%ntypat) 480 iblock_quadrupoles = ddb_lw%get_quadrupoles(ddb_hdr%ddb_version,lwsym, 33, qdrp_cart) 481 call ddb_lw%free() 482 end if 483 endif 484 485 ! The default value is 1. Here we set the flags to zero if Q* is not available. 486 if (iblock_quadrupoles == 0) then 487 dtset%dipquad = 0 488 dtset%quadquad = 0 489 end if 490 491 if (my_rank == master) then 492 if (iblock_quadrupoles == 0) then 493 call wrtout(units, sjoin("- Cannot find quadrupole tensor in DDB file:", ddb_filepath)) 494 call wrtout(units, " Values initialized with zeros.") 495 else 496 call wrtout(units, sjoin("- Found quadrupole tensor in DDB file:", ddb_filepath)) 497 end if 498 end if 499 500 call ifc%init(cryst, ddb, & 501 dtset%brav, dtset%asr, dtset%symdynmat, dtset%dipdip, dtset%rfmeth, & 502 dtset%ddb_ngqpt, ddb_nqshift, ddb_qshifts, dielt, zeff, & 503 qdrp_cart, nsphere0, dtset%rifcsph, prtsrlr0, dtset%enunit, comm, & 504 dipquad=dtset%dipquad, quadquad=dtset%quadquad) 505 506 ABI_FREE(ddb_qshifts) 507 if (my_rank == master) then 508 call ifc%print(unit=std_out) 509 !call ifc%print(unit=ab_out) 510 end if 511 512 ! Output phonon band structure (requires qpath) 513 if (dtset%prtphbands /= 0) call ifc_mkphbs(ifc, cryst, dtset, dtfil%filnam_ds(4), comm) 514 515 if (dtset%prtphdos == 1) then 516 call wrtout(std_out, " Computing Phonon DOS. Use prtphdos 0 to disable this part.") 517 wminmax = zero 518 do 519 call phdos%init(cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, dtset%ph_smear, dtset%ph_ngqpt, & 520 dtset%ph_nqshift, dtset%ph_qshift, "", wminmax, count_wminmax, comm) 521 if (all(count_wminmax == 0)) exit 522 wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05; wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05 523 call phdos%free() 524 write(msg, "(a, 2f8.5)") "Initial frequency mesh not large enough. Recomputing PHDOS with wmin, wmax: ",wminmax 525 call wrtout(std_out, msg) 526 end do 527 528 if (my_rank == master) then 529 if (dtset%prtvol > 0) then 530 ! Disabled by default because it's slow and we use netcdf that is much better. 531 path = strcat(dtfil%filnam_ds(4), "_PHDOS") 532 call wrtout(units, sjoin("- Writing phonon DOS to file:", path)) 533 call phdos%print(path) 534 end if 535 536 path = strcat(dtfil%filnam_ds(4), "_PHDOS.nc") 537 call wrtout(units, sjoin("- Writing phonon DOS to netcdf file:", path)) 538 NCF_CHECK_MSG(nctk_open_create(ncid, path, xmpi_comm_self), sjoin("Creating PHDOS.nc file:", path)) 539 NCF_CHECK(cryst%ncwrite(ncid)) 540 call phdos%ncwrite(ncid) 541 NCF_CHECK(nf90_close(ncid)) 542 end if 543 call phdos%free() 544 end if ! prtphdos 545 546 if (dtset%prtbltztrp == 1 .and. my_rank == master) then 547 call ifc%outphbtrap(cryst, dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, dtfil%filnam_ds(4)) 548 ! BoltzTraP output files in GENEric format 549 call ebands_prtbltztrp(ebands, cryst, dtfil%filnam_ds(4)) 550 end if 551 552 ! Output phonon isosurface in Xcrysden format. 553 if (dtset%prtphsurf == 1) then 554 path = strcat(dtfil%filnam_ds(4), "_PH.bxsf") 555 call wrtout(units, sjoin("- Writing phonon frequencies in Xcrysden format to file:", path)) 556 call ifc%printbxsf(cryst, dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, path, comm) 557 end if 558 559 call cwtime_report(" eph%ifc:", cpu, wall, gflops) 560 561 ! Initialize the object used to read DeltaVscf (required if eph_task /= 0) 562 if (use_dvdb) then 563 dvdb = dvdb_new(dvdb_filepath, comm) 564 565 ! DVDB cryst comes from DPPT --> no time-reversal if q /= 0 566 ! Change the value so that we use the same as the GS part. 567 dvdb%cryst%timrev = cryst%timrev 568 if (cryst%compare(dvdb%cryst, header=" Comparing WFK crystal with DVDB crystal") /= 0) then 569 ABI_ERROR("Crystal structure from WFK and DVDB do not agree! Check messages above!") 570 end if 571 if (dtset%prtvol > 10) dvdb%debug = .True. 572 573 ! This to symmetrize the DFPT potentials. 574 dvdb%symv1 = dtset%symv1scf 575 576 ! Copy brav variable 577 dvdb%brav = dtset%brav 578 579 ! Select algorithm for generating the list of R-points and the weigths used to compute W(r,R) 580 dvdb%rspace_cell = dtset%dvdb_rspace_cell 581 582 !call dvdb%load_ddb(dtset%prtvol, comm, ddb=ddb) 583 584 ! Set qdamp from frohl_params 585 dvdb%qdamp = dtset%dvdb_qdamp 586 587 ! Set quadrupoles 588 dvdb%qstar = qdrp_cart; if (iblock_quadrupoles /= 0) dvdb%has_quadrupoles = .True. 589 590 ! Set dielectric tensor, BECS and associated flags. 591 ! This flag activates automatically the treatment of the long-range term in the Fourier interpolation 592 ! of the DFPT potentials except when dvdb_add_lr == 0 593 dvdb%add_lr = dtset%dvdb_add_lr 594 if (iblock_dielt /= 0) then 595 dvdb%has_dielt = .True.; dvdb%dielt = dielt 596 end if 597 if (iblock_dielt_zeff /= 0) then 598 dvdb%has_zeff = .True.; dvdb%zeff = zeff; dvdb%zeff_raw = zeff_raw 599 end if 600 if (.not. dvdb%has_dielt .or. .not. (dvdb%has_zeff .or. dvdb%has_quadrupoles)) then 601 if (dvdb%add_lr /= 0) then 602 dvdb%add_lr = 0 603 ABI_WARNING("Setting dvdb_add_lr to 0. Long-range term won't be substracted in Fourier interpolation.") 604 end if 605 end if 606 607 if (dvdb%add_lr == 2) then 608 if (dvdb%has_quadrupoles) then 609 call wrtout(std_out, "dvdb_add_lr == 2 --> Quadrupoles are set to zero and won't be used in the interpolation") 610 end if 611 dvdb%has_quadrupoles = .False. 612 dvdb%qstar = zero 613 end if 614 615 if (my_rank == master) then 616 call dvdb%print() 617 call dvdb%list_perts([-1, -1, -1], npert_miss) 618 ABI_CHECK(npert_miss == 0, sjoin(itoa(npert_miss), "independent perturbation(s) are missing in the DVDB file!")) 619 end if 620 end if 621 622 call pawfgr_init(pawfgr, dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, & 623 gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=cryst%gmet, k0=k0) 624 625 call print_ngfft(ngfftc, header='Coarse FFT mesh used for the wavefunctions') 626 call print_ngfft(ngfftf, header='Dense FFT mesh used for densities and potentials') 627 628 ! Fake MPI_type for the sequential part. 629 call initmpi_seq(mpi_enreg) 630 call init_distribfft_seq(mpi_enreg%distribfft, 'c', ngfftc(2), ngfftc(3), 'all') 631 call init_distribfft_seq(mpi_enreg%distribfft, 'f', ngfftf(2), ngfftf(3), 'all') 632 633 ! I am not sure yet the EFMAS file will be needed as soon as eph_frohlichm/=0. To be decided later. 634 if (dtset%eph_frohlichm /= 0) then 635 NCF_CHECK(nctk_open_read(ncid, dtfil%fnameabi_efmas, xmpi_comm_self)) 636 call efmas_ncread(efmasdeg, efmasval, kpt_efmas, ncid) 637 NCF_CHECK(nf90_close(ncid)) 638 end if 639 640 ! =========================================== 641 ! === Open and read pseudopotential files === 642 ! =========================================== 643 call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm) 644 645 ! TODO: Make sure that all subdrivers work with useylm == 1 646 ABI_CHECK(dtset%useylm == 0, "useylm != 0 not implemented/tested") 647 648 ! Relase nkpt-based arrays in dtset to decrease memory requirement if dense sampling. 649 ! EPH routines should not access them after this point. 650 if (all(dtset%eph_task /= [6, 10])) call dtset%free_nkpt_arrays() 651 652 ! ==================================================== 653 ! === This is the real EPH stuff once all is ready === 654 ! ==================================================== 655 656 select case (dtset%eph_task) 657 case (0) 658 ! This is just to access the DDB post-processing tools 659 continue 660 661 case (1) 662 ! Compute phonon linewidths in metals. 663 call eph_phgamma(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 664 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 665 666 case (2, -2) 667 ! Compute e-ph matrix elements. 668 call eph_gkk(wfk0_path, wfq_path, dtfil, ngfftc, ngfftf, dtset, cryst, ebands, ebands_kq, dvdb, ifc, & 669 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 670 671 case (3) 672 ! Compute phonon self-energy. 673 call eph_phpi(wfk0_path, wfq_path, dtfil, ngfftc, ngfftf, dtset, cryst, ebands, ebands_kq, dvdb, ifc, & 674 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 675 676 case (4, -4) 677 ! Compute electron self-energy (phonon contribution) 678 call sigmaph(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, ebands, dvdb, ifc, wfk0_hdr, & 679 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 680 681 ! Compute transport properties in the RTA/IBTE only if sigma_erange has been used 682 if (dtset%eph_task == -4 .and. any(abs(dtset%sigma_erange) > zero)) then 683 if (dtset%ibte_prep > 0) then 684 call ibte_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) ! Solve IBTE 685 else 686 call rta_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) ! Compute RTA 687 end if 688 end if 689 690 case (5, -5) 691 ! Interpolate the phonon potential. 692 call dvdb%interpolate_and_write(dtset, dtfil%fnameabo_dvdb, ngfftc, ngfftf, cryst, & 693 ifc%ngqpt, ifc%nqshft, ifc%qshft, comm) 694 695 case (6) 696 ! Estimate zero-point renormalization and temperature-dependent electronic structure using the Frohlich model 697 if (my_rank == master) then 698 call frohlichmodel_zpr(frohlich, cryst, dtset, efmasdeg, efmasval, ifc) 699 end if 700 701 case (7) 702 ! Compute phonon-limited RTA from SIGEPH file. 703 call rta_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) 704 705 case (8) 706 ! Solve IBTE from SIGEPH file. 707 call ibte_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) 708 709 case (9) 710 ! Compute cumulant from SIGEPH.nc file. 711 call cumulant_driver(dtfil, dtset, ebands, cryst, comm) 712 713 case (10) 714 ! Estimate polaron effective mass in the triply-degenerate VB or CB cubic case 715 if (my_rank == master) then 716 call frohlichmodel_polaronmass(frohlich, cryst, dtset, efmasdeg, & 717 efmasval, ifc) 718 end if 719 720 case (11) 721 ! Write e-ph matrix elements to GSTORE file. 722 if (dtfil%filgstorein /= ABI_NOFILE) then 723 call wrtout(units, sjoin(" Restarting GSTORE computation from:", dtfil%filgstorein)) 724 call gstore%from_ncpath(dtfil%filgstorein, 1, dtset, cryst, ebands, ifc, comm) 725 else 726 path = strcat(dtfil%filnam_ds(4), "_GSTORE.nc") 727 call gstore%init(path, dtset, wfk0_hdr, cryst, ebands, ifc, comm) 728 call gstore%compute(wfk0_path, ngfftc, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 729 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 730 end if 731 call gstore%free() 732 733 case (12, -12) 734 ! Migdal-Eliashberg equations (isotropic/anisotropic case) 735 call gstore%from_ncpath(dtfil%filgstorein, with_cplex1, dtset, cryst, ebands, ifc, comm) 736 if (dtset%eph_task == -12) call migdal_eliashberg_iso(gstore, dtset, dtfil) 737 !if (dtset%eph_task == +12) call migdal_eliashberg_aniso(gstore, dtset, dtfil) 738 call gstore%free() 739 740 !case (13) 741 ! Variational polaron equations 742 !if (dtfil%filgstorein /= ABI_NOFILE) then 743 ! call wrtout(units, sjoin(" Computing variational polaron from pre-existent GSTORE file:", dtfil%filgstorein)) 744 ! call gstore%from_ncpath(dtfil%filgstorein, with_cplex2, dtset, cryst, ebands, ifc, comm) 745 !else 746 ! path = strcat(dtfil%filnam_ds(4), "_GSTORE.nc") 747 ! call wrtout(units, sjoin(" Computing GSTORE file:", dtfil%filgstorein)) 748 ! Customize input vars for this task. 749 ! dtset%gstore_qzone = "ibz"; dtset%gstore_kzone = "bz"; dtset%gstore_cplex = 2; dtset%gstore_with_vk = 1 750 ! call gstore%init(path, dtset, wfk0_hdr, cryst, ebands, ifc, comm) 751 ! call gstore%compute(wfk0_path, ngfftc, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 752 ! pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 753 ! call gstore%free() 754 ! call gstore%from_ncpath(path, with_cplex2, dtset, cryst, ebands, ifc, comm) 755 !end if 756 !call variational_polaron(gstore, dtset, dtfil) 757 !call gstore%free() 758 759 case (14) 760 ! Molecular Berry Curvature 761 if (dtfil%filgstorein /= ABI_NOFILE) then 762 call wrtout(units, sjoin(" Computing Berry curvature from pre-existent GSTORE file:", dtfil%filgstorein)) 763 call gstore%from_ncpath(dtfil%filgstorein, with_cplex2, dtset, cryst, ebands, ifc, comm) 764 else 765 path = strcat(dtfil%filnam_ds(4), "_GSTORE.nc") 766 call wrtout(units, sjoin(" Computing GSTORE file:", dtfil%filgstorein, "for Berry curvature from scracth")) 767 ! Customize input vars for this eph_task. 768 dtset%gstore_qzone = "ibz"; dtset%gstore_kzone = "bz"; dtset%gstore_cplex = 2; dtset%gstore_with_vk = 1 769 call gstore%init(path, dtset, wfk0_hdr, cryst, ebands, ifc, comm) 770 call gstore%compute(wfk0_path, ngfftc, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 771 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 772 call gstore%free() 773 call gstore%from_ncpath(path, with_cplex2, dtset, cryst, ebands, ifc, comm) 774 end if 775 776 call berry_curvature(gstore, dtset, dtfil) 777 call gstore%free() 778 779 case (15, -15) 780 ! Write average of DFPT potentials to file. 781 if (nprocs > 1) then 782 ABI_WARNING("eph_task in [15, -15] (average of DFPT potentials) does not support nprocs > 1. Running in sequential.") 783 end if 784 dvdb%comm = xmpi_comm_self 785 if (my_rank == master) then 786 call dvdb%open_read(ngfftf, xmpi_comm_self) 787 !call ephtk_set_pertables(cryst%natom, my_npert, pert_table, my_pinfo, comm) 788 !call dvdb%set_pert_distrib(sigma%comm_pert, sigma%my_pinfo, sigma%pert_table) 789 call dvdb%write_v1qavg(dtset, strcat(dtfil%filnam_ds(4), "_V1QAVG.nc")) 790 end if 791 792 case (-16, 16) 793 if (nprocs > 1) then 794 ABI_WARNING("eph_task in [16, -16] (test_phrotation) does not support nprocs > 1. Running in sequential.") 795 end if 796 797 qptopt = ebands%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt 798 call test_phrotation(ifc, cryst, qptopt, dtset%ph_ngqpt, comm) 799 800 dvdb%comm = xmpi_comm_self 801 if (my_rank == master) then 802 call dvdb%open_read(ngfftf, xmpi_comm_self) 803 804 ! Compute \delta V_{q,nu)(r) and dump results to netcdf file. 805 call ncwrite_v1qnu(dvdb, dtset, ifc, strcat(dtfil%filnam_ds(4), "_V1QNU.nc")) 806 end if 807 808 case default 809 ABI_ERROR(sjoin("Unsupported value of eph_task:", itoa(dtset%eph_task))) 810 end select 811 812 !===================== 813 !==== Free memory ==== 814 !===================== 815 call cryst%free() 816 call dvdb%free() 817 call ddb%free() 818 call ddb_hdr%free() 819 call ifc%free() 820 call wfk0_hdr%free() 821 call ebands_free(ebands) 822 call ebands_free(ebands_kq) 823 call pawfgr_destroy(pawfgr) 824 call destroy_mpi_enreg(mpi_enreg) 825 826 if (allocated(efmasdeg)) call efmasdeg_free_array(efmasdeg) 827 if (allocated(efmasval)) call efmasval_free_array(efmasval) 828 ABI_SFREE(kpt_efmas) 829 830 ! Deallocation for PAW. 831 if (dtset%usepaw == 1) then 832 !call pawrhoij_free(pawrhoij) 833 !ABI_FREE(pawrhoij) 834 !call pawfgrtab_free(pawfgrtab) 835 !ABI_FREE(pawfgrtab) 836 !call paw_ij_free(paw_ij) 837 !ABI_FREE(paw_ij) 838 !call paw_an_free(paw_an) 839 !ABI_FREE(paw_an) 840 end if 841 842 DBG_EXIT('COLL') 843 844 end subroutine eph