1!------------------------------------------------------------------------
2PROGRAM open_grid
3  !------------------------------------------------------------------------
4  !
5  USE kinds, ONLY : DP
6  USE io_global,  ONLY : stdout, ionode, ionode_id
7  USE mp_global,  ONLY : mp_startup
8  USE mp_images,  ONLY : intra_image_comm
9  USE mp_pools,   ONLY : npool
10  USE mp,         ONLY : mp_bcast
11  USE cell_base,  ONLY : at, bg, tpiba2, alat
12  USE klist,      ONLY : nks, nkstot, xk, wk, igk_k, ngk, qnorm
13  USE io_files,   ONLY : prefix, tmp_dir, nwordwfc, iunwfc, diropn
14  USE noncollin_module,   ONLY : noncolin, m_loc, angle1, angle2, nspin_lsda
15  USE spin_orb,           ONLY : domag
16  USE control_flags,      ONLY : gamma_only
17  USE environment,        ONLY : environment_start, environment_end
18  USE ions_base,          ONLY : nat, tau, ityp
19  USE symm_base,          ONLY : nrot, nsym, s, t_rev, fft_fact, find_sym
20  USE parameters,         ONLY : npk
21  USE exx_base,           ONLY : nq1,nq2,nq3, xkq_collect, &
22                                 nkqs, exx_mp_init, index_xk, exx_grid_init
23  USE exx,                ONLY : exxbuff, exxinit, use_ace, ecutfock
24  USE gvecw,              ONLY : ecutwfc, gcutw
25  USE gvect,              ONLY : g, ngm
26  USE funct,              ONLY : dft_force_hybrid
27  USE wvfct,              ONLY : nbnd, npwx, g2kin, et, wg
28  USE wavefunctions, ONLY : evc
29  USE buffers,            ONLY : save_buffer, open_buffer, close_buffer
30  USE scf,                ONLY : rho
31  USE lsda_mod,           ONLY : nspin, isk, lsda, starting_magnetization
32  USE io_rho_xml,         ONLY : write_scf
33  USE noncollin_module,   ONLY : nspin_mag, npol
34  USE fft_interfaces,     ONLY : fwfft
35  !
36  USE fft_base,           ONLY : dffts
37  USE control_flags,      ONLY : gamma_only, io_level
38  USE start_k, ONLY : init_start_k
39  USE extfield,           ONLY : gate
40  !
41  IMPLICIT NONE
42  !
43  CHARACTER(LEN=256), EXTERNAL :: trimcheck
44  !
45  INTEGER :: ios, ik, ibnd, ik_idx, ik_idx_kpt, ik_idx_exx, is, na
46  CHARACTER(len=4) :: spin_component
47  CHARACTER(len=256) :: outdir
48  INTEGER :: k1, k2, k3
49  LOGICAL :: exst, opnd, exst_mem, magnetic_sym
50  REAL(DP),ALLOCATABLE :: et0(:,:), wg0(:,:), yk(:,:), wk0(:)
51  INTEGER, EXTERNAL  :: n_plane_waves
52  COMPLEX(DP),ALLOCATABLE :: psic(:), evx(:,:)
53  !
54  LOGICAL           :: use_ace_back, exx_status_back
55  REAL(DP)          :: ecutfock_back
56  INTEGER           :: nq_back(3)
57
58  ! these are in wannier module.....-> integer :: ispinw, ikstart, ikstop, iknum
59  NAMELIST / inputpp / outdir, prefix !, nq
60  !
61  ! initialise environment
62  !
63#if defined(__MPI)
64  CALL mp_startup ( )
65#endif
66  !! not sure if this should be called also in 'library' mode or not !!
67  CALL environment_start ( 'OPEN_GRID' )
68  !
69  ios = 0
70  IF(ionode) THEN
71     !
72     ! Check to see if we are reading from a file
73     !
74     CALL input_from_file()
75     !
76     !   set default values for variables in namelist
77     !
78     CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir )
79     IF ( trim( outdir ) == ' ' ) outdir = './'
80     prefix = ' '
81     !nq = 0
82     !
83     READ (5, inputpp, iostat=ios)
84     tmp_dir = trimcheck(outdir)
85  ENDIF
86  !
87  !
88  CALL mp_bcast(outdir,ionode_id, intra_image_comm)
89  CALL mp_bcast(tmp_dir,ionode_id, intra_image_comm)
90  CALL mp_bcast(prefix,ionode_id, intra_image_comm)
91  !
92  WRITE(stdout,*)
93  WRITE(stdout,*) ' Reading nscf_save data'
94  CALL read_file
95  !print*, "initial", nks, nkstot
96  CALL open_buffer(iunwfc, 'wfc', nwordwfc, io_level, exst_mem, exst)
97  !
98  WRITE(stdout,*)
99  IF ( npool > 1 .and. nspin_mag>1) CALL errore( 'open_grid', &
100      'pools+spin not implemented, run this tool with only one pool', npool )
101  !
102  IF(gamma_only) CALL errore("open_grid", &
103      "not implemented, and pointless, for gamma-only",1)
104  !
105  ! Store some variables related to exx to be put back before writing
106  use_ace_back = use_ace
107  ecutfock_back = ecutfock
108  nq_back = (/ nq1, nq2, nq3 /)
109  exx_status_back = .true.
110  CALL dft_force_hybrid(exx_status_back)
111
112  magnetic_sym = noncolin .AND. domag
113  ALLOCATE(m_loc(3,nat))
114  IF (noncolin.and.domag) THEN
115     DO na = 1, nat
116        m_loc(1,na) = starting_magnetization(ityp(na)) * &
117                      SIN( angle1(ityp(na)) ) * COS( angle2(ityp(na)) )
118        m_loc(2,na) = starting_magnetization(ityp(na)) * &
119                      SIN( angle1(ityp(na)) ) * SIN( angle2(ityp(na)) )
120        m_loc(3,na) = starting_magnetization(ityp(na)) * &
121                      COS( angle1(ityp(na)) )
122     ENDDO
123  ENDIF
124  CALL find_sym ( nat, tau, ityp, magnetic_sym, m_loc, gate )
125
126  nq1 = -1
127  nq2 = -1
128  nq3 = -1
129  ecutfock = 4*ecutwfc
130  use_ace = .false.
131
132  CALL exx_grid_init()
133  CALL exx_mp_init()
134  !
135  CALL exxinit(.false.)
136  qnorm = 0._dp
137  !
138  CALL close_buffer(iunwfc,'keep')
139  !
140  ALLOCATE(et0(nbnd,nks), wg0(nbnd,nks), wk0(nks))
141  et0 = et
142  wg0 = wg
143  wk0 = wk(1:nks)
144  !
145  nks = nkqs
146  nkstot = nks
147  ! Temporary order of xk points, they will be resorted later (if nspin>1)
148  xk(:,1:nks) = xkq_collect(:,1:nks)
149  wk(1:nks) = 1._dp/DBLE(nks) !/DBLE(nspin_mag)
150  npwx = n_plane_waves(gcutw, nks, xk, g, ngm)
151  IF (nspin==2) THEN
152    isk(1:nks/2) = 1
153    isk(nks/2+1:nks) = 2
154  ENDIF
155  !
156  DEALLOCATE(igk_k, ngk, et, wg, g2kin)
157  ALLOCATE(igk_k(npwx,nks), ngk(nks), g2kin(npwx))
158  ALLOCATE(et(nbnd,nks), wg(nbnd,nks))
159  !
160  DEALLOCATE(evc)
161  ALLOCATE(evc(npwx*npol,nbnd))
162  !
163  prefix = TRIM(prefix)//"_open"
164  nwordwfc = nbnd * npwx * npol
165  CALL open_buffer(iunwfc, 'wfc', nwordwfc, +1, exst_mem, exst)
166  !
167  ! Write everything again with the new prefix
168  CALL write_scf(rho, nspin)
169  !
170  ALLOCATE(psic(dffts%nnr), evx(npol*npwx, nbnd))
171
172  DO ik = 1, nks
173    ik_idx_kpt = ik !+ (is-1)*(nks/nspin_mag) !(ik-1)*nspin_mag + is
174    ik_idx_exx = ik !+ (is-1)*(nks/nspin_mag)
175    xk(:,ik_idx_kpt) = xkq_collect(:,ik_idx_exx)
176!    wk(ik_idx_kpt) = 1._dp/nkqs*DBLE(nspin_mag)
177    !
178    CALL gk_sort (xk(:,ik_idx_kpt), ngm, g, ecutwfc / tpiba2, &
179                  ngk(ik_idx_kpt), igk_k(:,ik_idx_kpt), g2kin)
180    evx(:,:) = 0._dp
181    DO ibnd = 1, nbnd
182      psic(1:dffts%nnr) = exxbuff(1:dffts%nnr,ibnd,ik_idx_exx)
183      CALL fwfft('Wave', psic, dffts)
184      evx(1:ngk(ik_idx_kpt),ibnd) = psic(dffts%nl(igk_k(1:ngk(ik_idx_kpt),ik_idx_kpt)))
185      !
186      IF(noncolin)THEN
187        psic(1:dffts%nnr) = exxbuff(dffts%nnr+1:2*dffts%nnr,ibnd,ik_idx_exx)
188        CALL fwfft('Wave', psic, dffts)
189        evx(npwx+1:npwx+ngk(ik_idx_kpt),ibnd) = &
190                                  psic(dffts%nl(igk_k(1:ngk(ik_idx_kpt),ik_idx_kpt)))
191      ENDIF
192      !
193    ENDDO
194    !
195    CALL save_buffer( evx, nwordwfc, iunwfc, ik_idx_kpt )
196    et(:,ik_idx_kpt) = et0(:,index_xk(ik_idx_kpt))
197    wg(:,ik_idx_kpt) = wg0(:,index_xk(ik_idx_kpt))/wk0(index_xk(ik_idx_kpt))*wk(ik_idx_kpt)
198  !ENDDO
199  ENDDO
200  DEALLOCATE(psic, et0, wg0)
201  !
202  ! reconstruct input variables
203  k1 = 0
204  k2 = 0
205  k3 = 0
206  CALL init_start_k(nq1,nq2,nq3, k1, k2, k3, "automatic",nks/nspin_lsda, xk, wk)
207  !
208  ! Restore EXX variables
209  use_ace = use_ace_back
210  ecutfock = ecutfock_back
211  nq1 = nq_back(1)
212  nq2 = nq_back(2)
213  nq3 = nq_back(3)
214  CALL dft_force_hybrid(exx_status_back)
215  !
216  CALL punch('all')
217  !
218  ALLOCATE(yk(3,nks))
219  yk = xk
220  CALL cryst_to_cart(nks, yk, at, -1)
221  WRITE(stdout,'(5x,a)') "Grid of q-points"
222  WRITE(stdout,'(5x,a,3i4)') "Dimensions:", nq1, nq2, nq3
223  WRITE(stdout,'(5x,a,3i4)') "Shift:     ", k1,k2,k3
224  WRITE(stdout,'(5x,a)') "List to be put in the .win file of wannier90: &
225                          &(already in crystal/fractionary coordinates):"
226
227  DO ik = 1, nks/nspin_lsda
228    WRITE(stdout,'(3f21.15,3x,f13.10)') yk(:,ik), wk(ik)
229  ENDDO
230  DEALLOCATE(yk)
231  !
232  CALL environment_end ( 'OPEN_GRID' )
233  WRITE( stdout, *  )
234  CALL stop_pp
235  !
236END PROGRAM open_grid
237!
238