Line data Source code
1 : !!****m* ABINIT/m_wfk_analyze
2 : !! NAME
3 : !! m_wfk_analyze
4 : !!
5 : !! FUNCTION
6 : !! Post-processing tools for WFK file
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_wfk_analyze
23 :
24 : use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use m_abicore
27 : use m_xmpi
28 : use m_errors
29 : use m_hdr
30 : use m_crystal
31 : use m_ebands
32 : use m_nctk
33 : use m_wfd
34 : use m_dtset
35 : use m_dtfil
36 : use m_distribfft
37 :
38 : use m_io_tools, only : iomode_from_fname, get_unit
39 : use defs_datatypes, only : pseudopotential_type
40 : use defs_abitypes, only : mpi_type
41 : use m_time, only : timab
42 : use m_fstrings, only : strcat, sjoin, itoa, ftoa, ltoa, ktoa
43 : use m_fftcore, only : print_ngfft
44 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, init_mpi_enreg
45 : use m_esymm, only : esymm_t, esymm_free
46 : use m_ddk, only : ddkstore_t
47 : use m_ksdiago, only : psbands_t
48 : use m_pawang, only : pawang_type
49 : use m_pawrad, only : pawrad_type
50 : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
51 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
52 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
53 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init, pawfgrtab_print
54 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_inquire_dim
55 : use m_pawdij, only : pawdij, symdij
56 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
57 : use m_paw_sphharm, only : setsym_ylm
58 : use m_paw_init, only : pawinit, paw_gencond
59 : use m_paw_nhat, only : nhatgrid
60 : use m_paw_tools, only : chkpawovlp
61 : use m_paw_correlations,only : pawpuxinit
62 : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
63 : use m_classify_bands, only : classify_bands, dmats_t
64 : use m_kpts, only : kpts_ibz_from_kptrlatt
65 : use m_pspini, only : pspini
66 : use m_sigtk, only : sigtk_kpts_in_erange
67 : use m_iowf, only : prtkbff
68 : use m_wfd_wannier, only : wfd_run_wannier
69 : use m_wfk, only : wfk_to_bz, wfk_t, wfk_read_eigenvalues, wfk_check_symtab
70 :
71 : implicit none
72 :
73 : private
74 : !!***
75 :
76 : public :: wfk_analyze
77 : !!***
78 :
79 : contains
80 : !!***
81 :
82 : !!****f* ABINIT/wfk_analyze
83 : !! NAME
84 : !! wfk_analyze
85 : !!
86 : !! FUNCTION
87 : !! Main routine implementing postprocessing tools for the WFK file.
88 : !! Main differences wrt cut3d:
89 : !!
90 : !! - MPI support.
91 : !! - No interactive prompt.
92 : !! - Run the analysis once, store all the important results in netcdf files
93 : !! and use python tools to analyze data
94 : !!
95 : !! INPUTS
96 : !! acell(3)=Length scales of primitive translations (bohr)
97 : !! codvsn=Code version
98 : !! dtfil<datafiles_type>=Variables related to files.
99 : !! dtset<dataset_type>=All input variables for this dataset.
100 : !! pawang<pawang_type)>=PAW angular mesh and related data.
101 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
102 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
103 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
104 : !! Before entering the first time in the routine, a significant part of psps has been initialized :
105 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
106 : !! and the arrays dimensioned to npsp. All the remaining components of psps are to be initialized in
107 : !! the call to pspini. The next time the code enters bethe_salpeter, psps might be identical to the
108 : !! one of the previous dtset, in which case, no reinitialisation is scheduled in pspini.F90.
109 : !! rprim(3,3)=Dimensionless real space primitive translations.
110 : !! xred(3,natom)=Reduced atomic coordinates.
111 : !!
112 : !! NOTES
113 : !!
114 : !! ON THE USE OF FFT GRIDS:
115 : !! =================
116 : !! In case of PAW:
117 : !! ---------------
118 : !! Two FFT grids are used:
119 : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
120 : !! It is defined by nfft, ngfft, mgfft, ...
121 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
122 : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
123 : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
124 : !! In case of norm-conserving:
125 : !! ---------------------------
126 : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
127 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
128 : !!
129 : !! SOURCE
130 :
131 12 : subroutine wfk_analyze(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred)
132 :
133 : !Arguments ------------------------------------
134 : !scalars
135 : character(len=8),intent(in) :: codvsn
136 : type(datafiles_type),intent(in) :: dtfil
137 : type(dataset_type),intent(inout) :: dtset
138 : type(pawang_type),intent(inout) :: pawang
139 : type(pseudopotential_type),intent(inout) :: psps
140 : !arrays
141 : real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom)
142 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
143 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
144 :
145 : !Local variables ------------------------------
146 : !scalars
147 : integer,parameter :: master = 0, formeig0 = 0
148 : integer :: comm,nprocs,my_rank,mgfftf,nfftf !,nfftf_tot
149 : integer :: optcut,optgr0,optgr1,optgr2,optrad,psp_gencond,ii
150 : !integer :: option,option_test,option_dij,optrhoij
151 : integer :: band,ik_ibz,spin,nband_k, islice, ib, mpw, mcg, nb, npw_k ! first_band,last_band,
152 : integer :: ierr,usexcnhat, sc_mode, nspinor, nsto
153 : integer :: cplex,cplex_dij,cplex_rhoij,ndij,nspden_rhoij,gnt_option
154 : real(dp),parameter :: spinmagntarget=-99.99_dp
155 : real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff,gsqcut_shp, gs_fermie
156 : !real(dp) :: cpu,wall,gflops
157 : !real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,nelect,norm,oldefermi
158 : character(len=500) :: msg
159 : character(len=fnlen) :: wfk0_path, outwfk_path
160 : logical :: call_pawinit !, use_paw_aeur
161 36 : type(hdr_type) :: wfk0_hdr, hdr_bz, out_hdr
162 1248 : type(crystal_t) :: cryst, cryst_dtset
163 24 : type(ebands_t) :: ebands, ebands_bz
164 : type(pawfgr_type) :: pawfgr
165 : !type(paw_dmft_type) :: paw_dmft
166 12 : type(mpi_type) :: mpi_enreg
167 12 : type(wfd_t) :: wfd
168 12 : type(ddkstore_t) :: ds
169 96 : type(wfk_t) :: in_wfk, out_wfk
170 12 : type(dmats_t) :: dmats
171 : !type(dataset_type) :: my_dtset
172 : !arrays
173 24 : integer :: ngfftc(18),ngfftf(18), units(2), band_block(2), bstart, brange_spin(2, dtset%nsppol)
174 : integer :: nkibz_full, nkbz_full, ikbz, nstar_fail
175 12 : integer,allocatable :: l_size_atm(:), kg_k(:,:)
176 24 : real(dp),allocatable :: wtk_full(:), kibz_full(:,:), kbz_full(:,:)
177 : real(dp),parameter :: k0(3)=zero
178 12 : real(dp),pointer :: gs_eigen(:,:,:)
179 12 : real(dp),allocatable :: eig_k(:), occ_k(:), thetas(:) !, out_cg(:,:), work(:,:,:,:), allcg_k(:,:)
180 12 : real(dp),allocatable,target :: cg_k(:,:)
181 12 : complex(gwp),allocatable :: ur_ae(:)
182 12 : complex(dp),pointer :: cg_k_cplx(:,:)
183 12 : complex(dp),allocatable :: ps_ug(:,:)
184 12 : logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:)
185 : real(dp) :: tsec(2)
186 12 : type(Pawrhoij_type),allocatable :: pawrhoij(:)
187 12 : type(pawfgrtab_type),allocatable :: pawfgrtab(:)
188 : !type(paw_ij_type),allocatable :: paw_ij(:)
189 : !type(paw_an_type),allocatable :: paw_an(:)
190 : !type(esymm_t),allocatable :: esymm(:,:)
191 12 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
192 12 : type(psbands_t),allocatable :: psb_ks(:,:)
193 : !************************************************************************
194 :
195 : DBG_ENTER('COLL')
196 :
197 : ! abirules!
198 : if (.False.) write(std_out,*)acell,codvsn,rprim,xred
199 36 : units = [std_out, ab_out]
200 :
201 12 : comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
202 :
203 12 : wfk0_path = dtfil%fnamewffk
204 12 : if (my_rank == master) then
205 : ! Accept WFK file in Fortran or netcdf format.
206 12 : if (nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then
207 0 : ABI_ERROR(sjoin("Cannot find GS WFK file:", ch10, msg))
208 : end if
209 : end if
210 12 : call xmpi_bcast(wfk0_path, master, comm, ierr)
211 12 : call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path))
212 :
213 : !call cwtime(cpu,wall,gflops,"start")
214 :
215 : ! Construct crystal and ebands from the GS WFK file.
216 12 : call wfk_read_eigenvalues(wfk0_path, gs_eigen, wfk0_hdr, comm) !,gs_occ)
217 12 : call wfk0_hdr%vs_dtset(dtset)
218 12 : nspinor = dtset%nspinor
219 :
220 : ! Get fermie from the GS calculation.
221 : ! NB: It might understimate the real Fermi level, especially if the den was computed on a shifted k-mesh
222 : ! at present it's only used to implement pseudobands
223 12 : gs_fermie = wfk0_hdr%fermie
224 :
225 12 : cryst = wfk0_hdr%get_crystal()
226 12 : call cryst%print(header="Crystal structure from WFK file")
227 :
228 : ! Compare structure with the one computed from input file.
229 12 : cryst_dtset = dtset%get_crystal(1)
230 12 : if (cryst%compare(cryst_dtset, header=" Comparing WFK crystal with crystal from dtset") /= 0) then
231 0 : ABI_ERROR("Crystal structure from WFK and dataser do not agree! Check messages above!")
232 : end if
233 12 : call cryst_dtset%free()
234 :
235 119 : call ebands%from_hdr(wfk0_hdr, maxval(wfk0_hdr%nband), gs_eigen)
236 :
237 : !call ebands%update_occ(spinmagntarget)
238 24 : call ebands%print([std_out], header="Ground state energies", prtvol=dtset%prtvol)
239 12 : ABI_FREE(gs_eigen)
240 :
241 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
242 12 : gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0)
243 :
244 24 : call print_ngfft([std_out], ngfftc, header='Coarse FFT mesh used for the wavefunctions')
245 24 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh used for densities and potentials')
246 :
247 : ! Fake MPI_type for the sequential part.
248 12 : call initmpi_seq(mpi_enreg)
249 12 : call mpi_enreg%distribfft%init_seq('c',ngfftc(2),ngfftc(3),'all')
250 12 : call mpi_enreg%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
251 :
252 : ! ===========================================
253 : ! === Open and read pseudopotential files ===
254 : ! ===========================================
255 12 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm)
256 :
257 : ! ============================
258 : ! ==== PAW initialization ====
259 : ! ============================
260 12 : if (dtset%usepaw == 1) then
261 1 : call chkpawovlp(cryst%natom,cryst%ntypat,dtset%pawovlp,pawtab,cryst%rmet,cryst%typat,cryst%xred)
262 :
263 1 : cplex_dij=nspinor; cplex=1; ndij=1
264 :
265 5 : ABI_MALLOC(pawrhoij,(cryst%natom))
266 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
267 1 : nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
268 1 : call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,dtset%nsppol,cryst%typat,pawtab=pawtab)
269 :
270 : ! Initialize values for several basic arrays
271 1 : gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
272 :
273 : ! Test if we have to call pawinit
274 1 : call paw_gencond(dtset,gnt_option,"test",call_pawinit)
275 :
276 1 : if (psp_gencond==1 .or. call_pawinit) then
277 0 : call timab(553,1,tsec)
278 0 : gsqcut_shp = two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
279 : call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,dtset%pawlcutd,dtset%pawlmix,&
280 : psps%mpsang,dtset%pawnphi,cryst%nsym,dtset%pawntheta,pawang,Pawrad,&
281 0 : dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
282 0 : call timab(553,2,tsec)
283 :
284 : ! Update internal values
285 0 : call paw_gencond(dtset,gnt_option,"save",call_pawinit)
286 :
287 : else
288 1 : if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij =2
289 1 : if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla=2
290 : end if
291 :
292 2 : psps%n1xccc=MAXVAL(pawtab(1:cryst%ntypat)%usetcore)
293 :
294 : ! Initialize optional flags in pawtab to zero
295 : ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
296 2 : pawtab(:)%has_nabla = 0
297 2 : pawtab(:)%usepawu = 0
298 2 : pawtab(:)%useexexch = 0
299 2 : pawtab(:)%exchmix =zero
300 2 : pawtab(:)%lamb_shielding =zero
301 :
302 1 : call setsym_ylm(cryst%gprimd,pawang%l_max-1,cryst%nsym,dtset%pawprtvol,cryst%rprimd,cryst%symrec,pawang%zarot)
303 :
304 : ! Initialize and compute data for DFT+U
305 : !paw_dmft%use_dmft=dtset%usedmft
306 : !call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
307 : ! .false.,dtset%jpawu,dtset%lexexch,dtset%lpawu,cryst%ntypat,pawang,dtset%pawprtvol,&
308 : ! Pawrad,pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
309 : !ABI_CHECK(paw_dmft%use_dmft==0,"DMFT not available")
310 : !call destroy_sc_dmft(paw_dmft)
311 :
312 1 : if (my_rank == master) call pawtab_print(pawtab, unit=std_out)
313 :
314 : ! Get Pawrhoij from the header of the WFK file.
315 1 : call pawrhoij_copy(wfk0_hdr%pawrhoij,pawrhoij)
316 :
317 : ! Variables/arrays related to the fine FFT grid.
318 5 : ABI_MALLOC(pawfgrtab,(cryst%natom))
319 1 : call pawtab_get_lsize(pawtab,l_size_atm,cryst%natom,cryst%typat)
320 : cplex=1
321 1 : call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat)
322 1 : ABI_FREE(l_size_atm)
323 :
324 2 : usexcnhat=maxval(pawtab(:)%usexcnhat)
325 : ! 0 if Vloc in atomic data is Vbare (Blochl s formulation)
326 : ! 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
327 1 : call wrtout(std_out,sjoin("using usexcnhat= ",itoa(usexcnhat)))
328 : !
329 : ! Identify parts of the rectangular grid where the density has to be calculated ===
330 : !optcut=0; optgr0=dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-dtset%pawstgylm
331 : !if (dtset%xclevel==2 .and. usexcnhat>0) optgr1=dtset%pawstgylm
332 1 : optcut=1; optgr0=1; optgr1=1; optgr2=1; optrad=1
333 :
334 : call nhatgrid(cryst%atindx1,cryst%gmet,cryst%natom,cryst%natom,cryst%nattyp,ngfftf,cryst%ntypat,&
335 1 : optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,cryst%rprimd,cryst%typat,cryst%ucvol,cryst%xred)
336 :
337 3 : call pawfgrtab_print(pawfgrtab,cryst%natom,unit=std_out,prtvol=dtset%pawprtvol)
338 :
339 : !ABI_MALLOC(ks_nhat,(nfftf,dtset%nspden))
340 : !ks_nhat=zero
341 : else
342 11 : ABI_MALLOC(pawfgrtab,(0))
343 : end if !End of PAW Initialization
344 :
345 14 : select case (dtset%wfk_task)
346 :
347 : case (WFK_TASK_FULLBZ, WFK_TASK_OPTICS_FULLBZ)
348 : ! Read wfk0_path and build WFK in full BZ.
349 2 : if (my_rank == master) then
350 2 : outwfk_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) outwfk_path = nctk_ncify(outwfk_path)
351 2 : call wfk_to_bz(wfk0_path, dtset, psps, pawtab, outwfk_path, hdr_bz, ebands_bz)
352 2 : call ebands_bz%free()
353 :
354 : ! Write KB form factors.
355 2 : if (dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then
356 1 : call prtkbff(outwfk_path, hdr_bz, psps, dtset%prtvol)
357 : end if
358 2 : call hdr_bz%free()
359 : end if
360 2 : call xmpi_barrier(comm)
361 :
362 2 : if (dtset%wfk_task == WFK_TASK_OPTICS_FULLBZ) then
363 : ! Calculate the DDK matrix elements from the WFK file in the full BZ.
364 : ! This is needed for computing non-linear properties in optics as symmetries are not
365 : ! implemented correctly.
366 : ds%only_diago = .False.
367 0 : call ds%compute_ddk(outwfk_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
368 0 : call ds%free()
369 : end if
370 :
371 : case (WFK_TASK_KPTS_ERANGE)
372 2 : call sigtk_kpts_in_erange(dtset, cryst, ebands, psps, pawtab, dtfil%filnam_ds(4), comm)
373 :
374 : case (WFK_TASK_DDK, WFK_TASK_DDK_DIAGO)
375 : ! Calculate the DDK matrix elements from the WFK file
376 3 : ds%only_diago = .False.; if (dtset%wfk_task == WFK_TASK_DDK_DIAGO) ds%only_diago = .True.
377 3 : call ds%compute_ddk(wfk0_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
378 3 : call ds%free()
379 :
380 : case (WFK_TASK_EINTERP)
381 : ! Band structure interpolation from eigenvalues computed on the k-mesh.
382 0 : call ebands%interpolate_kpath(dtset, cryst, [0, 0], dtfil%filnam_ds(4), comm)
383 :
384 : case (WFK_TASK_CHECK_SYMTAB)
385 0 : call wfk_check_symtab(wfk0_path, comm)
386 :
387 : case (WFK_TASK_CLASSIFY)
388 : ! Band classification.
389 :
390 : ! New version
391 : ! Compute the mixing matrices D^{k}(S) from the wavefunctions stored in wfd_t.
392 0 : do spin=1,dtset%nsppol
393 0 : brange_spin(:,spin) = [1, dtset%mband]
394 : end do
395 0 : call dmats%init(wfk0_path, dtset, cryst, brange_spin, ngfftf, pawtab, psps, comm)
396 :
397 0 : if (my_rank == master) then
398 0 : call dmats%check([std_out], dtset%prtvol)
399 0 : call dmats%classify(dtset%prtvol)
400 :
401 : ! Independently validate the group-conjugation D-matrix reconstruction (dmats%check_star)
402 : ! by testing EVERY k-point in the full BZ mesh, not just the IBZ points dmats was built
403 : ! from: each full-BZ k-point is the symmetry-star image of some IBZ k-point, so this
404 : ! exercises dmats_get_star_dmats's multable/toinv composition logic (and, in particular,
405 : ! its still-unverified two-step analytic phase formula) across the whole mesh.
406 : call kpts_ibz_from_kptrlatt(cryst, ebands%kptrlatt, ebands%kptopt, ebands%nshiftk, ebands%shiftk, &
407 0 : nkibz_full, kibz_full, wtk_full, nkbz_full, kbz_full)
408 :
409 0 : nstar_fail = 0
410 0 : do spin=1,dtset%nsppol
411 0 : do ikbz=1,nkbz_full
412 0 : call dmats%check_star(spin, kbz_full(:,ikbz), [std_out], dtset%prtvol, ierr)
413 0 : if (ierr /= 0) then
414 0 : nstar_fail = nstar_fail + 1
415 : call wrtout(units, sjoin("check_star FAILED for spin:", itoa(spin), &
416 0 : ", kbz:", ktoa(kbz_full(:,ikbz)), ", ierr:", itoa(ierr)))
417 : end if
418 : end do
419 : end do
420 : call wrtout(units, sjoin("check_star: tested", itoa(nkbz_full * dtset%nsppol), &
421 0 : "(k,spin) points in the full BZ, failures:", itoa(nstar_fail)))
422 :
423 0 : ABI_FREE(kibz_full)
424 0 : ABI_FREE(wtk_full)
425 0 : ABI_FREE(kbz_full)
426 : end if
427 0 : call dmats%free()
428 :
429 : #if 0
430 : ! old implementation
431 : call read_wfd()
432 :
433 : ABI_MALLOC(esymm,(wfd%nkibz,wfd%nsppol))
434 : use_paw_aeur=.False. ! should pass ngfftf but the dense mesh is not forced to be symmetric
435 :
436 : do spin=1,wfd%nsppol
437 : do ik_ibz=1,wfd%nkibz
438 : first_band = 1
439 : last_band = wfd%nband(ik_ibz,spin)
440 : call classify_bands(wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,wfd%ngfft,&
441 : cryst,ebands,pawtab,pawrad,pawang,psps,dtset%tolsym,esymm(ik_ibz,spin))
442 : end do
443 : end do
444 :
445 : call esymm_free(esymm)
446 : ABI_FREE(esymm)
447 : #endif
448 :
449 : !case (WFK_TASK_UR)
450 : ! ! plot KSS wavefunctions. Change bks_mask to select particular states.
451 : ! ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol))
452 : ! bks_mask=.False.; bks_mask(1:4,1,1)=.True.
453 : ! call wfd%plot_ur(Cryst,Psps,Pawtab,Pawrad,ngfftf,bks_mask)
454 : ! ABI_FREE(bks_mask)
455 :
456 : case (WFK_TASK_PSEUDOBANDS)
457 0 : if (my_rank /= master) goto 100 ! NO MPI parallelism here
458 :
459 : ! out_hdr is the header the STO_WFK file.
460 0 : call wfk0_hdr%copy(out_hdr)
461 :
462 : ! Pre-compute slices for all k-points and spin so that we know the new number of bands in STO_WFK file.
463 0 : ABI_MALLOC(psb_ks, (ebands%nkpt, ebands%nsppol))
464 0 : do spin=1,ebands%nsppol
465 0 : do ik_ibz=1,ebands%nkpt
466 0 : associate (psb => psb_ks(ik_ibz, spin))
467 0 : nband_k = ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
468 0 : call psb%init(dtset, nband_k, ebands%eig(:, ik_ibz, spin), gs_fermie)
469 : ! Change the number of bands to account for pseudo bands.
470 : !print *, "nb_tot:", psb%nb_tot
471 0 : out_hdr%nband(ik_ibz + (spin-1)*ebands%nkpt) = psb%nb_tot
472 : end associate
473 : end do
474 : end do
475 :
476 : ! Compute new value of bantot
477 : ! TODO: Have to change all arrays in outhdr_hdr depending on nband_ks
478 0 : out_hdr%bantot = sum(out_hdr%nband)
479 0 : out_hdr%mband = maxval(out_hdr%nband)
480 0 : ABI_RECALLOC(out_hdr%occ, (out_hdr%bantot))
481 :
482 0 : outwfk_path = strcat(dtfil%filnam_ds(4), "_STO_WFK")
483 0 : if (dtset%iomode == IO_MODE_ETSF) outwfk_path = nctk_ncify(outwfk_path)
484 0 : call out_wfk%open_write(out_hdr, outwfk_path, formeig0, dtset%iomode, get_unit(), xmpi_comm_self)
485 :
486 0 : call in_wfk%open_read(wfk0_path, formeig0, iomode_from_fname(wfk0_path), get_unit(), xmpi_comm_self)
487 :
488 : ! The output arrays eig_k and occ_k contain the *full* set of eigenvalues and occupation
489 : ! factors stored in the file and are dimensioned with mband.
490 0 : ABI_MALLOC(eig_k, (wfk0_hdr%mband))
491 0 : ABI_MALLOC(occ_k, (wfk0_hdr%mband))
492 0 : mpw = maxval(wfk0_hdr%npwarr)
493 0 : ABI_MALLOC(kg_k, (3, mpw))
494 : !print *, "mpw:", mpw
495 0 : sc_mode = xmpio_single
496 :
497 0 : do spin=1,ebands%nsppol
498 0 : do ik_ibz=1,ebands%nkpt
499 0 : associate (psb => psb_ks(ik_ibz, spin))
500 : ! Read and write protected states.
501 0 : band_block = [1, psb%nb_protected]
502 0 : nb = band_block(2) - band_block(1) + 1
503 0 : npw_k = wfk0_hdr%npwarr(ik_ibz)
504 0 : mcg = npw_k * nspinor * nb
505 0 : ABI_MALLOC(cg_k, (2, mcg))
506 :
507 : call in_wfk%read_band_block(band_block, ik_ibz, spin, sc_mode, &
508 0 : kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
509 :
510 : !call wrtout(std_out, sjoin(" About to write islice:", itoa(0), "with band block:", ltoa(band_block)))
511 0 : eig_k(1:psb%nb_tot) = psb%ps_eig(:) ! Change eigenvalues to account for pseudo bands
512 : !occ_k = ???
513 : call out_wfk%write_band_block(band_block, ik_ibz, spin, sc_mode, &
514 0 : kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
515 0 : ABI_FREE(cg_k)
516 :
517 : ! ========================================
518 : ! Build pseudobands and write them to disk
519 : ! ========================================
520 0 : bstart = psb%nb_protected + 1
521 0 : do islice=1,psb%nslices
522 0 : band_block = psb%subspace(1:2, islice)
523 0 : nb = band_block(2) - band_block(1) + 1
524 0 : mcg = npw_k * nspinor * nb
525 0 : ABI_MALLOC(cg_k, (2, mcg))
526 :
527 0 : call in_wfk%read_band_block(band_block, ik_ibz, spin, sc_mode, cg_k=cg_k)
528 0 : call c_f_pointer(c_loc(cg_k), cg_k_cplx, [npw_k*nspinor, nb])
529 :
530 : ! Allocate pseudobands.
531 0 : nsto = psb%subspace(3, islice)
532 0 : ABI_CALLOC(ps_ug, (npw_k*nspinor, nsto))
533 0 : ABI_MALLOC(thetas, (nb))
534 :
535 0 : if (nsto == 1) then
536 : ! Use KS state.
537 0 : ps_ug = cg_k_cplx
538 : else
539 : ! Multiply by random phases.
540 0 : do ii=1,nsto
541 0 : call random_number(thetas)
542 0 : do ib=1,nb
543 0 : ps_ug(:,ii) = ps_ug(:,ii) + cg_k_cplx(:,ib) * exp(j_dpc*two_pi*thetas(ib)) / sqrt(one * nsto)
544 : end do
545 : end do
546 : end if
547 :
548 0 : band_block = [bstart, bstart + psb%subspace(3, islice) - 1]
549 : !call wrtout(std_out, sjoin(" About to write islice:", itoa(islice), "with band block:", ltoa(band_block)))
550 0 : call out_wfk%write_band_block(band_block, ik_ibz, spin, sc_mode, cg_k=cg_k)
551 : !kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
552 0 : bstart = bstart + psb%subspace(3, islice)
553 :
554 0 : ABI_FREE(thetas)
555 0 : ABI_FREE(cg_k)
556 0 : ABI_FREE(ps_ug)
557 : end do ! islice
558 :
559 0 : call psb%free()
560 : end associate
561 : end do ! ik_ibz
562 : end do ! spin
563 :
564 0 : ABI_FREE(eig_k)
565 0 : ABI_FREE(occ_k)
566 0 : ABI_FREE(kg_k)
567 0 : call in_wfk%close(); call out_wfk%close(); call out_hdr%free()
568 :
569 : ! DEBUG section. Try to read the output WFK file.
570 : !call out_wfk%open_read(outwfk_path, formeig0, iomode_from_fname(outwfk_path), get_unit(), xmpi_comm_self)
571 : !do spin=1,ebands%nsppol
572 : ! do ik_ibz=1,ebands%nkpt
573 : ! nband_k = out_wfk%hdr%nband(ik_ibz + (spin-1)*ebands%nkpt)
574 : ! band_block = [1, nband_k]
575 : ! call out_wfk%read_band_block(band_block, ik_ibz, spin, sc_mode, &
576 : ! kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
577 : ! end do
578 : !end do
579 : !call out_wfk%close()
580 :
581 : case (WFK_TASK_PAW_AEPSI)
582 : ! Compute AE PAW wavefunction in real space on the dense FFT mesh.
583 0 : call read_wfd()
584 :
585 0 : ABI_CHECK(wfd%usepaw == 1, "Not a PAW run")
586 0 : ABI_MALLOC(paw_onsite, (cryst%natom))
587 : call paw_pwaves_lmn_init(paw_onsite,cryst%natom,cryst%natom,cryst%ntypat, &
588 0 : cryst%rprimd,cryst%xcart,pawtab,pawrad,pawfgrtab)
589 :
590 : ! Use dense FFT mesh
591 0 : call wfd%change_ngfft(cryst,psps,ngfftf)
592 0 : band = 1; spin = 1; ik_ibz = 1
593 :
594 0 : ABI_MALLOC(ur_ae, (wfd%nfft*wfd%nspinor))
595 0 : call wfd%paw_get_aeur(band,ik_ibz,spin,cryst,paw_onsite,psps,pawtab,pawfgrtab,ur_ae)
596 0 : ABI_FREE(ur_ae)
597 :
598 0 : call paw_pwaves_lmn_free(paw_onsite)
599 0 : ABI_FREE(paw_onsite)
600 :
601 : case (WFK_TASK_WANNIER)
602 : ! Construct Wannier functions.
603 :
604 : ! This part was implemented by gmatteo to debug GaAs with a 8x8x8 k-mesh.
605 :
606 : !if (wfk0_hdr%kptopt == 1) then
607 : ! ! Generate WFK in the full BZ (only master works here)
608 : ! outwfk_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) outwfk_path = nctk_ncify(outwfk_path)
609 : ! if (my_rank == master) then
610 : ! call wrtout(units, sjoin("- Generating WFK file with kpoints in the full BZ and istwfk == 1", outwfk_path))
611 : ! call wfk_to_bz(wfk0_path, dtset, psps, pawtab, outwfk_path, hdr_bz, ebands_bz)
612 : ! call ebands_bz%free(); call hdr_bz%free()
613 : ! end if
614 : ! call xmpi_barrier(comm)
615 : ! my_dtset = dtset%copy()
616 : ! ebands_bz = wfk_read_ebands(outwfk_path, comm, hdr_bz)
617 : ! call hdr_transfer_nkpt_arrays(hdr_bz, my_dtset)
618 : ! my_dtset%kptopt = hdr_bz%kptopt
619 : ! call hdr_bz%vs_dtset(my_dtset)
620 : ! call wfd_run_wannier__(outwfk_path, my_dtset, ebands_bz, hdr_bz)
621 : ! call ebands_bz%free(); call hdr_bz%free(); call my_dtset%free()
622 :
623 : !else
624 5 : call wfk0_hdr%vs_dtset(dtset)
625 5 : call wfd_run_wannier__(wfk0_path, dtset, ebands, wfk0_hdr)
626 : !end if
627 :
628 : case default
629 12 : ABI_ERROR(sjoin("Wrong wfk_task:", itoa(dtset%wfk_task)))
630 : end select
631 :
632 : 100 continue
633 :
634 : ! Free memory
635 12 : call cryst%free(); call ebands%free(); call wfd%free(); call destroy_mpi_enreg(mpi_enreg); call wfk0_hdr%free()
636 12 : call pawfgr_destroy(pawfgr)
637 :
638 : ! Deallocation for PAW.
639 12 : if (dtset%usepaw==1) then
640 1 : call pawrhoij_free(pawrhoij)
641 3 : ABI_FREE(pawrhoij)
642 1 : call pawfgrtab_free(pawfgrtab)
643 : !call paw_ij_free(paw_ij)
644 : !ABI_FREE(paw_ij)
645 : !call paw_an_free(paw_an)
646 : !ABI_FREE(paw_an)
647 : end if
648 26 : ABI_FREE(pawfgrtab)
649 :
650 : DBG_EXIT('COLL')
651 :
652 : contains
653 : !!***
654 :
655 : !!****f* wfk_analyze/read_wfd
656 : !! NAME
657 : !! read_wfd
658 : !!
659 : !! FUNCTION
660 : !! Initialize the wavefunction descriptor from file.
661 : !!
662 : !! SOURCE
663 :
664 0 : subroutine read_wfd()
665 :
666 0 : ABI_MALLOC(keep_ur, (ebands%mband, ebands%nkpt, ebands%nsppol))
667 0 : ABI_MALLOC(bks_mask, (ebands%mband, ebands%nkpt, ebands%nsppol))
668 0 : keep_ur = .False.; bks_mask = .True.
669 :
670 : call wfd%init(cryst,pawtab,psps,keep_ur,ebands%mband,ebands%nband,ebands%nkpt,dtset%nsppol,bks_mask,&
671 : dtset%nspden,dtset%nspinor,ecut_eff,dtset%ecutsm,dtset%dilatmx,wfk0_hdr%istwfk,ebands%kptns,ngfftc,&
672 0 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
673 :
674 0 : ABI_FREE(keep_ur)
675 0 : ABI_FREE(bks_mask)
676 :
677 0 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
678 :
679 0 : end subroutine read_wfd
680 :
681 5 : subroutine wfd_run_wannier__(wfk_filepath, dtset_, ebands_, hdr_)
682 :
683 : type(dataset_type),intent(in) :: dtset_
684 : character(len=*),intent(in) :: wfk_filepath
685 : type(ebands_t),intent(in) :: ebands_
686 : type(hdr_type),intent(in) :: hdr_
687 :
688 25 : ABI_MALLOC(keep_ur, (ebands_%mband, ebands_%nkpt, ebands_%nsppol))
689 20 : ABI_MALLOC(bks_mask, (ebands_%mband, ebands_%nkpt, ebands_%nsppol))
690 941 : keep_ur = .False.; bks_mask = .True.
691 :
692 : ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts)
693 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
694 : !wfk0_hdr%istwfk = 1; ebands%istwfk = 1; dtset%istwfk = 1
695 :
696 : call wfd%init(cryst, pawtab, psps, keep_ur, ebands_%mband, ebands_%nband, ebands_%nkpt, dtset_%nsppol, bks_mask, &
697 : dtset_%nspden, dtset_%nspinor, ecut_eff, dtset_%ecutsm, dtset_%dilatmx, hdr_%istwfk, ebands_%kptns, ngfftc, &
698 5 : dtset_%nloalg, dtset_%prtvol, dtset_%pawprtvol, comm)
699 :
700 5 : ABI_FREE(keep_ur)
701 5 : ABI_FREE(bks_mask)
702 5 : call wfd%read_wfk(wfk_filepath, iomode_from_fname(wfk_filepath))
703 :
704 : call wfd_run_wannier(cryst=cryst, ebands=ebands_, hdr=hdr_, mpi_enreg=mpi_enreg, &
705 : ngfftc=ngfftc, ngfftf=ngfftf, wfd=wfd, dtset=dtset_, dtfil=dtfil, &
706 5 : pawang=pawang, pawrad=pawrad, pawtab=pawtab, psps=psps)
707 :
708 5 : end subroutine wfd_run_wannier__
709 :
710 : end subroutine wfk_analyze
711 : !!***
712 :
713 : subroutine hdr_transfer_nkpt_arrays(hdr, dtset)
714 :
715 : use m_copy, only : alloc_copy
716 :
717 : class(hdr_type),intent(in) :: hdr
718 : type(dataset_type),intent(inout) :: dtset
719 :
720 : ABI_SFREE(dtset%istwfk)
721 : ABI_SFREE(dtset%kpt)
722 : ABI_SFREE(dtset%kptns)
723 : ABI_SFREE(dtset%occ_orig)
724 : ABI_SFREE(dtset%wtk)
725 : ABI_SFREE(dtset%kptns_hf) ! Free HF k-points as well.
726 : ABI_SFREE(dtset%nband)
727 :
728 : dtset%nkpt = hdr%nkpt
729 : call alloc_copy(hdr%istwfk, dtset%istwfk)
730 : call alloc_copy(hdr%nband, dtset%nband)
731 : call alloc_copy(hdr%kptns, dtset%kpt)
732 : call alloc_copy(hdr%kptns, dtset%kptns)
733 : !call alloc_copy(hdr%occ, dtset%occ_orig(:,1)
734 : call alloc_copy(hdr%wtk, dtset%wtk)
735 : call alloc_copy(hdr%kptns, dtset%kptns_hf)
736 :
737 : end subroutine hdr_transfer_nkpt_arrays
738 :
739 : end module m_wfk_analyze
740 : !!***
|
