1!===============================================================================
2!
3! Routines:
4!
5! (1) input()   Originally By ?         Last Modified 10/5/2009 (gsm)
6!
7! Sets up various data structures. Reads in and distributes wavefunctions.
8! Stores wavefunctions in structures distributed in memory.
9!
10!===============================================================================
11
12#include "f_defs.h"
13
14module input_m
15
16  use global_m
17  use checkbz_m
18  use createpools_m
19  use eqpcor_m
20  use fftw_m
21  use fullbz_m
22  use input_utils_m
23  use misc_m
24  use read_rho_vxc_m
25  use scissors_m
26  use sort_m
27  use wfn_rho_vxc_io_m
28  use io_utils_m
29  use epsread_hdf5_m
30#ifdef HDF5
31  use hdf5
32#endif
33  use timing_m, only: timing => sigma_timing
34
35  implicit none
36
37  private
38  public :: input
39
40contains
41
42subroutine input(crys,gvec,syms,kp,wpg,sig,wfnk,iunit_c,iunit_k,fnc,fnk,wfnkqmpi,wfnkmpi, wfnk_phonq, wfnk_phonq_mpi, ep_read_in)
43  type (crystal), intent(out) :: crys
44  type (gspace), intent(out) :: gvec
45  type (symmetry), intent(out) :: syms
46  type (kpoints), intent(out) :: kp
47  type (wpgen), intent(out) :: wpg
48  type (siginfo), intent(inout) :: sig ! ZL modifies
49  type (wfnkstates), intent(out) :: wfnk, wfnk_phonq ! ZL adds wfnk_phonq
50  integer, intent(in) :: iunit_c, iunit_k
51  character*20, intent(in) :: fnc, fnk
52  type (wfnkqmpiinfo), intent(out) :: wfnkqmpi
53  type (wfnkmpiinfo), intent(out) :: wfnkmpi, wfnk_phonq_mpi ! ZL adds wfnk_phonq_mpi
54  ! ZL: EP temporary variables
55  type (crystal) :: crys_tmp
56  type (gspace) :: gvec_tmp
57  type (symmetry) :: syms_tmp
58  type (kpoints) :: kp_tmp
59
60
61  character :: fncor*32, tmpfn*16
62  character :: tmpstr*100,tmpstr1*16,tmpstr2*16
63  character :: filenameh5*80
64  integer :: i,ierr,ii,j,jj,k,kk
65  integer :: ikn, irk, nq
66  integer :: ncoul,ncoulb,ncouls,nmtx,nmtx_l,nmtx_col
67  integer, allocatable :: isrt(:), isrti(:)
68  integer :: ipe
69  integer :: Nrod,Nplane,Nfft(3),dNfft(3),dkmax(3),nmpinode
70  integer :: ntband,npools,ndiag_max,noffdiag_max,ntband_max
71  real(DP) :: scalarsize,qk(3),vcell,qtot,omega_plasma,rhog0
72  real(DP) :: mem,fmem,rmem,rmem2,scale,dscale
73  integer, allocatable :: gvecktmp(:,:)
74  logical :: do_ch_sum, skip_checkbz
75  type(grid) :: gr
76  real(DP) :: del
77  integer :: g0(3), i_looper ! ZL: to map EP sig%indk_phonq_g0
78
79  character(len=3) :: sheader
80  integer :: iflavor, freq_dep, error
81
82  integer :: ng_eps, nq_eps, nmtxmax_eps, qgrid_eps(3), nfreq_eps, nfreq_imag_eps
83  real(DP) :: ecuts
84
85  logical :: file_exists
86
87  ! ZL: electron phonon (EP) control
88  logical, intent(in), optional :: ep_read_in
89  logical :: ep_read
90  logical :: do_phonq ! ZL: if do_phonq = .false., wfnk_phonq and wfnk_phonq_mpi will not
91                     ! be used, nor allocated or evaluated
92
93  PUSH_SUB(input)
94  ! ZL: set up EP
95  do_phonq = .false.
96  ep_read = .false.
97
98!---------------
99! Read sigma.inp
100!  write(6,*) 'before inread'
101  if (.not.ep_read) then
102    call inread(sig) ! ZL: in inread(), if nphonq .eq. 1, program exits.
103                     ! ZL: so after inread(), simply use sig%elph
104  else
105    if (peinf%inode.eq.0) write(6,*) 'Skip inread sigma for dWFN.'
106  endif
107!  write(6,*) 'after inread'
108  ! ZL: conditions for generating wfnk_phonq, wfnk_phonq_mpi
109
110!------------------------
111! Initialize HDF5
112#ifdef HDF5
113  if(sig%use_hdf5 .or. sig%wfn_hdf5) call h5open_f(error)
114#endif
115
116  if(sig%nkn .lt. 1) then
117    if(peinf%inode.eq.0) write(0,*) 'sig%nkn < 1'
118    call die('sig%nkn')
119  endif
120
121!---------------------------------
122! Determine the available memory
123
124  call procmem(mem,nmpinode)
125
126  if(peinf%inode.eq.0) then
127    write(6,998) mem/1024.0d0**2
128  endif
129998 format(1x,'Memory available: ',f0.1,' MB per PE')
130
131  fmem=mem/8.0d0
132
133!-------------------------------
134! (gsm) Estimate the required memory
135
136  if(peinf%inode.eq.0) then
137
138! determine number of frequencies and number of q-points
139    if (sig%freq_dep.eq.-1) then
140      nq=sig%nq
141      sig%nFreq=1 ! ZL: for EP, this is fine reevaluate
142    endif
143
144    if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1.or.sig%freq_dep.eq.2.or.sig%freq_dep.eq.3) then
145
146#ifdef HDF5
147      if(sig%use_hdf5) then
148        filenameh5 = 'eps0mat.h5'
149        call read_eps_grid_sizes_hdf5(ng_eps, nq_eps, ecuts, nfreq_eps, &
150          nfreq_imag_eps, nmtxmax_eps, qgrid_eps, freq_dep, filenameh5)
151        write(6,*)
152        write(6,'(1x,a)') 'Read from eps0mat.h5:'
153        write(6,'(1x,a,i0)') '- Number of G-vectors in the charge-density grid: ', ng_eps
154        write(6,'(1x,a,i0)') '- Max. number of G-vectors for dielectric matrices: ', nmtxmax_eps
155        write(6,'(1x,a,f0.3)') '- Cutoff of the dielectric matrix (Ry): ', ecuts
156        write(6,'(1x,a,i0)') '- Number of imaginary frequencies: ', nfreq_imag_eps
157        write(6,'(1x,a,i0)') '- Total number of frequencies: ', nfreq_eps
158        write(6,'(1x,a,i0)') '- Frequency dependence: ', freq_dep
159        write(6,'(1x,a,i0)') '- Number of q-points: ', nq_eps
160        write(6,'(1x,a,3(i0,1x))') '- Q grid: ', qgrid_eps
161        write(6,*)
162        sig%nFreq = nfreq_eps
163        sig%nfreq_imag = nfreq_imag_eps
164        nq = nq_eps
165
166        INQUIRE(FILE="epsmat.h5", EXIST=file_exists)
167        if (file_exists) then
168          filenameh5 = 'epsmat.h5'
169          call read_eps_grid_sizes_hdf5(ng_eps, nq_eps, ecuts, nfreq_eps, nfreq_imag_eps, &
170            nmtxmax_eps, qgrid_eps, freq_dep, filenameh5)
171          write(6,'(1x,a)') 'Read from epsmat.h5:'
172          write(6,'(1x,a,i0)') '- Number of G-vectors in the charge-density grid: ', ng_eps
173          write(6,'(1x,a,i0)') '- Max. number of G-vectors for dielectric matrices: ', nmtxmax_eps
174          write(6,'(1x,a,f0.3)') '- Cutoff of the dielectric matrix: ', ecuts
175          write(6,'(1x,a,i0)') '- Number of imaginary frequencies: ', nfreq_imag_eps
176          write(6,'(1x,a,i0)') '- Total number of frequencies: ', nfreq_eps
177          write(6,'(1x,a,i0)') '- Frequency dependence: ', freq_dep
178          write(6,'(1x,a,i0)') '- Number of q-points: ', nq_eps
179          write(6,'(1x,a,3(i0,1x))') '- Q grid: ', qgrid_eps
180          write(6,*)
181          nq = nq + nq_eps
182        endif
183      else
184#endif
185
186        call open_file(10,file='eps0mat',form='unformatted',status='old')
187        read(10)
188        read(10) i, sig%nFreq
189        read(10)
190        read(10)
191        read(10)
192        read(10)
193        read(10) ecuts
194        read(10) nq
195        call close_file(10)
196        call open_file(11,file='epsmat',form='unformatted',status='old',iostat=ierr)
197        if (ierr.eq.0) then
198          read(11)
199          read(11)
200          read(11)
201          read(11)
202          read(11)
203          read(11)
204          read(11)
205          read(11) nq_eps
206          nq = nq + nq_eps
207          call close_file(11)
208        endif
209
210#ifdef HDF5
211    endif
212#endif
213
214    endif
215
216  endif
217
218! FHJ: Read header of WFN file (and WFNq, if requred), perform consistency
219! checks, and figure out the number of bands and ncrit.
220  sheader = 'WFN'
221  iflavor = 0
222  if(sig%wfn_hdf5) then
223  else
224    if(peinf%inode == 0) then
225      if(.not.ep_read) then
226        write(6,'(a)') ' Reading header of WFN_inner'
227        call open_file(25,file='WFN_inner',form='unformatted',status='old')
228      else
229        ! ZL: open dWFN
230        write(6,'(a)') ' Reading header of dWFN'
231        call open_file(25,file='dWFN',form='unformatted',status='old')
232      endif
233    endif
234    call read_binary_header_type(25, sheader, iflavor, kp, gvec, syms, crys)
235  endif
236  call check_trunc_kpts(sig%icutv, kp) ! This does not affect EP
237
238  if (sig%ecutb<TOL_ZERO) then
239    sig%ecutb = kp%ecutwfc
240  endif
241  if (sig%freq_dep/=-1 .and. sig%ecuts<TOL_ZERO) then
242#ifdef MPI
243    call MPI_Bcast(ecuts, 1, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr)
244#endif
245    sig%ecuts = ecuts
246  endif
247
248  if (sig%ntband==0 .and. (.not.ep_read)) sig%ntband = kp%mnband - 1
249
250! DVF: Here we re-define the number of bands in the dynamical self energy
251! summation to exclude the deep core states. This ensures that the pools will be
252! set up properly. This achieves most of what is needed to exclude core states.
253! There is a little more work in the i/o routines to make sure you`re getting the
254! higher valence states.
255! ZL: we disable this operation if doing elph to avoid possible mismaps of band indices
256  if ((sig%ncore_excl.ne.0) .and. sig%elph) then
257    call die('Electron-Phonon calculation is not consistent with sig%ncore_excl.ne.0 .')
258  endif
259  sig%ntband=sig%ntband-sig%ncore_excl
260
261  SAFE_ALLOCATE(kp%elda, (sig%ntband,kp%nrk,kp%nspin))
262  kp%el(:,:,:) = kp%el(:,:,:) - sig%avgpot / ryd
263! DVF : add sig%ncore_excl to make sure we are getting the right energies. kp%el is setup in
264! wfn_rho_vxc_io.f90 routine (see .p.f for something comprehensible) and includes core states.
265  kp%elda(1:sig%ntband, 1:kp%nrk, 1:kp%nspin)=ryd*kp%el(1+sig%ncore_excl:sig%ntband+sig%ncore_excl, 1:kp%nrk, 1:kp%nspin)
266  call scissors_shift(kp, sig%scis, sig%spl_tck)
267  ! If QP corrections requested, read the corrected QP energies from file (in eV)
268  if(sig%eqp_corrections .and. (.not.ep_read)) then
269    fncor='eqp.dat'
270! DVF : add sig%ncore_excl to make sure we are getting the right energies. The eqp.dat file
271! will have band indices referenced to the total number of states, including the core states.
272    call eqpcor(fncor,peinf%inode,peinf%npes,kp,1+sig%ncore_excl,sig%ntband+sig%ncore_excl,0,0,kp%el,kp%el,kp%el,1,0)
273    ! note: want in Ry since conversion occurs below
274  endif
275  if (peinf%inode==0 .and. (.not.ep_read)) then
276    call find_efermi(sig%rfermi, sig%efermi, sig%efermi_input, kp, sig%ntband+sig%ncore_excl, 1+sig%ncore_excl, &
277      "unshifted grid", should_search=.true., should_update=.true., write7=.false.)
278  endif
279
280  if (sig%nvband==0 .and. (.not.ep_read)) then
281    sig%nvband = minval(kp%ifmax)
282    sig%ncrit = maxval(kp%ifmax) - minval(kp%ifmax)
283  endif
284
285! DVF: Here we re-define the number of bands in the bare exchange
286! summation to exclude the deep core states.
287! ZL: this has been disabled already for EP
288  sig%nvband=sig%nvband-sig%ncore_excl
289
290#ifdef MPI
291    call MPI_Bcast(sig%efermi, 1, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr)
292    call MPI_Bcast(sig%nvband, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
293    call MPI_Bcast(sig%ncrit, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
294#endif
295
296! DVF: again, this is to exclude the core states.
297  if ((sig%freq_dep==0 .and. sig%exact_ch==1) .or. sig%freq_dep==-1) then
298    if (sig%ntband + sig%ncore_excl > max(maxval(sig%diag(1:sig%ndiag)), sig%nvband + sig%ncore_excl + sig%ncrit)) then
299      ! Do not waste time and memory with the unoccupied bands which will not be used
300      sig%ntband = max(maxval(sig%diag(1:sig%ndiag))-sig%ncore_excl, sig%nvband + sig%ncrit)
301      if (peinf%inode==0 .and. sig%freq_dep==0)  then
302        write(0,'(a)') "WARNING: Resetting number_bands to number of bands actually needed for exact static CH."
303      else if (peinf%inode==0) then
304        write(0,'(a)') "WARNING: Resetting number_bands to number of bands actually needed for exchange."
305      endif
306    endif
307  endif
308
309  if(peinf%inode == 0) then
310    if(.not.ep_read) then
311      call logit('Reading WFN_inner -- gvec info')
312    else
313      call logit('Reading dWFN -- gvec info')
314    endif
315  endif
316  SAFE_ALLOCATE(gvec%components, (3, gvec%ng))
317  if(sig%wfn_hdf5) then
318  else
319      call read_binary_gvectors(25, gvec%ng, gvec%ng, gvec%components)
320  endif
321
322!-----------------------
323! sort G-vectors with respect to their kinetic energy
324  SAFE_ALLOCATE(isrt, (gvec%ng))
325  SAFE_ALLOCATE(isrti, (gvec%ng))
326  SAFE_ALLOCATE(gvec%ekin, (gvec%ng))
327  call kinetic_energies(gvec, crys%bdot, gvec%ekin)
328  call sortrx(gvec%ng, gvec%ekin, isrt, gvec = gvec%components)
329  ncouls = gcutoff(gvec%ng, gvec%ekin, isrt, sig%ecuts)
330  ncoulb = gcutoff(gvec%ng, gvec%ekin, isrt, sig%ecutb)
331  SAFE_DEALLOCATE_P(gvec%ekin)
332
333  if (peinf%inode==0) then
334    write(6,*)
335    write(6,'(a)') ' Calculation parameters:'
336    write(6,'(1x,a,f0.2)') '- Cutoff of the bare Coulomb interaction (Ry): ', sig%ecutb
337    write(6,'(1x,a,f0.2)') '- Cutoff of the screened Coulomb interaction (Ry): ', sig%ecuts
338    write(6,'(1x,a,i0)') '- Number of G-vectors up to the bare int. cutoff: ', ncoulb
339    write(6,'(1x,a,i0)') '- Number of G-vectors up to the screened int. cutoff: ', ncouls
340    write(6,'(1x,a,i0)') '- Total number of bands in the calculation: ', sig%ntband
341    write(6,'(1x,a,i0)') '- Number of fully occupied valence bands: ', sig%nvband
342    write(6,'(1x,a,i0)') '- Number of partially occ. conduction bands: ', sig%ncrit
343    write(6,*)
344  endif
345
346  if(sig%ecuts > sig%ecutb) then
347    call die("The screened_coulomb_cutoff cannot be greater than the bare_coulomb_cutoff.", &
348      only_root_writes = .true.)
349  endif
350  if(sig%ntband > kp%mnband) then
351    if(.not.ep_read) then
352      call die("number_bands is larger than are available in WFN_inner", only_root_writes = .true.)
353    else
354      call die("number_bands is larger than are available in dWFN", only_root_writes = .true.)
355    endif
356  endif
357  do ii = 1, sig%ndiag
358    if (sig%diag(ii)>sig%ntband) then
359      write(0,'(a,3(i0,1x))') 'ERROR: the band_index_max cannot be greater than number_bands: ', &
360        ii, sig%diag(ii), sig%ntband
361      call die('band_index_max cannot be greater than number_bands', only_root_writes = .true.)
362    endif
363  enddo
364
365! FHJ: do vanilla WFN reading stuff
366
367  SAFE_ALLOCATE(gvecktmp, (3,gvec%ng))
368  gvecktmp(:,:)=gvec%components(:,:)
369  do i=1,gvec%ng
370    gvec%components(:,i)=gvecktmp(:,isrt(i))
371    isrti(isrt(i)) = i
372  enddo
373  SAFE_DEALLOCATE(gvecktmp)
374
375  call gvec_index(gvec)
376
377!------------------------------------------------------------------------------
378! Pools, distribution and memory estimation
379!------------------------------------------------------------------------------
380
381  if(peinf%inode == 0) then
382    ncoul=max(ncouls,ncoulb)
383    nmtx=ncouls
384
385! divide bands over processors (this is repeated below)
386    ntband=min(sig%ntband,kp%mnband)
387    if (peinf%npools .le. 0 .or. peinf%npools .gt. peinf%npes) then
388      call createpools(sig%ndiag,ntband,peinf%npes,npools,ndiag_max,ntband_max)
389    else
390      npools = peinf%npools
391      if (mod(sig%ndiag,npools).eq.0) then
392        ndiag_max=sig%ndiag/npools
393      else
394        ndiag_max=sig%ndiag/npools+1
395      endif
396      if (mod(ntband,peinf%npes/npools).eq.0) then
397        ntband_max=ntband/(peinf%npes/npools)
398      else
399        ntband_max=ntband/(peinf%npes/npools)+1
400      endif
401    endif
402    if (sig%noffdiag.gt.0) then
403      if (mod(sig%noffdiag,npools).eq.0) then
404        noffdiag_max=sig%noffdiag/npools
405      else
406        noffdiag_max=sig%noffdiag/npools+1
407      endif
408    endif
409
410    nmtx_l=int(sqrt(dble(nmtx)**2/dble(peinf%npes/npools)))
411    nmtx_col=int(dble(nmtx)/dble(peinf%npes/npools))+1
412
413    scalarsize = sizeof_scalar()
414
415! required memory
416    rmem=0.0d0
417! arrays eps and epstemp in program main
418    if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1) then
419      rmem=rmem+(dble(nmtx_l)**2+dble(nmtx))*scalarsize
420      rmem=rmem+(dble(nmtx_col)*nmtx)*scalarsize
421    endif
422! wtilde_array, imagloop_ig, imaginary_flag in mtxel_cor
423    if (sig%freq_dep.eq.1) then
424      rmem=rmem+(dble(nmtx_col)*nmtx)*16
425      rmem=rmem+(dble(nmtx_col)*nmtx)*4
426      rmem=rmem+(dble(nmtx_col)*nmtx)*1
427    endif
428! arrays epsR, epsRtemp, epsA, and epsAtemp in program main
429    if (sig%freq_dep.eq.2 .or.sig%freq_dep.eq.3) then
430      rmem=rmem+(dble(nmtx_l)**2+dble(nmtx)) &
431        *dble(sig%nFreq)*2.0d0*scalarsize
432    endif
433! array epsmpi%eps in subroutine epscopy
434    if (sig%freq_dep.eq.0.or.sig%freq_dep.eq.1) then
435      rmem=rmem+(dble(nmtx_col)*nmtx)*dble(nq)*scalarsize
436    endif
437! arrays epsmpi%epsR and epsmpi%epsA in subroutine epscopy
438    if (sig%freq_dep.eq.2 .or.sig%freq_dep.eq.3) then
439      rmem=rmem+(nmtx_col*nmtx)*dble(nq) &
440        *dble(sig%nFreq)*2.0d0*scalarsize
441    endif
442! array aqs in program main
443    rmem=rmem+dble(ntband_max)*dble(ncoul)*scalarsize
444! array aqsaug in program main
445    if (sig%noffdiag.gt.0) then
446      rmem=rmem+dble(ntband_max)*dble(ncoul) &
447        *dble(sig%ndiag)*dble(sig%nspin)*scalarsize
448    endif
449! array aqsch in program main
450    if (sig%exact_ch.eq.1) then
451      rmem=rmem+dble(ncoul)*scalarsize
452    endif
453! arrays aqsaugchd and aqsaugcho in program main
454    if (sig%exact_ch.eq.1.and.nq.gt.1) then
455      rmem=rmem+dble(ncoul)*dble(ndiag_max) &
456        *dble(sig%nspin)*scalarsize
457      if (sig%noffdiag.gt.0) then
458        rmem=rmem+dble(ncoul)*dble(noffdiag_max)* &
459          dble(sig%nspin)*scalarsize
460      endif
461    endif
462! array wfnk%zk in subroutine input
463    rmem=rmem+dble(ndiag_max)*dble(kp%ngkmax)* &
464      dble(sig%nspin*kp%nspinor)*scalarsize
465! array wfnkq%zkq in subroutine genwf
466    rmem=rmem+dble(ntband_max)*dble(kp%ngkmax) &
467      *dble(sig%nspin*kp%nspinor)*scalarsize
468! arrays zc in subroutine input and zin in subroutine genwf
469    rmem=rmem+dble(kp%ngkmax)*dble(kp%nspin*kp%nspinor)*scalarsize
470! array wfnkoff%zk in program main and subroutines mtxel_vxc and mtxel_cor
471    if ((sig%exact_ch.eq.1.and.sig%noffdiag.gt.0).or. &
472      (.not.sig%use_vxcdat .and. .not.sig%sigma_correction .and. .not. sig%is_EXX)) then
473      rmem=rmem+2.0d0*dble(kp%ngkmax)*scalarsize
474    endif
475! array gvec%index_vec in input
476    rmem=rmem+dble(gvec%nFFTgridpts)*4.0d0
477! arrays fftbox1 and fftbox2 in subroutines mtxel and mtxel_occ
478    call setup_FFT_sizes(gvec%FFTgrid,Nfft,scale)
479    rmem=rmem+dble(Nfft(1)*Nfft(2)*Nfft(3))*32.0d0
480! arrays wfnkqmpi%isort and wfnkqmpi%cg in subroutine input
481    rmem=rmem+dble(kp%ngkmax)*dble(kp%nrk)*4.0d0+ &
482      dble(kp%ngkmax)*dble(ntband_max)*dble(sig%nspin*kp%nspinor)* &
483      dble(kp%nrk)*dble(scalarsize)
484! arrays wfnkmpi%isort and wfnkmpi%cg in subroutine input
485    rmem=rmem+dble(kp%ngkmax)*dble(sig%nkn)*4.0d0+ &
486      dble((ndiag_max*kp%ngkmax)/(peinf%npes/npools))* &
487      dble(sig%nspin*kp%nspinor)*dble(sig%nkn)* &
488      dble(scalarsize)
489
490989 format(1x,'Memory required for execution: ',f0.1,' MB per PE')
491    write(6,989) rmem/1024.0d0**2
492
493! random numbers
494    rmem=0.0D0
495    if (sig%icutv/=TRUNC_BOX) then
496! arrays ran, qran, and qran2
497! (ran is deallocated before qran2 is allocated)
498      rmem=rmem+6.0D0*dble(nmc)*8.0D0
499    endif
500! various truncation schemes
501    rmem2=0.0d0
502! cell wire truncation
503    if (sig%icutv==TRUNC_WIRE) then
504      dkmax(1) = gvec%FFTgrid(1) * n_in_wire
505      dkmax(2) = gvec%FFTgrid(2) * n_in_wire
506      dkmax(3) = 1
507      call setup_FFT_sizes(dkmax,dNfft,dscale)
508! array fftbox_2D
509      rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*16.0d0
510! array inv_indx
511      rmem2=rmem2+dble(Nfft(1))*dble(Nfft(2))*dble(Nfft(3))* &
512        4.0d0
513! array qran
514      rmem2=rmem2+3.0D0*dble(nmc)*8.0D0
515    endif
516! cell box truncation (parallel version only)
517    if (sig%icutv==TRUNC_BOX) then
518      dkmax(1) = gvec%FFTgrid(1) * n_in_box
519      dkmax(2) = gvec%FFTgrid(2) * n_in_box
520      dkmax(3) = gvec%FFTgrid(3) * n_in_box
521      call setup_FFT_sizes(dkmax,dNfft,dscale)
522      if (mod(dNfft(3),peinf%npes) == 0) then
523        Nplane = dNfft(3)/peinf%npes
524      else
525        Nplane = dNfft(3)/peinf%npes+1
526      endif
527      if (mod(dNfft(1)*dNfft(2),peinf%npes) == 0) then
528        Nrod = (dNfft(1)*dNfft(2))/peinf%npes
529      else
530        Nrod = (dNfft(1)*dNfft(2))/peinf%npes+1
531      endif
532! array fftbox_2D
533      rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*dble(Nplane)* &
534        16.0d0
535! array fftbox_1D
536      rmem2=rmem2+dble(dNfft(3))*dble(Nrod)*16.0d0
537! array dummy
538!          rmem2=rmem2+dble(dNfft(1))*dble(dNfft(2))*16.0d0
539! arrays dummy1 and dummy2
540      rmem2=rmem2+dble(Nrod)*dble(peinf%npes+1)*16.0d0
541! array inv_indx
542      rmem2=rmem2+dble(Nfft(1))*dble(Nfft(2))*dble(Nfft(3))* &
543        4.0d0
544    endif
545    if (rmem2 .gt. rmem) rmem = rmem2
546988 format(1x,'Memory required for vcoul: ',f0.1,' MB per PE')
547    write(6,988) rmem/1024.0d0**2
548    write(6,*)
549  endif
550
551  if(peinf%inode == 0) then
552!----------------------------------------------
553! Compute cell volume from wave number metric
554    call get_volume(vcell,crys%bdot)
555    if (abs(crys%celvol-vcell).gt.TOL_Small) then
556      call die('volume mismatch')
557    endif
558
559! (gsm) check consistency of spin indices
560
561    do k=1,sig%nspin
562      if (sig%spin_index(k).lt.1.or.sig%spin_index(k).gt.kp%nspin) &
563        call die('inconsistent spin indices')
564    enddo
565
566    if(sig%ntband.gt.kp%mnband) then
567      write(tmpstr1,660) sig%ntband
568      write(tmpstr2,660) kp%mnband
569      write(0,666) TRUNC(tmpstr1), TRUNC(tmpstr2)
570      if(.not.ep_read) then
571        call die('More bands specified in sigma.inp than available in WFN_inner.')
572      else
573        call die('More bands specified in sigma.inp than available in dWFN.')
574      endif
575660   format(i16)
576666   format(1x,'The total number of bands (',a,') specified in sigma.inp',/, &
577        3x,'is larger than the number of bands (',a,') available in WFN_inner.',/)
578    endif
579    do_ch_sum = .not.((sig%freq_dep == -1) .or. (sig%freq_dep == 0 .and. sig%exact_ch == 1))
580    if(sig%ntband.eq.kp%mnband) then
581      if(.not.ep_read) then
582        call die("You must provide one more band in WFN_inner than used in sigma.inp number_bands in order to assess degeneracy.")
583      else
584        call die("You must provide one more band in dWFN than used in sigma.inp number_bands in order to assess degeneracy.")
585      endif
586    endif
587
588!--------------------------------------------------------------------
589! SIB:  Find the k-points in sig%kpt in the list kp%rk (die if not found).
590! sig%indkn holds the index of the k-points in sig%kpt in kp%rk, i.e.
591! kp%rk(sig%indkn(ikn))=sig%kpt(ikn)
592
593    SAFE_ALLOCATE(sig%indkn, (sig%nkn))
594    do ikn=1,sig%nkn
595      sig%indkn(ikn)=0
596      qk(:)=sig%kpt(:,ikn)
597      do irk=1,kp%nrk
598        if(all(abs(kp%rk(1:3,irk)-qk(1:3)) .lt. TOL_Small)) sig%indkn(ikn)=irk
599      enddo
600
601      ! ZL: print information
602      if(do_phonq .and. peinf%inode.eq.0) write(6,*) "sig%indkn(",ikn,") =", sig%indkn(ikn)
603
604      if(sig%indkn(ikn) .eq. 0) then
605        if(.not.ep_read) then
606          write(0,'(a,3f10.5,a)') 'Could not find the k-point ', (qk(i),i=1,3), ' among those read from WFN_inner :'
607        else
608          write(0,'(a,3f10.5,a)') 'Could not find the k-point ', (qk(i),i=1,3), ' among those read from dWFN :'
609        endif
610        write(0,'(3f10.5)') ((kp%rk(i,irk),i=1,3),irk=1,kp%nrk)
611        call die('k-point in sigma.inp k_points block not available.')
612      endif
613    enddo
614
615    ! ZL: EP, when reading regular WFN_inner, NOT dWFN, we map k+phonq points
616    if(do_phonq .and. (.not. ep_read) .and. (sig%nphonq .eq. 1)) then  ! ZL: only nphonq=1 supported for now
617      SAFE_ALLOCATE(sig%indk_phonq, (sig%nkn * sig%nphonq))
618      SAFE_ALLOCATE(sig%indk_phonq_g0, (3, sig%nkn * sig%nphonq))
619      do ikn=1,sig%nkn*sig%nphonq
620        sig%indk_phonq(ikn)=0
621        qk(:)=sig%k_phonq(:,ikn)
622        irk_loop_g0: do irk=1,kp%nrk
623          do i_looper=1,3
624            ! ZL: this idea is taken from genwf_mpi.f90
625            !     the most interesting part is that
626            !     it uses the fact that g0(3) is INTEGER
627            !     therefore it will automatically drop the fractional part
628            !     and keep only the integer part
629            !     This is also guaranteed by adding TOL_small so that 0.999999
630            !     is indeed 1
631            !     Then we just need to compare if del is close to integer
632            ! NOTE: gmap uses a different idea, it maps between old and new wfn
633            !       but here we map to master G-list
634            del = qk(i_looper) - kp%rk(i_looper,irk)
635            if (del .ge. 0.0d0) g0(i_looper) = del + TOL_Small
636            if (del .lt. 0.0d0) g0(i_looper) = del - TOL_Small
637            if (abs(del-g0(i_looper)) .gt. TOL_Small) cycle irk_loop_g0
638          enddo
639          sig%indk_phonq(ikn)=irk
640          sig%indk_phonq_g0(:,ikn)=g0(:)
641          exit irk_loop_g0
642        enddo irk_loop_g0
643
644        ! ZL: print information
645        if(peinf%inode .eq. 0) then
646          write(6,*) "sig%indk_phonq(",ikn,") =", sig%indk_phonq(ikn)
647          write(6,*) "sig%indk_phonq_g0(",ikn,") =", sig%indk_phonq_g0(:,ikn)
648        endif
649
650        if(sig%indk_phonq(ikn) .eq. 0) then
651          write(0,'(a,3f10.5,a)') 'Could not find the k_phonq-point ', (qk(i),i=1,3), ' among those in WFN_inner :'
652          write(0,'(3f10.5)') ((kp%rk(i,irk),i=1,3),irk=1,kp%nrk)
653          call die('k_phonq-point from sigma.inp k_points+phonq not available.')
654        endif
655      enddo
656    endif ! EP k+phonq mapping
657
658
659    if(do_ch_sum .and. (.not.ep_read)) then
660      if(any(abs(kp%el(sig%ntband+sig%ncore_excl, 1:kp%nrk, 1:kp%nspin) - &
661                 kp%el(sig%ntband+sig%ncore_excl+1, 1:kp%nrk, 1:kp%nspin)) .lt. sig%tol/RYD)) then
662        if(sig%degeneracy_check_override) then
663          write(0,'(a)') &
664            "WARNING: Selected number of bands for CH sum (number_bands) breaks degenerate subspace. " // &
665            "Run degeneracy_check.x for allowable numbers."
666          write(0,*)
667        else
668          write(0,'(a)') &
669            "Run degeneracy_check.x for allowable numbers, or use keyword " // &
670            "degeneracy_check_override to run anyway (at your peril!)."
671          call die("Selected number of bands for CH sum (number_bands) breaks degenerate subspace.")
672        endif
673      endif
674    endif
675  endif
676
677!----------------------------------------------------------------
678! Check consistency
679
680  if (kp%mnband < maxval(sig%diag) .and. peinf%inode == 0) then
681    write(0,*) 'The highest requested band is ', maxval(sig%diag),' but WFN_inner contains only ', kp%mnband,' bands.'
682    call die('too few bands')
683  endif
684
685  ! qgridsym has no meaning when there is only one q-point (Gamma)
686  ! setting this to false means warnings about its applicability will not be triggered
687  if(all(abs(kp%rk(1:3,1:kp%nrk)) < TOL_Zero)) sig%qgridsym = .false.
688
689  if(peinf%inode == 0) then
690    ! ZL: assess_degeneracies will allocate kp%degeneracy
691    ! kp%el is still in Ry here, but sig%tol is in eV.
692    call assess_degeneracies(kp, kp%el(sig%ntband+sig%ncore_excl+1, :, :), sig%ntband, sig%efermi, sig%tol/RYD, sig = sig, &
693                             ncore_excl=sig%ncore_excl)
694    if(.not.ep_read) then
695      call calc_qtot(kp, crys%celvol, sig%efermi, qtot, omega_plasma, write7=.false.)
696    endif
697  endif
698
699  if ( peinf%inode == 0 ) call timing%start(timing%fullbz)
700  gr%nr = kp%nrk
701  SAFE_ALLOCATE(gr%r, (3, gr%nr))
702  gr%r = kp%rk
703  call fullbz(crys,syms,gr,syms%ntran,skip_checkbz,wigner_seitz=.false.,paranoid=.true.)
704  if(.not.ep_read) then
705    tmpfn='WFN_inner'
706  else
707    tmpfn='dWFN'
708  endif
709  if (.not. skip_checkbz) then
710    call checkbz(gr%nf,gr%f,kp%kgrid,kp%shift,crys%bdot, &
711      tmpfn,'k',.false.,sig%freplacebz,sig%fwritebz)
712  endif
713  call dealloc_grid(gr)
714  if ( peinf%inode == 0 ) call timing%stop(timing%fullbz)
715
716  ! For discussion of how q-symmetry may (and may not) be used with degenerate states,
717  ! see Hybertsen and Louie, Phys. Rev. B 35, 5585 (1987) Appendix A
718  if(sig%qgridsym .and. sig%noffdiag > 0) then
719    if(peinf%inode == 0) then
720      write(0,'(a)') "WARNING: Cannot calculate offdiagonal elements unless no_symmetries_q_grid is set."
721      write(0,'(a)') "This flag is being reset to enable the calculation."
722    endif
723    sig%qgridsym = .false.
724  endif
725
726  ! ZL: for EP, turn off q-symmetry. Interestingly, this coincides with off-diag
727  if(sig%qgridsym .and. ep_read) then
728    if(peinf%inode == 0) then
729      write(0,'(a)') "WARNING: Cannot calculate electron phonon unless no_symmetries_q_grid is set."
730      write(0,'(a)') "This flag is being reset to enable the calculation."
731    endif
732    sig%qgridsym = .false.
733  endif
734
735  if(peinf%inode == 0) then
736    if(sig%qgridsym) then
737      write(6,'(1x,a)') 'Q-grid symmetries are being used.'
738    else
739      write(6,'(1x,a)') 'Q-grid symmetries are not being used.'
740    endif
741  endif
742
743  kp%el(:,:,:)=kp%el(:,:,:)*ryd
744
745 !----------------------------------------------------------------
746 ! Distribute data
747
748#ifdef MPI
749  if(sig%qgridsym .or. sig%noffdiag > 0 .or. sig%ncrit > 0) then
750    if(peinf%inode.ne.0) then
751      SAFE_ALLOCATE(kp%degeneracy, (sig%ntband, kp%nrk, kp%nspin))
752    endif
753    call MPI_Bcast(kp%degeneracy(1,1,1), sig%ntband * kp%nrk * kp%nspin, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
754  endif
755
756  if(peinf%inode.ne.0) then
757    SAFE_ALLOCATE(sig%indkn, (sig%nkn))
758    if(do_phonq) then
759      SAFE_ALLOCATE(sig%indk_phonq, (sig%nkn*sig%nphonq))
760      SAFE_ALLOCATE(sig%indk_phonq_g0, (3, sig%nkn*sig%nphonq))
761    endif
762  endif
763  call MPI_Bcast(sig%indkn, sig%nkn, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
764  if(do_phonq) then
765    call MPI_Bcast(sig%indk_phonq, sig%nkn*sig%nphonq, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
766    call MPI_Bcast(sig%indk_phonq_g0, 3*sig%nkn*sig%nphonq, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
767  endif
768#endif
769
770!----------------------------------------------------------------
771! CHP: Set the init_freq_eval relative to the Fermi energy in case
772!      of full frequency. This part should NOT be put before the
773!      Fermi energy is set.
774
775  if ((sig%freq_dep==2 .or. (sig%freq_dep==1.and.sig%fdf==-3)) .and. sig%freq_grid_shift==1) then
776    ! FHJ: only do this for FE-based freq. grids
777    sig%freqevalmin = sig%freqevalmin + sig%efermi
778  endif
779
780!--------------------------------------------------------------------
781! Read in the exchange-correlation potential and store in array vxc
782
783  if(.not.sig%use_vxcdat .and. .not.sig%sigma_correction .and. .not. sig%is_EXX) then
784
785    call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 1)
786
787  endif ! not using vxc.dat
788
789  if(.not. sig%use_vxc2dat) then
790    ! This is for hybrid functional like calculation (one shot)
791    if (sig%freq_dep.eq.-1 .and. ((1.0d0 - sig%coulomb_mod%long_range_frac_fock > TOL_SMALL) .or. &
792                  (1.0d0 - sig%coulomb_mod%short_range_frac_fock > TOL_SMALL)))  &
793      call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 2)
794
795  endif ! not using vxc2.dat
796
797  ! ZL: read dVXC when reading dWFN
798  if (sig%elph .and. ep_read) then
799    call read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, 3)
800  endif
801
802!--------------------------------------------------------------------
803! Read in the charge density and store in array rho (formerly known as CD95)
804! CD95 Ref: http://www.nature.com/nature/journal/v471/n7337/full/nature09897.html
805
806  if(sig%freq_dep.eq.1) then
807    call timing%start(timing%input_read)
808
809    call read_rho(wpg, gvec, kp, syms, crys, isrti, isrt, 'WFN_inner') ! ZL: only read from WFN_inner, NOT dWFN, for RHO
810
811    rhog0 = sum(wpg%nelec(1:kp%nspin))
812    if(peinf%inode == 0 .and. abs(rhog0 - qtot) > TOL_Small) then
813      write(0,'(a,g12.6,a,g12.6)') 'WARNING: RHO has total charge ', rhog0, ' but WFN_inner occupations give charge ', qtot
814    endif
815
816    call timing%stop(timing%input_read)
817  endif ! sig%freq_dep
818
819!------------------------------------------------------------------------
820! Divide bands over processors
821!
822! sig%ntband           number of total (valence and conduction) bands
823! peinf%npools         number of parallel sigma calculations
824! peinf%ndiag_max      maximum number of diag sigma calculations
825! peinf%noffdiag_max   maximum number of offdiag sigma calculations
826! peinf%ntband_max     maximum number of total bands per node
827!
828! ntband_node          number of total bands per current node
829! nvband_node          number of valence bands per current node
830! peinf%indext(itb)    indices of total bands belonging to current node
831!
832! peinf%index_diag     band index for diag sigma calculation
833! peinf%flag_diag      flag for storing diag sigma calculation
834! peinf%index_offdiag  band index for offdiag sigma calculation
835! peinf%flag_offdiag   flag for storing offdiag sigma calculation
836!
837! flags are needed in case of uneven distribution, nodes still do the
838! calculation because they share epsilon and wavefunctions and need to
839! participate in global communications, but the result is not stored
840
841  if (peinf%npools .le. 0 .or. peinf%npools .gt. peinf%npes) then
842    call createpools(sig%ndiag,sig%ntband,peinf%npes,npools,ndiag_max,ntband_max)
843    peinf%npools = npools
844    peinf%ndiag_max = ndiag_max
845    peinf%ntband_max = ntband_max
846  else
847    if (mod(sig%ndiag,peinf%npools).eq.0) then
848      peinf%ndiag_max=sig%ndiag/peinf%npools
849    else
850      peinf%ndiag_max=sig%ndiag/peinf%npools+1
851    endif
852
853    if (mod(sig%ntband,peinf%npes/peinf%npools).eq.0) then
854      peinf%ntband_max=sig%ntband/(peinf%npes/peinf%npools)
855    else
856      peinf%ntband_max=sig%ntband/(peinf%npes/peinf%npools)+1
857    endif
858  endif
859
860! JRD Create Separate Comms for Each Pool
861
862  peinf%npes_pool = peinf%npes/peinf%npools
863  peinf%my_pool=peinf%inode/peinf%npes_pool
864  peinf%pool_rank = mod(peinf%inode,peinf%npes_pool)
865
866#ifdef MPI
867!  call MPI_Comm_Group(MPI_COMM_WORLD,orig_group,mpierr)
868!
869!  SAFE_ALLOCATE(ranks,(peinf%npes))
870!  do ii=1,peinf%npes
871!    ranks(ii)=ii-1
872!  enddo
873!
874!  if (peinf%my_pool .ne. -1) then
875!    call MPI_Group_incl(orig_group, peinf%npes_pool, &
876!      ranks(peinf%my_pool*peinf%npes_pool+1:(peinf%my_pool+1)*peinf%npes_pool), peinf%pool_group,mpierr)
877!  else
878!! We need a group to send comm_create below. If you are not in a legit pool, we just use the first group
879!    call MPI_Group_incl(orig_group, peinf%npes_pool, &
880!      ranks(0*peinf%npes_pool+1:(0+1)*peinf%npes_pool), peinf%pool_group,mpierr)
881!  endif
882!
883!! All procs in original comm need to do this
884!
885!  call MPI_Comm_create(MPI_COMM_WORLD,peinf%pool_group,peinf%pool_comm,mpierr)
886
887  call MPI_Comm_split(MPI_COMM_WORLD, peinf%my_pool, peinf%pool_rank, peinf%pool_comm, mpierr)
888
889  if (peinf%my_pool .gt. (peinf%npools-1)) then
890    peinf%my_pool = -1
891  endif
892
893!  !write(6,*) peinf%inode,"Created groups comms"
894
895!  SAFE_DEALLOCATE(ranks)
896#endif
897
898  if (sig%noffdiag.gt.0) then
899    if (mod(sig%noffdiag,peinf%npools).eq.0) then
900      peinf%noffdiag_max=sig%noffdiag/peinf%npools
901    else
902      peinf%noffdiag_max=sig%noffdiag/peinf%npools+1
903    endif
904  endif
905
906  SAFE_ALLOCATE(peinf%index_diag, (peinf%ndiag_max))
907  SAFE_ALLOCATE(peinf%flag_diag, (peinf%ndiag_max))
908  peinf%index_diag=1
909  peinf%flag_diag=.false.
910  do ii=1,peinf%ndiag_max*peinf%npools
911    jj=(ii-1)/peinf%npools+1 ! which number in the pool it is
912    kk=mod(ii-1,peinf%npools) ! which pool it is in
913    if (peinf%inode/(peinf%npes/peinf%npools).eq.kk) then
914      if (ii.le.sig%ndiag) then
915        peinf%index_diag(jj)=ii
916        peinf%flag_diag(jj)=.true.
917      else
918        peinf%index_diag(jj)=1
919        peinf%flag_diag(jj)=.false.
920      endif
921    endif
922  enddo
923
924  if (sig%noffdiag.gt.0) then
925    SAFE_ALLOCATE(peinf%index_offdiag, (peinf%noffdiag_max))
926    SAFE_ALLOCATE(peinf%flag_offdiag, (peinf%noffdiag_max))
927    peinf%index_offdiag=1
928    peinf%flag_offdiag=.false.
929    do ii=1,peinf%noffdiag_max*peinf%npools
930      jj=(ii-1)/peinf%npools+1
931      kk=mod(ii-1,peinf%npools)
932      if (peinf%inode/(peinf%npes/peinf%npools).eq.kk) then
933        if (ii.le.sig%noffdiag) then
934          peinf%index_offdiag(jj)=ii
935          peinf%flag_offdiag(jj)=.true.
936        else
937          peinf%index_offdiag(jj)=1
938          peinf%flag_offdiag(jj)=.false.
939        endif
940      endif
941    enddo
942  endif
943
944  SAFE_ALLOCATE(peinf%indext, (peinf%ntband_max))
945
946!JRD Now is over procs per pool
947
948  SAFE_ALLOCATE(peinf%indext_dist, (peinf%ntband_max,peinf%npes_pool))
949  SAFE_ALLOCATE(peinf%ntband_dist, (peinf%npes_pool))
950
951  peinf%ntband_node=0
952  peinf%ntband_dist=0
953  peinf%nvband_node=0
954  peinf%indext=0
955  peinf%indext_dist=0
956
957  do ii=1,sig%ntband
958      ! FHJ: The bands are assigned to the processors in a sequential fashion.
959      ! This is more efficient for the parallel IO. Eg: for 2  bands/processor,
960      ! v1->PE0, v2->PE0, v3->PE1, etc.
961      ipe = (ii-1)/peinf%ntband_max
962      if (peinf%pool_rank.eq.ipe) then
963        peinf%ntband_node=peinf%ntband_node+1
964        if(ii.le.(sig%nvband+sig%ncrit)) then
965          peinf%nvband_node=peinf%nvband_node+1
966        endif
967        peinf%indext(peinf%ntband_node)=ii
968      endif
969      peinf%ntband_dist(ipe+1)=peinf%ntband_dist(ipe+1)+1
970      peinf%indext_dist(peinf%ntband_dist(ipe+1),ipe+1)=ii
971  enddo
972#ifdef MPI
973  if (peinf%verb_debug) then
974    call logit('Begin distrib report')
975    FLUSH(0)
976    FLUSH(6)
977    call MPI_Barrier(MPI_COMM_WORLD, mpierr)
978    if (peinf%inode==0) then
979      do ipe=1,peinf%npes_pool
980        write(6,'(/,a,i0,2x,a,i0,a)') '### pool_ipe=',ipe,' ; num bands=',peinf%ntband_dist(ipe),' ###'
981        do ii=1,peinf%ntband_dist(ipe)
982          write(6,'(a,i0,a,i0)') 'local:',ii,' -> global:', peinf%indext_dist(ii,ipe)
983        enddo
984      enddo
985      write(6,*)
986    endif
987    FLUSH(6)
988    call MPI_Barrier(MPI_COMM_WORLD, mpierr)
989    call logit('End distrib report')
990  endif
991#endif
992
993!------------------------------------------------------------------------
994! Report distribution of bands over processors
995
996  if (peinf%inode==0) then
997    write(6,*)
998    write(6,'(1x,a)') 'Parallelization report:'
999701 format(1x,'- Using ',i0,' processor(s), ' ,i0,' pool(s), ',i0,' processor(s) per pool.')
1000    write(6,701) peinf%npes, peinf%npools, &
1001      peinf%npes/peinf%npools
1002702 format(1x,' - Note: distribution is not ideal; ',i0,' processor(s) is/are idle.')
1003    if (mod(peinf%npes,peinf%npools).ne.0) &
1004      write(6,702) peinf%npes-(peinf%npes/peinf%npools)* &
1005      peinf%npools
1006    if (mod(sig%ndiag,peinf%npools).eq.0) then
1007703 format(1x,'- Each pool is computing ',i0,' diagonal sigma matrix element(s).')
1008      write(6,703) peinf%ndiag_max
1009    else
1010704 format(1x,'- Each pool is computing ',i0,' to ',i0,' diagonal sigma matrix element(s).')
1011      write(6,704) peinf%ndiag_max-1, peinf%ndiag_max
1012705 format(1x,'- Note: distribution is not ideal because the number of diagonal sigma',/,&
1013      3x,'matrix elements (',i0,') is not divisible by the number of pools (',i0,').')
1014      write(6,705) sig%ndiag, peinf%npools
1015    endif
1016    if (mod(sig%noffdiag,peinf%npools).eq.0) then
1017706 format(1x,'- Each pool is computing ',i0,' off-diagonal sigma matrix element(s).')
1018      write(6,706) peinf%noffdiag_max
1019    else
1020707 format(1x,'- Each pool is computing ',i0,' to ',i0,' off-diagonal sigma matrix element(s).')
1021      write(6,707) peinf%noffdiag_max-1, peinf%noffdiag_max
1022708 format(1x,'- Note: distribution is not ideal because the number of off-diagonal sigma',/,&
1023      3x,'matrix elements (',i0,') is not divisible by the number of pools (',i0,').')
1024      write(6,708) sig%noffdiag, peinf%npools
1025    endif
1026    if (mod(sig%ntband,peinf%npes/peinf%npools).eq.0) then
1027709 format(1x,'- Each processor is holding ',i0, ' band(s).')
1028      write(6,709) peinf%ntband_max
1029    else
1030710 format(1x,'- Each pool is holding ',i0,' to ',i0,' band(s).')
1031      write(6,710) minval(peinf%ntband_dist), maxval(peinf%ntband_dist)
1032711 format(1x,'- Note: distribution is not ideal because the total number of bands',/,&
1033      3x,'(',i0,') is not divisible by the number of processors per pool (',i0,').')
1034      write(6,711) sig%ntband, peinf%npes/peinf%npools
1035    endif
1036    write(6,*)
1037  endif
1038
1039!-----------------------------------------------------------
1040!
1041!     LOOP OVER K-POINT GRID AND READ IN WAVEFUNCTIONS
1042!
1043!-----------------------------------------------------------
1044
1045  if(sig%wfn_hdf5) then
1046  else
1047    ! ZL: we call read_wavefunctions() only once in input()
1048    if(.not.ep_read) then
1049      if(.not.do_phonq) then
1050        call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi)
1051      else
1052        call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi,&
1053                                ep_read_in=.false., do_phonq_in=do_phonq, wfnk_phonq=wfnk_phonq, wfnk_phonq_mpi=wfnk_phonq_mpi)
1054      endif
1055    else
1056      call read_wavefunctions(kp, gvec, sig, wfnk, iunit_c, iunit_k, fnc, fnk, wfnkqmpi, wfnkmpi, ep_read_in = ep_read)
1057    endif
1058    if(peinf%inode.eq.0) then
1059      call close_file(25)
1060    endif
1061  endif
1062
1063  if(peinf%inode.eq.0) then
1064! Write out information about self-energy calculation
1065    call scissors_write(6, sig%scis)
1066    call scissors_write(6, sig%scis_outer, "outer")
1067    write(6,'(a)')
1068
1069    if(sig%freq_dep.eq.1) then
1070      do k=1,sig%nspin
1071        write(6,160) sig%spin_index(k),wpg%nelec &
1072          (sig%spin_index(k)),sqrt(wpg%wpsq(sig%spin_index(k)))
1073160     format(/,1x,'Data for sum rule:',/,&
1074          1x,'- rho(0,',i0,') = ',f0.6,' electrons',/,&
1075          1x,'- wp = ',f0.6,' eV',/)
1076      enddo
1077    endif ! sig%freq_dep
1078    write(6,'(1x,a,i0)') 'Number of bands to compute diagonal self-energy matrix elements: ', sig%ndiag
1079    write(6,'(1x,a)') 'Bands:'
1080    write(6,'(1x,"- ",i0)') sig%diag(:)
1081    write(6,'(1x,a,i0)') 'Number of off-diagonal bands to compute self-energy matrix elements: ', sig%noffdiag
1082    if (sig%noffdiag>0) then
1083      write(6,'(1x,a)') 'Off-diagonal bands:'
1084      do k=1,sig%noffdiag
1085174 format(1x,'- offmap(:, ',i6,') =',3i6)
1086        write(6,174) k, (sig%offmap(k,ii), ii = 1, 3)
1087      enddo
1088    endif
1089    write(6,*)
1090
1091  endif ! node 0
1092
1093  SAFE_DEALLOCATE(isrt)
1094  SAFE_DEALLOCATE(isrti)
1095  SAFE_DEALLOCATE_P(kp%w)
1096  SAFE_DEALLOCATE_P(kp%el)
1097  SAFE_DEALLOCATE_P(kp%elda)
1098
1099  POP_SUB(input)
1100
1101  return
1102
1103end subroutine input
1104
1105
1106subroutine read_wavefunctions(kp,gvec,sig,wfnk,iunit_c,iunit_k,fnc,fnk,wfnkqmpi,wfnkmpi, &
1107                              ep_read_in,do_phonq_in,wfnk_phonq,wfnk_phonq_mpi)
1108  type (kpoints), intent(in) :: kp
1109  type (gspace), intent(in) :: gvec
1110  type (siginfo), intent(in) :: sig
1111  type (wfnkstates), intent(out) :: wfnk
1112  integer, intent(in) :: iunit_c, iunit_k
1113  character*20, intent(in) :: fnc, fnk
1114  type (wfnkqmpiinfo), intent(out) :: wfnkqmpi
1115  type (wfnkmpiinfo), intent(out) :: wfnkmpi
1116  ! ZL: for EP
1117  type (wfnkstates), optional, intent(out) :: wfnk_phonq  ! similar to wfnk
1118  type (wfnkmpiinfo), optional, intent(out) :: wfnk_phonq_mpi  ! similar to wfnkmpi
1119
1120  character :: tmpstr*100,tmpstr1*16,tmpstr2*16
1121  integer :: i,ii,j,k
1122  integer :: ikn, irk
1123  integer :: istore,nstore,inum,g1,g2,iknstore, istore_ep, iknstore_ep ! ZL adds istore_ep, iknstore_ep
1124  integer :: ndvmax, ndvmax_l
1125  integer, allocatable :: isort(:)
1126  real(DP) :: qk(3)
1127  SCALAR, allocatable :: zc(:,:)
1128  logical :: dont_read
1129  type(gspace) :: gvec_kpt
1130
1131  type(progress_info) :: prog_info !< a user-friendly progress report
1132
1133  ! ZL: electron phonon (EP)
1134  logical, intent(in), optional :: ep_read_in, do_phonq_in
1135  logical :: ep_read, do_phonq
1136  integer :: induse ! index in use, for indkn or indk_phonq
1137  integer :: tmpngk
1138
1139  PUSH_SUB(read_wavefunctions)
1140
1141  ! ZL: set up EP
1142  ep_read = .false.
1143  if(present(ep_read_in)) then
1144    ep_read = ep_read_in
1145  endif
1146
1147  do_phonq = .false.
1148  if(present(do_phonq_in)) then
1149    do_phonq = do_phonq_in
1150  endif
1151
1152  ! ZL: although it should never happen, in case someone mis-uses, double check here
1153  if(do_phonq .and. ep_read) then
1154    call die("Doing k+phonq points (designed for outer) is NOT allowed for dWFN (only used as inner).", &
1155              only_root_writes = .true.)
1156  endif
1157
1158  if (do_phonq) then
1159    if ((.not.present(wfnk_phonq)) .or. (.not.present(wfnk_phonq_mpi))) then
1160      call die("Doing k_phonq wfn needs wfklkrq and wfnk_phonq_mpi as arguments!", only_root_writes = .true.)
1161    endif
1162  else
1163    if (present(wfnk_phonq) .or. present(wfnk_phonq_mpi)) then
1164      call die("You are not doing k_phonq wfn, then do not pass wfnk_phonq or wfnk_phonq_mpi.", only_root_writes=.true.)
1165    endif
1166  endif
1167  ! ZL: done EP setup
1168
1169  SAFE_ALLOCATE(isort, (gvec%ng))
1170
1171  SAFE_ALLOCATE(wfnkqmpi%nkptotal, (kp%nrk))
1172  SAFE_ALLOCATE(wfnkqmpi%isort, (kp%ngkmax,kp%nrk))
1173  SAFE_ALLOCATE(wfnkqmpi%band_index, (peinf%ntband_max,kp%nrk))
1174  SAFE_ALLOCATE(wfnkqmpi%qk, (3,kp%nrk))
1175  SAFE_ALLOCATE(wfnkqmpi%el, (sig%ntband,sig%nspin,kp%nrk))
1176  ! ZL: ntband_max is total number of bands distributed to each processor
1177  !     wfnkqmpi%cg is distributed over bands
1178  SAFE_ALLOCATE(wfnkqmpi%cg, (kp%ngkmax,peinf%ntband_max,sig%nspin*kp%nspinor,kp%nrk))
1179
1180  if (sig%nkn.gt.1) then  ! ZL: same logic for EP
1181    ! it is using ngkmax, so same for EP, wfnkmpi%cg is distributed over ndv
1182    ndvmax=peinf%ndiag_max*kp%ngkmax
1183    if (mod(ndvmax,peinf%npes/peinf%npools).eq.0) then
1184      ndvmax_l=ndvmax/(peinf%npes/peinf%npools)
1185    else
1186      ndvmax_l=ndvmax/(peinf%npes/peinf%npools)+1
1187    endif
1188    SAFE_ALLOCATE(wfnkmpi%nkptotal, (sig%nkn))
1189    SAFE_ALLOCATE(wfnkmpi%isort, (kp%ngkmax,sig%nkn))
1190    SAFE_ALLOCATE(wfnkmpi%qk, (3,sig%nkn))
1191    SAFE_ALLOCATE(wfnkmpi%el, (sig%ntband,sig%nspin,sig%nkn))
1192    SAFE_ALLOCATE(wfnkmpi%elda, (sig%ntband,sig%nspin,sig%nkn))
1193    SAFE_ALLOCATE(wfnkmpi%cg, (ndvmax_l,sig%nspin*kp%nspinor,sig%nkn))
1194
1195    ! ZL: allocate wfnk_phonq_mpi
1196    !     if sig%nkn=1, then no mpi, only use wfnk_phonq which stores k+phonq
1197    if (do_phonq) then
1198      ! ZL: the size should in principle be sig%nkn*sig%nphonq,
1199      !     since we consider only nphonq=1, we omit it here
1200      SAFE_ALLOCATE(wfnk_phonq_mpi%nkptotal, (sig%nkn))
1201      SAFE_ALLOCATE(wfnk_phonq_mpi%isort, (kp%ngkmax,sig%nkn))
1202      SAFE_ALLOCATE(wfnk_phonq_mpi%qk, (3,sig%nkn))
1203      SAFE_ALLOCATE(wfnk_phonq_mpi%el, (sig%ntband,sig%nspin,sig%nkn))
1204      SAFE_ALLOCATE(wfnk_phonq_mpi%elda, (sig%ntband,sig%nspin,sig%nkn))
1205      SAFE_ALLOCATE(wfnk_phonq_mpi%cg, (ndvmax_l,sig%nspin*kp%nspinor,sig%nkn))
1206    endif
1207  endif
1208
1209!-----------------------------------
1210! Read in wavefunction information
1211
1212  SAFE_ALLOCATE(wfnk%isrtk, (gvec%ng))
1213  SAFE_ALLOCATE(wfnk%ek, (sig%ntband,sig%nspin))
1214  SAFE_ALLOCATE(wfnk%elda, (sig%ntband,sig%nspin))
1215
1216  ! ZL: allocate wfnk_phonq
1217  if (do_phonq) then
1218    SAFE_ALLOCATE(wfnk_phonq%isrtk, (gvec%ng))
1219    SAFE_ALLOCATE(wfnk_phonq%ek, (sig%ntband,sig%nspin))
1220    SAFE_ALLOCATE(wfnk_phonq%elda, (sig%ntband,sig%nspin))
1221  endif
1222
1223  if(.not.ep_read) then
1224    call progress_init(prog_info, 'reading wavefunctions (WFN_inner)', 'state', kp%nrk*sig%ntband)
1225  else
1226    call progress_init(prog_info, 'reading wavefunctions (dWFN)', 'state', kp%nrk*sig%ntband)
1227  endif
1228
1229  ! ZL: outmost loop, irk over kp%nrk
1230  do irk=1,kp%nrk
1231    if(.not.ep_read) then
1232      write(tmpstr,*) 'Reading WFN_inner -> cond/val wfns irk=',irk
1233    else
1234      write(tmpstr,*) 'Reading dWFN -> cond/val wfns irk=',irk
1235    endif
1236    call logit(tmpstr)
1237    qk(:)=kp%rk(:,irk)
1238
1239!----------------------------
1240! Read in and sort gvectors
1241    SAFE_ALLOCATE(gvec_kpt%components, (3, kp%ngk(irk)))
1242    if ( peinf%inode == 0 ) call timing%start(timing%input_read)
1243    call read_binary_gvectors(25, kp%ngk(irk), kp%ngk(irk), gvec_kpt%components)
1244    if ( peinf%inode == 0 ) call timing%stop(timing%input_read)
1245
1246    do i = 1, kp%ngk(irk)
1247      call findvector(isort(i), gvec_kpt%components(:, i), gvec)
1248      if (isort(i) .eq. 0) then
1249        write(0,*) 'could not find gvec', kp%ngk(irk), i, gvec_kpt%components(1:3, i)
1250        call die('findvector')
1251      endif
1252    enddo
1253    SAFE_DEALLOCATE_P(gvec_kpt%components)
1254
1255!--------------------------------------------------------
1256! Determine if Sigma must be computed for this k-point.
1257! If so, store the bands and wavefunctions on file iunit_k.
1258! If there is only one k-point, store directly in wfnk.
1259
1260    istore=0
1261    if(do_phonq) istore_ep=0 ! ZL: initialize
1262    do ikn=1,sig%nkn ! ZL: here we take advantage that we implement for one phonq point for now
1263                     !     otherwize, need another outer loop: do ikn=1,sig%nkn*sig%nphonq
1264
1265      if(sig%indkn(ikn).eq.irk) then  ! ZL: KEY if-statement: here it judges whether store wfnkmpi
1266        istore=1  ! ZL: for a given irk point, only one ind(ikn) index can match
1267        iknstore=ikn ! ZL: and this ind(ikn) is recorded as iknstore
1268
1269        wfnk%nkpt=kp%ngk(irk)
1270        wfnk%ndv=peinf%ndiag_max*kp%ngk(irk)
1271        wfnk%k(:)=qk(:)
1272        wfnk%isrtk(:)=isort(:)
1273        do k=1,sig%nspin
1274          wfnk%ek(1:sig%ntband,k)= &
1275            kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k))
1276          wfnk%elda(1:sig%ntband,k)= &
1277            kp%elda(1:sig%ntband,irk,sig%spin_index(k))
1278        enddo
1279        SAFE_ALLOCATE(wfnk%zk, (wfnk%ndv,sig%nspin*kp%nspinor))
1280        wfnk%zk=ZERO
1281      endif
1282
1283      ! ZL: for do_phonq
1284      if(do_phonq) then
1285        if (sig%indk_phonq(ikn).eq.irk) then
1286          istore_ep=1
1287          iknstore_ep=ikn
1288
1289          wfnk_phonq%nkpt=kp%ngk(irk)
1290          wfnk_phonq%ndv=peinf%ndiag_max*kp%ngk(irk)
1291          wfnk_phonq%k(:)=qk(:) + sig%indk_phonq_g0(:,ikn)
1292          wfnk_phonq%isrtk(:)=isort(:)
1293
1294          call shift_phonq_g0(sig%indk_phonq_g0(:,ikn), wfnk_phonq, gvec)
1295
1296          do k=1,sig%nspin
1297            wfnk_phonq%ek(1:sig%ntband,k)= &
1298              kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k))
1299            wfnk_phonq%elda(1:sig%ntband,k)= &
1300              kp%elda(1:sig%ntband,irk,sig%spin_index(k))
1301          enddo
1302          SAFE_ALLOCATE(wfnk_phonq%zk, (wfnk_phonq%ndv,sig%nspin*kp%nspinor))
1303          wfnk_phonq%zk=ZERO
1304        endif
1305      endif
1306
1307    enddo ! ikn
1308
1309    ! ZL: nkptotal is ngk
1310    wfnkqmpi%nkptotal(irk) = kp%ngk(irk)
1311    wfnkqmpi%isort(1:kp%ngk(irk),irk) = isort(1:kp%ngk(irk))
1312    do k = 1, sig%nspin
1313      wfnkqmpi%el(1:sig%ntband,k,irk) = &
1314        kp%el(1+ sig%ncore_excl:sig%ntband+ sig%ncore_excl,irk,sig%spin_index(k))
1315    enddo
1316    wfnkqmpi%qk(1:3,irk) = qk(1:3)
1317
1318!-------------------------------------------------------------------------
1319! SIB:  Read wave functions from file WFN_inner (ZL: or dWFN) (unit=25) and have
1320! the appropriate processor write it to iunit_c after checking norm.
1321! The wavefunctions for bands where Sigma matrix elements are
1322! requested are stored in wfnk%zk for later writing to unit iunit_k.
1323! If band index is greater than sig%ntband, we will actually
1324! not do anything with this band (see code below).
1325!  We still have to read it though, in order
1326! to advance the file to get to the data for the next k-point.
1327
1328    SAFE_ALLOCATE(zc, (kp%ngk(irk), kp%nspin*kp%nspinor))
1329
1330    inum=0
1331    do i=1,kp%mnband
1332
1333      if ( peinf%inode == 0 ) call timing%start(timing%input_read)
1334
1335! DVF: don`t read deep core states.
1336      dont_read = (i > sig%ntband+sig%ncore_excl .or. i <= sig%ncore_excl)
1337
1338      call read_binary_data(25, kp%ngk(irk), kp%ngk(irk), kp%nspin*kp%nspinor, zc, dont_read = dont_read)
1339
1340      if ( peinf%inode == 0 ) call timing%stop(timing%input_read)
1341
1342      if (.not.dont_read) then
1343        call progress_step(prog_info, sig%ntband*(irk-1) + i)
1344        if(peinf%inode.eq.0) then
1345          do k=1,sig%nspin
1346            if(.not.ep_read) then
1347              call checknorm('WFN_inner',i,irk,kp%ngk(irk),k,kp%nspinor,zc(:,:))
1348            endif
1349          enddo
1350        endif
1351
1352        nstore=0  ! ZL: for bands
1353        do j=1,peinf%ndiag_max
1354          if (.not.peinf%flag_diag(j)) cycle
1355          if (i==sig%diag(peinf%index_diag(j))) nstore=j
1356        enddo
1357
1358        ! ZL: if this is the kpoint in outer, save in wfnk
1359        if((istore.eq.1).and.(nstore.ne.0)) then
1360          do k=1,sig%nspin*kp%nspinor
1361            do j=1,kp%ngk(irk)
1362              wfnk%zk((nstore-1)*kp%ngk(irk)+j,k) = zc(j,sig%spin_index(k))
1363            enddo
1364          enddo
1365        endif
1366
1367        ! ZL: for do_phonq
1368        if(do_phonq.and.(istore_ep.eq.1).and.(nstore.ne.0)) then
1369          do k=1,sig%nspin*kp%nspinor
1370            do j=1,kp%ngk(irk)
1371              wfnk_phonq%zk((nstore-1)*kp%ngk(irk)+j,k) = zc(j,sig%spin_index(k))
1372            enddo
1373          enddo
1374        endif
1375
1376! DVF : i is indexed including the deep core states (since it runs over all the bands
1377! in the wavefunction), while the indexing arrays
1378! are not. We have to subtract sig%ncore_excl from i to get the right states.
1379        j=0  ! ZL: index to decide if save this band
1380        do ii=1,peinf%ntband_node
1381          if (peinf%indext(ii)==i-sig%ncore_excl) j=1
1382        enddo
1383
1384        if ( peinf%inode == 0 ) call timing%start(timing%input_write)
1385
1386        if (j.eq.1) then
1387          inum=inum+1
1388          wfnkqmpi%band_index(inum,irk)=i-sig%ncore_excl
1389          do k=1,sig%nspin*kp%nspinor
1390            wfnkqmpi%cg(1:kp%ngk(irk),inum,k,irk)= &
1391              zc(1:kp%ngk(irk),sig%spin_index(k))
1392          enddo
1393        endif
1394
1395        if ( peinf%inode == 0 ) call timing%stop(timing%input_write)
1396      else
1397        ! FHJ: the following lines were introduced in r6294 and are supposed to
1398        ! be a shortcut if we are past the last band of the last k-point. However,
1399        ! in light of a previous bug (#223), this feature is commented out for now.
1400        !! FHJ: shortcut if this is past the last band of the last k-point
1401        !if (irk==kp%nrk) exit
1402      endif
1403
1404    enddo ! i (loop over bands) ZL: kp%mnband
1405    SAFE_DEALLOCATE(zc)
1406
1407    ! ZL: this is the needed kpoint for wfnkmpi, copy from wfnk to wfnkmpi
1408    if((istore.eq.1).and.(sig%nkn.gt.1)) then  ! ZL: here it stores wfnkmpi
1409
1410      if ( peinf%inode == 0 ) call timing%start(timing%input_write)
1411
1412      ikn=iknstore
1413      wfnkmpi%nkptotal(ikn)=kp%ngk(irk)
1414      wfnkmpi%isort(1:kp%ngk(irk),ikn)=wfnk%isrtk(1:kp%ngk(irk))
1415      wfnkmpi%qk(1:3,ikn)=qk(1:3)
1416      wfnkmpi%el(1:sig%ntband,1:sig%nspin,ikn)= &
1417        wfnk%ek(1:sig%ntband,1:sig%nspin)
1418      wfnkmpi%elda(1:sig%ntband,1:sig%nspin,ikn)= &
1419        wfnk%elda(1:sig%ntband,1:sig%nspin)
1420      do k=1,sig%nspin*kp%nspinor
1421#ifdef MPI
1422        i=mod(peinf%inode,peinf%npes/peinf%npools)
1423        if (mod(wfnk%ndv,peinf%npes/peinf%npools).eq.0) then
1424          j=wfnk%ndv/(peinf%npes/peinf%npools)
1425        else
1426          j=wfnk%ndv/(peinf%npes/peinf%npools)+1
1427        endif
1428        g1=1+i*j
1429        g2=min(j+i*j,wfnk%ndv)
1430        if (g2.ge.g1) then
1431          wfnkmpi%cg(1:g2-g1+1,k,ikn)=wfnk%zk(g1:g2,k)
1432        endif ! g2.ge.g1
1433#else
1434        wfnkmpi%cg(1:wfnk%ndv,k,ikn)=wfnk%zk(1:wfnk%ndv,k)
1435#endif
1436      enddo
1437
1438      if ( peinf%inode == 0 ) call timing%stop(timing%input_write)
1439
1440      ! ZL: if sig%nkn.gt.1, clean wfnk%zk, otherwise (only one kpoint in outer) we keep it
1441      SAFE_DEALLOCATE_P(wfnk%zk)
1442    endif
1443
1444    ! ZL: for do_phonq
1445    if(do_phonq.and.(istore_ep.eq.1).and.(sig%nkn.gt.1)) then
1446
1447      if ( peinf%inode == 0 ) call timing%start(timing%input_write)
1448
1449      ikn=iknstore_ep
1450      wfnk_phonq_mpi%nkptotal(ikn)=kp%ngk(irk)
1451      wfnk_phonq_mpi%isort(1:kp%ngk(irk),ikn)=wfnk_phonq%isrtk(1:kp%ngk(irk))
1452      wfnk_phonq_mpi%qk(1:3,ikn)=wfnk_phonq%k(:)
1453      wfnk_phonq_mpi%el(1:sig%ntband,1:sig%nspin,ikn)= &
1454        wfnk_phonq%ek(1:sig%ntband,1:sig%nspin)
1455      wfnk_phonq_mpi%elda(1:sig%ntband,1:sig%nspin,ikn)= &
1456        wfnk_phonq%elda(1:sig%ntband,1:sig%nspin)
1457      do k=1,sig%nspin*kp%nspinor
1458#ifdef MPI
1459        i=mod(peinf%inode,peinf%npes/peinf%npools)
1460        if (mod(wfnk_phonq%ndv,peinf%npes/peinf%npools).eq.0) then
1461          j=wfnk_phonq%ndv/(peinf%npes/peinf%npools)
1462        else
1463          j=wfnk_phonq%ndv/(peinf%npes/peinf%npools)+1
1464        endif
1465        g1=1+i*j
1466        g2=min(j+i*j,wfnk_phonq%ndv)
1467        if (g2.ge.g1) then
1468          wfnk_phonq_mpi%cg(1:g2-g1+1,k,ikn)=wfnk_phonq%zk(g1:g2,k)
1469        endif ! g2.ge.g1
1470#else
1471        wfnk_phonq_mpi%cg(1:wfnk_phonq%ndv,k,ikn)=wfnk_phonq%zk(1:wfnk_phonq%ndv,k)
1472#endif
1473      enddo
1474
1475      if ( peinf%inode == 0 ) call timing%stop(timing%input_write)
1476
1477      ! ZL: if sig%nkn.gt.1, remove wfnk%zk, otherwise (only one kpoint in outer) we keep it
1478      SAFE_DEALLOCATE_P(wfnk_phonq%zk)
1479    endif
1480
1481  enddo ! irk (loop over k-points)
1482  call progress_free(prog_info)
1483
1484  SAFE_DEALLOCATE(isort)
1485
1486  PUSH_SUB(read_wavefunctions)
1487
1488end subroutine read_wavefunctions
1489
1490
1491end module input_m
1492
1493
1494