TABLE OF CONTENTS


ABINIT/m_wfk_analyze [ Modules ]

[ Top ] [ Modules ]

NAME

  m_wfk_analyze

FUNCTION

  Post-processing tools for WFK file

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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 defs_basis
25  use m_abicore
26  use m_xmpi
27  use m_errors
28  use m_hdr
29  use m_crystal
30  use m_ebands
31  use m_nctk
32  use m_wfk
33  use m_wfd
34  use m_dtset
35  use m_dtfil
36  use m_distribfft
37 
38  use defs_datatypes,    only : pseudopotential_type, ebands_t
39  use defs_abitypes,     only : mpi_type
40  use m_time,            only : timab
41  use m_fstrings,        only : strcat, sjoin, itoa, ftoa
42  use m_fftcore,         only : print_ngfft
43  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
44  use m_esymm,           only : esymm_t, esymm_free
45  use m_ddk,             only : ddkstore_t
46  use m_pawang,          only : pawang_type
47  use m_pawrad,          only : pawrad_type
48  use m_pawtab,          only : pawtab_type, pawtab_print, pawtab_get_lsize
49  use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
50  use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
51  use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init, pawfgrtab_print
52  use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_inquire_dim
53  use m_pawdij,          only : pawdij, symdij
54  use m_pawfgr,          only : pawfgr_type, pawfgr_init, pawfgr_destroy
55  use m_paw_sphharm,     only : setsym_ylm
56  use m_paw_init,        only : pawinit, paw_gencond
57  use m_paw_nhat,        only : nhatgrid
58  use m_paw_tools,       only : chkpawovlp
59  use m_paw_correlations,only : pawpuxinit
60  use m_paw_pwaves_lmn,  only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
61  use m_classify_bands,  only : classify_bands
62  use m_pspini,          only : pspini
63  use m_sigtk,           only : sigtk_kpts_in_erange
64  use m_iowf,            only : prtkbff
65 
66  implicit none
67 
68  private

ABINIT/wfk_analyze [ Functions ]

[ Top ] [ Functions ]

NAME

  wfk_analyze

FUNCTION

 Main routine implementing postprocessing tools for the WFK file.
 Main differences wrt cut3d:

   - MPI support.
   - No interactive prompt.
   - Run the analysis once, store all the important results in netcdf files
     and use python tools to analyze data

INPUTS

 acell(3)=Length scales of primitive translations (bohr)
 codvsn=Code version
 dtfil<datafiles_type>=Variables related to files.
 dtset<dataset_type>=All input variables for this dataset.
 pawang<pawang_type)>=PAW angular mesh and related data.
 pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
 pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
 psps<pseudopotential_type>=Variables related to pseudopotentials.
   Before entering the first time in the routine, a significant part of psps has been initialized :
   the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
   and the arrays dimensioned to npsp. All the remaining components of psps are to be initialized in
   the call to pspini. The next time the code enters bethe_salpeter, psps might be identical to the
   one of the previous dtset, in which case, no reinitialisation is scheduled in pspini.F90.
 rprim(3,3)=Dimensionless real space primitive translations.
 xred(3,natom)=Reduced atomic coordinates.

NOTES

 ON THE USE OF FFT GRIDS:
 =================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.

SOURCE

126 subroutine wfk_analyze(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred)
127 
128 !Arguments ------------------------------------
129 !scalars
130  character(len=8),intent(in) :: codvsn
131  type(datafiles_type),intent(in) :: dtfil
132  type(dataset_type),intent(in) :: dtset
133  type(pawang_type),intent(inout) :: pawang
134  type(pseudopotential_type),intent(inout) :: psps
135 !arrays
136  real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom)
137  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
138  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
139 
140 !Local variables ------------------------------
141 !scalars
142  integer,parameter :: master=0
143  integer :: comm,nprocs,my_rank,mgfftf,nfftf !,nfftf_tot
144  integer :: optcut,optgr0,optgr1,optgr2,optrad,psp_gencond !,ii
145  !integer :: option,option_test,option_dij,optrhoij
146  integer :: band,ik_ibz,spin,first_band,last_band
147  integer :: ierr,usexcnhat
148  integer :: cplex,cplex_dij,cplex_rhoij,ndij,nspden_rhoij,gnt_option
149  real(dp),parameter :: spinmagntarget=-99.99_dp
150  real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff,gsqcut_shp
151  !real(dp) :: cpu,wall,gflops
152  !real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,nelect,norm,oldefermi
153  character(len=500) :: msg
154  character(len=fnlen) :: wfk0_path,wfkfull_path
155  logical :: call_pawinit, use_paw_aeur
156  type(hdr_type) :: wfk0_hdr, hdr_kfull
157  type(crystal_t) :: cryst
158  type(ebands_t) :: ebands
159  type(pawfgr_type) :: pawfgr
160  !type(paw_dmft_type) :: paw_dmft
161  type(mpi_type) :: mpi_enreg
162  type(wfd_t) :: wfd
163  type(ddkstore_t) :: ds
164 !arrays
165  integer :: ngfftc(18),ngfftf(18)
166  integer,allocatable :: l_size_atm(:)
167  real(dp),parameter :: k0(3)=zero
168  !real(dp) :: nelect_per_spin(dtset%nsppol),n0(dtset%nsppol)
169  real(dp),pointer :: gs_eigen(:,:,:)
170  complex(gwpc),allocatable :: ur_ae(:)
171  logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:)
172  real(dp) :: tsec(2)
173  type(Pawrhoij_type),allocatable :: pawrhoij(:)
174  type(pawfgrtab_type),allocatable :: pawfgrtab(:)
175  !type(paw_ij_type),allocatable :: paw_ij(:)
176  !type(paw_an_type),allocatable :: paw_an(:)
177  type(esymm_t),allocatable :: esymm(:,:)
178  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
179 
180 !************************************************************************
181 
182  DBG_ENTER('COLL')
183 
184  ! abirules!
185  if (.False.) write(std_out,*)acell,codvsn,rprim,xred
186 
187  comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
188 
189  wfk0_path = dtfil%fnamewffk
190  if (my_rank == master) then
191    ! Accept WFK file in Fortran or netcdf format.
192    if (nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then
193      ABI_ERROR(sjoin("Cannot find GS WFK file:", ch10, msg))
194    end if
195  end if
196  call xmpi_bcast(wfk0_path, master, comm, ierr)
197  call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path))
198 
199  !call cwtime(cpu,wall,gflops,"start")
200 
201  ! Costruct crystal and ebands from the GS WFK file.
202  call wfk_read_eigenvalues(wfk0_path, gs_eigen, wfk0_hdr, comm) !,gs_occ)
203  call wfk0_hdr%vs_dtset(dtset)
204 
205  cryst = wfk0_hdr%get_crystal()
206  call cryst%print(header="Crystal structure from WFK file")
207 
208  ebands = ebands_from_hdr(wfk0_hdr, maxval(wfk0_hdr%nband), gs_eigen)
209 
210  !call ebands_update_occ(ebands, spinmagntarget)
211  call ebands_print(ebands,header="Ground state energies", prtvol=dtset%prtvol)
212  ABI_FREE(gs_eigen)
213 
214  call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,&
215                   gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0)
216 
217  call print_ngfft(ngfftc, header='Coarse FFT mesh used for the wavefunctions')
218  call print_ngfft(ngfftf, header='Dense FFT mesh used for densities and potentials')
219 
220  ! Fake MPI_type for the sequential part.
221  call initmpi_seq(mpi_enreg)
222  call init_distribfft_seq(mpi_enreg%distribfft,'c',ngfftc(2),ngfftc(3),'all')
223  call init_distribfft_seq(mpi_enreg%distribfft,'f',ngfftf(2),ngfftf(3),'all')
224 
225  ! ===========================================
226  ! === Open and read pseudopotential files ===
227  ! ===========================================
228  call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm)
229 
230  ! ============================
231  ! ==== PAW initialization ====
232  ! ============================
233  if (dtset%usepaw == 1) then
234    call chkpawovlp(cryst%natom,cryst%ntypat,dtset%pawovlp,pawtab,cryst%rmet,cryst%typat,cryst%xred)
235 
236    cplex_dij=dtset%nspinor; cplex=1; ndij=1
237 
238    ABI_MALLOC(pawrhoij,(cryst%natom))
239    call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
240                              nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
241    call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,cryst%typat,pawtab=pawtab)
242 
243    ! Initialize values for several basic arrays
244    gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
245 
246    ! Test if we have to call pawinit
247    call paw_gencond(dtset,gnt_option,"test",call_pawinit)
248 
249    if (psp_gencond==1 .or. call_pawinit) then
250      call timab(553,1,tsec)
251      gsqcut_shp = two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
252      call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,dtset%pawlcutd,dtset%pawlmix,&
253                   psps%mpsang,dtset%pawnphi,cryst%nsym,dtset%pawntheta,pawang,Pawrad,&
254                   dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
255      call timab(553,2,tsec)
256 
257      ! Update internal values
258      call paw_gencond(dtset,gnt_option,"save",call_pawinit)
259 
260    else
261      if (pawtab(1)%has_kij  ==1) pawtab(1:cryst%ntypat)%has_kij  =2
262      if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla=2
263    end if
264 
265    psps%n1xccc=MAXVAL(pawtab(1:cryst%ntypat)%usetcore)
266 
267    ! Initialize optional flags in pawtab to zero
268    ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed)
269    pawtab(:)%has_nabla = 0
270    pawtab(:)%usepawu   = 0
271    pawtab(:)%useexexch = 0
272    pawtab(:)%exchmix   =zero
273    pawtab(:)%lamb_shielding   =zero
274 
275    call setsym_ylm(cryst%gprimd,pawang%l_max-1,cryst%nsym,dtset%pawprtvol,cryst%rprimd,cryst%symrec,pawang%zarot)
276 
277    ! Initialize and compute data for DFT+U
278    !paw_dmft%use_dmft=dtset%usedmft
279    !call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
280    !    .false.,dtset%jpawu,dtset%lexexch,dtset%lpawu,cryst%ntypat,pawang,dtset%pawprtvol,&
281    !    Pawrad,pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
282    !ABI_CHECK(paw_dmft%use_dmft==0,"DMFT not available")
283    !call destroy_sc_dmft(paw_dmft)
284 
285    if (my_rank == master) call pawtab_print(pawtab, unit=std_out)
286 
287    ! Get Pawrhoij from the header of the WFK file.
288    call pawrhoij_copy(wfk0_hdr%pawrhoij,pawrhoij)
289 
290    ! Variables/arrays related to the fine FFT grid.
291    ABI_MALLOC(pawfgrtab,(cryst%natom))
292    call pawtab_get_lsize(pawtab,l_size_atm,cryst%natom,cryst%typat)
293    cplex=1
294    call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat)
295    ABI_FREE(l_size_atm)
296 
297    usexcnhat=maxval(pawtab(:)%usexcnhat)
298    ! 0 if Vloc in atomic data is Vbare    (Blochl s formulation)
299    ! 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
300    call wrtout(std_out,sjoin("using usexcnhat= ",itoa(usexcnhat)))
301    !
302    ! Identify parts of the rectangular grid where the density has to be calculated ===
303    !optcut=0; optgr0=dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-dtset%pawstgylm
304    !if (dtset%xclevel==2 .and. usexcnhat>0) optgr1=dtset%pawstgylm
305    optcut=1; optgr0=1; optgr1=1; optgr2=1; optrad=1
306 
307    call nhatgrid(cryst%atindx1,cryst%gmet,cryst%natom,cryst%natom,cryst%nattyp,ngfftf,cryst%ntypat,&
308                 optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,cryst%rprimd,cryst%typat,cryst%ucvol,cryst%xred)
309 
310    call pawfgrtab_print(pawfgrtab,cryst%natom,unit=std_out,prtvol=dtset%pawprtvol)
311 
312    !ABI_MALLOC(ks_nhat,(nfftf,dtset%nspden))
313    !ks_nhat=zero
314  else
315    ABI_MALLOC(pawfgrtab,(0))
316  end if !End of PAW Initialization
317 
318  select case (dtset%wfk_task)
319 
320  case (WFK_TASK_FULLBZ, WFK_TASK_OPTICS_FULLBZ)
321    ! Read wfk0_path and build WFK in full BZ.
322    if (my_rank == master) then
323      wfkfull_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path)
324      call wfk_tofullbz(wfk0_path, dtset, psps, pawtab, wfkfull_path, hdr_kfull)
325 
326      ! Write KB form factors.
327      if (dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then
328        call prtkbff(wfkfull_path, hdr_kfull, psps, dtset%prtvol)
329      end if
330      call hdr_kfull%free()
331    end if
332    call xmpi_barrier(comm)
333 
334    if (dtset%wfk_task == WFK_TASK_OPTICS_FULLBZ) then
335      ! Calculate the DDK matrix elements from the WFK file in the full BZ.
336      ! This is needed fo computing non-linear properties in optics as symmetries are not
337      ! implemented correctly.
338      ds%only_diago = .False.
339      call ds%compute_ddk(wfkfull_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
340      call ds%free()
341    end if
342 
343  case (WFK_TASK_KPTS_ERANGE)
344    call sigtk_kpts_in_erange(dtset, cryst, ebands, psps, pawtab, dtfil%filnam_ds(4), comm)
345 
346  !case (WFK_TASK_LDOS)
347 
348  case (WFK_TASK_DDK, WFK_TASK_DDK_DIAGO)
349    ! Calculate the DDK matrix elements from the WFK file
350    ds%only_diago = .False.; if (dtset%wfk_task == WFK_TASK_DDK_DIAGO) ds%only_diago = .True.
351    call ds%compute_ddk(wfk0_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm)
352    call ds%free()
353 
354  case (WFK_TASK_EINTERP)
355    ! Band structure interpolation from eigenvalues computed on the k-mesh.
356    call ebands_interpolate_kpath(ebands, dtset, cryst, [0, 0], dtfil%filnam_ds(4), comm)
357 
358  case (WFK_TASK_CHECK_SYMTAB)
359    call wfk_check_symtab(wfk0_path, comm)
360 
361  case (WFK_TASK_CLASSIFY)
362    ! Band classification.
363    call read_wfd()
364 
365    ABI_MALLOC(esymm,(wfd%nkibz,wfd%nsppol))
366    use_paw_aeur=.False. ! should pass ngfftf but the dense mesh is not forced to be symmetric
367 
368    do spin=1,wfd%nsppol
369      do ik_ibz=1,wfd%nkibz
370        first_band = 1
371        last_band  = wfd%nband(ik_ibz,spin)
372        call classify_bands(wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,wfd%ngfft,&
373        cryst,ebands,pawtab,pawrad,pawang,psps,dtset%tolsym,esymm(ik_ibz,spin))
374      end do
375    end do
376 
377    call esymm_free(esymm)
378    ABI_FREE(esymm)
379 
380  !case (WFK_TASK_UR)
381  !  ! plot KSS wavefunctions. Change bks_mask to select particular states.
382  !  ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol))
383  !  bks_mask=.False.; bks_mask(1:4,1,1)=.True.
384  !  call wfd%plot_ur(Cryst,Psps,Pawtab,Pawrad,ngfftf,bks_mask)
385  !  ABI_FREE(bks_mask)
386 
387  case (WFK_TASK_PAW_AEPSI)
388    ! Compute AE PAW wavefunction in real space on the dense FFT mesh.
389    call read_wfd()
390 
391    ABI_CHECK(wfd%usepaw == 1, "Not a PAW run")
392    ABI_MALLOC(paw_onsite, (cryst%natom))
393    call paw_pwaves_lmn_init(paw_onsite,cryst%natom,cryst%natom,cryst%ntypat,&
394    cryst%rprimd,cryst%xcart,pawtab,pawrad,pawfgrtab)
395 
396    ! Use dense FFT mesh
397    call wfd%change_ngfft(cryst,psps,ngfftf)
398    band = 1; spin = 1; ik_ibz = 1
399 
400    ABI_MALLOC(ur_ae, (wfd%nfft*wfd%nspinor))
401    call wfd%paw_get_aeur(band,ik_ibz,spin,cryst,paw_onsite,psps,pawtab,pawfgrtab,ur_ae)
402    ABI_FREE(ur_ae)
403 
404    call paw_pwaves_lmn_free(paw_onsite)
405    ABI_FREE(paw_onsite)
406 
407  !case ("paw_aeden")
408 
409  case default
410    ABI_ERROR(sjoin("Wrong task:", itoa(dtset%wfk_task)))
411  end select
412 
413  ! Free memory
414  call cryst%free()
415  call ebands_free(ebands)
416  call wfd%free()
417  call pawfgr_destroy(pawfgr)
418  call destroy_mpi_enreg(mpi_enreg)
419  call wfk0_hdr%free()
420 
421  ! Deallocation for PAW.
422  if (dtset%usepaw==1) then
423    call pawrhoij_free(pawrhoij)
424    ABI_FREE(pawrhoij)
425    call pawfgrtab_free(pawfgrtab)
426    !call paw_ij_free(paw_ij)
427    !ABI_FREE(paw_ij)
428    !call paw_an_free(paw_an)
429    !ABI_FREE(paw_an)
430  end if
431  ABI_FREE(pawfgrtab)
432 
433  DBG_EXIT('COLL')
434 
435  contains

wfk_analyze/read_wfd [ Functions ]

[ Top ] [ wfk_analyze ] [ Functions ]

NAME

  read_wfd

FUNCTION

  Initialize the wavefunction descriptor

SOURCE

447 subroutine read_wfd()
448 
449 ! *************************************************************************
450 
451    ABI_MALLOC(keep_ur, (ebands%mband, ebands%nkpt, ebands%nsppol))
452    ABI_MALLOC(bks_mask, (ebands%mband, ebands%nkpt, ebands%nsppol))
453    keep_ur = .False.; bks_mask = .True.
454 
455    call wfd_init(wfd,cryst,pawtab,psps,keep_ur,ebands%mband,ebands%nband,ebands%nkpt,dtset%nsppol,bks_mask,&
456      dtset%nspden,dtset%nspinor,ecut_eff,dtset%ecutsm,dtset%dilatmx,wfk0_hdr%istwfk,ebands%kptns,ngfftc,&
457      dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
458 
459    ABI_FREE(keep_ur)
460    ABI_FREE(bks_mask)
461 
462    call wfd%read_wfk(wfk0_path,IO_MODE_MPI)
463    !call wfd%test_ortho(cryst, pawtab)
464 
465  end subroutine read_wfd
466 
467 end subroutine wfk_analyze