TABLE OF CONTENTS


ABINIT/m_ddk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddk

FUNCTION

  Objects and methods to extract data from DDK files.
  The DDK files are binary (soon also netcdf) files with Hamiltonian derivatives
  wrt k, and the corresponding wave functions

COPYRIGHT

 Copyright (C) 2016-2018 ABINIT group (MJV)
 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.

PARENTS

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_ddk
29 
30  use defs_basis
31  use m_abicore
32  use m_errors
33  use m_xmpi
34  use m_nctk
35  use m_hdr
36  use m_kptrank
37  use m_fstab
38  use m_wfk
39 #ifdef HAVE_NETCDF
40  use netcdf
41 #endif
42 
43  use m_fstrings,      only : sjoin, itoa, endswith
44  use m_symtk,         only : matr3inv
45  use m_io_tools,      only : iomode_from_fname
46  use defs_abitypes,   only : hdr_type
47  use m_geometry,      only : mkradim
48  use m_crystal,       only : crystal_t, crystal_free
49  use m_crystal_io,    only : crystal_from_hdr
50  use defs_datatypes,  only : ebands_t
51 
52  implicit none
53 
54  private

m_ddk/ddk_free [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_free

FUNCTION

 Close the file and release the memory allocated.

PARENTS

      eph

CHILDREN

      wrtout

SOURCE

794 subroutine ddk_free(ddk)
795 
796 
797 !This section has been created automatically by the script Abilint (TD).
798 !Do not modify the following lines by hand.
799 #undef ABI_FUNC
800 #define ABI_FUNC 'ddk_free'
801 !End of the abilint section
802 
803  implicit none
804 
805 !Arguments ------------------------------------
806 !scalars
807  type(ddk_t),intent(inout) :: ddk
808 
809 !************************************************************************
810 
811  ! integer arrays
812 
813  ! real arrays
814  if (allocated(ddk%velocity)) then
815    ABI_DEALLOCATE(ddk%velocity)
816  end if
817  if (allocated(ddk%velocity_fsavg)) then
818    ABI_DEALLOCATE(ddk%velocity_fsavg)
819  end if
820 
821  ! types
822  call crystal_free(ddk%cryst)
823 
824 end subroutine ddk_free

m_ddk/ddk_fs_average_veloc [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_fs_average_veloc

FUNCTION

  Perform Fermi surface average of velocity squared then square rooted, print and store in ddk object

INPUTS

   ddk = object with electron band velocities
   fstab(ddk%nsppol)=Tables with the correspondence between points of the Fermi surface (FS)
     and the k-points in the IBZ
   comm=MPI communicator

PARENTS

      m_phgamma

CHILDREN

      wrtout

SOURCE

707 subroutine ddk_fs_average_veloc(ddk, ebands, fstab, sigmas)
708 
709 
710 !This section has been created automatically by the script Abilint (TD).
711 !Do not modify the following lines by hand.
712 #undef ABI_FUNC
713 #define ABI_FUNC 'ddk_fs_average_veloc'
714 !End of the abilint section
715 
716  implicit none
717 
718 !Arguments ------------------------------------
719 !scalars
720 !integer,intent(in) :: comm  ! could distribute this over k in the future
721  real(dp),intent(in) :: sigmas(:)
722  type(ebands_t),intent(in) :: ebands
723  type(ddk_t),intent(inout) :: ddk
724  type(fstab_t),target,intent(in) :: fstab(ddk%nsppol)
725 
726 !Local variables-------------------------------
727 !scalars
728  integer :: idir, ikfs, isppol, ik_ibz, iene
729  integer :: iband
730  integer :: mnb, nband_k
731  integer :: nsig
732  type(fstab_t), pointer :: fs
733 !arrays
734  real(dp), allocatable :: wtk(:,:)
735 
736 !************************************************************************
737 
738  ddk%nene = fstab(1)%nene
739  ABI_MALLOC(ddk%velocity_fsavg, (3,ddk%nene,ddk%nsppol))
740  ddk%velocity_fsavg = zero
741 
742  nsig = size(sigmas, dim=1)
743  mnb = 1
744  do isppol=1,ddk%nsppol
745    fs => fstab(isppol)
746    mnb = max(mnb, maxval(fs%bstcnt_ibz(2, :)))
747  end do
748  ABI_MALLOC(wtk, (nsig,mnb))
749 
750  do isppol=1,ddk%nsppol
751    fs => fstab(isppol)
752    do iene = 1, fs%nene
753      do ikfs=1,fs%nkfs
754        ik_ibz = fs%istg0(1,ikfs)
755        nband_k = fs%bstcnt_ibz(2, ik_ibz)
756        call fstab_weights_ibz(fs, ebands, ik_ibz, isppol, sigmas, wtk, iene)
757 
758        do idir = 1,3
759          do iband = 1, nband_k
760            ddk%velocity_fsavg(idir, iene, isppol) = ddk%velocity_fsavg(idir, iene, isppol) + &
761 &             wtk(1,iband) * ddk%velocity(idir, iband, ikfs, isppol)**2
762 !&             fs%tetra_wtk_ene(iband,ik_ibz,iene) * ddk%velocity(idir, iband, ikfs, isppol)**2
763          end do
764        end do ! idir
765      end do ! ikfs
766    end do ! iene
767   ! sqrt is element wise on purpose
768    ddk%velocity_fsavg(:,:,isppol) = sqrt(ddk%velocity_fsavg(:,:,isppol)) / dble(fs%nkfs)
769  end do ! isppol
770 
771  ABI_DEALLOCATE(wtk)
772 
773 end subroutine ddk_fs_average_veloc

m_ddk/ddk_init [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_init

FUNCTION

  Initialize the object from file. This is a COLLECTIVE procedure that must be called
  by each process in the MPI communicator comm.

INPUTS

   paths=3 Filenames (could be either DFPT WFK files of DKK.nc files.
   comm=MPI communicator.

PARENTS

      eph

CHILDREN

      wrtout

SOURCE

178 subroutine ddk_init(ddk, paths, comm)
179 
180 
181 !This section has been created automatically by the script Abilint (TD).
182 !Do not modify the following lines by hand.
183 #undef ABI_FUNC
184 #define ABI_FUNC 'ddk_init'
185 !End of the abilint section
186 
187  implicit none
188 
189 !Arguments ------------------------------------
190 !scalars
191  character(len=*),intent(in) :: paths(3)
192  integer,intent(in) :: comm
193  type(ddk_t),intent(inout) :: ddk
194 
195 !Local variables-------------------------------
196 !scalars
197  integer,parameter :: timrev2=2,fform2=2
198  integer :: my_rank, restart, restartpaw, ii
199 !arrays
200  integer :: fforms(3)
201  type(hdr_type) :: hdrs(3)
202 
203 !************************************************************************
204 
205  my_rank = xmpi_comm_rank(comm)
206  ddk%paths = paths; ddk%comm = comm
207  ddk%iomode = iomode_from_fname(paths(1))
208 
209  ! In this calls everything is broadcast properly to the whole comm
210  do ii=1,3
211    ddk%use_ncddk(ii) = endswith(paths(ii), "_EVK.nc")
212    call hdr_read_from_fname(hdrs(ii), paths(ii), fforms(ii), comm)
213    if (ddk%debug) call hdr_echo(hdrs(ii), fforms(ii), 4, unit=std_out)
214    ! check that 2 headers are compatible
215    if (ii > 1) call hdr_check(fform2, fform2, hdrs(ii-1), hdrs(ii), 'COLL', restart, restartpaw)
216  end do
217 
218  ddk%nsppol = hdrs(1)%nsppol
219  ddk%nspinor = hdrs(1)%nspinor
220  ddk%usepaw = hdrs(1)%usepaw
221  ABI_CHECK(ddk%usepaw == 0, "PAW not yet supported")
222 
223  ! Init crystal_t
224  call crystal_from_hdr(ddk%cryst, hdrs(1), timrev2)
225 
226  ! Compute rprim, and gprimd. Used for slow FFT q--r if multiple shifts
227  call mkradim(ddk%acell,ddk%rprim,ddk%cryst%rprimd)
228  call matr3inv(ddk%rprim,ddk%gprim)
229 
230  do ii=1,3
231    call hdr_free(hdrs(ii))
232  end do
233 
234 end subroutine ddk_init

m_ddk/ddk_print [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_print

FUNCTION

  Print info on the object.

INPUTS

 [unit]=the unit number for output
 [prtvol]=verbosity level
 [mode_paral]=either "COLL" or "PERS"

OUTPUT

  Only printing.

PARENTS

CHILDREN

      wrtout

SOURCE

851 subroutine ddk_print(ddk, header, unit, prtvol, mode_paral)
852 
853 
854 !This section has been created automatically by the script Abilint (TD).
855 !Do not modify the following lines by hand.
856 #undef ABI_FUNC
857 #define ABI_FUNC 'ddk_print'
858 !End of the abilint section
859 
860  implicit none
861 
862 !Arguments ------------------------------------
863 !scalars
864  integer,optional,intent(in) :: prtvol,unit
865  character(len=4),optional,intent(in) :: mode_paral
866  character(len=*),optional,intent(in) :: header
867  type(ddk_t),intent(in) :: ddk
868 
869 !Local variables-------------------------------
870 !scalars
871  integer :: my_unt,my_prtvol
872  character(len=4) :: my_mode
873  character(len=500) :: msg
874 
875 ! *************************************************************************
876 
877  my_unt =std_out; if (PRESENT(unit)) my_unt   =unit
878  my_prtvol=0    ; if (PRESENT(prtvol)) my_prtvol=prtvol
879  my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
880 
881  msg=' ==== Info on the ddk% object ==== '
882  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
883  call wrtout(my_unt,msg,my_mode)
884 
885  write(std_out,*)"Number of FS bands: ",ddk%maxnb
886  write(std_out,*)"Number of FS k-points: ",ddk%nkfs
887  write(std_out,*)"Number of spin channels: ",ddk%nsppol
888  write(std_out,*)"Paths to files: "
889  write(std_out,*)"  ", trim(ddk%paths(1))
890  write(std_out,*)"  ", trim(ddk%paths(2))
891  write(std_out,*)"  ", trim(ddk%paths(3))
892 
893 end subroutine ddk_print

m_ddk/ddk_read_fsvelocities [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  ddk_read_fsvelocities

FUNCTION

  Read FS velocities from DDK files. Returned in ddk%velocity

INPUTS

   fstab(ddk%nsppol)=Tables with the correspondence between points of the Fermi surface (FS)
     and the k-points in the IBZ
   comm=MPI communicator

PARENTS

      m_phgamma

CHILDREN

      wrtout

SOURCE

565 subroutine ddk_read_fsvelocities(ddk, fstab, comm)
566 
567 
568 !This section has been created automatically by the script Abilint (TD).
569 !Do not modify the following lines by hand.
570 #undef ABI_FUNC
571 #define ABI_FUNC 'ddk_read_fsvelocities'
572 !End of the abilint section
573 
574  implicit none
575 
576 !Arguments ------------------------------------
577 !scalars
578  integer,intent(in) :: comm
579  type(ddk_t),intent(inout) :: ddk
580  type(fstab_t),target,intent(in) :: fstab(ddk%nsppol)
581 
582 !Local variables-------------------------------
583 !scalars
584  integer :: idir, ikfs, isppol, ik_ibz, ii
585  integer :: symrankkpt, ikpt_ddk,iband, bd2tot_index
586  integer :: bstart_k, nband_k, nband_in, vdim
587 #ifdef HAVE_NETCDF
588  integer :: ncid, varid, nc_fform
589 #endif
590  type(hdr_type) :: hdr1
591  type(kptrank_type) :: kptrank_t
592  type(fstab_t), pointer :: fs
593  character(len=500) :: msg
594 !arrays
595  real(dp), allocatable :: eigen1(:)
596  real(dp), allocatable :: velocityp(:,:)
597 
598 !************************************************************************
599 
600  if (ddk%rw_mode /= ddk_NOMODE) then
601    MSG_ERROR("ddk should be in ddk_NOMODE before open_read is called.")
602  end if
603  ddk%rw_mode = ddk_READMODE
604 
605  ddk%maxnb = maxval(fstab(:)%maxnb)
606  ddk%nkfs = maxval(fstab(:)%nkfs)
607  ABI_MALLOC(ddk%velocity, (3,ddk%maxnb,ddk%nkfs,ddk%nsppol))
608 
609  call wrtout(std_out, sjoin('Read DDK FILES with iomode=', itoa(ddk%iomode)), 'COLL')
610  do idir = 1,3
611    ! Open the files. All procs in comm receive hdr1 and eigen1
612    if (ddk%use_ncddk(idir)) then
613 #ifdef HAVE_NETCDF
614      NCF_CHECK(nctk_open_read(ncid, ddk%paths(ii), comm))
615      call hdr_ncread(hdr1, ncid, nc_fform)
616      varid = nctk_idname(ncid, "h1_matrix_elements")
617      ABI_MALLOC(eigen1, (2*hdr1%mband*hdr1%mband*hdr1%nkpt*ddk%nsppol))
618      !NCF_CHECK(nctk_set_collective(ncid, varid))
619      NCF_CHECK(nf90_get_var(ncid, varid, eigen1, count=[2, hdr1%mband, hdr1%mband, hdr1%nkpt, hdr1%nsppol]))
620      NCF_CHECK(nf90_close(ncid))
621 #else
622      MSG_ERROR("Netcdf not available!")
623 #endif
624    else
625      call wfk_read_h1mat(ddk%paths(idir), eigen1, hdr1, comm)
626    end if
627    nband_in = maxval(hdr1%nband)
628 
629    ! need correspondence hash between the DDK and the fs k-points
630    call mkkptrank (hdr1%kptns,hdr1%nkpt,kptrank_t)
631    do isppol=1,ddk%nsppol
632      fs => fstab(isppol)
633      do ikfs=1,fs%nkfs
634        ik_ibz = fs%istg0(1,ikfs)
635        call get_rank_1kpt (fs%kpts(:,ikfs),symrankkpt, kptrank_t)
636        ikpt_ddk = kptrank_t%invrank(symrankkpt)
637        if (ikpt_ddk == -1) then
638          write(msg, "(3a)")&
639            "Error in correspondence between ddk and fsk kpoint sets",ch10,&
640            "kpt sets in fsk and ddk files must agree."
641          MSG_ERROR(msg)
642        end if
643 
644        bd2tot_index=2*nband_in**2*(ikpt_ddk-1) + 2*nband_in**2*hdr1%nkpt*(isppol-1)
645        bstart_k = fs%bstcnt_ibz(1, ik_ibz)
646        nband_k = fs%bstcnt_ibz(2, ik_ibz)
647        ! first derivative eigenvalues for k-point. Diagonal of eigen1 is real -> only use that part
648        do iband = bstart_k, bstart_k+nband_k-1
649          ddk%velocity(idir, iband-bstart_k+1, ikfs, isppol)=eigen1(bd2tot_index + 2*nband_in*(iband-1) + iband)
650        end do
651      end do
652    end do
653 
654    ABI_FREE(eigen1)
655    call destroy_kptrank(kptrank_t)
656    call hdr_free(hdr1)
657  end do ! idir
658 
659  ! process the eigenvalues(1): rotate to cartesian and divide by 2 pi
660  ! use DGEMM better here on whole matrix and then reshape?
661  vdim = ddk%maxnb*ddk%nkfs*ddk%nsppol
662  ABI_MALLOC(velocityp, (3,vdim))
663  velocityp = zero
664  call dgemm('n','n',3,vdim,3,one,ddk%cryst%rprimd,3,ddk%velocity,3,zero,velocityp,3)
665 
666 ! do isppol = 1, ddk%nsppol
667 !   do ikpt = 1, ddk%nkfs
668 !     do iband = 1, ddk%maxnb
669 !       tmpveloc = ddk%velocity (:, iband, ikpt, isppol)
670 !       call dgemm('n','n',3,3,3,one,ddk%cryst%rprimd,3,tmpveloc,3,zero,tmpveloc2,3)
671 !       ddk%velocity (:, iband, ikpt, isppol) = tmpveloc2
672 !     end do
673 !   end do
674 ! end do
675 
676  ddk%velocity = reshape (velocityp, [3,ddk%maxnb,ddk%nkfs,ddk%nsppol])
677 
678  ABI_FREE(velocityp)
679  call destroy_kptrank (kptrank_t)
680 
681 end subroutine ddk_read_fsvelocities

m_ddk/ddk_t [ Types ]

[ Top ] [ m_ddk ] [ Types ]

NAME

  ddk_t

FUNCTION

  object containing ddk derivatives ([H,r] proportional to band velocities)

SOURCE

 72  type,public :: ddk_t
 73 
 74   !integer :: fh(3)
 75   ! file handler
 76   !  Fortran unit number if iomode==IO_MODE_FORTRAN
 77   !  MPI file handler if iomode==IO_MODE_MPI
 78 
 79   integer :: comm
 80   ! MPI communicator used for IO.
 81 
 82   !integer :: version
 83   ! File format version read from file.
 84 
 85   integer :: iomode=IO_MODE_FORTRAN
 86   ! Method used to access the DDK file:
 87   !   IO_MODE_FORTRAN for usual Fortran IO routines
 88   !   IO_MODE_MPI if MPI/IO routines.
 89   !   IO_MODE_ETSF netcdf files in etsf format
 90 
 91   integer :: rw_mode = DDK_NOMODE
 92    ! (Read|Write) mode
 93 
 94   !integer :: current_fpos
 95   ! The current position of the file pointer used for sequential access with Fortran-IO
 96   !  FPOS_EOF signals the end of file
 97 
 98   integer :: nsppol
 99    ! Number of spin polarizations.
100 
101   integer :: nspinor
102    ! Number of spinor components.
103 
104   integer :: nene
105     ! Number of energy points we may need the fsavg at
106 
107   integer :: nkfs
108     ! Number of k-points on the Fermi-surface (FS-BZ).
109 
110   integer :: maxnb
111    ! Max number of bands on the FS.
112    ! TODO: Maybe maxnbfs
113 
114   integer :: usepaw
115    ! 1 if PAW calculation, 0 otherwise
116 
117   integer :: prtvol=0
118    ! Verbosity level
119 
120   logical :: debug=.False.
121    ! Debug flag
122 
123   character(len=fnlen) :: paths(3) = ABI_NOFILE
124    ! File name
125 
126   real(dp) :: acell(3),rprim(3,3),gprim(3,3)
127    ! TODO: Are these really needed?
128 
129   real(dp), allocatable :: velocity (:,:,:,:)
130    ! (3,maxnb,nkfs,nsppol)
131    ! velocity on the FS in cartesian coordinates.
132 
133   real(dp), allocatable :: velocity_fsavg (:,:,:)
134    ! (3,nene,nsppol)
135    ! velocity on the FS in cartesian coordinates.
136 
137   logical :: use_ncddk(3)
138    ! True if we are readin DDK matrix elements from EVK.nc instead of WFK file
139 
140   type(crystal_t) :: cryst
141    ! Crystal structure read from file
142 
143  end type ddk_t
144 
145  public :: ddk_init              ! Initialize the object.
146  public :: ddk_read_fsvelocities ! Read FS velocities from file.
147  public :: ddk_fs_average_veloc  ! find FS average of velocity squared
148  public :: ddk_free              ! Close the file and release the memory allocated.
149  public :: ddk_print             ! output values
150  public :: eph_ddk               ! calculate eph ddk

m_ddk/eph_ddk [ Functions ]

[ Top ] [ m_ddk ] [ Functions ]

NAME

  eph_ddk

FUNCTION

  Calculate the DDK matrix elements using the commutator formulation.

INPUTS

PARENTS

      wfk_analyse

CHILDREN

SOURCE

258 subroutine eph_ddk(wfk_path,dtfil,dtset,&
259                    psps,pawtab,inclvkb,ngfftc,mpi_enreg,comm)
260 
261  use defs_basis
262  use defs_datatypes
263  use defs_abitypes
264  use m_profiling_abi
265  use m_xmpi
266  use m_errors
267  use m_wfk
268  use m_wfd
269 
270  use m_ebands,          only : ebands_ncwrite
271  use m_time,            only : cwtime, sec2str
272  use m_vkbr,            only : vkbr_t, nc_ihr_comm, vkbr_init, vkbr_free
273  use m_fstrings,        only : strcat, sjoin, itoa, ftoa, ktoa
274  use m_io_tools,        only : iomode_from_fname, get_unit
275  use m_cgtools,         only : dotprod_g
276  use m_fftcore,         only : get_kg, kpgsph, sphere
277  use m_crystal,         only : crystal_t
278  use m_crystal_io,      only : crystal_ncwrite
279  use m_pawtab,          only : pawtab_type
280 
281 !This section has been created automatically by the script Abilint (TD).
282 !Do not modify the following lines by hand.
283 #undef ABI_FUNC
284 #define ABI_FUNC 'eph_ddk'
285 !End of the abilint section
286 
287  implicit none
288 
289 !Arguments ------------------------------------
290 !scalars
291  character(len=*),intent(in) :: wfk_path
292  integer,intent(in) :: comm
293  type(datafiles_type),intent(in) :: dtfil
294  type(dataset_type),intent(in) :: dtset
295  type(wfk_t),target :: in_wfk
296  type(wfd_t),target :: in_wfd
297  type(vkbr_t) :: vkbr
298  type(ebands_t) :: ebands
299  type(crystal_t) :: cryst
300  type(hdr_type) :: hdr_tmp
301  type(pseudopotential_type),intent(in) :: psps
302  type(mpi_type),intent(inout) :: mpi_enreg
303  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
304 
305 !Local variables ------------------------------
306 !scalars
307  integer,parameter :: dummy_npw=0, formeig0=0, paral_kgb=0, master=0
308  logical,parameter :: force_istwfk1=.True.
309  integer :: iomode, mband, nbcalc, nsppol, ib_v, ib_c, inclvkb, dummy_gvec(3,dummy_npw)
310  integer :: mpw, spin, nspinor, nkpt, nband_k, npw_k
311  integer :: in_iomode, ii, ik, bandmin, bandmax, istwf_k
312  integer :: my_rank, nproc, ierr
313 #ifdef HAVE_NETCDF
314  integer :: ncerr,ncid
315 #endif
316 !arrays
317  integer,intent(in) :: ngfftc(18)
318  logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
319  integer,allocatable :: task_distrib(:,:,:,:)
320  integer,allocatable :: nband(:,:)
321  integer,allocatable :: kg_k(:,:)
322  character(len=500) :: msg
323  character(len=fnlen) :: fname
324  real(dp) :: cpu,wall,gflops,ecut
325  real(dp) :: kbz(3)
326  real(dp),allocatable :: dipoles(:,:,:,:,:,:)
327  complex(gwpc),allocatable :: ihrc(:,:)
328  complex(dp)               :: vg(3), vr(3)
329  complex(gwpc),allocatable :: ug_c(:),ug_v(:)
330 
331 !************************************************************************
332 
333  write(msg, '(2a)') "Computation of electron-photon coupling matrix elements (ddk)", ch10
334  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
335  call wrtout(std_out, msg, "COLL", do_flush=.True.)
336 
337  if (psps%usepaw == 1) then
338    MSG_ERROR("PAW not implemented")
339  end if
340 
341 #ifndef HAVE_NETCDF
342   MSG_ERROR("The matrix elements are only written in NETCDF format")
343 #endif
344 
345  ! paralelism
346  my_rank = xmpi_comm_rank(comm)
347  nproc = xmpi_comm_size(comm)
348 
349  ! Open input file, extract dimensions and allocate workspace arrays.
350  in_iomode = iomode_from_fname(wfk_path)
351  call wfk_open_read(in_wfk,wfk_path,formeig0,in_iomode,get_unit(),xmpi_comm_self)
352 
353  !read crystal
354  call crystal_from_hdr(cryst, in_wfk%hdr, 2)
355 
356  !read ebands
357  ebands = wfk_read_ebands(wfk_path,comm)
358 
359  mpw     = maxval(in_wfk%Hdr%npwarr)
360  nkpt    = in_wfk%nkpt
361  nsppol  = in_wfk%nsppol
362  nspinor = in_wfk%nspinor
363  mband   = in_wfk%mband
364  ecut    = in_wfk%hdr%ecut
365 
366  !TODO: hardcoded for now but should be an arugment
367  bandmin = 1
368  bandmax = mband
369  nbcalc  = bandmax-bandmin
370 
371  ABI_MALLOC(ug_c,    (mpw*nspinor))
372  ABI_MALLOC(ug_v,    (mpw*nspinor))
373  ABI_MALLOC(kg_k,    (3,mpw))
374  ABI_CALLOC(dipoles, (3,2,mband,mband,nkpt,nsppol))
375  ABI_MALLOC(ihrc,    (3, nspinor**2))
376 
377  ABI_MALLOC(nband,   (nkpt, nsppol))
378  ABI_MALLOC(keep_ur, (mband, nkpt, nsppol))
379  ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
380 
381  write(std_out,*) 'inclvkb: ', inclvkb
382  write(std_out,*) 'nkpoints:', nkpt
383  write(std_out,*) 'nbands:  ', mband
384  write(std_out,*) 'spin:    ', nsppol
385  write(std_out,*) 'spinor:  ', nspinor
386  write(std_out,*) 'ngfft:   ', in_wfk%hdr%ngfft
387  write(std_out,*) 'mpw:     ', mpw
388  write(std_out,*) 'ecut:    ', ecut
389 
390  !create distribution of the wavefunctions mask
391  keep_ur = .false.
392  bks_mask = .false.
393  nband = mband
394 
395  ! Distribute the k-points and bands over the processors
396  ABI_MALLOC(task_distrib,(bandmin:bandmax,bandmin:bandmax,nkpt,nsppol))
397  call xmpi_distab(nproc,task_distrib)
398 
399  ! create bks_mask to load the wavefunctions
400  do spin=1,nsppol ! Loop over spins
401    do ik=1,nkpt ! Loop over kpoints
402      do ib_v=bandmin,bandmax ! Loop over v bands
403        do ib_c=bandmin,bandmax ! Loop over c bands
404          if (task_distrib(ib_c,ib_v,ik,spin) == my_rank) then
405            bks_mask(ib_v,ik,spin) = .true.
406            bks_mask(ib_c,ik,spin) = .true.
407          end if
408        end do
409      end do
410    end do
411  end do
412 
413  !initialize distributed wavefunctions object
414  call wfd_init(in_wfd,cryst,pawtab,psps,keep_ur,paral_kgb,dummy_npw,mband,nband,nkpt,nsppol,&
415    bks_mask,dtset%nspden,nspinor,dtset%ecutsm,dtset%dilatmx,ebands%istwfk,ebands%kptns,&
416    ngfftc,dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut)
417 
418  ABI_FREE(bks_mask)
419  ABI_FREE(keep_ur)
420  ABI_FREE(nband)
421 
422  call wfd_print(in_wfd,header="Wavefunctions on the k-points grid",mode_paral='PERS')
423 
424  !Read Wavefunctions
425  iomode = iomode_from_fname(wfk_path)
426  call wfd_read_wfk(in_wfd,wfk_path,iomode)
427 
428 do spin=1,nsppol ! Loop over spins
429 
430    do ik=1,nkpt ! Loop over kpoints
431      ! Only do a subset a k-points
432      if (all(task_distrib(bandmin:bandmax,bandmin:bandmax,ik,spin) /= my_rank)) cycle
433      call cwtime(cpu,wall,gflops,"start")
434 
435      nband_k  = in_wfk%nband(ik,spin)
436      istwf_k  = in_wfk%hdr%istwfk(ik)
437      kbz      = in_wfk%hdr%kptns(:,ik)
438      npw_k    = in_wfk%hdr%npwarr(ik)
439 
440      ! Read WF
441      kg_k(:,1:npw_k) = in_wfd%kdata(ik)%kg_k
442 
443      ! Allocate KB form factors
444      if (inclvkb/=0) then ! Prepare term i <n,k|[Vnl,r]|n"k>
445        call vkbr_init(vkbr,cryst,psps,inclvkb,istwf_k,npw_k,kbz,kg_k)
446      end if
447 
448      ! Loop over bands
449      do ib_v=bandmin,bandmax
450        if (all(task_distrib(:,ib_v,ik,spin) /= my_rank)) cycle
451        ug_v(1:npw_k*nspinor) = in_wfd%wave(ib_v,ik,spin)%ug
452 
453        ! Loop over bands
454        do ib_c=ib_v,bandmax
455          if (task_distrib(ib_c,ib_v,ik,spin) /= my_rank) cycle
456          ug_c(1:npw_k*nspinor) = in_wfd%wave(ib_c,ik,spin)%ug
457 
458          ! Calculate matrix elements of i[H,r] for NC pseudopotentials.
459          ihrc = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,&
460                             kbz,ug_c,ug_v,kg_k)
461 
462          ! HM: 24/07/2018
463          ! Transform dipoles to be consistent with results from DFPT
464          ! Perturbations with DFPT are along the reciprocal lattice vectors
465          ! Perturbations with Commutator are along real space lattice vectors
466          ! dot(A, DFPT) = X
467          ! dot(B, COMM) = X
468          ! B = 2 pi (A^{-1})^T =>
469          ! dot(B^T B,COMM) = 2 pi DFPT
470          vr = (2*pi)*(2*pi)*sum(ihrc(:,:),2)
471          vg(1) = dot_product(Cryst%gmet(1,:),vr)
472          vg(2) = dot_product(Cryst%gmet(2,:),vr)
473          vg(3) = dot_product(Cryst%gmet(3,:),vr)
474 
475          ! Save matrix elements of i*r in the IBZ
476          dipoles(:,1,ib_c,ib_v,ik,spin) = real(vg)
477          dipoles(:,1,ib_v,ib_c,ik,spin) = real(vg) ! Hermitian conjugate
478          if (ib_v == ib_c) then
479             dipoles(:,2,ib_c,ib_v,ik,spin) = 0
480             dipoles(:,2,ib_v,ib_c,ik,spin) = 0
481          else
482             dipoles(:,2,ib_c,ib_v,ik,spin) =  aimag(vg)
483             dipoles(:,2,ib_v,ib_c,ik,spin) = -aimag(vg) ! Hermitian conjugate
484          end if
485        end do
486      end do
487 
488      ! Free KB form factors
489      call vkbr_free(vkbr)
490 
491      ! loop over k-points
492      call cwtime(cpu,wall,gflops,"stop")
493      write(msg,'(2(a,i0),2(a,f8.2))')"k-point [",ik,"/",nkpt,"] completed. cpu:",cpu,", wall:",wall
494      call wrtout(std_out, msg, do_flush=.True.)
495 
496    end do ! loop over k-points
497  end do ! loop over spin
498 
499  ABI_FREE(ug_c)
500  ABI_FREE(ug_v)
501  ABI_FREE(kg_k)
502  ABI_FREE(ihrc)
503 
504  ! Gather the k-points computed by all processes
505  call xmpi_sum_master(dipoles,master,comm,ierr)
506 
507  !write the matrix elements
508 #ifdef HAVE_NETCDF
509    ! Output DDK file in netcdf format.
510    if (my_rank == master) then
511 
512      ! Have to build hdr on k-grid with info about perturbation.
513      call hdr_copy(in_wfk%hdr, hdr_tmp)
514      hdr_tmp%qptn = [0,0,0]
515 
516      do ii=1,3
517          fname = strcat(dtfil%filnam_ds(4), '_', itoa(ii), "_EVK.nc")
518          NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
519          hdr_tmp%pertcase = (cryst%natom*3)+ii
520          NCF_CHECK(hdr_ncwrite(hdr_tmp, ncid, 43, nc_define=.True.))
521          NCF_CHECK(crystal_ncwrite(cryst, ncid))
522          NCF_CHECK(ebands_ncwrite(ebands, ncid))
523          ncerr = nctk_def_arrays(ncid, [ &
524            nctkarr_t('h1_matrix_elements', "dp", &
525             "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
526          NCF_CHECK(ncerr)
527          NCF_CHECK(nctk_set_datamode(ncid))
528          ncerr = nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), dipoles(ii,:,:,:,:,:) )
529          NCF_CHECK(ncerr)
530          NCF_CHECK(nf90_close(ncid))
531      end do
532      call hdr_free(hdr_tmp)
533 
534    end if
535 #endif
536 
537  ABI_FREE(dipoles)
538 
539 end subroutine eph_ddk