TABLE OF CONTENTS
ABINIT/eph [ Functions ]
NAME
eph
FUNCTION
Main routine to compute electron phonon coupling matrix elements and calculate related properties - superconductin Tc, phonon linewidths, electronic renormalization due to phonons and temperature effects...
COPYRIGHT
Copyright (C) 2009-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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.
PARENTS
driver
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.
CHILDREN
crystal_free,crystal_from_hdr,crystal_print,cwtime,ddb_free ddb_from_file,ddk_free,ddk_init,destroy_mpi_enreg,dvdb_free,dvdb_init dvdb_interpolate_and_write,dvdb_list_perts,dvdb_print,ebands_free ebands_print,ebands_prtbltztrp,ebands_set_fermie,ebands_set_scheme ebands_update_occ,ebands_write,edos_free,edos_print,edos_write,eph_gkk eph_phgamma,eph_phpi,hdr_free,hdr_vs_dtset,ifc_free,ifc_init,ifc_mkphbs ifc_outphbtrap,ifc_print,ifc_printbxsf,ifc_test_phinterp init_distribfft_seq,initmpi_seq,mkphdos,pawfgr_destroy,pawfgr_init phdos_free,phdos_ncwrite,phdos_print,phdos_print_thermo,print_ngfft pspini,sigmaph,wfk_read_eigenvalues,wrtout,xmpi_bcast,xmpi_end
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 subroutine eph(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,xred) 76 77 use defs_basis 78 use defs_datatypes 79 use defs_abitypes 80 use m_profiling_abi 81 use m_xmpi 82 use m_xomp 83 use m_errors 84 use m_hdr 85 use m_crystal 86 use m_crystal_io 87 use m_ebands 88 use m_ddk 89 use m_ddb 90 use m_dvdb 91 use m_ifc 92 use m_phonons 93 use m_nctk 94 use m_wfk 95 #ifdef HAVE_NETCDF 96 use netcdf 97 #endif 98 99 use m_io_tools, only : file_exists 100 use m_time, only : cwtime 101 use m_fstrings, only : strcat, sjoin, ftoa, itoa 102 use m_fftcore, only : print_ngfft 103 use m_mpinfo, only : destroy_mpi_enreg 104 use m_pawang, only : pawang_type 105 use m_pawrad, only : pawrad_type 106 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 107 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 108 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 109 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init 110 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, symrhoij 111 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 112 use m_phgamma, only : eph_phgamma 113 use m_gkk, only : eph_gkk, ncwrite_v1qnu 114 use m_phpi, only : eph_phpi 115 use m_sigmaph, only : sigmaph 116 117 !This section has been created automatically by the script Abilint (TD). 118 !Do not modify the following lines by hand. 119 #undef ABI_FUNC 120 #define ABI_FUNC 'eph' 121 use interfaces_14_hidewrite 122 use interfaces_51_manage_mpi 123 use interfaces_56_io_mpi 124 use interfaces_64_psp 125 !End of the abilint section 126 127 implicit none 128 129 !Arguments ------------------------------------ 130 !scalars 131 character(len=6),intent(in) :: codvsn 132 type(datafiles_type),intent(in) :: dtfil 133 type(dataset_type),intent(in) :: dtset 134 type(pawang_type),intent(inout) :: pawang 135 type(pseudopotential_type),intent(inout) :: psps 136 !arrays 137 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom) 138 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 139 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 140 141 !Local variables ------------------------------ 142 !scalars 143 integer,parameter :: master=0,level40=40,natifc0=0,timrev2=2,selectz0=0 144 integer,parameter :: brav1=-1 ! WARNING. This choice is only to insure backwards compatibility with the tests, 145 !while eph is developed. Actually, should be switched to brav1=1 as soon as possible ... 146 integer,parameter :: nsphere0=0,prtsrlr0=0 147 integer :: ii,comm,nprocs,my_rank,psp_gencond,mgfftf,nfftf !,nfftf_tot 148 integer :: iblock,ddb_nqshift,ierr 149 integer :: omp_ncpus, work_size, nks_per_proc 150 real(dp):: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem !,ug_mem,ur_mem,cprj_mem 151 #ifdef HAVE_NETCDF 152 integer :: ncid,ncerr 153 #endif 154 real(dp),parameter :: rifcsph0=zero 155 real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff 156 real(dp) :: cpu,wall,gflops 157 logical :: use_wfk,use_wfq,use_dvdb 158 character(len=500) :: msg 159 character(len=fnlen) :: wfk0_path,wfq_path,ddb_path,dvdb_path,path 160 character(len=fnlen) :: ddk_path(3) 161 type(hdr_type) :: wfk0_hdr, wfq_hdr 162 type(crystal_t) :: cryst,cryst_ddb 163 type(ebands_t) :: ebands, ebands_kq 164 type(ddb_type) :: ddb 165 type(dvdb_t) :: dvdb 166 type(ddk_t) :: ddk 167 type(ifc_type) :: ifc 168 type(pawfgr_type) :: pawfgr 169 type(mpi_type) :: mpi_enreg 170 type(phonon_dos_type) :: phdos 171 !arrays 172 integer :: ngfftc(18),ngfftf(18) 173 integer,allocatable :: dummy_atifc(:) 174 real(dp),parameter :: k0(3)=zero 175 real(dp) :: dielt(3,3),zeff(3,3,dtset%natom) 176 real(dp),pointer :: gs_eigen(:,:,:) !,gs_occ(:,:,:) 177 real(dp),allocatable :: ddb_qshifts(:,:) 178 !real(dp) :: tsec(2) 179 !type(pawfgrtab_type),allocatable :: pawfgrtab(:) 180 !type(paw_ij_type),allocatable :: paw_ij(:) 181 !type(paw_an_type),allocatable :: paw_an(:) 182 183 !************************************************************************ 184 185 ! This part performs the initialization of basic objects used to perform e-ph calculations i.e: 186 ! 187 ! 1) crystal structure `cryst` 188 ! 2) Ground state band energies: `ebands` 189 ! 3) Interatomic force constants: `ifc` 190 ! 4) DVDB database with the dvscf potentials 191 ! 5) Pseudos and PAW basic objects. 192 ! 193 ! Once we have these objects, we can call specialized routines for e-ph calculations. 194 ! Notes: 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 specialized routines that will employ different MPI distribution schemes. 198 199 DBG_ENTER('COLL') 200 201 if (psps%usepaw == 1) then 202 MSG_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 211 ! Initialize filenames 212 wfk0_path = dtfil%fnamewffk 213 wfq_path = dtfil%fnamewffq 214 ddb_path = dtfil%filddbsin 215 dvdb_path = dtfil%filddbsin; ii=len_trim(dvdb_path); dvdb_path(ii-2:ii+1) = "DVDB" 216 use_wfk = (dtset%eph_task /= 5) 217 use_wfq = (dtset%irdwfq/=0 .or. dtset%getwfq/=0) 218 use_dvdb = (dtset%eph_task /= 0) 219 220 ddk_path(1) = strcat(dtfil%fnamewffddk, itoa(3*dtset%natom+1)) 221 ddk_path(2) = strcat(dtfil%fnamewffddk, itoa(3*dtset%natom+2)) 222 ddk_path(3) = strcat(dtfil%fnamewffddk, itoa(3*dtset%natom+3)) 223 224 if (my_rank == master) then 225 if (.not. file_exists(ddb_path)) MSG_ERROR(sjoin("Cannot find DDB file:", ddb_path)) 226 if (use_dvdb .and. .not. file_exists(dvdb_path)) MSG_ERROR(sjoin("Cannot find DVDB file:", dvdb_path)) 227 228 ! Accept WFK file in Fortran or netcdf format. 229 if (use_wfk .and. nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then 230 MSG_ERROR(sjoin("Cannot find GS WFK file:", wfk0_path, msg)) 231 end if 232 ! WFQ file 233 if (use_wfq) then 234 if (nctk_try_fort_or_ncfile(wfq_path, msg) /= 0) then 235 MSG_ERROR(sjoin("Cannot find GS WFQ file:", wfq_path, msg)) 236 end if 237 end if 238 239 if (dtset%eph_transport > 0) then 240 do ii=1,3 241 if (nctk_try_fort_or_ncfile(ddk_path(ii), msg) /= 0) then 242 MSG_ERROR(sjoin("Cannot find DDK file:", ddk_path(ii), msg)) 243 end if 244 end do 245 end if 246 247 end if ! master 248 249 ! Broadcast filenames (needed because they might have been changed if we are using netcdf files) 250 if (use_wfk) then 251 call xmpi_bcast(wfk0_path,master,comm,ierr) 252 call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path)) 253 end if 254 if (use_wfq) then 255 call xmpi_bcast(wfq_path,master,comm,ierr) 256 call wrtout(ab_out, sjoin("- Reading GS states from WFQ file:", wfq_path) ) 257 end if 258 call wrtout(ab_out, sjoin("- Reading DDB from file:", ddb_path)) 259 if (use_dvdb) call wrtout(ab_out, sjoin("- Reading DVDB from file:", dvdb_path)) 260 if (dtset%eph_transport > 0) then 261 call xmpi_bcast(ddk_path,master,comm,ierr) 262 call wrtout(ab_out, sjoin("- Reading DDK x from file:", ddk_path(1))) 263 call wrtout(ab_out, sjoin("- Reading DDK y from file:", ddk_path(2))) 264 call wrtout(ab_out, sjoin("- Reading DDK z from file:", ddk_path(3))) 265 ! Read header in DDK files and init basic dimensions. 266 ! subdrivers will use ddk to get the matrix elements from file. 267 call ddk_init(ddk, ddk_path, comm) 268 ! TODO: Should perform consistency check 269 !call hdr_vs_dtset(ddk_hdr(ii), dtset) 270 end if 271 272 ! autoparal section 273 ! TODO: This just to activate autoparal in abipy. Lot of things should be improved. 274 if (dtset%max_ncpus /=0) then 275 write(ab_out,'(a)')"--- !Autoparal" 276 write(ab_out,"(a)")"# Autoparal section for EPH runs" 277 write(ab_out,"(a)") "info:" 278 write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal 279 write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 280 write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt 281 write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol 282 write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor 283 write(ab_out,"(a,i0)")" mband: ",dtset%mband 284 write(ab_out,"(a,i0)")" eph_task: ",dtset%eph_task 285 286 work_size = dtset%nkpt * dtset%nsppol 287 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI. 288 nonscal_mem = zero 289 max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp 290 291 ! List of configurations. 292 ! Assuming an OpenMP implementation with perfect speedup! 293 write(ab_out,"(a)")"configurations:" 294 295 do ii=1,dtset%max_ncpus 296 nks_per_proc = work_size / ii 297 nks_per_proc = nks_per_proc + mod(work_size, ii) 298 eff = (one * work_size) / (ii * nks_per_proc) 299 ! Add the non-scalable part and increase by 10% to account for other datastructures. 300 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 301 302 do omp_ncpus=1,1 !xomp_get_max_threads() 303 write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 304 write(ab_out,"(a,i0)")" mpi_ncpus: ",ii 305 write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus 306 write(ab_out,"(a,f12.9)")" efficiency: ",eff 307 write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 308 end do 309 end do 310 write(ab_out,'(a)')"..." 311 MSG_ERROR_NODUMP("aborting now") 312 end if 313 314 call cwtime(cpu,wall,gflops,"start") 315 316 ! Construct crystal and ebands from the GS WFK file. 317 if (use_wfk) then 318 call wfk_read_eigenvalues(wfk0_path,gs_eigen,wfk0_hdr,comm) !,gs_occ) 319 call hdr_vs_dtset(wfk0_hdr,dtset) 320 321 call crystal_from_hdr(cryst,wfk0_hdr,timrev2) 322 call crystal_print(cryst,header="crystal structure from WFK file") 323 324 ebands = ebands_from_hdr(wfk0_hdr,maxval(wfk0_hdr%nband),gs_eigen) 325 call hdr_free(wfk0_hdr) 326 ABI_FREE(gs_eigen) 327 end if 328 329 ! Read WFQ and construct ebands on the shifted grid. 330 if (use_wfq) then 331 call wfk_read_eigenvalues(wfq_path,gs_eigen,wfq_hdr,comm) !,gs_occ) 332 ! GKA TODO: Have to construct a header with the proper set of q-shifted k-points then compare against file. 333 !call hdr_vs_dtset(wfq_hdr,dtset) 334 ebands_kq = ebands_from_hdr(wfq_hdr,maxval(wfq_hdr%nband),gs_eigen) 335 call hdr_free(wfq_hdr) 336 ABI_FREE(gs_eigen) 337 end if 338 339 ! Here we change the GS bands (fermi level, scissors operator ...) 340 ! All the modifications to ebands should be done here. 341 if (use_wfk) then 342 343 if (dtset%occopt /= ebands%occopt .or. abs(dtset%tsmear - ebands%tsmear) > tol12) then 344 write(msg,"(2a,2(a,i0,a,f14.6,a))")& 345 " Changing occupation scheme as input occopt and tsmear differ from those read from WFK file.",ch10,& 346 " From WFK file: occopt = ",ebands%occopt,", tsmear = ",ebands%tsmear,ch10,& 347 " From input: occopt = ",dtset%occopt,", tsmear = ",dtset%tsmear,ch10 348 call wrtout(ab_out,msg) 349 call ebands_set_scheme(ebands,dtset%occopt,dtset%tsmear,dtset%spinmagntarget,dtset%prtvol) 350 if (use_wfq) then 351 call ebands_set_scheme(ebands_kq,dtset%occopt,dtset%tsmear,dtset%spinmagntarget,dtset%prtvol) 352 end if 353 end if 354 355 if (dtset%eph_fermie /= zero) then ! default value of eph_fermie is zero hence no tolerance is used! 356 ABI_CHECK(abs(dtset%eph_extrael) <= tol12, "eph_fermie and eph_extrael are mutually exclusive") 357 call wrtout(ab_out, sjoin(" Fermi level set by the user at:",ftoa(dtset%eph_fermie))) 358 call ebands_set_fermie(ebands, dtset%eph_fermie, msg) 359 call wrtout(ab_out,msg) 360 if (use_wfq) then 361 call ebands_set_fermie(ebands_kq, dtset%eph_fermie, msg) 362 call wrtout(ab_out,msg) 363 end if 364 365 else if (abs(dtset%eph_extrael) > tol12) then 366 NOT_IMPLEMENTED_ERROR() 367 ! TODO: Be careful with the trick used in elphon for passing the concentration 368 !call ebands_set_nelect(ebands, dtset%eph_extrael, dtset%spinmagntarget, msg) 369 !call wrtout(ab_out,msg) 370 !if (use_wfq) then 371 ! call ebands_set_nelect(ebands_kq, dtset%eph_extrael, dtset%spinmagntarget, msg) 372 ! call wrtout(ab_out,msg) 373 !end if 374 end if 375 376 ! Recompute occupations. This is needed if WFK files have been produced in a NSCF run 377 ! since occ are set to zero, and fermie is taken from the previous density. 378 call ebands_update_occ(ebands, dtset%spinmagntarget, prtvol=dtset%prtvol) 379 call ebands_print(ebands,header="Ground state energies",prtvol=dtset%prtvol) 380 if (use_wfq) then 381 call ebands_update_occ(ebands_kq, dtset%spinmagntarget, prtvol=dtset%prtvol) 382 call ebands_print(ebands_kq,header="Ground state energies (K+Q)", prtvol=dtset%prtvol) 383 end if 384 385 end if 386 387 call cwtime(cpu,wall,gflops,"stop") 388 write(msg,'(2(a,f8.2))')"eph%init: cpu: ",cpu,", wall: ",wall 389 call wrtout(std_out, msg, do_flush=.True.) 390 call cwtime(cpu,wall,gflops,"start") 391 392 ! ======================================= 393 ! Output useful info on electronic bands 394 ! ======================================= 395 if (my_rank == master) then 396 ! Fermi Surface 397 if (dtset%prtfsurf /= 0) then 398 path = strcat(dtfil%filnam_ds(4), "_BXSF") 399 call wrtout(ab_out, sjoin("- Writing Fermi surface to file:", path)) 400 if (ebands_write_bxsf(ebands,cryst,path) /= 0) then 401 msg = "Cannot produce file for Fermi surface, check log file for more info" 402 MSG_WARNING(msg) 403 call wrtout(ab_out,msg) 404 end if 405 end if 406 407 ! Nesting factor (requires qpath) 408 if (dtset%prtnest /= 0 .and. dtset%ph_nqpath > 0) then 409 path = strcat(dtfil%filnam_ds(4), "_NEST") 410 call wrtout(ab_out, sjoin("- Writing nesting factor to file:", path)) 411 if (ebands_write_nesting(ebands,cryst,path,dtset%prtnest,& 412 dtset%tsmear,dtset%fermie_nest,dtset%ph_qpath(:,1:dtset%ph_nqpath),msg) /= 0) then 413 MSG_WARNING(msg) 414 call wrtout(ab_out,msg) 415 end if 416 end if 417 418 if (use_wfk) then 419 call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4)) 420 end if 421 end if 422 423 !if (.False.) then 424 !!if (.True.) then 425 ! !call ebands_set_interpolator(ebands, cryst, bstart, bcount, mode, espline_ords, eskw_ratio, comm) 426 ! call ebands_test_interpolator(ebands, dtset, cryst, dtfil%filnam_ds(4), comm) 427 ! MSG_ERROR("interpolation done") 428 !end if 429 430 call cwtime(cpu,wall,gflops,"stop") 431 write(msg,'(2(a,f8.2))')"eph%edos: cpu:",cpu,", wall: ",wall 432 call wrtout(std_out, msg, do_flush=.True.) 433 call cwtime(cpu,wall,gflops,"start") 434 435 ! Read the DDB file. 436 ABI_CALLOC(dummy_atifc, (dtset%natom)) 437 438 if (use_wfk) then 439 call ddb_from_file(ddb,ddb_path,brav1,dtset%natom,natifc0,dummy_atifc,cryst_ddb,comm, prtvol=dtset%prtvol) 440 call crystal_free(cryst_ddb) 441 else 442 call ddb_from_file(ddb,ddb_path,brav1,dtset%natom,natifc0,dummy_atifc,cryst,comm, prtvol=dtset%prtvol) 443 end if 444 ABI_FREE(dummy_atifc) 445 446 ddb_nqshift = 1 447 ABI_CALLOC(ddb_qshifts, (3,ddb_nqshift)) 448 449 ! Set the q-shift for the DDB 450 ddb_qshifts(:,1) = dtset%ddb_shiftq(:) 451 452 ! Get Dielectric Tensor and Effective Charges 453 ! (initialized to one_3D and zero if the derivatives are not available in the DDB file) 454 iblock = ddb_get_dielt_zeff(ddb,cryst,dtset%rfmeth,dtset%chneut,selectz0,dielt,zeff) 455 if (my_rank == master) then 456 if (iblock == 0) then 457 call wrtout(ab_out, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_path)) 458 call wrtout(ab_out, "Values initialized with zeros") 459 else 460 call wrtout(ab_out, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_path)) 461 end if 462 end if 463 464 ! Build the inter-atomic force constants. 465 ! WARNING : brav1 has been set to -1 at the initialisation of eph.F90, see the message there. Should be turned to 1 as soon as possible. 466 call ifc_init(ifc,cryst,ddb,& 467 brav1,dtset%asr,dtset%symdynmat,dtset%dipdip,dtset%rfmeth,dtset%ddb_ngqpt,ddb_nqshift,ddb_qshifts,dielt,zeff,& 468 nsphere0,rifcsph0,prtsrlr0,dtset%enunit,comm) 469 ABI_FREE(ddb_qshifts) 470 call ifc_print(ifc, unit=std_out) 471 472 ! Test B-spline interpolation of phonons 473 !if (.True.) then 474 if (.False.) then 475 call ifc_test_phinterp(ifc, cryst, [8,8,8], 1, [zero,zero,zero], [3,3,3], comm) 476 !call ifc_set_interpolator(ifc, cryst, nustart, nucount, mode, phspline_ords, phskw_ratio, comm) 477 !call ifc_test_intepolator(ifc, dtset, dtfil, comm) 478 call xmpi_end() 479 end if 480 481 ! Output phonon band structure (requires qpath) 482 if (dtset%prtphbands /= 0) call ifc_mkphbs(ifc, cryst, dtset, dtfil%filnam_ds(4), comm) 483 484 if (dtset%prtphdos == 1) then 485 486 ! Phonon Density of States. 487 ! FIXME: mkphdos expects qshift(3) instead of qshift(3, nqshift) 488 ! TODO: Parallelize this routine. 489 call mkphdos(phdos,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,dtset%ph_smear,dtset%ph_ngqpt,dtset%ph_qshift,comm) 490 491 if (my_rank == master) then 492 path = strcat(dtfil%filnam_ds(4), "_PHDOS") 493 call wrtout(ab_out, sjoin("- Writing phonon DOS to file:", path)) 494 call phdos_print(phdos, path) 495 496 !call phdos_print_debye(phdos, crystal%ucvol) 497 498 !TODO: do we want to pass the temper etc... from anaddb_dtset into the full dtset for abinit? 499 ! Otherwise just leave these defaults. 500 !MG: 1) Disabled for the time being because of SIGFPE in v8[41] 501 ! 2) I've addeded a new abinit variable (tmesh) to specifiy the list of temperatures. 502 path = strcat(dtfil%filnam_ds(4), "_MSQD_T") 503 !MG: Disabled for the time being because of SIGFPE in v8[41] 504 !call phdos_print_msqd(phdos, path, 1000, one, one) 505 path = strcat(dtfil%filnam_ds(4), "_THERMO") 506 call phdos_print_thermo(PHdos, path, 1000, zero, one) 507 508 #ifdef HAVE_NETCDF 509 path = strcat(dtfil%filnam_ds(4), "_PHDOS.nc") 510 ncerr = nctk_open_create(ncid, path, xmpi_comm_self) 511 NCF_CHECK_MSG(ncerr, sjoin("Creating PHDOS.nc file:", path)) 512 NCF_CHECK(crystal_ncwrite(cryst, ncid)) 513 call phdos_ncwrite(phdos, ncid) 514 NCF_CHECK(nf90_close(ncid)) 515 #endif 516 end if 517 call phdos_free(phdos) 518 end if 519 520 if (dtset%prtbltztrp == 1 .and. my_rank == master) then 521 call ifc_outphbtrap(ifc,cryst,dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,dtfil%filnam_ds(4)) 522 523 ! BoltzTraP output files in GENEric format 524 call ebands_prtbltztrp(ebands, cryst, dtfil%filnam_ds(4)) 525 end if 526 527 ! Output phonon isosurface in Xcrysden format. 528 if (dtset%prtphsurf == 1) then 529 path = strcat(dtfil%filnam_ds(4), "_PH.bxsf") 530 call wrtout(ab_out, sjoin("- Writing phonon frequencies in Xcrysden format to file:", path)) 531 call ifc_printbxsf(ifc, cryst, dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, path, comm) 532 end if 533 534 call cwtime(cpu,wall,gflops,"stop") 535 write(msg,'(2(a,f8.2))')"eph%ifc: cpu:",cpu,", wall: ",wall 536 call wrtout(std_out, msg, do_flush=.True.) 537 call cwtime(cpu,wall,gflops,"start") 538 539 ! Initialize the object used to read DeltaVscf (required if eph_task /= 0) 540 if (use_dvdb) then 541 call dvdb_init(dvdb, dvdb_path, comm) 542 if (my_rank == master) then 543 call dvdb_print(dvdb) 544 call dvdb_list_perts(dvdb, [-1,-1,-1], unit=ab_out) 545 end if 546 547 if (iblock /= 0) then 548 dvdb%dielt = dielt 549 dvdb%zeff = zeff 550 dvdb%has_dielt_zeff = .True. 551 end if 552 553 ! Compute \delta V_{q,nu)(r) and dump results to netcdf file. 554 if (.False. .and. my_rank == master) then 555 call ncwrite_v1qnu(dvdb, cryst, ifc, dvdb%nqpt, dvdb%qpts, dtset%prtvol, strcat(dtfil%filnam_ds(4), "_V1QNU.nc")) 556 end if 557 end if 558 559 ! TODO Recheck getng, should use same trick as that used in screening and sigma. 560 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 561 gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0) 562 563 call print_ngfft(ngfftc,header='Coarse FFT mesh used for the wavefunctions') 564 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 565 566 ! Fake MPI_type for the sequential part. 567 call initmpi_seq(mpi_enreg) 568 call init_distribfft_seq(mpi_enreg%distribfft,'c',ngfftc(2),ngfftc(3),'all') 569 call init_distribfft_seq(mpi_enreg%distribfft,'f',ngfftf(2),ngfftf(3),'all') 570 571 ! =========================================== 572 ! === Open and read pseudopotential files === 573 ! =========================================== 574 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,& 575 & pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm) 576 577 ! ==================================================== 578 ! === This is the real epc stuff once all is ready === 579 ! ==================================================== 580 ! TODO: decide whether to make several driver functions. 581 ! before that, however, need to encapsulate most of the functionalities in eph_phgamma 582 ! otherwise there will be tons of duplicated code 583 584 ! TODO: Make sure that all subdrivers work with useylm == 1 585 ABI_CHECK(dtset%useylm == 0, "useylm != 0 not implemented/tested") 586 select case (dtset%eph_task) 587 case (0) 588 continue 589 590 case (1) 591 ! Compute phonon linewidths in metals. 592 call eph_phgamma(wfk0_path,dtfil,ngfftc,ngfftf,dtset,cryst,ebands,dvdb,ddk,ifc,& 593 pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm) 594 595 case (2) 596 ! Compute electron-phonon matrix elements 597 call eph_gkk(wfk0_path,wfq_path,dtfil,ngfftc,ngfftf,dtset,cryst,ebands,ebands_kq,dvdb,ifc,& 598 pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm) 599 600 case (3) 601 ! Compute phonon self-energy 602 call eph_phpi(wfk0_path,wfq_path,dtfil,ngfftc,ngfftf,dtset,cryst,ebands,ebands_kq,dvdb,ifc,& 603 pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm) 604 605 case (4) 606 ! Compute electron self-energy (phonon contribution) 607 call sigmaph(wfk0_path,dtfil,ngfftc,ngfftf,dtset,cryst,ebands,dvdb,ifc,& 608 pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm) 609 610 case (5) 611 ! Interpolate the phonon potential 612 call dvdb_interpolate_and_write(dtfil,ngfftc,ngfftf,cryst,dvdb,& 613 & ifc%ngqpt,ifc%nqshft,ifc%qshft, & 614 & dtset%eph_ngqpt_fine,dtset%qptopt,mpi_enreg,comm) 615 616 case default 617 MSG_ERROR(sjoin("Unsupported value of eph_task:", itoa(dtset%eph_task))) 618 end select 619 620 !===================== 621 !==== Free memory ==== 622 !===================== 623 call crystal_free(cryst) 624 call dvdb_free(dvdb) 625 call ddb_free(ddb) 626 call ddk_free(ddk) 627 call ifc_free(ifc) 628 if (use_wfk) call ebands_free(ebands) 629 if (use_wfq) call ebands_free(ebands_kq) 630 call pawfgr_destroy(pawfgr) 631 call destroy_mpi_enreg(mpi_enreg) 632 633 ! Deallocation for PAW. 634 if (dtset%usepaw==1) then 635 !call pawrhoij_free(pawrhoij) 636 !ABI_DT_FREE(pawrhoij) 637 !call pawfgrtab_free(pawfgrtab) 638 !ABI_DT_FREE(pawfgrtab) 639 !call paw_ij_free(paw_ij) 640 !ABI_DT_FREE(paw_ij) 641 !call paw_an_free(paw_an) 642 !ABI_DT_FREE(paw_an) 643 end if 644 645 DBG_EXIT('COLL') 646 647 end subroutine eph