TABLE OF CONTENTS


ABINIT/m_eph_driver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_eph_driver

FUNCTION

   Driver for EPH calculations

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_eph_driver
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_pspini,          only : pspini
34 
35  implicit none
36 
37  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 - superconductin 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.

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

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