1!-----------------------------------------------------------------------
2subroutine read_export (pp_file,kunit,uspp_spsi, ascii, single_file, raw)
3  !-----------------------------------------------------------------------
4  !
5  use kinds,          ONLY : DP
6  use pwcom
7  use control_flags,  ONLY : gamma_only
8  use becmod,         ONLY : bec_type, becp, calbec, &
9                             allocate_bec_type, deallocate_bec_type
10!  use symme,          ONLY : nsym, s, invsym, sname, irt, ftau
11!  use symme,          ONLY : nsym, s, invsym, irt, ftau
12!  use char,           ONLY : sname
13! sname is in symme which is now outside pwcom
14  use  uspp,          ONLY : nkb, vkb
15  use wavefunctions,  ONLY : evc
16  use io_files,       ONLY : prefix, iunwfc, nwordwfc, iunsat, nwordatwfc
17  use io_files,       ONLY : pseudo_dir, psfile
18  use io_global,      ONLY : ionode, stdout
19  USE ions_base,      ONLY : atm, nat, ityp, tau, nsp
20  use mp_pools,       ONLY : nproc_pool, my_pool_id, intra_pool_comm, inter_pool_comm
21  USE mp_world, ONLY : world_comm, nproc
22  use mp,             ONLY : mp_sum, mp_max
23  use ldaU,           ONLY :  lda_plus_u
24  USE gvecw,     ONLY : ecutwfc
25  USE cell_base, ONLY : bg
26  USE gvect,                ONLY : ngm, ngm_g,g, gstart, ig_l2g, mill
27!  USE pwcom, ONLY : igk_k
28  USE cell_base,            ONLY : tpiba2
29
30  implicit none
31
32  CHARACTER(5), PARAMETER :: fmt_name="QEXPT"
33  CHARACTER(5), PARAMETER :: fmt_version="1.1.0"
34
35  integer, intent(in) :: kunit
36  character(80), intent(in) :: pp_file
37  logical, intent(in) :: uspp_spsi, ascii, single_file, raw
38
39  integer :: i, j, k, ig, ik, ibnd, na, ngg,ig_, ierr
40  integer, allocatable :: kisort(:)
41  real(DP) :: xyz(3), tmp(3)
42  integer :: npool, nkbl, nkl, nkr, npwx_g
43  integer :: ike, iks, npw_g, ispin, local_pw
44  integer, allocatable :: ngk_g( : )
45  integer, allocatable :: itmp_g( :, : )
46  real(DP),allocatable :: rtmp_g( :, : )
47  real(DP),allocatable :: rtmp_gg( : )
48  integer, allocatable :: itmp1( : )
49  integer, allocatable :: igwk( :, : )
50  integer, allocatable :: l2g_new( : )
51  integer, allocatable :: igk_l2g( :, : )
52
53
54  real(DP) :: wfc_scal
55  logical :: twf0, twfm
56  complex(DP), allocatable :: sevc (:,:)
57
58  CHARACTER(LEN=256) :: outdir
59
60  write(stdout,*) "nkstot=", nkstot
61
62  IF( nkstot > 0 ) THEN
63
64     IF( ( kunit < 1 ) .OR. ( MOD( nkstot, kunit ) /= 0 ) ) &
65       CALL errore( ' write_export ',' wrong kunit ', 1 )
66
67     IF( ( nproc_pool > nproc ) .OR. ( MOD( nproc, nproc_pool ) /= 0 ) ) &
68       CALL errore( ' write_export ',' nproc_pool ', 1 )
69
70     !  find out the number of pools
71     npool = nproc / nproc_pool
72
73     !  find out number of k points blocks
74     nkbl = nkstot / kunit
75
76     !  k points per pool
77     nkl = kunit * ( nkbl / npool )
78
79     !  find out the reminder
80     nkr = ( nkstot - nkl * npool ) / kunit
81
82     !  Assign the reminder to the first nkr pools
83     IF( my_pool_id < nkr ) nkl = nkl + kunit
84
85     !  find out the index of the first k point in this pool
86     iks = nkl * my_pool_id + 1
87     IF( my_pool_id >= nkr ) iks = iks + nkr * kunit
88
89     !  find out the index of the last k point in this pool
90     ike = iks + nkl - 1
91
92  END IF
93
94  write(stdout,*) "after first init"
95
96  ! find out the global number of G vectors: ngm_g
97  ngm_g = ngm
98  call mp_sum( ngm_g , intra_pool_comm)
99
100  ! collect all G vectors across processors within the pools
101  ! and compute their modules
102  !
103  allocate( itmp_g( 3, ngm_g ) )
104  allocate( rtmp_g( 3, ngm_g ) )
105  allocate( rtmp_gg( ngm_g ) )
106
107  itmp_g = 0
108  do  ig = 1, ngm
109    itmp_g( 1, ig_l2g( ig ) ) = mill(1, ig )
110    itmp_g( 2, ig_l2g( ig ) ) = mill(2, ig )
111    itmp_g( 3, ig_l2g( ig ) ) = mill(3, ig )
112  end do
113  call mp_sum( itmp_g , intra_pool_comm)
114  !
115  ! here we are in crystal units
116  rtmp_g(1:3,1:ngm_g) = REAL( itmp_g(1:3,1:ngm_g) )
117  !
118  ! go to cartesian units (tpiba)
119  call cryst_to_cart( ngm_g, rtmp_g, bg , 1 )
120  !
121  ! compute squared moduli
122  do  ig = 1, ngm_g
123     rtmp_gg(ig) = rtmp_g(1,ig)**2 + rtmp_g(2,ig)**2 + rtmp_g(3,ig)**2
124  enddo
125  deallocate( rtmp_g )
126
127  ! build the G+k array indexes
128  allocate ( igk_l2g ( npwx, nks ) )
129  allocate ( kisort( npwx ) )
130  do ik = 1, nks
131     kisort = 0
132     npw = npwx
133     call gk_sort (xk (1, ik+iks-1), ngm, g, ecutwfc / tpiba2, npw, kisort(1), g2kin)
134     !
135     ! mapping between local and global G vector index, for this kpoint
136     !
137     DO ig = 1, npw
138        !
139        igk_l2g(ig,ik) = ig_l2g( kisort(ig) )
140        !
141     END DO
142     !
143     igk_l2g( npw+1 : npwx, ik ) = 0
144     !
145     ngk (ik) = npw
146  end do
147  deallocate (kisort)
148
149  ! compute the global number of G+k vectors for each k point
150  allocate( ngk_g( nkstot ) )
151  ngk_g = 0
152  ngk_g( iks:ike ) = ngk( 1:nks )
153  CALL mp_sum( ngk_g,world_comm )
154
155  ! compute the Maximum G vector index among all G+k and processors
156  npw_g = MAXVAL( igk_l2g(:,:) )
157  CALL mp_max( npw_g,world_comm )
158
159  ! compute the Maximum number of G vector among all k points
160  npwx_g = MAXVAL( ngk_g( 1:nkstot ) )
161
162  deallocate(rtmp_gg)
163
164  allocate( igwk( npwx_g,nkstot ) )
165
166  write(stdout,*) "after g stuff"
167
168! wfc grids
169
170  DO ik = 1, nkstot
171    igwk(:,ik) = 0
172    !
173    ALLOCATE( itmp1( npw_g ), STAT= ierr )
174    IF ( ierr/=0 ) CALL errore('pw_export','allocating itmp1', ABS(ierr) )
175    itmp1 = 0
176    !
177    IF( ik >= iks .AND. ik <= ike ) THEN
178      DO  ig = 1, ngk( ik-iks+1 )
179        itmp1( igk_l2g( ig, ik-iks+1 ) ) = igk_l2g( ig, ik-iks+1 )
180      END DO
181    END IF
182    !
183    CALL mp_sum( itmp1 ,world_comm)
184    !
185    ngg = 0
186    DO  ig = 1, npw_g
187      IF( itmp1( ig ) == ig ) THEN
188        ngg = ngg + 1
189        igwk( ngg , ik) = ig
190      END IF
191    END DO
192    IF( ngg /= ngk_g( ik ) ) THEN
193      WRITE( stdout,*) ' ik, ngg, ngk_g = ', ik, ngg, ngk_g( ik )
194    END IF
195    !
196    DEALLOCATE( itmp1 )
197    !
198  ENDDO
199  !
200  deallocate( itmp_g )
201
202  write(stdout,*)"after wfc waves"
203
204#ifdef __MPI
205  call poolrecover (et, nbnd, nkstot, nks)
206#endif
207
208  wfc_scal = 1.0d0
209  twf0 = .true.
210  twfm = .false.
211
212  do ik = 1, nkstot
213     local_pw = 0
214     IF( (ik >= iks) .AND. (ik <= ike) ) THEN
215
216       call davcio (evc, 2*nwordwfc, iunwfc, (ik-iks+1), - 1)
217       !IF ( lda_plus_u ) CALL davcio( swfcatom, nwordatwfc, iunsat, (ik-iks+1), -1 )
218       local_pw = ngk(ik-iks+1)
219
220     ENDIF
221
222
223     allocate(l2g_new(local_pw))
224
225     l2g_new = 0
226     do ig = 1, local_pw
227       ngg = igk_l2g(ig,ik-iks+1)
228       do ig_ = 1, ngk_g(ik)
229         if(ngg == igwk(ig_,ik)) then
230           l2g_new(ig) = ig_
231           exit
232         end if
233       end do
234     end do
235
236
237     ispin = isk( ik )
238     !  WRITE(0,*) ' ### ', ik,nkstot,iks,ike,kunit,nproc,nproc_pool
239     deallocate(l2g_new)
240  end do
241  !
242
243  write(stdout,*) "after davcio"
244
245  ! If specified and if USPP are used the wfcs S_psi are written
246  ! | spsi_nk > = \hat S | psi_nk >
247  ! where S is the overlap operator of US PP
248  !
249  IF ( uspp_spsi .AND. nkb > 0 ) THEN
250
251       ALLOCATE( sevc(npwx,nbnd), STAT=ierr )
252       IF (ierr/=0) CALL errore( ' read_export ',' Unable to allocate SEVC ', ABS(ierr) )
253
254       CALL allocate_bec_type (nkb,nbnd,becp)
255
256       do ik = 1, nkstot
257
258           local_pw = 0
259           IF( (ik >= iks) .AND. (ik <= ike) ) THEN
260
261               CALL gk_sort (xk (1, ik+iks-1), ngm, g, ecutwfc / tpiba2, npw, igk_k(1,ik), g2kin)
262               CALL davcio (evc, 2*nwordwfc, iunwfc, (ik-iks+1), - 1)
263
264               CALL init_us_2(npw, igk_k(1,ik), xk(1, ik), vkb)
265               local_pw = ngk(ik-iks+1)
266
267               IF ( gamma_only ) THEN
268                  if(nkb>0) CALL calbec ( ngk_g(ik), vkb, evc, becp )
269               ELSE
270                  CALL calbec ( npw, vkb, evc, becp )
271               ENDIF
272               CALL s_psi(npwx, npw, nbnd, evc, sevc)
273           ENDIF
274
275           ALLOCATE(l2g_new(local_pw))
276
277           l2g_new = 0
278           DO ig = 1, local_pw
279             ngg = igk_l2g(ig,ik-iks+1)
280             DO ig_ = 1, ngk_g(ik)
281               IF(ngg == igwk(ig_,ik)) THEN
282                 l2g_new(ig) = ig_
283                 EXIT
284               ENDIF
285             ENDDO
286           ENDDO
287
288           ispin = isk( ik )
289           DEALLOCATE(l2g_new)
290       ENDDO
291
292       DEALLOCATE( sevc, STAT=ierr )
293       IF ( ierr/= 0 ) CALL errore('read_export','Unable to deallocate SEVC',ABS(ierr))
294       CALL deallocate_bec_type ( becp )
295  ENDIF
296
297  DEALLOCATE( igk_l2g )
298  DEALLOCATE( igwk )
299  DEALLOCATE ( ngk_g )
300
301end subroutine read_export
302