1!===============================================================================
2!
3! Routines:
4!
5! (1) input_outer()       Originally By gsm       Last Modified 10/5/2009 (gsm)
6!
7! Adapted from subroutine input. Reads in and distributes outer
8! wavefunctions. Stores wavefunctions in structures distributed in memory.
9!
10!===============================================================================
11
12#include "f_defs.h"
13
14subroutine input_outer(crys,gvec,syms,kp,sig,wfnk,iunit_k,fnk,wfnkmpi)
15
16  use global_m
17  use eqpcor_m
18  use input_utils_m
19  use misc_m
20  use scissors_m
21  use wfn_rho_vxc_io_m
22  use io_utils_m
23  implicit none
24
25  type (crystal), intent(in) :: crys
26  type (symmetry), intent(in) :: syms
27  type (gspace), intent(in) :: gvec
28  type (kpoints), intent(in) :: kp
29  type (siginfo), intent(inout) :: sig
30  type (wfnkstates), intent(inout) :: wfnk !< do not lose the allocations from input
31  integer, intent(in) :: iunit_k
32  character*20, intent(in) :: fnk
33  type (wfnkmpiinfo), intent(inout) :: wfnkmpi !< do not lose the allocations from input
34
35  character :: fncor*32
36  integer :: i,j,k,ikn,irk,nbnmin,nbnmax, &
37    iouter,istore,nstore,nstore_sum,g1,g2,iknstore, &
38    ntband_outer
39  integer, allocatable :: isort(:)
40  real(DP) :: qk(3)
41  SCALAR, allocatable :: zc(:,:)
42  logical :: dont_read
43
44  character(len=3) :: sheader
45  integer :: iflavor
46  type(gspace) :: gvec_outer, gvec_kpt
47  type(kpoints) :: kp_outer
48  type(crystal) :: crys_outer
49  type(symmetry) :: syms_outer
50  logical, allocatable :: kpt_outer_required(:)
51  type(progress_info) :: prog_info !< a user-friendly progress report
52
53  PUSH_SUB(input_outer)
54
55  call logit('Opening WFN_outer')
56
57!-------------------------------------------
58! Try to open WFN_outer file
59! Return to caller if error is encountered
60
61  if(peinf%inode.eq.0) call open_file(25,file='WFN_outer',form='unformatted',status='old',iostat=iouter)
62#ifdef MPI
63  call MPI_Bcast(iouter,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)
64#endif
65
66  sig%wfn_outer_present = (iouter == 0)
67  if (iouter.ne.0) then
68    POP_SUB(input_outer)
69    return
70  endif ! iouter
71
72  sheader = 'WFN'
73  iflavor = 0
74  call read_binary_header_type(25, sheader, iflavor, kp_outer, gvec_outer, syms_outer, crys_outer, warn = .false.)
75  call check_trunc_kpts(sig%icutv, kp_outer)
76
77  call check_header('WFN_inner', kp, gvec, syms, crys, 'WFN_outer', kp_outer, gvec_outer, syms_outer, crys_outer, is_wfn = .true.)
78
79  ! DAS: there is no fundamental reason for the condition below, but seems to be assumed in the code
80  ! and it would require some rewriting to make that kp_outer need only contain the kpts in sig%kpt (the ones in the input file).
81  if(kp_outer%nrk /= kp%nrk) call die('different number of k-points in WFN_outer')
82  if(kp_outer%mnband < maxval(sig%diag)) call die('not enough bands in WFN_outer')
83
84  ntband_outer = min(kp_outer%mnband, sig%ntband)
85
86  kp_outer%el(:,:,:) = kp_outer%el(:,:,:) - sig%avgpot_outer / ryd
87
88  SAFE_ALLOCATE(kp_outer%elda, (ntband_outer,kp_outer%nrk,kp_outer%nspin))
89  kp_outer%elda(1:ntband_outer,1:kp_outer%nrk,1:kp_outer%nspin)=ryd*kp_outer%el(1:ntband_outer,1:kp_outer%nrk,1:kp_outer%nspin)
90  call scissors_shift(kp_outer, sig%scis_outer, sig%spl_tck_outer)
91  kp_outer%el(:,:,:)=ryd*kp_outer%el(:,:,:)
92
93!---------------------------------------------------------------
94! If quasi-particle corrections requested, read in qp energies
95  if(sig%eqp_outer_corrections) then
96    fncor='eqp_outer.dat'
97    nbnmin=sig%ntband+1
98    nbnmax=0
99    do i=1,sig%ndiag
100      if (sig%diag(i).lt.nbnmin) nbnmin=sig%diag(i)
101      if (sig%diag(i).gt.nbnmax) nbnmax=sig%diag(i)
102    enddo
103
104    ! FHJ: Determine which k-points from WFN_outer/eqp_outer.dat we need.
105    SAFE_ALLOCATE(kpt_outer_required, (kp_outer%nrk))
106    kpt_outer_required(:) = .false.
107    do ikn = 1, sig%nkn
108      kpt_outer_required(sig%indkn(ikn)) = .true.
109    enddo
110    call eqpcor(fncor,peinf%inode,peinf%npes,kp_outer,nbnmin,nbnmax,0,0,&
111      kp_outer%el,kp_outer%el,kp_outer%el,0,0,kpt_required=kpt_outer_required)
112    SAFE_DEALLOCATE(kpt_outer_required)
113  endif
114
115!-------------------
116! Read in g-vectors
117
118  SAFE_ALLOCATE(gvec_outer%components, (3, gvec%ng))
119  call read_binary_gvectors(25, gvec%ng, gvec%ng, gvec_outer%components)
120  SAFE_DEALLOCATE_P(gvec_outer%components)
121
122  SAFE_ALLOCATE(isort, (gvec%ng))
123
124!---------------------
125! Loop over k-points
126
127  call progress_init(prog_info, 'reading wavefunctions (WFN_outer)', 'state', kp_outer%nrk*ntband_outer)
128  do irk=1,kp_outer%nrk
129    qk(:)=kp_outer%rk(:,irk)
130
131!--------------------------------------------
132! Read in and sort g-vectors for k-point irk
133
134    SAFE_ALLOCATE(gvec_kpt%components, (3, kp_outer%ngk(irk)))
135    call read_binary_gvectors(25, kp_outer%ngk(irk), kp_outer%ngk(irk), gvec_kpt%components)
136
137    do i = 1, kp_outer%ngk(irk)
138      call findvector(isort(i), gvec_kpt%components(:, i), gvec)
139      if (isort(i) == 0) call die('outer wfn: could not find G-vector')
140    enddo
141
142    SAFE_DEALLOCATE_P(gvec_kpt%components)
143
144!--------------------------------------------------------
145! Determine if Sigma must be computed for this k-point.
146! If so, store the bands and wavefunctions on file iunit_k.
147! If there is only one k-point, store directly in wfnk.
148
149    istore=0
150    do ikn=1,sig%nkn
151      if(sig%indkn(ikn).eq.irk) then
152        istore=1
153        iknstore=ikn
154        wfnk%nkpt=kp_outer%ngk(irk)
155        wfnk%ndv=peinf%ndiag_max*kp_outer%ngk(irk)
156        wfnk%k(:)=qk(:)
157        wfnk%isrtk(:)=isort(:)
158        do k=1,sig%nspin
159          wfnk%ek(1:ntband_outer,k)= &
160            kp_outer%el(1:ntband_outer,irk,sig%spin_index(k))
161          wfnk%elda(1:ntband_outer,k)= &
162            kp_outer%elda(1:ntband_outer,irk,sig%spin_index(k))
163        enddo
164        SAFE_ALLOCATE(wfnk%zk, (wfnk%ndv,sig%nspin*kp_outer%nspinor))
165        wfnk%zk=ZERO
166      endif
167    enddo ! ikn
168
169!---------------------------------------------------------------
170! Read the wavefunctions, check the norm, and store in wfnk%zk
171
172    SAFE_ALLOCATE(zc, (kp_outer%ngk(irk),kp_outer%nspin*kp_outer%nspinor))
173    do i=1,kp_outer%mnband
174
175! Check whether i-th band is needed for diagonal Sigma matrix elements.
176! We don`t need to check for off-diagonal Sigma matrix elements separately
177! because whenever off-diagonals are requested the corresponding diagonals
178! are also included in the calculation as enforced in Sigma/inread.f90.
179
180      if (i.le.ntband_outer) then
181        call progress_step(prog_info, ntband_outer*(irk-1) + i)
182        nstore=0
183        do j=1,peinf%ndiag_max
184          if (.not.peinf%flag_diag(j)) cycle
185          if (i==sig%diag(peinf%index_diag(j))) nstore=j
186        enddo
187#ifdef MPI
188        call MPI_Allreduce(nstore,nstore_sum,1, &
189          MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,mpierr)
190#else
191        nstore_sum=nstore
192#endif
193      else
194        nstore_sum=0
195      endif
196
197      dont_read = (istore == 0 .OR. nstore_sum == 0)
198      call read_binary_data(25, kp_outer%ngk(irk), kp_outer%ngk(irk), kp_outer%nspin*kp_outer%nspinor, zc, dont_read = dont_read)
199
200      if(.NOT. dont_read) then
201        if(peinf%inode.eq.0) then
202          do k=1,sig%nspin
203            call checknorm('WFN_outer',i,irk,kp_outer%ngk(irk),k,kp_outer%nspinor,zc(:,:))
204          enddo
205        endif
206
207        if((istore.eq.1).and.(nstore.ne.0)) then
208          do k=1,sig%nspin*kp_outer%nspinor
209            do j=1,kp_outer%ngk(irk)
210              wfnk%zk((nstore-1)*kp_outer%ngk(irk)+j,k) = zc(j,sig%spin_index(k))
211            enddo
212          enddo
213        endif
214      else
215        ! FHJ: the following lines were introduced in r6294 and are supposed to
216        ! be a shortcut if we are past the last band of the last k-point. However,
217        ! in light of a previous bug (#223), this feature is commented out for now.
218        !! FHJ: shortcut if this is past the last band of the last k-point
219        !if (irk==kp_outer%nrk) exit
220      endif ! .NOT. dont_read
221    enddo ! i (loop over bands)
222    SAFE_DEALLOCATE(zc)
223
224!---------------------------------------------------------
225! Write the wavefunctions stored in wfnk%zk to file iunit_k
226
227    if((istore.eq.1).and.(sig%nkn.gt.1)) then
228      ikn=iknstore
229      wfnkmpi%nkptotal(ikn)=kp_outer%ngk(irk)
230      wfnkmpi%isort(1:kp_outer%ngk(irk),ikn)=wfnk%isrtk(1:kp_outer%ngk(irk))
231      wfnkmpi%qk(1:3,ikn)=qk(1:3)
232      wfnkmpi%el(1:ntband_outer,1:sig%nspin,ikn)=wfnk%ek(1:ntband_outer,1:sig%nspin)
233      wfnkmpi%elda(1:ntband_outer,1:sig%nspin,ikn)=wfnk%elda(1:ntband_outer,1:sig%nspin)
234      do k=1,sig%nspin*kp_outer%nspinor
235#ifdef MPI
236        i=mod(peinf%inode,peinf%npes/peinf%npools)
237        if (mod(wfnk%ndv,peinf%npes/peinf%npools).eq.0) then
238          j=wfnk%ndv/(peinf%npes/peinf%npools)
239        else
240          j=wfnk%ndv/(peinf%npes/peinf%npools)+1
241        endif
242        g1=1+i*j
243        g2=min(j+i*j,wfnk%ndv)
244        if (g2.ge.g1) then
245          wfnkmpi%cg(1:g2-g1+1,k,ikn)=wfnk%zk(g1:g2,k)
246        endif ! g2.ge.g1
247#else
248        wfnkmpi%cg(1:wfnk%ndv,k,ikn)=wfnk%zk(1:wfnk%ndv,k)
249#endif
250      enddo
251      SAFE_DEALLOCATE_P(wfnk%zk)
252    endif
253
254  enddo ! irk (loop over k-points)
255  call progress_free(prog_info)
256
257  if(peinf%inode.eq.0) then
258    call close_file(25) ! WFN_outer
259  endif
260
261  SAFE_DEALLOCATE(isort)
262  SAFE_DEALLOCATE_P(kp_outer%w)
263  SAFE_DEALLOCATE_P(kp_outer%el)
264  SAFE_DEALLOCATE_P(kp_outer%elda)
265
266  POP_SUB(input_outer)
267
268  return
269
270end subroutine input_outer
271