TABLE OF CONTENTS
ABINIT/wfk_analyze [ 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
COPYRIGHT
Copyright (C) 2009-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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.
PARENTS
driver
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.
CHILDREN
wfd_init,wfd_read_wfk,wfd_test_ortho
SOURCE
64 #if defined HAVE_CONFIG_H 65 #include "config.h" 66 #endif 67 68 #include "abi_common.h" 69 70 subroutine wfk_analyze(acell,codvsn,dtfil,dtset,pawang,pawrad,pawtab,psps,rprim,xred) 71 72 use defs_basis 73 use defs_datatypes 74 use defs_abitypes 75 use m_profiling_abi 76 use m_xmpi 77 use m_errors 78 use m_hdr 79 use m_crystal 80 use m_crystal_io 81 use m_ebands 82 use m_tetrahedron 83 use m_nctk 84 use m_wfk 85 use m_wfd 86 87 !use m_time, only : cwtime 88 use m_fstrings, only : strcat, sjoin, itoa, ftoa 89 use m_fftcore, only : print_ngfft 90 use m_kpts, only : tetra_from_kptrlatt 91 use m_bz_mesh, only : kpath_t, kpath_new, kpath_free 92 use m_mpinfo, only : destroy_mpi_enreg 93 use m_esymm, only : esymm_t, esymm_free, esymm_failed 94 use m_pawang, only : pawang_type 95 use m_pawrad, only : pawrad_type 96 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 97 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 98 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 99 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init, pawfgrtab_print 100 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_get_nspden, symrhoij 101 use m_pawdij, only : pawdij, symdij 102 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 103 !use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 104 !use m_paw_dmft, only : paw_dmft_type 105 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 106 107 !This section has been created automatically by the script Abilint (TD). 108 !Do not modify the following lines by hand. 109 #undef ABI_FUNC 110 #define ABI_FUNC 'wfk_analyze' 111 use interfaces_14_hidewrite 112 use interfaces_18_timing 113 use interfaces_51_manage_mpi 114 use interfaces_56_io_mpi 115 use interfaces_64_psp 116 use interfaces_65_paw 117 use interfaces_69_wfdesc 118 !End of the abilint section 119 120 implicit none 121 122 !Arguments ------------------------------------ 123 !scalars 124 character(len=6),intent(in) :: codvsn 125 type(datafiles_type),intent(in) :: dtfil 126 type(dataset_type),intent(in) :: dtset 127 type(pawang_type),intent(inout) :: pawang 128 type(pseudopotential_type),intent(inout) :: psps 129 !arrays 130 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom) 131 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 132 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 133 134 !Local variables ------------------------------ 135 !scalars 136 integer,parameter :: master=0,level40=40,brav1=1,timrev2=2,dummy_npw=1 137 integer :: comm,nprocs,my_rank,mgfftf,nfftf !,nfftf_tot 138 integer :: optcut,optgr0,optgr1,optgr2,optrad,psp_gencond !,ii 139 !integer :: option,option_test,option_dij,optrhoij 140 integer :: band,ik_ibz,spin,first_band,last_band 141 integer :: ierr,usexcnhat !,edos_intmeth 142 integer :: cplex_dij,cplex,ndij,nspden_rhoij,gnt_option 143 real(dp),parameter :: spinmagntarget=-99.99_dp 144 real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff,gsqcut_shp 145 !real(dp) :: edos_step,edos_broad 146 !real(dp) :: cpu,wall,gflops 147 !real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,nelect,norm,oldefermi 148 character(len=500) :: msg 149 character(len=fnlen) :: wfk0_path,wfkfull_path 150 logical :: call_pawinit, use_paw_aeur 151 type(hdr_type) :: wfk0_hdr 152 type(crystal_t) :: cryst 153 type(ebands_t) :: ebands 154 type(t_tetrahedron) :: tetra 155 !type(edos_t) :: edos 156 type(pawfgr_type) :: pawfgr 157 !type(paw_dmft_type) :: paw_dmft 158 type(mpi_type) :: mpi_enreg 159 type(wfd_t) :: wfd 160 !arrays 161 integer :: dummy_gvec(3,dummy_npw),ngfftc(18),ngfftf(18) 162 integer,allocatable :: l_size_atm(:) 163 real(dp),parameter :: k0(3)=zero 164 !real(dp) :: nelect_per_spin(dtset%nsppol),n0(dtset%nsppol) 165 real(dp),pointer :: gs_eigen(:,:,:) 166 complex(gwpc),allocatable :: ur_ae(:) 167 logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:) 168 real(dp) :: tsec(2) 169 type(Pawrhoij_type),allocatable :: pawrhoij(:) 170 type(pawfgrtab_type),allocatable :: pawfgrtab(:) 171 !type(paw_ij_type),allocatable :: paw_ij(:) 172 !type(paw_an_type),allocatable :: paw_an(:) 173 type(esymm_t),allocatable :: esymm(:,:) 174 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 175 176 !SHIRLEY 177 integer :: min_bsize,nbounds,ndivsm,intp_mband 178 real(dp) :: sh_coverage 179 !integer,allocatable :: intp_nband(:,:) 180 real(dp) :: ewin(2,dtset%nsppol) 181 real(dp),allocatable :: kptbounds(:,:) 182 !type(ebands_t) :: obands 183 type(kpath_t) :: kpath 184 character(len=fnlen) :: sh_fname 185 !type(wfd_t) :: owsh 186 !END SHIRLEY 187 188 !************************************************************************ 189 190 DBG_ENTER('COLL') 191 192 ! abirules! 193 if (.False.) write(std_out,*)acell,codvsn,rprim,xred 194 195 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 196 197 wfk0_path = dtfil%fnamewffk 198 if (my_rank == master) then 199 ! Accept WFK file in Fortran or netcdf format. 200 if (nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then 201 MSG_ERROR(sjoin("Cannot find GS WFK file:", ch10, msg)) 202 end if 203 end if 204 call xmpi_bcast(wfk0_path,master,comm,ierr) 205 call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path) ) 206 207 !call cwtime(cpu,wall,gflops,"start") 208 209 ! Costruct crystal and ebands from the GS WFK file. 210 call wfk_read_eigenvalues(wfk0_path,gs_eigen,wfk0_hdr,comm) !,gs_occ) 211 call hdr_vs_dtset(wfk0_hdr,dtset) 212 213 call crystal_from_hdr(cryst,wfk0_hdr,timrev2) 214 call crystal_print(cryst,header="crystal structure from WFK file") 215 216 ebands = ebands_from_hdr(wfk0_hdr,maxval(wfk0_hdr%nband),gs_eigen) 217 218 ! TODO: 219 ! Make sure everything is OK if WFK comes from a NSCF run since occ are set to zero 220 ! fermie is set to 0 if nscf! 221 222 ! Here we change the GS bands (fermi level, scissors operator ...) 223 ! All the modifications to ebands should be done here. 224 225 !if (dtset%occopt /= ebands%occopt .or. abs(dtset%tsmear - ebands%tsmear) > tol12) then 226 ! write(msg,"(2a,2(a,i0,a,f14.6,a))")& 227 ! " Changing occupation scheme as input occopt and tsmear differ from those read from WFK file.",ch10,& 228 ! " From WFK file: occopt = ",ebands%occopt,", tsmear = ",ebands%tsmear,ch10,& 229 ! " From input: occopt = ",dtset%occopt,", tsmear = ",dtset%tsmear,ch10 230 ! call wrtout(ab_out,msg) 231 ! call ebands_set_scheme(ebands,dtset%occopt,dtset%tsmear,spinmagntarget,dtset%prtvol) 232 !end if 233 234 !if (dtset%eph_fermie /= zero) then ! default value of eph_fermie is zero hence no tolerance is used! 235 ! ABI_CHECK(abs(dtset%eph_extrael) <= tol12, "eph_fermie and eph_extrael are mutually exclusive") 236 ! call wrtout(ab_out, sjoin(" Fermi level set by the user at:",ftoa(dtset%eph_fermie))) 237 ! call ebands_set_fermie(ebands, dtset%eph_fermie, msg) 238 ! call wrtout(ab_out,msg) 239 240 !else if (abs(dtset%eph_extrael) > tol12) then 241 ! NOT_IMPLEMENTED_ERROR() 242 ! ! TODO: Be careful with the trick used in elphon for passing the concentration 243 ! !call ebands_set_nelect(ebands, dtset%eph_extrael, spinmagntarget, msg) 244 ! call wrtout(ab_out,msg) 245 !end if 246 247 !call ebands_update_occ(ebands, spinmagntarget) 248 call ebands_print(ebands,header="Ground state energies",prtvol=dtset%prtvol) 249 ABI_FREE(gs_eigen) 250 251 ! Compute electron DOS. 252 ! TODO: Optimize this part. Really slow if tetra and lots of points 253 ! Could just do DOS around efermi 254 !edos_intmeth = 2; if (dtset%prtdos == 1) edos_intmeth = 1 255 !edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear 256 !edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha 257 !call edos_init(edos,ebands,cryst,edos_intmeth,edos_step,edos_broad,comm,ierr) 258 !ABI_CHECK(ierr==0, "Error in edos_init, see message above.") 259 260 ! Store DOS per spin channels 261 !n0(:) = edos%gef(1:edos%nsppol) 262 !if (my_rank == master) then 263 ! path = strcat(dtfil%filnam_ds(4), "_EDOS") 264 ! call edos_write(edos, path) 265 ! !call edos_print(edos) 266 ! write(ab_out,"(a)")sjoin("- Writing electron DOS to file:", path) 267 ! write(ab_out,'(a,es16.8,a)')' Fermi level: ',edos%mesh(edos%ief)*Ha_eV," [eV]" 268 ! write(ab_out,"(a,es16.8)")" Total electron DOS in states/eV : ",edos%gef(0) / Ha_eV 269 ! if (ebands%nsppol == 2) then 270 ! write(ab_out,"(a,es16.8)")" Spin up: ",edos%gef(1) / Ha_eV 271 ! write(ab_out,"(a,es16.8)")" Spin down:",edos%gef(2) / Ha_eV 272 ! end if 273 !end if 274 !call edos_free(edos) 275 276 ! Output useful info on the electronic bands. 277 ! Fermi Surface 278 !if (dtset%prtfsurf /= 0 .and. my_rank == master) then 279 ! path = strcat(dtfil%filnam_ds(4), "_BXSF") 280 ! if (ebands_write_bxsf(ebands,cryst,path) /= 0) then 281 ! MSG_WARNING("Cannot produce file for Fermi surface, check log file for more info") 282 ! end if 283 !end if 284 285 ! TODO Recheck getng, should use same trick as that used in screening and sigma. 286 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 287 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0) 288 289 call print_ngfft(ngfftc,header='Coarse FFT mesh used for the wavefunctions') 290 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 291 292 ! Fake MPI_type for the sequential part. 293 call initmpi_seq(mpi_enreg) 294 call init_distribfft_seq(mpi_enreg%distribfft,'c',ngfftc(2),ngfftc(3),'all') 295 call init_distribfft_seq(mpi_enreg%distribfft,'f',ngfftf(2),ngfftf(3),'all') 296 297 ! =========================================== 298 ! === Open and read pseudopotential files === 299 ! =========================================== 300 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,& 301 & pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm) 302 303 ! ============================ 304 ! ==== PAW initialization ==== 305 ! ============================ 306 if (dtset%usepaw == 1) then 307 call chkpawovlp(cryst%natom,cryst%ntypat,dtset%pawovlp,pawtab,cryst%rmet,cryst%typat,cryst%xred) 308 309 cplex_dij=dtset%nspinor; cplex=1; ndij=1 310 311 ABI_DT_MALLOC(pawrhoij,(cryst%natom)) 312 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 313 call pawrhoij_alloc(pawrhoij,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,dtset%nsppol,cryst%typat,pawtab=pawtab) 314 315 ! Initialize values for several basic arrays 316 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 317 318 ! Test if we have to call pawinit 319 call paw_gencond(dtset,gnt_option,"test",call_pawinit) 320 321 if (psp_gencond==1 .or. call_pawinit) then 322 call timab(553,1,tsec) 323 gsqcut_shp = two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 324 call pawinit(gnt_option,gsqcut_shp,zero,dtset%pawlcutd,dtset%pawlmix,& 325 & psps%mpsang,dtset%pawnphi,cryst%nsym,dtset%pawntheta,pawang,Pawrad,& 326 & dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero) 327 call timab(553,2,tsec) 328 329 ! Update internal values 330 call paw_gencond(dtset,gnt_option,"save",call_pawinit) 331 332 else 333 if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij =2 334 if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla=2 335 end if 336 337 psps%n1xccc=MAXVAL(pawtab(1:cryst%ntypat)%usetcore) 338 339 ! Initialize optional flags in pawtab to zero 340 ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed) 341 pawtab(:)%has_nabla = 0 342 pawtab(:)%usepawu = 0 343 pawtab(:)%useexexch = 0 344 pawtab(:)%exchmix =zero 345 346 call setsymrhoij(cryst%gprimd,pawang%l_max-1,cryst%nsym,dtset%pawprtvol,cryst%rprimd,cryst%symrec,pawang%zarot) 347 348 ! Initialize and compute data for LDA+U 349 !paw_dmft%use_dmft=dtset%usedmft 350 !if (dtset%usepawu>0.or.dtset%useexexch>0) then 351 ! call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 352 ! dtset%jpawu,dtset%lexexch,dtset%lpawu,cryst%ntypat,pawang,dtset%pawprtvol,& 353 ! Pawrad,pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu) 354 !end if 355 !ABI_CHECK(paw_dmft%use_dmft==0,"DMFT not available") 356 !call destroy_sc_dmft(paw_dmft) 357 358 if (my_rank == master) call pawtab_print(pawtab, unit=std_out) 359 360 ! Get Pawrhoij from the header of the WFK file. 361 call pawrhoij_copy(wfk0_hdr%pawrhoij,pawrhoij) 362 363 ! Variables/arrays related to the fine FFT grid === 364 ABI_DT_MALLOC(pawfgrtab,(cryst%natom)) 365 call pawtab_get_lsize(pawtab,l_size_atm,cryst%natom,cryst%typat) 366 cplex=1 367 call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat) 368 ABI_FREE(l_size_atm) 369 370 usexcnhat=maxval(pawtab(:)%usexcnhat) 371 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 372 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 373 call wrtout(std_out,sjoin("using usexcnhat= ",itoa(usexcnhat))) 374 ! 375 ! Identify parts of the rectangular grid where the density has to be calculated === 376 !optcut=0; optgr0=dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-dtset%pawstgylm 377 !if (dtset%xclevel==2 .and. usexcnhat>0) optgr1=dtset%pawstgylm 378 optcut=1; optgr0=1; optgr1=1; optgr2=1; optrad=1 379 380 call nhatgrid(cryst%atindx1,cryst%gmet,cryst%natom,cryst%natom,cryst%nattyp,ngfftf,cryst%ntypat,& 381 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,cryst%rprimd,cryst%typat,cryst%ucvol,cryst%xred) 382 383 call pawfgrtab_print(pawfgrtab,cryst%natom,unit=std_out,prtvol=dtset%pawprtvol) 384 385 !ABI_MALLOC(ks_nhat,(nfftf,dtset%nspden)) 386 !ks_nhat=zero 387 else 388 ABI_DT_MALLOC(pawfgrtab,(0)) 389 end if !End of PAW Initialization 390 391 select case (dtset%wfk_task) 392 393 case (WFK_TASK_FULLBZ) 394 ! Read wfk0_path and build WFK in full BZ. 395 if (my_rank == master) then 396 397 wfkfull_path = dtfil%fnameabo_wfk 398 if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path) 399 call wfk_tofullbz(wfk0_path, dtset, psps, pawtab, wfkfull_path) 400 401 ! Write tetrahedron tables. 402 tetra = tetra_from_kptrlatt(cryst, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, & 403 dtset%shiftk, dtset%nkpt, dtset%kptns, msg, ierr) 404 if (ierr == 0) then 405 call tetra_write(tetra, dtset%nkpt, dtset%kptns, strcat(dtfil%filnam_ds(4), "_TETRA")) 406 else 407 MSG_WARNING(sjoin("Cannot produce TETRA file", ch10, msg)) 408 end if 409 410 call destroy_tetra(tetra) 411 end if 412 call xmpi_barrier(comm) 413 414 !case ("pjdos") 415 416 case (WFK_TASK_CLASSIFY) 417 ! Band classification. 418 call read_wfd() 419 420 ABI_DT_MALLOC(esymm,(wfd%nkibz,wfd%nsppol)) 421 use_paw_aeur=.False. ! should pass ngfftf but the dense mesh is not forced to be symmetric 422 423 do spin=1,wfd%nsppol 424 do ik_ibz=1,wfd%nkibz 425 first_band = 1 426 last_band = wfd%nband(ik_ibz,spin) 427 call classify_bands(wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,wfd%ngfft,& 428 cryst,ebands,pawtab,pawrad,pawang,psps,dtset%tolsym,esymm(ik_ibz,spin)) 429 end do 430 end do 431 432 call esymm_free(esymm) 433 ABI_DT_FREE(esymm) 434 435 !case (WFK_TASK_UR) 436 ! ! plot KSS wavefunctions. Change bks_mask to select particular states. 437 ! ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 438 ! bks_mask=.False.; bks_mask(1:4,1,1)=.True. 439 ! call wfd_plot_ur(Wfd,Cryst,Psps,Pawtab,Pawrad,ngfftf,bks_mask) 440 ! ABI_FREE(bks_mask) 441 442 case (WFK_TASK_PAW_AEPSI) 443 ! Compute AE PAW wavefunction in real space on the dense FFT mesh. 444 call read_wfd() 445 446 ABI_CHECK(wfd%usepaw == 1, "Not a PAW run") 447 ABI_DT_MALLOC(paw_onsite, (cryst%natom)) 448 call paw_pwaves_lmn_init(paw_onsite,cryst%natom,cryst%natom,cryst%ntypat,& 449 cryst%rprimd,cryst%xcart,pawtab,pawrad,pawfgrtab) 450 451 ! Use dense FFT mesh 452 call wfd_change_ngfft(wfd,cryst,psps,ngfftf) 453 band = 1; spin = 1; ik_ibz = 1 454 455 ABI_MALLOC(ur_ae, (wfd%nfft*wfd%nspinor)) 456 call wfd_paw_get_aeur(wfd,band,ik_ibz,spin,cryst,paw_onsite,psps,pawtab,pawfgrtab,ur_ae) 457 ABI_FREE(ur_ae) 458 459 call paw_pwaves_lmn_free(paw_onsite) 460 ABI_DT_FREE(paw_onsite) 461 462 !case ("paw_aeden") 463 !case ("outqmc") 464 !subroutine outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs) 465 466 case (WFK_TASK_SHIRLEY) 467 ! Work in progress. 468 469 call read_wfd() 470 471 ! Energy window 472 !ABI_CHECK(dtset%gwpara==1,"Must use gwpara==1 for shirley") 473 ewin(1,:) = smallest_real; ewin(2,:) = greatest_real 474 if (dtset%userrc > dtset%userrb) then 475 ewin(1,:) = dtset%userrb 476 ewin(2,:) = dtset%userrc 477 write(std_out,*)"Using energy window",ewin*Ha_eV," [eV]" 478 end if 479 480 ! Select the list k-points and bands for the interpolation. 481 nbounds = 7; ndivsm=20 482 ABI_MALLOC(kptbounds, (3,nbounds)) 483 484 kptbounds = reshape([half,zero,zero, & 485 zero,zero,zero, & 486 zero,half,half, & 487 half,zero,zero, & 488 half,zero,half, & 489 zero,zero,half, & 490 zero,zero,zero], [3,nbounds]) 491 492 kptbounds = reshape([zero, zero, zero, & !Gamma 493 half, zero, zero, & !M 494 one/3, one/3, zero, & !K 495 one/3, one/3, half, & !H 496 half, zero, half, & !L 497 zero, zero, half, & !A 498 zero, zero, zero], [3,nbounds]) 499 500 kpath = kpath_new(kptbounds,cryst%gprimd,ndivsm) 501 ABI_FREE(kptbounds) 502 503 intp_mband=dtset%nband(1); min_bsize=intp_mband; sh_coverage=dtset%userra 504 505 sh_fname = strcat(dtfil%filnam_ds(4),"_SHBANDS_GSR.nc") 506 #if 0 507 if (dtset%userie==-667) then 508 ! Interpolate bands directly 509 call wrtout(std_out," Calling shirley_bands","COLL") 510 511 ! Bands on path 512 call shirley_bands(wfd,dtset,cryst,kmesh,ebands,psps,pawtab,pawang,pawrad,pawfgr,pawrhoij,paw_ij,& 513 sh_coverage,ewin,ngfftc,ngfftf,nfftf,ks_vtrial,intp_mband,kpath%points,[(zero, ii=1,kpath%npts)],sh_fname,obands,owsh) 514 515 else if (dtset%userie==-666) then 516 ! Use energy window, then interpolate bands. 517 call wrtout(std_out," Calling shirley_window","COLL") 518 519 call shirley_window(wfd,dtset,cryst,kmesh,ebands,psps,pawtab,pawang,pawrad,pawfgr,pawrhoij,paw_ij,& 520 sh_coverage,ewin,ngfftc,ngfftf,nfftf,ks_vtrial,owsh) 521 522 ABI_MALLOC(intp_nband, (kpath%npts,Wfd%nsppol)) 523 intp_nband=intp_mband 524 525 ! Energies only 526 call shirley_interp(owsh,"N",dtset,cryst,psps,pawtab,pawfgr,pawang,pawrad,& 527 pawrhoij,paw_ij,ngfftc,ngfftf,nfftf,ks_vtrial,& 528 intp_nband,intp_mband,kpath%npts,kpath%points,sh_fname,obands) 529 530 ABI_FREE(intp_nband) 531 end if 532 533 ! Write Shirley wavefunctions to file. 534 !call wfd_mkheader(oWfd,Crystal,oBands,Psps,Pawtab,KS_Pawrhoij,iHdr,pawecutdg,ngfftdg,shHdr) 535 !call wfd_write_wfk(oWsh,shHdr,oBands,"shirley_WFK") 536 !call hdr_clean(shHdr) 537 538 ! Free memory 539 call kpath_free(kpath) 540 call ebands_free(obands) 541 call wfd_free(owsh) 542 #endif 543 544 case default 545 MSG_ERROR(sjoin("Wrong task:",itoa(dtset%wfk_task))) 546 end select 547 548 !===================== 549 !==== Free memory ==== 550 !===================== 551 call crystal_free(cryst) 552 call ebands_free(ebands) 553 call wfd_free(wfd) 554 call pawfgr_destroy(pawfgr) 555 call destroy_mpi_enreg(mpi_enreg) 556 call hdr_free(wfk0_hdr) 557 558 ! Deallocation for PAW. 559 if (dtset%usepaw==1) then 560 call pawrhoij_free(pawrhoij) 561 ABI_DT_FREE(pawrhoij) 562 call pawfgrtab_free(pawfgrtab) 563 !call paw_ij_free(paw_ij) 564 !ABI_DT_FREE(paw_ij) 565 !call paw_an_free(paw_an) 566 !ABI_DT_FREE(paw_an) 567 end if 568 ABI_DT_FREE(pawfgrtab) 569 570 DBG_EXIT('COLL') 571 572 contains
wfk_analyze/read_wfd [ Functions ]
[ Top ] [ wfk_analyze ] [ Functions ]
NAME
read_wfd
FUNCTION
Initialize the wavefunction descriptor
PARENTS
wfk_analyze
CHILDREN
wfd_init,wfd_read_wfk,wfd_test_ortho
SOURCE
590 subroutine read_wfd() 591 592 593 !This section has been created automatically by the script Abilint (TD). 594 !Do not modify the following lines by hand. 595 #undef ABI_FUNC 596 #define ABI_FUNC 'read_wfd' 597 !End of the abilint section 598 599 implicit none 600 601 !Arguments ------------------------------------ 602 603 !Local variables------------------------------- 604 605 ! ************************************************************************* 606 607 ABI_MALLOC(keep_ur, (ebands%mband, ebands%nkpt, ebands%nsppol)) 608 ABI_MALLOC(bks_mask, (ebands%mband, ebands%nkpt, ebands%nsppol)) 609 keep_ur = .False.; bks_mask = .True. 610 611 call wfd_init(wfd,cryst,pawtab,psps,keep_ur,dtset%paral_kgb,dummy_npw,& 612 ebands%mband,ebands%nband,ebands%nkpt,dtset%nsppol,bks_mask,& 613 dtset%nspden,dtset%nspinor,dtset%ecutsm,dtset%dilatmx,wfk0_hdr%istwfk,ebands%kptns,ngfftc,& 614 dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut_eff) 615 616 ABI_FREE(keep_ur) 617 ABI_FREE(bks_mask) 618 619 call wfd_read_wfk(wfd,wfk0_path,IO_MODE_MPI) 620 call wfd_test_ortho(wfd,cryst,pawtab) 621 622 end subroutine read_wfd 623 624 end subroutine wfk_analyze