Branch data 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, ktoa
43 : : use m_fftcore, only : print_ngfft
44 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, init_mpi_enreg
45 : : use m_ddk, only : ddkstore_t
46 : : use m_ksdiago, only : psbands_t
47 : : use m_pawang, only : pawang_type
48 : : use m_pawrad, only : pawrad_type
49 : : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
50 : : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init, pawfgrtab_print
51 : : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_inquire_dim
52 : : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
53 : : use m_paw_sphharm, only : setsym_ylm
54 : : use m_paw_init, only : pawinit, paw_gencond
55 : : use m_paw_nhat, only : nhatgrid
56 : : use m_paw_tools, only : chkpawovlp
57 : : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
58 : : use m_classify_bands, only : dmats_t
59 : : use m_kpts, only : kpts_ibz_from_kptrlatt
60 : : use m_pspini, only : pspini
61 : : use m_sigtk, only : sigtk_kpts_in_erange
62 : : use m_iowf, only : prtkbff
63 : : use m_wfd_wannier, only : wfd_run_wannier
64 : : use m_wfk, only : wfk_to_bz, wfk_t, wfk_read_eigenvalues, wfk_check_symtab
65 : :
66 : : implicit none
67 : :
68 : : private
69 : : !!***
70 : :
71 : : public :: wfk_analyze
72 : : !!***
73 : :
74 : : contains
75 : : !!***
76 : :
77 : : !!****f* ABINIT/wfk_analyze
78 : : !! NAME
79 : : !! wfk_analyze
80 : : !!
81 : : !! FUNCTION
82 : : !! Driver for optdriver 8: read the WFK file given by getwfk/irdwfk (dtfil%fnamewffk) and run
83 : : !! the post-processing tool selected by wfk_task (WFK_TASK_* in defs_basis, parsed by str2wfktask).
84 : : !! Unlike cut3d, the tools support MPI, have no interactive prompt and write results to netcdf files.
85 : : !!
86 : : !! The crystal structure of the WFK must agree with the one of the input file.
87 : : !! Pseudopotentials, PAW data and FFT meshes are initialized for all tasks.
88 : : !!
89 : : !! wfk_task Output MPI Notes
90 : : !! ----------------------------------------------------------------------------------------------------
91 : : !! wfk_fullbz WFK in the full BZ (+KB master only istwfk > 1 converted to full G-sphere,
92 : : !! form factors if prtkbff 1) prtkbff needs iomode 3, NC only
93 : : !! wfk_optics_fullbz WFK in full BZ + EVK files parallel DDK for non-linear optics in optic
94 : : !! wfk_ddk EVK files (3 directions) parallel same k-points as the WFK
95 : : !! wfk_ddk_diago EVK files, diagonal only parallel
96 : : !! wfk_einterp band structure (SKW) comm einterp variables
97 : : !! wfk_kpts_erange KERANGE.nc comm sigma_ngkpt, sigma_erange, einterp
98 : : !! wannier Wannier90 files untested NC only if prtwant /= 0 (chkinp);
99 : : !! mpi_enreg from initmpi_seq (me=0)
100 : : !! classify log only parallel init developer: D-matrices, irreps, check_star
101 : : !! check_symtab log only master only developer: symmetry tables vs WFK
102 : : !! paw_aepsi none - developer stub
103 : : !! pseudobands STO_WFK master only WIP, rejected by chkinp
104 : : !!
105 : : !! INPUTS
106 : : !! acell(3), codvsn, rprim(3,3), xred(3,natom): not used (the structure is read from the WFK file).
107 : : !! dtfil<datafiles_type>=Variables related to files.
108 : : !! dtset<dataset_type>=All input variables for this dataset.
109 : : !! pawang<pawang_type)>=PAW angular mesh and related data.
110 : : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
111 : : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
112 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials, completed by pspini.
113 : : !!
114 : : !! NOTES
115 : : !!
116 : : !! ON THE USE OF FFT GRIDS:
117 : : !! =================
118 : : !! In case of PAW:
119 : : !! ---------------
120 : : !! Two FFT grids are used:
121 : : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
122 : : !! It is defined by nfft, ngfft, mgfft, ...
123 : : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
124 : : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
125 : : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
126 : : !! In case of norm-conserving:
127 : : !! ---------------------------
128 : : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
129 : : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
130 : : !!
131 : : !! SOURCE
132 : :
133 : 16 : subroutine wfk_analyze(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred)
134 : :
135 : : !Arguments ------------------------------------
136 : : !scalars
137 : : character(len=8),intent(in) :: codvsn
138 : : type(datafiles_type),intent(in) :: dtfil
139 : : type(dataset_type),intent(inout) :: dtset
140 : : type(pawang_type),intent(inout) :: pawang
141 : : type(pseudopotential_type),intent(inout) :: psps
142 : : !arrays
143 : : real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom)
144 : : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
145 : : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
146 : :
147 : : !Local variables ------------------------------
148 : : !scalars
149 : : integer,parameter :: master = 0, formeig0 = 0
150 : : integer :: comm,nprocs,my_rank,mgfftf,nfftf !,nfftf_tot
151 : : integer :: optcut,optgr0,optgr1,optgr2,optrad,psp_gencond,ii
152 : : integer :: band,ik_ibz,spin,nband_k, islice, ib, mpw, mcg, nb, npw_k ! first_band,last_band,
153 : : integer :: ierr,usexcnhat, sc_mode, nspinor, nsto
154 : : integer :: cplex,cplex_dij,cplex_rhoij,ndij,nspden_rhoij,gnt_option
155 : : real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff,gsqcut_shp, gs_fermie
156 : : character(len=500) :: msg
157 : : character(len=fnlen) :: wfk0_path, outwfk_path
158 : : logical :: call_pawinit
159 : 48 : type(hdr_type) :: wfk0_hdr, hdr_bz, out_hdr
160 [ + + + + : 1664 : type(crystal_t) :: cryst, cryst_dtset
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + ]
161 : 32 : type(ebands_t) :: ebands, ebands_bz
162 : : type(pawfgr_type) :: pawfgr
163 : 16 : type(mpi_type) :: mpi_enreg
164 : 16 : type(wfd_t) :: wfd
165 : 16 : type(ddkstore_t) :: ds
166 [ + + + + ]: 128 : type(wfk_t) :: in_wfk, out_wfk
167 : 16 : type(dmats_t) :: dmats
168 : : !arrays
169 : 32 : integer :: ngfftc(18),ngfftf(18), units(2), band_block(2), bstart, brange_spin(2, dtset%nsppol)
170 : : integer :: nkibz_full, nkbz_full, ikbz, nstar_fail
171 : 16 : integer,allocatable :: l_size_atm(:), kg_k(:,:)
172 : 32 : real(dp),allocatable :: wtk_full(:), kibz_full(:,:), kbz_full(:,:)
173 : : real(dp),parameter :: k0(3)=zero
174 : 16 : real(dp),pointer :: gs_eigen(:,:,:)
175 : 16 : real(dp),allocatable :: eig_k(:), occ_k(:), thetas(:) !, out_cg(:,:), work(:,:,:,:), allcg_k(:,:)
176 : 16 : real(dp),allocatable,target :: cg_k(:,:)
177 : 16 : complex(gwp),allocatable :: ur_ae(:)
178 : 16 : complex(dp),pointer :: cg_k_cplx(:,:)
179 : 16 : complex(dp),allocatable :: ps_ug(:,:)
180 : 16 : logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:)
181 : : real(dp) :: tsec(2)
182 : 16 : type(Pawrhoij_type),allocatable :: pawrhoij(:)
183 : 16 : type(pawfgrtab_type),allocatable :: pawfgrtab(:)
184 : 16 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
185 : 16 : type(psbands_t),allocatable :: psb_ks(:,:)
186 : : !************************************************************************
187 : :
188 : : DBG_ENTER('COLL')
189 : :
190 : : ! abirules!
191 : : if (.False.) write(std_out,*)acell,codvsn,rprim,xred
192 [ + + ]: 48 : units = [std_out, ab_out]
193 : :
194 : 16 : comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
195 : :
196 : 16 : wfk0_path = dtfil%fnamewffk
197 [ + - ]: 16 : if (my_rank == master) then
198 : : ! Accept WFK file in Fortran or netcdf format.
199 [ - + ]: 16 : if (nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then
200 : 0 : ABI_ERROR(sjoin("Cannot find GS WFK file:", ch10, msg))
201 : : end if
202 : : end if
203 : 16 : call xmpi_bcast(wfk0_path, master, comm, ierr)
204 : 16 : call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path))
205 : :
206 : : !call cwtime(cpu,wall,gflops,"start")
207 : :
208 : : ! Construct crystal and ebands from the GS WFK file.
209 : 16 : call wfk_read_eigenvalues(wfk0_path, gs_eigen, wfk0_hdr, comm) !,gs_occ)
210 : 16 : call wfk0_hdr%vs_dtset(dtset)
211 : 16 : nspinor = dtset%nspinor
212 : :
213 : : ! Get fermie from the GS calculation.
214 : : ! NB: It might understimate the real Fermi level, especially if the den was computed on a shifted k-mesh
215 : : ! at present it's only used to implement pseudobands
216 : 16 : gs_fermie = wfk0_hdr%fermie
217 : :
218 [ - + - + : 16 : cryst = wfk0_hdr%get_crystal()
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
219 : 16 : call cryst%print(header="Crystal structure from WFK file")
220 : :
221 : : ! Compare structure with the one computed from input file.
222 : 16 : cryst_dtset = dtset%get_crystal(1)
223 [ - + ]: 16 : if (cryst%compare(cryst_dtset, header=" Comparing WFK crystal with crystal from dtset") /= 0) then
224 : 0 : ABI_ERROR("Crystal structure from WFK and dataset do not agree! Check messages above!")
225 : : end if
226 : 16 : call cryst_dtset%free()
227 : :
228 [ + + - + ]: 211 : call ebands%from_hdr(wfk0_hdr, maxval(wfk0_hdr%nband), gs_eigen)
229 : :
230 : : !call ebands%update_occ(spinmagntarget)
231 [ + + ]: 32 : call ebands%print([std_out], header="Ground state energies", prtvol=dtset%prtvol)
232 [ - + ]: 16 : ABI_FREE(gs_eigen)
233 : :
234 : : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
235 : 16 : gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0)
236 : :
237 [ + + ]: 32 : call print_ngfft([std_out], ngfftc, header='Coarse FFT mesh used for the wavefunctions')
238 [ + + ]: 32 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh used for densities and potentials')
239 : :
240 : : ! Fake MPI_type for the sequential part.
241 : 16 : call initmpi_seq(mpi_enreg)
242 : 16 : call mpi_enreg%distribfft%init_seq('c',ngfftc(2),ngfftc(3),'all')
243 : 16 : call mpi_enreg%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
244 : :
245 : : ! ===========================================
246 : : ! === Open and read pseudopotential files ===
247 : : ! ===========================================
248 : 16 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm)
249 : :
250 : : ! ============================
251 : : ! ==== PAW initialization ====
252 : : ! ============================
253 [ + + ]: 16 : if (dtset%usepaw == 1) then
254 : 1 : call chkpawovlp(cryst%natom,cryst%ntypat,dtset%pawovlp,pawtab,cryst%rmet,cryst%typat,cryst%xred)
255 : :
256 : 1 : cplex_dij=nspinor; cplex=1; ndij=1
257 : :
258 [ - + + - : 5 : ABI_MALLOC(pawrhoij,(cryst%natom))
+ - + - -
+ - + - +
+ + ]
259 : : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
260 : 1 : nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
261 : 1 : call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,dtset%nsppol,cryst%typat,pawtab=pawtab)
262 : :
263 : : ! Initialize values for several basic arrays
264 [ + - + - : 1 : gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
- + ]
265 : :
266 : : ! Test if we have to call pawinit
267 : 1 : call paw_gencond(dtset,gnt_option,"test",call_pawinit)
268 : :
269 [ + - - + ]: 1 : if (psp_gencond==1 .or. call_pawinit) then
270 : 0 : call timab(553,1,tsec)
271 : 0 : gsqcut_shp = two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
272 : : call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,dtset%pawlcutd,dtset%pawlmix,&
273 : : psps%mpsang,dtset%pawnphi,cryst%nsym,dtset%pawntheta,pawang,Pawrad,&
274 : 0 : dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
275 : 0 : call timab(553,2,tsec)
276 : :
277 : : ! Update internal values
278 : 0 : call paw_gencond(dtset,gnt_option,"save",call_pawinit)
279 : :
280 : : else
281 [ - + - - ]: 1 : if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij =2
282 [ - + - - ]: 1 : if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla=2
283 : : end if
284 : :
285 [ + + ]: 2 : psps%n1xccc=MAXVAL(pawtab(1:cryst%ntypat)%usetcore)
286 : :
287 : : ! Initialize optional flags in pawtab to zero
288 : : ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
289 [ + + ]: 2 : pawtab(:)%has_nabla = 0
290 [ + + ]: 2 : pawtab(:)%usepawu = 0
291 [ + + ]: 2 : pawtab(:)%useexexch = 0
292 [ + + ]: 2 : pawtab(:)%exchmix =zero
293 [ + + ]: 2 : pawtab(:)%lamb_shielding =zero
294 : :
295 : 1 : call setsym_ylm(cryst%gprimd,pawang%l_max-1,cryst%nsym,dtset%pawprtvol,cryst%rprimd,cryst%symrec,pawang%zarot)
296 : :
297 : : ! Initialize and compute data for DFT+U
298 : : !paw_dmft%use_dmft=dtset%usedmft
299 : : !call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
300 : : ! .false.,dtset%jpawu,dtset%lexexch,dtset%lpawu,cryst%ntypat,pawang,dtset%pawprtvol,&
301 : : ! Pawrad,pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
302 : : !ABI_CHECK(paw_dmft%use_dmft==0,"DMFT not available")
303 : : !call destroy_sc_dmft(paw_dmft)
304 : :
305 [ + - ]: 1 : if (my_rank == master) call pawtab_print(pawtab, unit=std_out)
306 : :
307 : : ! Get Pawrhoij from the header of the WFK file.
308 : 1 : call pawrhoij_copy(wfk0_hdr%pawrhoij,pawrhoij)
309 : :
310 : : ! Variables/arrays related to the fine FFT grid.
311 [ - + + - : 5 : ABI_MALLOC(pawfgrtab,(cryst%natom))
+ - + - -
+ - + +
+ ]
312 : 1 : call pawtab_get_lsize(pawtab,l_size_atm,cryst%natom,cryst%typat)
313 : : cplex=1
314 : 1 : call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat)
315 [ - + ]: 1 : ABI_FREE(l_size_atm)
316 : :
317 [ + + ]: 2 : usexcnhat=maxval(pawtab(:)%usexcnhat)
318 : : ! 0 if Vloc in atomic data is Vbare (Blochl s formulation)
319 : : ! 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
320 : 1 : call wrtout(std_out,sjoin("using usexcnhat= ",itoa(usexcnhat)))
321 : : !
322 : : ! Identify parts of the rectangular grid where the density has to be calculated ===
323 : : !optcut=0; optgr0=dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-dtset%pawstgylm
324 : : !if (dtset%xclevel==2 .and. usexcnhat>0) optgr1=dtset%pawstgylm
325 : 1 : optcut=1; optgr0=1; optgr1=1; optgr2=1; optrad=1
326 : :
327 : : call nhatgrid(cryst%atindx1,cryst%gmet,cryst%natom,cryst%natom,cryst%nattyp,ngfftf,cryst%ntypat,&
328 : 1 : optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,cryst%rprimd,cryst%typat,cryst%ucvol,cryst%xred)
329 : :
330 : 3 : call pawfgrtab_print(pawfgrtab,cryst%natom,unit=std_out,prtvol=dtset%pawprtvol)
331 : :
332 : : !ABI_MALLOC(ks_nhat,(nfftf,dtset%nspden))
333 : : !ks_nhat=zero
334 : : else
335 [ - + ]: 15 : ABI_MALLOC(pawfgrtab,(0))
336 : : end if !End of PAW Initialization
337 : :
338 : 19 : select case (dtset%wfk_task)
339 : :
340 : : case (WFK_TASK_FULLBZ, WFK_TASK_OPTICS_FULLBZ)
341 : : ! Read wfk0_path and build WFK in full BZ (only master writes the file).
342 : : ! outwfk_path must be defined on all procs as compute_ddk reads it with comm.
343 [ + - ]: 3 : outwfk_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) outwfk_path = nctk_ncify(outwfk_path)
344 [ + - ]: 3 : if (my_rank == master) then
345 : 3 : call wfk_to_bz(wfk0_path, dtset, psps, pawtab, outwfk_path, hdr_bz, ebands_bz)
346 : 3 : call ebands_bz%free()
347 : :
348 : : ! Write KB form factors.
349 [ + + + - : 3 : if (dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then
+ - ]
350 : 1 : call prtkbff(outwfk_path, hdr_bz, psps, dtset%prtvol)
351 : : end if
352 : 3 : call hdr_bz%free()
353 : : end if
354 : 3 : call xmpi_barrier(comm)
355 : :
356 [ + + ]: 3 : if (dtset%wfk_task == WFK_TASK_OPTICS_FULLBZ) then
357 : : ! Calculate the DDK matrix elements from the WFK file in the full BZ.
358 : : ! This is needed for computing non-linear properties in optics as symmetries are not
359 : : ! implemented correctly.
360 : : ds%only_diago = .False.
361 : 1 : call ds%compute_ddk(outwfk_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
362 : 1 : call ds%free()
363 : : end if
364 : :
365 : : case (WFK_TASK_KPTS_ERANGE)
366 : 2 : call sigtk_kpts_in_erange(dtset, cryst, ebands, psps, pawtab, dtfil%filnam_ds(4), comm)
367 : :
368 : : case (WFK_TASK_DDK, WFK_TASK_DDK_DIAGO)
369 : : ! Calculate the DDK matrix elements from the WFK file
370 [ + + ]: 4 : ds%only_diago = .False.; if (dtset%wfk_task == WFK_TASK_DDK_DIAGO) ds%only_diago = .True.
371 : 4 : call ds%compute_ddk(wfk0_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
372 : 4 : call ds%free()
373 : :
374 : : case (WFK_TASK_EINTERP)
375 : : ! Band structure interpolation from eigenvalues computed on the k-mesh.
376 : 1 : call ebands%interpolate_kpath(dtset, cryst, [0, 0], dtfil%filnam_ds(4), comm)
377 : :
378 : : case (WFK_TASK_CHECK_SYMTAB)
379 : 1 : call wfk_check_symtab(wfk0_path, comm)
380 : :
381 : : case (WFK_TASK_CLASSIFY)
382 : : ! Band classification.
383 : :
384 : : ! New version
385 : : ! Compute the mixing matrices D^{k}(S) from the wavefunctions stored in wfd_t.
386 [ # # ]: 0 : do spin=1,dtset%nsppol
387 [ # # ]: 0 : brange_spin(:,spin) = [1, dtset%mband]
388 : : end do
389 : 0 : call dmats%init(wfk0_path, dtset, cryst, brange_spin, ngfftf, pawtab, psps, comm)
390 : :
391 [ # # ]: 0 : if (my_rank == master) then
392 [ # # ]: 0 : call dmats%check([std_out], dtset%prtvol)
393 : 0 : call dmats%classify(dtset%prtvol)
394 : :
395 : : ! Self-check of the D-matrices of the star (dmats%check_star) for every k-point in the full BZ,
396 : : ! not only the IBZ points used to build dmats: this tests the group-multiplication logic
397 : : ! of dmats_get_star_dmats (multable, toinv, phases) on the whole mesh.
398 : : call kpts_ibz_from_kptrlatt(cryst, ebands%kptrlatt, ebands%kptopt, ebands%nshiftk, ebands%shiftk, &
399 : 0 : nkibz_full, kibz_full, wtk_full, nkbz_full, kbz_full)
400 : :
401 : 0 : nstar_fail = 0
402 [ # # ]: 0 : do spin=1,dtset%nsppol
403 [ # # ]: 0 : do ikbz=1,nkbz_full
404 [ # # ]: 0 : call dmats%check_star(spin, kbz_full(:,ikbz), [std_out], dtset%prtvol, ierr)
405 [ # # ]: 0 : if (ierr /= 0) then
406 : 0 : nstar_fail = nstar_fail + 1
407 : : call wrtout(units, sjoin("check_star FAILED for spin:", itoa(spin), &
408 : 0 : ", kbz:", ktoa(kbz_full(:,ikbz)), ", ierr:", itoa(ierr)))
409 : : end if
410 : : end do
411 : : end do
412 : : call wrtout(units, sjoin("check_star: tested", itoa(nkbz_full * dtset%nsppol), &
413 : 0 : "(k,spin) points in the full BZ, failures:", itoa(nstar_fail)))
414 : :
415 [ # # ]: 0 : ABI_FREE(kibz_full)
416 [ # # ]: 0 : ABI_FREE(wtk_full)
417 [ # # ]: 0 : ABI_FREE(kbz_full)
418 : : end if
419 : 0 : call dmats%free()
420 : :
421 : : case (WFK_TASK_PSEUDOBANDS)
422 : : ! WIP: rejected by chkinp. psbands_init is not finished and the pseudobands (ps_ug) are not written yet.
423 [ # # ]: 0 : if (my_rank /= master) goto 100 ! NO MPI parallelism here
424 : :
425 : : ! out_hdr is the header the STO_WFK file.
426 : 0 : call wfk0_hdr%copy(out_hdr)
427 : :
428 : : ! Pre-compute slices for all k-points and spin so that we know the new number of bands in STO_WFK file.
429 [ # # # # : 0 : ABI_MALLOC(psb_ks, (ebands%nkpt, ebands%nsppol))
# # # # #
# # # # #
# # # # #
# ]
430 [ # # ]: 0 : do spin=1,ebands%nsppol
431 [ # # ]: 0 : do ik_ibz=1,ebands%nkpt
432 : 0 : associate (psb => psb_ks(ik_ibz, spin))
433 : 0 : nband_k = ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
434 : 0 : call psb%init(dtset, nband_k, ebands%eig(:, ik_ibz, spin), gs_fermie)
435 : : ! Change the number of bands to account for pseudo bands.
436 : : !print *, "nb_tot:", psb%nb_tot
437 : 0 : out_hdr%nband(ik_ibz + (spin-1)*ebands%nkpt) = psb%nb_tot
438 : : end associate
439 : : end do
440 : : end do
441 : :
442 : : ! Compute new value of bantot
443 : : ! TODO: Have to change all arrays in outhdr_hdr depending on nband_ks
444 [ # # ]: 0 : out_hdr%bantot = sum(out_hdr%nband)
445 [ # # ]: 0 : out_hdr%mband = maxval(out_hdr%nband)
446 [ # # # # : 0 : ABI_RECALLOC(out_hdr%occ, (out_hdr%bantot))
# # # # #
# # # # #
# # # # ]
447 : :
448 [ # # ]: 0 : outwfk_path = strcat(dtfil%filnam_ds(4), "_STO_WFK")
449 [ # # ]: 0 : if (dtset%iomode == IO_MODE_ETSF) outwfk_path = nctk_ncify(outwfk_path)
450 : 0 : call out_wfk%open_write(out_hdr, outwfk_path, formeig0, dtset%iomode, get_unit(), xmpi_comm_self)
451 : :
452 : 0 : call in_wfk%open_read(wfk0_path, formeig0, iomode_from_fname(wfk0_path), get_unit(), xmpi_comm_self)
453 : :
454 : : ! The output arrays eig_k and occ_k contain the *full* set of eigenvalues and occupation
455 : : ! factors stored in the file and are dimensioned with mband.
456 [ # # # # : 0 : ABI_MALLOC(eig_k, (wfk0_hdr%mband))
# # # # #
# # # ]
457 [ # # # # : 0 : ABI_MALLOC(occ_k, (wfk0_hdr%mband))
# # # # ]
458 [ # # ]: 0 : mpw = maxval(wfk0_hdr%npwarr)
459 [ # # # # : 0 : ABI_MALLOC(kg_k, (3, mpw))
# # # # #
# # # ]
460 : : !print *, "mpw:", mpw
461 : 0 : sc_mode = xmpio_single
462 : :
463 [ # # ]: 0 : do spin=1,ebands%nsppol
464 [ # # ]: 0 : do ik_ibz=1,ebands%nkpt
465 : 0 : associate (psb => psb_ks(ik_ibz, spin))
466 : : ! Read and write protected states.
467 [ # # ]: 0 : band_block = [1, psb%nb_protected]
468 : 0 : nb = band_block(2) - band_block(1) + 1
469 : 0 : npw_k = wfk0_hdr%npwarr(ik_ibz)
470 : 0 : mcg = npw_k * nspinor * nb
471 [ # # # # : 0 : ABI_MALLOC(cg_k, (2, mcg))
# # # # #
# # # #
# ]
472 : :
473 : : call in_wfk%read_band_block(band_block, ik_ibz, spin, sc_mode, &
474 : 0 : kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
475 : :
476 : : !call wrtout(std_out, sjoin(" About to write islice:", itoa(0), "with band block:", ltoa(band_block)))
477 [ # # ]: 0 : eig_k(1:psb%nb_tot) = psb%ps_eig(:) ! Change eigenvalues to account for pseudo bands
478 : : !occ_k = ???
479 : : call out_wfk%write_band_block(band_block, ik_ibz, spin, sc_mode, &
480 : 0 : kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
481 [ # # ]: 0 : ABI_FREE(cg_k)
482 : :
483 : : ! ========================================
484 : : ! Build pseudobands and write them to disk
485 : : ! ========================================
486 : 0 : bstart = psb%nb_protected + 1
487 [ # # ]: 0 : do islice=1,psb%nslices
488 [ # # ]: 0 : band_block = psb%subspace(1:2, islice)
489 : 0 : nb = band_block(2) - band_block(1) + 1
490 : 0 : mcg = npw_k * nspinor * nb
491 [ # # # # : 0 : ABI_MALLOC(cg_k, (2, mcg))
# # # # #
# # # #
# ]
492 : :
493 : 0 : call in_wfk%read_band_block(band_block, ik_ibz, spin, sc_mode, cg_k=cg_k)
494 [ # # ]: 0 : call c_f_pointer(c_loc(cg_k), cg_k_cplx, [npw_k*nspinor, nb])
495 : :
496 : : ! Allocate pseudobands.
497 : 0 : nsto = psb%subspace(3, islice)
498 [ # # # # : 0 : ABI_CALLOC(ps_ug, (npw_k*nspinor, nsto))
# # # # #
# # # # #
# # # # #
# ]
499 [ # # # # : 0 : ABI_MALLOC(thetas, (nb))
# # # # #
# # # #
# ]
500 : :
501 [ # # ]: 0 : if (nsto == 1) then
502 : : ! Use KS state.
503 [ # # # # : 0 : ps_ug = cg_k_cplx
# # # # #
# # # # #
# # ]
504 : : else
505 : : ! Multiply by random phases.
506 [ # # ]: 0 : do ii=1,nsto
507 : 0 : call random_number(thetas)
508 [ # # ]: 0 : do ib=1,nb
509 [ # # ]: 0 : ps_ug(:,ii) = ps_ug(:,ii) + cg_k_cplx(:,ib) * exp(j_dpc*two_pi*thetas(ib)) / sqrt(one * nsto)
510 : : end do
511 : : end do
512 : : end if
513 : :
514 [ # # ]: 0 : band_block = [bstart, bstart + psb%subspace(3, islice) - 1]
515 : : !call wrtout(std_out, sjoin(" About to write islice:", itoa(islice), "with band block:", ltoa(band_block)))
516 : 0 : call out_wfk%write_band_block(band_block, ik_ibz, spin, sc_mode, cg_k=cg_k)
517 : : !kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k)
518 : 0 : bstart = bstart + psb%subspace(3, islice)
519 : :
520 [ # # ]: 0 : ABI_FREE(thetas)
521 [ # # ]: 0 : ABI_FREE(cg_k)
522 [ # # ]: 0 : ABI_FREE(ps_ug)
523 : : end do ! islice
524 : :
525 : 0 : call psb%free()
526 : : end associate
527 : : end do ! ik_ibz
528 : : end do ! spin
529 : :
530 [ # # ]: 0 : ABI_FREE(eig_k)
531 [ # # ]: 0 : ABI_FREE(occ_k)
532 [ # # ]: 0 : ABI_FREE(kg_k)
533 [ # # # # : 0 : ABI_FREE(psb_ks)
# # # # #
# ]
534 : 0 : call in_wfk%close(); call out_wfk%close(); call out_hdr%free()
535 : :
536 : : case (WFK_TASK_PAW_AEPSI)
537 : : ! Compute AE PAW wavefunction in real space on the dense FFT mesh.
538 : : ! Developer stub: only (band=1, k=1, spin=1) is computed and the result is not written.
539 : 0 : call read_wfd()
540 : :
541 [ # # ]: 0 : ABI_CHECK(wfd%usepaw == 1, "Not a PAW run")
542 [ # # # # : 0 : ABI_MALLOC(paw_onsite, (cryst%natom))
# # # # #
# # # #
# ]
543 : : call paw_pwaves_lmn_init(paw_onsite,cryst%natom,cryst%natom,cryst%ntypat, &
544 : 0 : cryst%rprimd,cryst%xcart,pawtab,pawrad,pawfgrtab)
545 : :
546 : : ! Use dense FFT mesh
547 : 0 : call wfd%change_ngfft(cryst,psps,ngfftf)
548 : 0 : band = 1; spin = 1; ik_ibz = 1
549 : :
550 [ # # # # : 0 : ABI_MALLOC(ur_ae, (wfd%nfft*wfd%nspinor))
# # # # #
# # # ]
551 : 0 : call wfd%paw_get_aeur(band,ik_ibz,spin,cryst,paw_onsite,psps,pawtab,pawfgrtab,ur_ae)
552 [ # # ]: 0 : ABI_FREE(ur_ae)
553 : :
554 : 0 : call paw_pwaves_lmn_free(paw_onsite)
555 [ # # # # : 0 : ABI_FREE(paw_onsite)
# # # # #
# # # # #
# # ]
556 : :
557 : : case (WFK_TASK_WANNIER)
558 : : ! Construct Wannier functions. States in the full BZ are reconstructed by symmetry from the IBZ.
559 : 5 : call wfd_run_wannier__(wfk0_path, dtset, ebands, wfk0_hdr)
560 : :
561 : : case default
562 [ + + + + : 16 : ABI_ERROR(sjoin("Wrong wfk_task:", itoa(dtset%wfk_task)))
+ - - - +
- ]
563 : : end select
564 : :
565 : : 100 continue
566 : :
567 : : ! Free memory
568 : 16 : call cryst%free(); call ebands%free(); call wfd%free(); call destroy_mpi_enreg(mpi_enreg); call wfk0_hdr%free()
569 : 16 : call pawfgr_destroy(pawfgr)
570 : :
571 : : ! Deallocation for PAW.
572 [ + + ]: 16 : if (dtset%usepaw==1) then
573 : 1 : call pawrhoij_free(pawrhoij)
574 [ + - + + : 3 : ABI_FREE(pawrhoij)
- + - + -
+ - + - +
- + - + ]
575 : 1 : call pawfgrtab_free(pawfgrtab)
576 : : end if
577 [ + - + + : 34 : ABI_FREE(pawfgrtab)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - +
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - + - -
- - - - -
- - - - -
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
578 : :
579 : : DBG_EXIT('COLL')
580 : :
581 : : contains
582 : : !!***
583 : :
584 : : !!****f* wfk_analyze/read_wfd
585 : : !! NAME
586 : : !! read_wfd
587 : : !!
588 : : !! FUNCTION
589 : : !! Initialize the wavefunction descriptor from file.
590 : : !!
591 : : !! SOURCE
592 : :
593 : 0 : subroutine read_wfd()
594 : :
595 [ # # # # : 0 : ABI_MALLOC(keep_ur, (ebands%mband, ebands%nkpt, ebands%nsppol))
# # # # #
# # # # #
# # # # #
# # # ]
596 [ # # # # : 0 : ABI_MALLOC(bks_mask, (ebands%mband, ebands%nkpt, ebands%nsppol))
# # # # #
# # # # #
# # # # ]
597 [ # # # # : 0 : keep_ur = .False.; bks_mask = .True.
# # # # #
# # # ]
598 : :
599 : : call wfd%init(cryst,pawtab,psps,keep_ur,ebands%mband,ebands%nband,ebands%nkpt,dtset%nsppol,bks_mask,&
600 : : dtset%nspden,dtset%nspinor,ecut_eff,dtset%ecutsm,dtset%dilatmx,wfk0_hdr%istwfk,ebands%kptns,ngfftc,&
601 : 0 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
602 : :
603 [ # # ]: 0 : ABI_FREE(keep_ur)
604 [ # # ]: 0 : ABI_FREE(bks_mask)
605 : :
606 : 0 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
607 : :
608 : 0 : end subroutine read_wfd
609 : :
610 : 5 : subroutine wfd_run_wannier__(wfk_filepath, dtset_, ebands_, hdr_)
611 : :
612 : : type(dataset_type),intent(in) :: dtset_
613 : : character(len=*),intent(in) :: wfk_filepath
614 : : type(ebands_t),intent(in) :: ebands_
615 : : type(hdr_type),intent(in) :: hdr_
616 : :
617 [ - + + - : 25 : ABI_MALLOC(keep_ur, (ebands_%mband, ebands_%nkpt, ebands_%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
618 [ - + + - : 20 : ABI_MALLOC(bks_mask, (ebands_%mband, ebands_%nkpt, ebands_%nsppol))
- + + - -
+ + - - +
- + - + ]
619 [ + + + + : 941 : keep_ur = .False.; bks_mask = .True.
+ + + + +
+ + + ]
620 : :
621 : : ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts)
622 : : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
623 : : !wfk0_hdr%istwfk = 1; ebands%istwfk = 1; dtset%istwfk = 1
624 : :
625 : : call wfd%init(cryst, pawtab, psps, keep_ur, ebands_%mband, ebands_%nband, ebands_%nkpt, dtset_%nsppol, bks_mask, &
626 : : dtset_%nspden, dtset_%nspinor, ecut_eff, dtset_%ecutsm, dtset_%dilatmx, hdr_%istwfk, ebands_%kptns, ngfftc, &
627 : 5 : dtset_%nloalg, dtset_%prtvol, dtset_%pawprtvol, comm)
628 : :
629 [ - + ]: 5 : ABI_FREE(keep_ur)
630 [ - + ]: 5 : ABI_FREE(bks_mask)
631 : 5 : call wfd%read_wfk(wfk_filepath, iomode_from_fname(wfk_filepath))
632 : :
633 : : call wfd_run_wannier(cryst=cryst, ebands=ebands_, hdr=hdr_, mpi_enreg=mpi_enreg, &
634 : : ngfftc=ngfftc, ngfftf=ngfftf, wfd=wfd, dtset=dtset_, dtfil=dtfil, &
635 : 5 : pawang=pawang, pawrad=pawrad, pawtab=pawtab, psps=psps)
636 : :
637 : 5 : end subroutine wfd_run_wannier__
638 : :
639 : : end subroutine wfk_analyze
640 : : !!***
641 : :
642 : : end module m_wfk_analyze
643 : : !!***
|
