TABLE OF CONTENTS
- ABINIT/m_wfk
- m_wfk/fill_or_check
- m_wfk/kvars_t
- m_wfk/mpio_read_eigocc_k
- m_wfk/mpio_read_kg_k
- m_wfk/mpio_write_eigocc_k
- m_wfk/mpio_write_kg_k
- m_wfk/wfk_check_symtab
- m_wfk/wfk_check_wfkfile
- m_wfk/wfk_close
- m_wfk/wfk_compare
- m_wfk/wfk_create_wfkfile
- m_wfk/wfk_diff
- m_wfk/wfk_findk
- m_wfk/wfk_klist2mesh
- m_wfk/wfk_nc2fort
- m_wfk/wfk_ncdef_dims_vars
- m_wfk/wfk_open_read
- m_wfk/wfk_open_write
- m_wfk/wfk_print
- m_wfk/wfk_prof
- m_wfk/wfk_read_band_block
- m_wfk/wfk_read_bks
- m_wfk/wfk_read_bmask
- m_wfk/wfk_read_ebands
- m_wfk/wfk_read_eigenvalues
- m_wfk/wfk_read_eigk
- m_wfk/wfk_read_h1mat
- m_wfk/wfk_read_my_kptbands
- m_wfk/wfk_rewind
- m_wfk/wfk_seek
- m_wfk/wfk_t
- m_wfk/wfk_tofullbz
- m_wfk/wfk_update_f90ptr
- m_wfk/wfk_validate_ks
- m_wfk/wfk_write_band_block
- m_wfk/wfk_write_h1mat
- m_wfk/wfk_write_my_kptbands
- m_wfkfile/wfk_compute_offsets
- m_wfkfile/wfk_show_offsets
ABINIT/m_wfk [ Modules ]
NAME
m_wfk
FUNCTION
This module defines the wfk_t object providing a high-level API to perform common IO operations on the WFK file produced by ABINIT. The API wraps thee different formats/io-libraries: 1) binary Fortran files with sequential Fortran IO (read, write) 2) binary Fortran files with MPI-IO primitives (C-steam + Fortran records) 3) Netcdf files with parallel IO a.k.a HDF5 and emulate random access when binary Fortran files are used in *read-only* mode. See notes below for more info.
COPYRIGHT
Copyright (C) 2009-2024 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.
NOTES
1) The wfk_t object supports random access also when plain Fortran-IO is used. One can easily *read* the block of wavefunctions with a given (kpt,spin) by simply passing the appropriate indices (ik_ibz,spin) to the wfk_read_ routines. Note however that the same feature is not available in write mode when Fortran IO is used. In this case indeed one should access the block of wavefunctions according to their (kpt,spin) indices in order to write the correct record markers MPI-IO and NETCDF do not have such limitation. 2) MPI-IO read operations are done with file views even for contiguous data of the same type. I found, indeed, that mixing views with explicit offset calls causes wrong results unless the file is closed and re-open! Very strange since, according to the documentation, the two APIs do not interfere and can be mixed. Calls to MPI_FILE_SEEK to reset the pointer to the start of the file do not solve the problem. Don't know if it's a feature or a bug (the problem showed up with MPICH2, I haven't tested other MPI libraries)
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 module m_wfk 49 50 use defs_basis 51 use m_abicore 52 use m_errors 53 use m_dtset 54 #ifdef HAVE_MPI2 55 use mpi 56 #endif 57 use m_xmpi 58 use m_mpiotk 59 use m_hdr 60 use m_sort 61 use m_crystal 62 use m_pawtab 63 use m_ebands 64 use m_pawrhoij 65 use m_wffile 66 use m_nctk 67 use netcdf 68 use m_clib 69 use m_symkpt 70 71 use defs_abitypes, only : MPI_type 72 use defs_datatypes, only : pseudopotential_type, ebands_t 73 use defs_wvltypes, only : wvl_internal_type 74 use m_build_info, only : abinit_version 75 use m_geometry, only : metric 76 use m_time, only : cwtime, cwtime_report, asctime 77 use m_fstrings, only : sjoin, strcat, endswith, itoa, ktoa, ftoa 78 use m_io_tools, only : get_unit, mvrecord, iomode_from_fname, iomode2str, open_file, close_unit, delete_file, file_exists 79 use m_numeric_tools,only : mask2blocks, stats_t, stats_eval, wrap2_pmhalf 80 use m_cgtk, only : cgtk_rotate, cgtk_rotate_symrec 81 use m_fftcore, only : get_kg, ngfft_seq 82 use m_distribfft, only : init_distribfft_seq 83 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 84 use m_rwwf, only : rwwf 85 use m_kpts, only : listkk, kpts_timrev_from_kptopt 86 87 implicit none 88 89 private 90 91 #ifdef HAVE_MPI1 92 include 'mpif.h' 93 #endif 94 95 integer,private,parameter :: WFK_NOMODE = 0 96 integer,private,parameter :: WFK_READMODE = 1 97 integer,private,parameter :: WFK_WRITEMODE = 2
m_wfk/fill_or_check [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
fill_or_check
FUNCTION
INPUTS
SOURCE
5178 subroutine fill_or_check(task,Hdr,Kvars,ik_ibz,spin,formeig,kg_k,cg_k,eig_k,occ_k,ierr) 5179 5180 !Arguments ------------------------------------ 5181 !scalars 5182 integer,intent(in) :: ik_ibz,spin,formeig 5183 integer,intent(out) :: ierr 5184 character(len=*),intent(in) :: task 5185 !arrays 5186 integer,intent(in) :: kg_k(:,:) 5187 real(dp),intent(inout) :: cg_k(:,:),eig_k(:),occ_k(:) 5188 type(hdr_type),intent(in) :: Hdr 5189 type(kvars_t),intent(in) :: Kvars 5190 5191 !Local variables------------------------------- 5192 !scalars 5193 integer :: nkpt,nsppol,nspinor,nband_k,npw_k,band,ipw,kspad,ii,base,idx,mpw,eigsz 5194 character(len=500) :: msg 5195 !arrays 5196 integer,allocatable :: ref_kg_k(:,:) 5197 real(dp),allocatable :: ref_eig_k(:),ref_occ_k(:),ref_cg_k(:,:) 5198 !************************************************************************ 5199 5200 ierr = 0 5201 5202 nkpt = Hdr%nkpt 5203 nsppol = Hdr%nsppol 5204 nspinor = Hdr%nspinor 5205 5206 nband_k = Hdr%nband(ik_ibz + (spin-1)*nkpt) 5207 npw_k = Hdr%npwarr(ik_ibz) 5208 5209 ABI_MALLOC(ref_kg_k,(3,npw_k)) 5210 ABI_MALLOC(ref_eig_k,((2*nband_k)**formeig*nband_k)) 5211 ABI_MALLOC(ref_occ_k,(nband_k)) 5212 ABI_MALLOC(ref_cg_k,(2,npw_k*nspinor*nband_k)) 5213 5214 ref_kg_k = Kvars%kg_k 5215 5216 ! Pad values according to (k,s). 5217 kspad = (spin-1)*nkpt + (ik_ibz-1) * nband_k 5218 5219 if (formeig==0) then 5220 eigsz = nband_k 5221 do band=1,nband_k 5222 ref_eig_k(band) = half * (kspad + band) 5223 ref_occ_k(band) = two * (kspad + band) 5224 end do 5225 else if (formeig==1) then 5226 eigsz = 2*nband_k**2 5227 base=0 5228 do band=1,nband_k 5229 do ii=1,2*nband_k 5230 idx = base + ii 5231 ref_eig_k(idx) = ii*(kspad + band) 5232 end do 5233 base = base + 2*nband_k 5234 end do 5235 end if 5236 5237 mpw = npw_k*nspinor*nband_k 5238 do ipw=1,mpw 5239 ref_cg_k(1,ipw) = ipw + kspad 5240 ref_cg_k(2,ipw) = -ipw + kspad 5241 end do 5242 5243 SELECT CASE (task) 5244 CASE ("fill") 5245 cg_k(:,1:mpw) = ref_cg_k(:,1:mpw) 5246 if (formeig==0) then 5247 eig_k(1:nband_k) = ref_eig_k 5248 occ_k(1:nband_k) = ref_occ_k 5249 else 5250 eig_k(1:2*nband_k**2) = ref_eig_k 5251 end if 5252 5253 CASE ("check") 5254 5255 if (ANY( ABS(cg_k(:,1:mpw) - ref_cg_k) > zero)) then 5256 ierr = ierr + 1 5257 ABI_WARNING("Difference in cg_k") 5258 end if 5259 5260 if (ANY( ABS(kg_k - ref_kg_k) > zero)) then 5261 ierr = ierr + 2 5262 ABI_WARNING("Difference in kg_k") 5263 !write(std_out,*)"ref_kg_k",ref_kg_k 5264 !write(std_out,*)"kg_k",kg_k 5265 end if 5266 5267 if (ANY( ABS(eig_k(1:eigsz) - ref_eig_k) > zero)) then 5268 ierr = ierr + 4 5269 ABI_WARNING("Difference in eig_k") 5270 !write(std_out,*)"ref_eig_k",ref_eig_k 5271 !write(std_out,*)"eig_k",eig_k 5272 end if 5273 5274 if (formeig==0) then 5275 if (ANY( ABS(occ_k(1:nband_k) - ref_occ_k) > zero)) then 5276 ierr = ierr + 8 5277 ABI_WARNING("occ_k") 5278 !write(std_out,*)"ref_occ_k",ref_occ_k 5279 !write(std_out,*)"occ_k",occ_k 5280 end if 5281 end if 5282 5283 write(msg,"(a,3(i0,2x))")" (ik_ibz, spin, ierr) ",ik_ibz,spin,ierr 5284 if (ierr/=0) then 5285 ABI_WARNING(TRIM(msg)//": FAILED") 5286 else 5287 call wrtout(std_out,TRIM(msg)//": OK") 5288 end if 5289 5290 CASE DEFAULT 5291 ABI_ERROR("Wrong task") 5292 END SELECT 5293 5294 ABI_FREE(ref_kg_k) 5295 ABI_FREE(ref_eig_k) 5296 ABI_FREE(ref_occ_k) 5297 ABI_FREE(ref_cg_k) 5298 5299 end subroutine fill_or_check
m_wfk/kvars_t [ Types ]
NAME
FUNCTION
SOURCE
295 type,public :: kvars_t 296 integer,allocatable :: kg_k(:,:) 297 real(dp),pointer :: occ_k(:) => null() 298 real(dp),pointer :: eig_k(:) => null() 299 end type kvars_t 300 301 CONTAINS
m_wfk/mpio_read_eigocc_k [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
mpio_read_eigocc_k
FUNCTION
Helper functions to read the eigenvalues and the occupations with MPI-IO
INPUTS
fh
offset
nband_disk
formeig
sc_mode= MPI-IO option
xmpio_single ==> for reading by current proc.
xmpio_collective ==> for collective reading.
OUTPUTS
buffer(:)
mpierr=MPI error status.
SOURCE
4151 #ifdef HAVE_MPI_IO 4152 4153 subroutine mpio_read_eigocc_k(fh,offset,nband_disk,formeig,sc_mode,buffer,mpierr) 4154 4155 !Arguments ------------------------------------ 4156 !scalars 4157 integer,intent(in) :: fh,nband_disk,formeig,sc_mode 4158 integer(XMPI_OFFSET_KIND),intent(in) :: offset 4159 integer,intent(out) :: mpierr 4160 !arrays 4161 real(dp),pointer :: buffer(:) 4162 4163 !Local variables------------------------------- 4164 !scalars 4165 integer :: myfh,bufsz,gkk_type,eneocc_type 4166 integer(XMPI_OFFSET_KIND) :: my_offset,my_offpad !,fmarker 4167 !arrays 4168 integer :: sizes(2),subsizes(2),starts(2) 4169 4170 !************************************************************************ 4171 4172 ! Workaround for XLF 4173 myfh = fh 4174 4175 SELECT CASE (formeig) 4176 CASE (0) 4177 ! 4178 ! Read both eig and occ in buffer 4179 bufsz = 2*nband_disk 4180 my_offset = offset 4181 ABI_MALLOC(buffer, (bufsz)) 4182 4183 call MPI_TYPE_CONTIGUOUS(bufsz, MPI_DOUBLE_PRECISION, eneocc_type, mpierr) 4184 ABI_HANDLE_MPIERR(mpierr) 4185 4186 call MPI_TYPE_COMMIT(eneocc_type,mpierr) 4187 ABI_HANDLE_MPIERR(mpierr) 4188 4189 call MPI_FILE_SET_VIEW(myfh, my_offset, MPI_BYTE, eneocc_type, 'native', xmpio_info, mpierr) 4190 ABI_HANDLE_MPIERR(mpierr) 4191 4192 call MPI_TYPE_FREE(eneocc_type,mpierr) 4193 ABI_HANDLE_MPIERR(mpierr) 4194 4195 if (sc_mode==xmpio_collective) then 4196 call MPI_FILE_READ_ALL(myfh,buffer,bufsz,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr) 4197 else if (sc_mode==xmpio_single) then 4198 call MPI_FILE_READ(myfh,buffer,bufsz,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr) 4199 else 4200 ABI_ERROR("Wrong sc_mode") 4201 end if 4202 ABI_HANDLE_MPIERR(mpierr) 4203 4204 CASE (1) 4205 ! Read the (nband_k,nband_k) matrix with the (complex) GKK matrix elements. 4206 bufsz = (nband_disk**2) 4207 sizes = [nband_disk, nband_disk] 4208 subsizes = [nband_disk, nband_disk] 4209 starts = [1, 1] 4210 4211 ABI_MALLOC(buffer, (2*bufsz)) 4212 4213 !my_offset = offset - xmpio_bsize_frm 4214 !call xmpio_read_dp(myfh,my_offset,sc_mode,2*nband_disk,buffer,fmarker,mpierr) 4215 !write(std_out,*)buffer(1:2*nband_disk) 4216 !ABI_ERROR("Done") 4217 4218 call xmpio_create_fsubarray_2D(sizes,subsizes,starts,MPI_DOUBLE_COMPLEX,gkk_type,my_offpad,mpierr) 4219 ABI_HANDLE_MPIERR(mpierr) 4220 4221 ! TODO: Rationalize the offsets 4222 my_offset = offset + my_offpad - xmpio_bsize_frm 4223 4224 call MPI_FILE_SET_VIEW(myfh,my_offset,MPI_BYTE,gkk_type,'native',xmpio_info,mpierr) 4225 ABI_HANDLE_MPIERR(mpierr) 4226 4227 call MPI_TYPE_FREE(gkk_type, mpierr) 4228 ABI_HANDLE_MPIERR(mpierr) 4229 4230 if (sc_mode==xmpio_collective) then 4231 call MPI_FILE_READ_ALL(myfh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 4232 else if (sc_mode==xmpio_single) then 4233 call MPI_FILE_READ(myfh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 4234 else 4235 ABI_ERROR("Wrong sc_mode") 4236 end if 4237 ABI_HANDLE_MPIERR(mpierr) 4238 4239 CASE DEFAULT 4240 ABI_ERROR("formeig not in [0,1]") 4241 END SELECT 4242 4243 end subroutine mpio_read_eigocc_k 4244 #endif
m_wfk/mpio_read_kg_k [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
mpio_read_kg_k
FUNCTION
Helper functions to read the G-vectors with MPI-IO
INPUTS
sc_mode= MPI-IO option
xmpio_single ==> for reading by current proc.
xmpio_collective ==> for collective reading.
OUTPUTS
kg_k=(3,npw_disk) = G-vectors
mpierr=MPI error status (error check is delegated to the caller)
SOURCE
4002 #ifdef HAVE_MPI_IO 4003 4004 subroutine mpio_read_kg_k(fh,offset,npw_disk,sc_mode,kg_k,mpierr) 4005 4006 !Arguments ------------------------------------ 4007 !scalars 4008 integer,intent(in) :: fh,npw_disk,sc_mode 4009 integer(XMPI_OFFSET_KIND),intent(in) :: offset 4010 integer,intent(out) :: mpierr 4011 !arrays 4012 integer,intent(out) :: kg_k(3,npw_disk) 4013 4014 !Local variables------------------------------- 4015 !scalars 4016 integer :: kg_k_type,ncount,myfh 4017 integer(XMPI_OFFSET_KIND) :: my_offset 4018 4019 !************************************************************************ 4020 4021 ! Workarounds for XLF 4022 myfh = fh 4023 ncount = 3*npw_disk 4024 my_offset = offset 4025 4026 call MPI_TYPE_CONTIGUOUS(ncount, MPI_INTEGER, kg_k_type, mpierr) 4027 ABI_HANDLE_MPIERR(mpierr) 4028 4029 call MPI_TYPE_COMMIT(kg_k_type,mpierr) 4030 ABI_HANDLE_MPIERR(mpierr) 4031 4032 call MPI_FILE_SET_VIEW(myfh,my_offset,MPI_BYTE,kg_k_type,'native',xmpio_info,mpierr) 4033 ABI_HANDLE_MPIERR(mpierr) 4034 4035 if (sc_mode==xmpio_collective) then 4036 call MPI_FILE_READ_ALL(myfh,kg_k,ncount,MPI_INTEGER,MPI_STATUS_IGNORE,mpierr) 4037 else if (sc_mode==xmpio_single) then 4038 !call MPI_File_seek(myfh, 0, MPI_SEEK_SET,mpierr) 4039 call MPI_FILE_READ(myfh,kg_k,ncount,MPI_INTEGER, MPI_STATUS_IGNORE,mpierr) 4040 else 4041 ABI_ERROR("Wrong sc_mode") 4042 end if 4043 4044 ABI_HANDLE_MPIERR(mpierr) 4045 4046 call MPI_TYPE_FREE(kg_k_type,mpierr) 4047 ABI_HANDLE_MPIERR(mpierr) 4048 4049 end subroutine mpio_read_kg_k 4050 #endif
m_wfk/mpio_write_eigocc_k [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
mpio_write_eigocc_k
FUNCTION
Helper functions to write the eigenvalues and the occupations with MPI-IO
INPUTS
fh
offset
nband_disk
formeig
sc_mode= MPI-IO option
xmpio_single ==> for writing by current proc.
xmpio_collective ==> for collective write.
OUTPUTS
buffer(:)
mpierr=MPI error status.
SOURCE
4271 #ifdef HAVE_MPI_IO 4272 4273 subroutine mpio_write_eigocc_k(fh,offset,nband_disk,formeig,sc_mode,buffer,mpierr) 4274 4275 !Arguments ------------------------------------ 4276 !scalars 4277 integer,intent(in) :: fh,nband_disk,formeig,sc_mode 4278 integer(XMPI_OFFSET_KIND),intent(in) :: offset 4279 integer,intent(out) :: mpierr 4280 !arrays 4281 real(dp),intent(in) :: buffer(:) 4282 4283 !Local variables------------------------------- 4284 !scalars 4285 integer :: bufsz,gkk_type,eneocc_type,myfh 4286 integer(XMPI_OFFSET_KIND) :: my_offset,my_offpad 4287 !arrays 4288 integer :: sizes(2),subsizes(2),starts(2) 4289 4290 !************************************************************************ 4291 4292 ! Workaround for XLF 4293 myfh = fh 4294 4295 SELECT CASE (formeig) 4296 CASE (0) 4297 ! 4298 ! write both eig and occ in buffer 4299 my_offset = offset 4300 4301 bufsz = 2*nband_disk 4302 ABI_CHECK(SIZE(buffer) >= bufsz, "buffer too small") 4303 4304 call MPI_TYPE_CONTIGUOUS(bufsz, MPI_DOUBLE_PRECISION, eneocc_type, mpierr) 4305 ABI_HANDLE_MPIERR(mpierr) 4306 4307 call MPI_TYPE_COMMIT(eneocc_type,mpierr) 4308 ABI_HANDLE_MPIERR(mpierr) 4309 4310 call MPI_FILE_SET_VIEW(myfh,my_offset,MPI_BYTE,eneocc_type,'native',xmpio_info,mpierr) 4311 ABI_HANDLE_MPIERR(mpierr) 4312 4313 call MPI_TYPE_FREE(eneocc_type,mpierr) 4314 ABI_HANDLE_MPIERR(mpierr) 4315 4316 if (sc_mode==xmpio_collective) then 4317 call MPI_FILE_WRITE_ALL(myfh,buffer,bufsz,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr) 4318 else if (sc_mode==xmpio_single) then 4319 call MPI_FILE_WRITE(myfh,buffer,bufsz,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr) 4320 else 4321 ABI_ERROR("Wrong sc_mode") 4322 end if 4323 ABI_HANDLE_MPIERR(mpierr) 4324 4325 CASE (1) 4326 !ABI_ERROR("formeig ==1 with MPI-IO not tested") 4327 ! write the (nband_k,nband_k) matrix with the (complex) GKK matrix elements. 4328 bufsz = (nband_disk**2) 4329 sizes = [nband_disk, nband_disk] 4330 subsizes = [nband_disk, nband_disk] 4331 starts = [1, 1] 4332 4333 ABI_CHECK(SIZE(buffer) >= bufsz, "buffer too small") 4334 4335 call xmpio_create_fsubarray_2D(sizes,subsizes,starts,MPI_DOUBLE_COMPLEX,gkk_type,my_offpad,mpierr) 4336 ABI_HANDLE_MPIERR(mpierr) 4337 4338 ! TODO: Rationalize the offsets 4339 my_offset = offset + my_offpad - xmpio_bsize_frm 4340 4341 call MPI_FILE_SET_VIEW(myfh,my_offset,MPI_BYTE,gkk_type,'native',xmpio_info,mpierr) 4342 ABI_HANDLE_MPIERR(mpierr) 4343 4344 call MPI_TYPE_FREE(gkk_type, mpierr) 4345 ABI_HANDLE_MPIERR(mpierr) 4346 4347 if (sc_mode==xmpio_collective) then 4348 call MPI_FILE_WRITE_ALL(myfh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 4349 else if (sc_mode==xmpio_single) then 4350 call MPI_FILE_WRITE(myfh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 4351 else 4352 ABI_ERROR("Wrong sc_mode") 4353 end if 4354 ABI_HANDLE_MPIERR(mpierr) 4355 4356 CASE DEFAULT 4357 ABI_ERROR("formeig not in [0,1]") 4358 END SELECT 4359 4360 end subroutine mpio_write_eigocc_k
m_wfk/mpio_write_kg_k [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
mpio_write_kg_k
FUNCTION
Helper function to write the G-vectors with MPI-IO
INPUTS
sc_mode= MPI-IO option
xmpio_single ==> for writing by current proc.
xmpio_collective ==> for collective write.
kg_k=(3,npw_disk) = G-vectors
OUTPUTS
mpierr=MPI error status (error check is delegated to the caller)
SOURCE
4073 #ifdef HAVE_MPI_IO 4074 4075 subroutine mpio_write_kg_k(fh,offset,npw_disk,sc_mode,kg_k,mpierr) 4076 4077 !Arguments ------------------------------------ 4078 !scalars 4079 integer,intent(in) :: fh,npw_disk,sc_mode 4080 integer(XMPI_OFFSET_KIND),intent(in) :: offset 4081 integer,intent(out) :: mpierr 4082 !arrays 4083 integer,intent(in) :: kg_k(3,npw_disk) 4084 4085 !Local variables------------------------------- 4086 !scalars 4087 integer :: myfh,kg_k_type,ncount 4088 integer(XMPI_OFFSET_KIND) :: my_offset 4089 4090 !************************************************************************ 4091 4092 DBG_ENTER("COLL") 4093 4094 ! Workarounds for XLF 4095 myfh = fh 4096 ncount = 3*npw_disk 4097 my_offset = offset 4098 4099 call MPI_TYPE_CONTIGUOUS(ncount, MPI_INTEGER, kg_k_type, mpierr) 4100 ABI_HANDLE_MPIERR(mpierr) 4101 4102 call MPI_TYPE_COMMIT(kg_k_type,mpierr) 4103 ABI_HANDLE_MPIERR(mpierr) 4104 4105 call MPI_FILE_SET_VIEW(myfh,my_offset,MPI_BYTE,kg_k_type,'native',xmpio_info,mpierr) 4106 ABI_HANDLE_MPIERR(mpierr) 4107 4108 call MPI_TYPE_FREE(kg_k_type,mpierr) 4109 ABI_HANDLE_MPIERR(mpierr) 4110 4111 if (sc_mode==xmpio_collective) then 4112 call MPI_FILE_WRITE_ALL(myfh,kg_k,ncount,MPI_INTEGER,MPI_STATUS_IGNORE,mpierr) 4113 else if (sc_mode==xmpio_single) then 4114 call MPI_FILE_WRITE(myfh,kg_k,ncount,MPI_INTEGER,MPI_STATUS_IGNORE,mpierr) 4115 else 4116 ABI_ERROR("Wrong sc_mode") 4117 end if 4118 4119 ABI_HANDLE_MPIERR(mpierr) 4120 4121 DBG_EXIT("COLL") 4122 4123 end subroutine mpio_write_kg_k 4124 #endif
m_wfk/wfk_check_symtab [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_check_symtab
FUNCTION
INPUTS
in_wfkpath = Input WFK file generated with kptopt 3
Only GS WFK files supported (formeig==0)
SOURCE
5761 subroutine wfk_check_symtab(in_wfkpath, comm) 5762 5763 use m_krank, only : krank_t, krank_new, krank_from_kptrlatt, get_ibz2bz, star_from_ibz_idx 5764 use m_kpts, only : kpts_ibz_from_kptrlatt, kpts_timrev_from_kptopt, kpts_map, kpts_map_print, kpts_pack_in_stars 5765 use m_cgtools, only : fxphas_and_cmp 5766 5767 !Arguments ------------------------------------ 5768 !scalars 5769 character(len=*),intent(in) :: in_wfkpath 5770 integer,intent(in) :: comm 5771 5772 !Local variables------------------------------- 5773 !scalars 5774 integer,parameter :: formeig0 = 0, master = 0, kptopt1 = 1 5775 integer :: spin, nband_k, mpw, mband, nspinor, ik_ibz, ik_bz !, ierr, ikf 5776 integer :: nsppol, iomode, npw_kf, npw_ki, istwf_kf, istwf_ki, ii, my_rank, ebands_timrev, nkibz, nkbz 5777 integer :: isym_k, trev_k, g0_k(3) 5778 logical :: isirr_k 5779 character(len=500) :: msg 5780 character(len=fnlen) :: my_inpath 5781 type(wfk_t) :: wfk 5782 type(crystal_t) :: cryst 5783 type(krank_t) :: krank_ibz 5784 type(ebands_t) :: ks_ebands 5785 !arrays 5786 integer :: work_ngfft(18), gmax(3), gmax_kf(3), gmax_ki(3) 5787 integer,allocatable :: symrec_kbz2ibz(:,:), symrel_kbz2ibz(:,:), symrec_ibz2bz(:), symrel_ibz2bz(:) 5788 integer,allocatable :: kg_kf(:,:), kg_ki(:,:) 5789 real(dp) :: ki(3), kf(3) 5790 real(dp),allocatable :: kibz(:,:), kbz(:,:), wtk(:), cg_kf(:,:), cg_ki(:,:), cg_symrel(:,:), cg_symrec(:,:) !, eig_k(:), occ_k(:) 5791 real(dp),allocatable :: work(:,:,:,:) 5792 5793 ! ************************************************************************* 5794 5795 my_rank = xmpi_comm_rank(comm); if (my_rank /= master) return 5796 5797 call wrtout(std_out, " In wfk_check_symtab") 5798 5799 ! Open WFK file with k-point list, extract dimensions and allocate workspace arrays. 5800 my_inpath = in_wfkpath 5801 if (nctk_try_fort_or_ncfile(my_inpath, msg) /= 0) then 5802 ABI_ERROR(msg) 5803 end if 5804 ks_ebands = wfk_read_ebands(my_inpath, xmpi_comm_self) 5805 ABI_CHECK_IEQ(ks_ebands%kptopt, 3, "kptopt should be 3") 5806 5807 iomode = iomode_from_fname(my_inpath) 5808 call wfk_open_read(wfk, my_inpath, formeig0, iomode, get_unit(), xmpi_comm_self) 5809 mband = wfk%mband; nsppol = wfk%nsppol; nspinor = wfk%nspinor 5810 5811 cryst = wfk%hdr%get_crystal() 5812 ebands_timrev = kpts_timrev_from_kptopt(kptopt1) 5813 !ebands_timrev = kpts_timrev_from_kptopt(ks_ebands%kptopt) 5814 5815 ! Get IBZ with kptopt1 5816 call kpts_ibz_from_kptrlatt(cryst, ks_ebands%kptrlatt, kptopt1, ks_ebands%nshiftk, ks_ebands%shiftk, & 5817 nkibz, kibz, wtk, nkbz, kbz) !, bz2ibz=bz2ibz) 5818 5819 ABI_CHECK(all(abs(ks_ebands%kptns - kbz) < tol12), "Wrong kbz!") 5820 5821 krank_ibz = krank_from_kptrlatt(nkibz, kibz, ks_ebands%kptrlatt, compute_invrank=.False.) 5822 5823 ! Build symmetry tables using the two conventions. 5824 5825 ABI_MALLOC(symrec_kbz2ibz, (6, nkbz)) 5826 if (kpts_map("symrec", ebands_timrev, cryst, krank_ibz, nkbz, kbz, symrec_kbz2ibz) /= 0) then 5827 ABI_ERROR("Cannot map kBZ to IBZ!") 5828 end if 5829 ! Index of IBZ k-point in the full BZ (used to access IBZ in the WFK) 5830 ABI_MALLOC(symrec_ibz2bz, (nkibz)) 5831 do ik_bz=1,nkbz 5832 ik_ibz = symrec_kbz2ibz(1,ik_bz); isym_k = symrec_kbz2ibz(2,ik_bz) 5833 trev_k = symrec_kbz2ibz(6,ik_bz); g0_k = symrec_kbz2ibz(3:5,ik_bz) 5834 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 5835 if (isirr_k) then 5836 !print *, "ik_bz, ik_ibz", ik_bz, ik_ibz 5837 symrec_ibz2bz(ik_ibz) = ik_bz 5838 end if 5839 end do 5840 5841 ABI_MALLOC(symrel_kbz2ibz, (6, nkbz)) 5842 if (kpts_map("symrel", ebands_timrev, cryst, krank_ibz, nkbz, kbz, symrel_kbz2ibz) /= 0) then 5843 ABI_ERROR("Cannot map kBZ to IBZ!") 5844 end if 5845 ! Index of IBZ k-point in the full BZ (used to access IBZ in the WFK) 5846 ABI_MALLOC(symrel_ibz2bz, (nkibz)) 5847 do ik_bz=1,nkbz 5848 ik_ibz = symrel_kbz2ibz(1,ik_bz); isym_k = symrel_kbz2ibz(2,ik_bz) 5849 trev_k = symrel_kbz2ibz(6,ik_bz); g0_k = symrel_kbz2ibz(3:5,ik_bz) 5850 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 5851 if (isirr_k) symrel_ibz2bz(ik_ibz) = ik_bz 5852 end do 5853 5854 ! Allocate workspace arrays for wavefunction block. 5855 mpw = maxval(ks_ebands%npwarr) 5856 ABI_MALLOC(kg_kf, (3, mpw)) 5857 ABI_MALLOC(kg_ki, (3, mpw)) 5858 ABI_MALLOC(cg_kf, (2, mpw * nspinor * mband)) 5859 ABI_MALLOC(cg_ki, (2, mpw * nspinor * mband)) 5860 ABI_MALLOC(cg_symrel, (2, mpw * nspinor * mband)) 5861 ABI_MALLOC(cg_symrec, (2, mpw * nspinor * mband)) 5862 !ABI_MALLOC(eig_k, ((2*mband)**wfk%formeig * mband) ) 5863 !ABI_MALLOC(occ_k, (mband)) 5864 5865 do spin=1,nsppol 5866 ! Note how we loop over the full BZ as this is what we have in the WFK file. 5867 do ik_bz=1,nkbz 5868 5869 nband_k = wfk%nband(ik_bz, spin) 5870 nband_k = min(4, nband_k) 5871 5872 ! Read wavefunctions at full k. 5873 npw_kf = wfk%hdr%npwarr(ik_bz) 5874 istwf_kf = wfk%hdr%istwfk(ik_bz) 5875 kf = ks_ebands%kptns(:, ik_bz) 5876 call wfk%read_band_block([1, nband_k], ik_bz, spin, xmpio_single, cg_k=cg_kf, kg_k=kg_kf) 5877 5878 ! ----------------------------------------------- 5879 ! Build ik_bz from ik_ibz using symrel convention 5880 ! ----------------------------------------------- 5881 5882 ik_ibz = symrel_kbz2ibz(1,ik_bz); isym_k = symrel_kbz2ibz(2,ik_bz) 5883 trev_k = symrel_kbz2ibz(6,ik_bz); g0_k = symrel_kbz2ibz(3:5,ik_bz) 5884 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 5885 5886 if (isirr_k) cycle 5887 5888 ii = symrel_ibz2bz(ik_ibz) 5889 ki = ks_ebands%kptns(:, ii) 5890 npw_ki = wfk%hdr%npwarr(ii) 5891 istwf_ki = wfk%hdr%istwfk(ii) 5892 !write(std_out, *), "kf: ", trim(ktoa(kf)), "istwf_kf:", istwf_kf 5893 !write(std_out, *), "ki: ", trim(ktoa(ki)), "istwf_ki:", istwf_ki 5894 !ABI_CHECK_IEQ(istwf_ki, istwf_kf, "istwf_ki /= istwf_kf") 5895 5896 call wfk%read_band_block([1, nband_k], ii, spin, xmpio_single, cg_k=cg_ki, kg_k=kg_ki) 5897 5898 ! FFT box must enclose the two spheres centered on kdisk and kf 5899 gmax_kf = maxval(abs(kg_kf(:, 1:npw_kf)), dim=2) 5900 gmax_ki = maxval(abs(kg_ki(:, 1:npw_ki)), dim=2) 5901 do ii=1,3 5902 gmax(ii) = max(gmax_kf(ii), gmax_ki(ii)) 5903 end do 5904 gmax = 2 * gmax + 1 5905 call ngfft_seq(work_ngfft, gmax) 5906 ABI_CALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6))) 5907 5908 call cgtk_rotate(cryst, ki, isym_k, trev_k, g0_k, nspinor, nband_k, & 5909 npw_ki, kg_ki, npw_kf, kg_kf, istwf_ki, istwf_kf, cg_ki, cg_symrel, work_ngfft, work) 5910 5911 ! Compare cg_kf with cg_symrel taking into account a possible gauge. 5912 if (.not. fxphas_and_cmp(npw_kf, nspinor, nband_k, istwf_kf, cg_kf, cg_symrel, ks_ebands%eig(:, ik_bz, spin), msg)) then 5913 call wrtout(std_out, msg) 5914 end if 5915 5916 ! ----------------------------------------------- 5917 ! Build ik_bz from ik_ibz using symrec convention 5918 ! ----------------------------------------------- 5919 ik_ibz = symrec_kbz2ibz(1,ik_bz); isym_k = symrec_kbz2ibz(2,ik_bz) 5920 trev_k = symrec_kbz2ibz(6,ik_bz); g0_k = symrec_kbz2ibz(3:5,ik_bz) 5921 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 5922 5923 ii = symrec_ibz2bz(ik_ibz) 5924 ki = ks_ebands%kptns(:, ii) 5925 npw_ki = wfk%hdr%npwarr(ii) 5926 istwf_ki = wfk%hdr%istwfk(ii) 5927 !write(std_out, *), "kf: ", trim(ktoa(kf)), "istwf_kf:", istwf_kf 5928 !write(std_out, *), "ki: ", trim(ktoa(ki)), "istwf_ki:", istwf_ki 5929 !ABI_CHECK_IEQ(istwf_ki, istwf_kf, "istwf_ki /= istwf_kf") 5930 5931 call wfk%read_band_block([1, nband_k], ii, spin, xmpio_single, cg_k=cg_ki, kg_k=kg_ki) 5932 5933 !call cgtk_rotate_symrec(cryst, ki, isym_k, trev_k, g0_k, nspinor, nband_k, & 5934 ! npw_ki, kg_ki, npw_kf, kg_kf, istwf_ki, istwf_kf, cg_ki, cg_symrec, work_ngfft, work) 5935 5936 ! Compare cg_kf with cg_symrel taking into account a possible gauge. 5937 !if (.not. fxphas_and_cmp(npw_kf, nspinor, nband_k, istwf_kf, cg_kf, cg_symrec, ks_ebands%eig(:, ik_bz, spin), msg)) then 5938 ! call wrtout(std_out, msg) 5939 !end if 5940 5941 ABI_FREE(work) 5942 end do 5943 end do 5944 5945 ! Free memory 5946 ABI_FREE(symrec_kbz2ibz) 5947 ABI_FREE(symrel_kbz2ibz) 5948 ABI_FREE(symrec_ibz2bz) 5949 ABI_FREE(symrel_ibz2bz) 5950 ABI_FREE(kibz) 5951 ABI_FREE(kbz) 5952 ABI_FREE(wtk) 5953 call krank_ibz%free() 5954 5955 ABI_FREE(kg_kf) 5956 ABI_FREE(cg_kf) 5957 ABI_FREE(kg_ki) 5958 ABI_FREE(cg_ki) 5959 ABI_FREE(cg_symrel) 5960 ABI_FREE(cg_symrec) 5961 ABI_SFREE(work) 5962 !ABI_FREE(eig_k) 5963 !ABI_FREE(occ_k) 5964 5965 call cryst%free() 5966 call ebands_free(ks_ebands) 5967 call wfk%close() 5968 5969 end subroutine wfk_check_symtab
m_wfk/wfk_check_wfkfile [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_check_wfkfile
FUNCTION
INPUTS
SOURCE
5042 subroutine wfk_check_wfkfile(wfk_fname,Hdr,iomode,method,formeig,Kvars,cwtimes,comm,ierr) 5043 5044 !Arguments ------------------------------------ 5045 !scalars 5046 integer,intent(in) :: iomode,formeig,comm,method 5047 integer,intent(out) :: ierr 5048 character(len=*),intent(in) :: wfk_fname 5049 !arrays 5050 real(dp),intent(out) :: cwtimes(2) 5051 type(hdr_type),intent(in) :: Hdr 5052 type(kvars_t),intent(in) :: Kvars(Hdr%nkpt) 5053 5054 !Local variables------------------------------- 5055 !scalars 5056 integer :: nkpt,nsppol,nspinor,ik_ibz,spin,funt,nband_k,npw_k,sc_mode 5057 integer :: my_ierr,restart,restartpaw,is,ik,ntests,test,mband 5058 real(dp) :: cpu,wall,gflops 5059 character(len=500) :: msg 5060 type(wfk_t) :: Wfk 5061 !arrays 5062 integer :: nband(Hdr%nkpt,Hdr%nsppol),spins(Hdr%nsppol),kindices(Hdr%nkpt) 5063 integer,allocatable :: kg_k(:,:) 5064 real(dp),allocatable :: cg_k(:,:),eig_k(:),occ_k(:) 5065 logical,allocatable :: bmask(:) 5066 5067 !************************************************************************ 5068 5069 !write(msg,"(3a,i2)")"Checking file: ",TRIM(wfk_fname),", with iomode = ",iomode 5070 !call wrtout(std_out,msg,"COLL") 5071 5072 ierr = 0 5073 cwtimes = zero 5074 5075 nband = RESHAPE(Hdr%nband, (/Hdr%nkpt,Hdr%nsppol/) ) 5076 nkpt = Hdr%nkpt 5077 nsppol = Hdr%nsppol 5078 nspinor = Hdr%nspinor 5079 5080 ! Open the file for writing. 5081 call cwtime(cpu,wall,gflops,"start") 5082 funt = get_unit() 5083 5084 call wfk_open_read(Wfk,wfk_fname,formeig,iomode,funt,comm) 5085 mband = Wfk%mband 5086 5087 call cwtime(cpu,wall,gflops,"stop") 5088 cwtimes = cwtimes + (/cpu,wall/) 5089 5090 ntests = 2 5091 5092 do test=1,ntests 5093 spins = (/(spin, spin=1,Hdr%nsppol)/) 5094 kindices = (/(ik_ibz, ik_ibz=1,Hdr%nkpt)/) 5095 5096 if (test==2) then ! Reverse the indices 5097 spins = (/(spin, spin=Hdr%nsppol,1,-1)/) 5098 kindices = (/(ik_ibz, ik_ibz=Hdr%nkpt,1,-1)/) 5099 end if 5100 ! 5101 do is=1,SIZE(spins) 5102 spin = spins(is) 5103 do ik=1,SIZE(kindices) 5104 ik_ibz = kindices(ik) 5105 5106 if (Wfk%debug) then 5107 call hdr_check(Wfk%fform,Wfk%fform,Hdr,Wfk%Hdr,"COLL",restart,restartpaw) 5108 end if 5109 5110 nband_k = nband(ik_ibz,spin) 5111 npw_k = Hdr%npwarr(ik_ibz) 5112 5113 ABI_MALLOC(kg_k, (3,npw_k)) 5114 ABI_MALLOC(cg_k, (2,npw_k*nspinor*nband_k)) 5115 ABI_MALLOC(eig_k, ((2*mband)**Wfk%formeig*mband) ) 5116 ABI_MALLOC(occ_k, (mband)) 5117 ! 5118 !sc_mode = xmpio_collective 5119 sc_mode = xmpio_single 5120 ABI_MALLOC(bmask, (mband)) 5121 bmask = .FALSE. 5122 bmask(1:nband_k) = .TRUE. 5123 5124 call cwtime(cpu,wall,gflops,"start") 5125 5126 if (method==0) then 5127 call wfk%read_band_block((/1,nband_k/),ik_ibz,spin,sc_mode,kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 5128 else if (method==1) then 5129 call wfk%read_bmask(bmask,ik_ibz,spin,sc_mode,kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 5130 else 5131 ABI_ERROR("Wrong method") 5132 end if 5133 5134 !call wfk%read_eigk(ik_ibz,spin,sc_mode,eig_k) 5135 !write(std_out,*)"eig_k",eig_k 5136 5137 call cwtime(cpu,wall,gflops,"stop") 5138 cwtimes = cwtimes + (/cpu,wall/) 5139 5140 ! Check the correctness of the reading. 5141 call fill_or_check("check",Hdr,Kvars(ik_ibz),ik_ibz,spin,formeig,kg_k,cg_k,eig_k,occ_k,my_ierr) 5142 5143 if (my_ierr/=0) then 5144 write(msg,"(a,i0)")"fill_or_check returned my_ierr: ",my_ierr 5145 ierr = my_ierr 5146 ABI_WARNING(msg) 5147 end if 5148 5149 ABI_FREE(kg_k) 5150 ABI_FREE(cg_k) 5151 ABI_FREE(eig_k) 5152 ABI_FREE(occ_k) 5153 5154 ABI_FREE(bmask) 5155 end do 5156 end do 5157 ! 5158 end do ! test 5159 5160 ! Close the file 5161 call wfk%close() 5162 5163 end subroutine wfk_check_wfkfile
m_wfk/wfk_close [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_close
FUNCTION
Close the wavefunction file handler and release the memory allocated Delete the file if `delete` is True. Default: False
SOURCE
671 subroutine wfk_close(Wfk, delete) 672 673 !Arguments ------------------------------------ 674 !scalars 675 class(wfk_t),intent(inout) :: Wfk 676 logical,optional,intent(in) :: delete 677 678 !Local variables------------------------------- 679 !scalars 680 integer :: ierr 681 !character(len=500) :: msg 682 #ifdef HAVE_MPI_IO 683 integer :: mpierr, nfrec 684 integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecords(:) 685 #endif 686 687 ! ************************************************************************* 688 689 DBG_ENTER("COLL") 690 691 ! Close the file only if it was open. 692 if (wfk%rw_mode /= WFK_NOMODE) then 693 Wfk%rw_mode = WFK_NOMODE 694 695 select case (Wfk%iomode) 696 case (IO_MODE_FORTRAN) 697 close(wfk%fh) 698 699 #ifdef HAVE_MPI_IO 700 case (IO_MODE_MPI) 701 call MPI_FILE_CLOSE(Wfk%fh, mpierr) 702 ABI_CHECK_MPI(mpierr, "FILE_CLOSE!") 703 704 if (Wfk%debug .and. Wfk%my_rank == Wfk%master) then 705 ! Check the fortran records. 706 call MPI_FILE_OPEN(xmpi_comm_self, Wfk%fname, MPI_MODE_RDONLY, xmpio_info, Wfk%fh, mpierr) 707 ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN") 708 call hdr_bsize_frecords(Wfk%Hdr,Wfk%formeig,nfrec,bsize_frecords) 709 call xmpio_check_frmarkers(Wfk%fh,Wfk%hdr_offset,xmpio_single,nfrec,bsize_frecords,ierr) 710 ABI_CHECK(ierr==0,"xmpio_check_frmarkers returned ierr!=0") 711 ABI_FREE(bsize_frecords) 712 call MPI_FILE_CLOSE(Wfk%fh,mpierr) 713 ABI_CHECK_MPI(mpierr, "FILE_CLOSE!") 714 end if 715 #endif 716 717 case (IO_MODE_ETSF) 718 NCF_CHECK(nf90_close(wfk%fh)) 719 720 case default 721 ABI_ERROR(sjoin('Wrong/unsupported value of iomode: ', itoa(Wfk%iomode))) 722 end select 723 end if 724 725 ! Free memory. 726 call Wfk%Hdr%free() 727 728 ABI_SFREE(Wfk%nband) 729 ABI_SFREE(Wfk%recn_ks) 730 ABI_SFREE(Wfk%offset_ks) 731 732 if (present(delete)) then 733 if (delete) call delete_file(wfk%fname, ierr) 734 end if 735 736 DBG_EXIT("COLL") 737 738 end subroutine wfk_close
m_wfk/wfk_compare [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_compare
FUNCTION
Test two wfk_t objects for consistency. Return non-zero value if test fails.
INPUTS
wfk1, wfk2 <class(wfk_t)> = WFK handlers to be compared
OUTPUT
ierr
SOURCE
1006 integer function wfk_compare(wfk1, wfk2) result(ierr) 1007 1008 !Arguments ------------------------------------ 1009 !scalars 1010 class(wfk_t),intent(in) :: wfk1, wfk2 1011 1012 !Local variables------------------------------- 1013 !scalars 1014 integer :: restart,restartpaw 1015 !character(len=500) :: msg 1016 1017 !************************************************************************ 1018 1019 ierr = 0 1020 1021 ierr=wfk1%hdr%compare(wfk2%hdr) 1022 1023 ! Test basic dimensions 1024 !if (wfk1%hdr%nsppol /= wfk2%hdr%nsppol) then 1025 ! ierr = ierr + 1; ABI_WARNING("Different nsppol") 1026 !end if 1027 !if (wfk1%hdr%nspinor /= wfk2%hdr%nspinor) then 1028 ! ierr = ierr + 1; ABI_WARNING("Different nspinor") 1029 !end if 1030 !if (wfk1%hdr%nspden /= wfk2%hdr%nspden) then 1031 ! ierr = ierr + 1; ABI_WARNING("Different nspden") 1032 !end if 1033 !if (wfk1%hdr%nkpt /= wfk2%hdr%nkpt) then 1034 ! ierr = ierr + 1; ABI_WARNING("Different nkpt") 1035 !end if 1036 if (wfk1%formeig /= wfk2%formeig) then 1037 ierr = ierr + 1; ABI_WARNING("Different formeig") 1038 end if 1039 !if (wfk1%hdr%usepaw /= wfk2%hdr%usepaw) then 1040 ! ierr = ierr + 1; ABI_WARNING("Different usepaw") 1041 !end if 1042 !if (wfk1%hdr%ntypat /= wfk2%hdr%ntypat) then 1043 ! ierr = ierr + 1; ABI_WARNING("Different ntypat") 1044 !end if 1045 !if (wfk1%hdr%natom /= wfk2%hdr%natom) then 1046 ! ierr = ierr + 1; ABI_WARNING("Different natom") 1047 !end if 1048 !if (wfk1%hdr%fform /= wfk2%hdr%fform) then 1049 ! ierr = ierr + 1; ABI_WARNING("Different fform") 1050 !end if 1051 1052 ! Return immediately if important dimensions are not equal. 1053 if (ierr /= 0) return 1054 1055 ! Test important arrays (rprimd is not tested) 1056 !if (any(wfk1%hdr%typat /= wfk2%hdr%typat)) then 1057 ! ierr = ierr + 1; ABI_WARNING("Different typat") 1058 !end if 1059 !if (any(wfk1%hdr%npwarr /= wfk2%hdr%npwarr)) then 1060 ! ierr = ierr + 1; ABI_WARNING("Different npwarr array") 1061 !end if 1062 if (any(wfk1%nband /= wfk2%nband)) then 1063 ierr = ierr + 1; ABI_WARNING("Different nband array") 1064 end if 1065 !if (any(abs(wfk1%hdr%kptns - wfk2%hdr%kptns) > tol6)) then 1066 ! ierr = ierr + 1; ABI_WARNING("Different kptns array") 1067 !end if 1068 1069 ! Call hdr_check to get a nice diff of the header but don't check restart and restartpaw. 1070 call hdr_check(wfk1%fform,wfk2%fform,wfk1%hdr,wfk2%hdr,"PERS",restart,restartpaw) 1071 1072 end function wfk_compare
m_wfk/wfk_create_wfkfile [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_create_wfkfile
FUNCTION
INPUTS
SOURCE
4942 subroutine wfk_create_wfkfile(wfk_fname,Hdr,iomode,formeig,Kvars,cwtimes,comm) 4943 4944 !Arguments ------------------------------------ 4945 !scalars 4946 integer,intent(in) :: iomode,formeig,comm 4947 character(len=*),intent(in) :: wfk_fname 4948 !arrays 4949 real(dp),intent(out) :: cwtimes(2) 4950 type(hdr_type),intent(in) :: Hdr 4951 type(kvars_t),target,intent(out) :: Kvars(Hdr%nkpt) 4952 4953 !Local variables------------------------------- 4954 !scalars 4955 integer :: nkpt,nsppol,nspinor,ierr,sc_mode 4956 integer :: ik_ibz,spin,funt,nband_k,npw_k,istwfk_k 4957 real(dp) :: cpu,wall,gflops,ucvol 4958 type(wfk_t) :: Wfk 4959 !arrays 4960 integer :: nband(Hdr%nkpt,Hdr%nsppol) 4961 integer,pointer :: kg_k(:,:) 4962 real(dp) :: kpoint(3),gmet(3,3),gprimd(3,3),rmet(3,3) 4963 real(dp),allocatable :: cg_k(:,:),eig_k(:),occ_k(:) 4964 4965 !************************************************************************ 4966 4967 cwtimes = zero 4968 4969 nband = RESHAPE(Hdr%nband, (/Hdr%nkpt,Hdr%nsppol/) ) 4970 nkpt = Hdr%nkpt 4971 nsppol = Hdr%nsppol 4972 nspinor = Hdr%nspinor 4973 4974 call metric(gmet,gprimd,dev_null,rmet,Hdr%rprimd,ucvol) 4975 4976 ! Generate the G-vectors from input Hdr%ecut. 4977 do ik_ibz=1,nkpt 4978 kpoint = Hdr%kptns(:,ik_ibz) 4979 istwfk_k = Hdr%istwfk(ik_ibz) 4980 call get_kg(kpoint,istwfk_k,Hdr%ecut,gmet,npw_k,Kvars(ik_ibz)%kg_k) 4981 ABI_CHECK(npw_k == Hdr%npwarr(ik_ibz),"npw_k != Hdr%npwarr(ik)") 4982 end do 4983 ! 4984 ! Open the file for writing. 4985 sc_mode = xmpio_collective 4986 4987 call cwtime(cpu,wall,gflops,"start") 4988 funt = get_unit() 4989 4990 call wfk%open_write(Hdr,wfk_fname,formeig,iomode,funt,comm,write_frm=.TRUE.) 4991 4992 call cwtime(cpu,wall,gflops,"stop") 4993 cwtimes = cwtimes + [cpu, wall] 4994 4995 do spin=1,nsppol 4996 do ik_ibz=1,nkpt 4997 4998 nband_k = nband(ik_ibz,spin) 4999 npw_k = Hdr%npwarr(ik_ibz) 5000 ABI_MALLOC(cg_k, (2,npw_k*nspinor*nband_k)) 5001 ABI_MALLOC(eig_k, ((2*Wfk%mband)**formeig*Wfk%mband) ) 5002 ABI_MALLOC(occ_k, (Wfk%mband)) 5003 5004 kg_k => Kvars(ik_ibz)%kg_k 5005 ! 5006 ! Fill cg_k, eig_k, occ_k using a deterministic algorithm so that 5007 ! we can check the correctness of the reading. 5008 call fill_or_check("fill",Hdr,Kvars(ik_ibz),ik_ibz,spin,formeig,kg_k,cg_k,eig_k,occ_k,ierr) 5009 ABI_CHECK(ierr==0,"filling") 5010 5011 call cwtime(cpu,wall,gflops,"start") 5012 5013 call wfk%write_band_block((/1,nband_k/),ik_ibz,spin,sc_mode,kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 5014 5015 call cwtime(cpu,wall,gflops,"stop") 5016 cwtimes = cwtimes + (/cpu,wall/) 5017 5018 ABI_FREE(cg_k) 5019 ABI_FREE(eig_k) 5020 ABI_FREE(occ_k) 5021 end do 5022 end do 5023 5024 ! Close the file 5025 call wfk%close() 5026 5027 end subroutine wfk_create_wfkfile
m_wfk/wfk_diff [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_diff
FUNCTION
Compare two WFK file for binary equality
INPUTS
SOURCE
5315 subroutine wfk_diff(fname1,fname2,formeig,comm,ierr) 5316 5317 !Arguments ------------------------------------ 5318 integer,intent(in) :: formeig,comm 5319 integer,intent(out) :: ierr 5320 character(len=*),intent(in) :: fname1,fname2 5321 5322 !Local variables------------------------------- 5323 !scalars 5324 integer,parameter :: master=0 5325 integer :: iomode1,iomode2,ik_ibz,spin,mband,nband_k 5326 integer :: npw_k,mcg,fform1,fform2,sc_mode,my_rank,nproc 5327 character(len=500) :: msg 5328 type(hdr_type) :: Hdr1,Hdr2 5329 type(wfk_t) :: Wfk1,Wfk2 5330 !arrays 5331 integer,allocatable :: kg1_k(:,:),kg2_k(:,:) 5332 real(dp),allocatable :: eig1_k(:),cg1_k(:,:),occ1_k(:) 5333 real(dp),allocatable :: eig2_k(:),cg2_k(:,:),occ2_k(:) 5334 5335 ! ************************************************************************* 5336 5337 call wrtout(std_out, "wfk_diff: comparing "//TRIM(fname1)//" "//TRIM(fname2)) 5338 5339 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 5340 sc_mode = xmpio_collective 5341 5342 call hdr_read_from_fname(Hdr1,fname1,fform1,comm) 5343 call hdr_read_from_fname(Hdr2,fname2,fform2,comm) 5344 5345 ABI_CHECK(fform1==fform2,"fform1 != fform2") 5346 ABI_CHECK(Hdr1%nsppol==Hdr2%nsppol,"nsppol1 != nsppol2") 5347 ABI_CHECK(Hdr1%nspinor==Hdr2%nspinor,"nspinor1 != nspinor2") 5348 ABI_CHECK(Hdr1%nkpt==Hdr2%nkpt,"nkpt1 != nkpt2") 5349 !call hdr_check(fform,fform0,hdr1,hdr2,"COLL",restart,restartpaw) 5350 5351 iomode1 = iomode_from_fname(fname1) 5352 iomode2 = iomode_from_fname(fname1) 5353 ABI_CHECK(iomode1==iomode2,"iomode1 != iomode2") 5354 !iomode1 = IO_MODE_FORTRAN 5355 !iomode2 = IO_MODE_MPI 5356 5357 call wfk_open_read(Wfk1,fname1,formeig,iomode1,get_unit(),comm) 5358 call wfk_open_read(Wfk2,fname2,formeig,iomode2,get_unit(),comm) 5359 5360 if (wfk1%compare(wfk2) /= 0) then 5361 ABI_ERROR("WFK files are not consistent. See above messages") 5362 end if 5363 5364 mband = Wfk1%mband 5365 ABI_CHECK(mband==Wfk2%mband,"different mband") 5366 ABI_CHECK(all(Wfk1%nband==Wfk2%nband),"different nband") 5367 ABI_CHECK(all(Wfk1%hdr%npwarr==Wfk2%hdr%npwarr),"different npwarr") 5368 5369 ierr = 0 5370 do spin=1,Wfk1%nsppol 5371 do ik_ibz=1,Wfk1%nkpt 5372 npw_k = Wfk1%Hdr%npwarr(ik_ibz) 5373 nband_k = Wfk1%nband(ik_ibz,spin) 5374 ABI_CHECK(npw_k ==Wfk2%Hdr%npwarr(ik_ibz),"different npw_k") 5375 ABI_CHECK(nband_k==Wfk2%nband(ik_ibz,spin),"different nband_k") 5376 5377 mcg = npw_k*Hdr1%nspinor*nband_k 5378 5379 ABI_MALLOC(eig1_k,((2*mband)**formeig*mband)) 5380 ABI_MALLOC(occ1_k,(mband)) 5381 ABI_MALLOC(kg1_k,(3,npw_k)) 5382 ABI_MALLOC_OR_DIE(cg1_k,(2,mcg), ierr) 5383 5384 ABI_MALLOC(eig2_k,((2*mband)**formeig*mband)) 5385 ABI_MALLOC(occ2_k,(mband)) 5386 ABI_MALLOC(kg2_k,(3,npw_k)) 5387 ABI_MALLOC_OR_DIE(cg2_k,(2,mcg), ierr) 5388 5389 ! Read the block of bands for this (k,s). 5390 call wfk1%read_band_block([1, nband_k],ik_ibz,spin,sc_mode,kg_k=kg1_k,eig_k=eig1_k,occ_k=occ1_k) !, cg_k=cg1_k, 5391 call wfk2%read_band_block([1, nband_k],ik_ibz,spin,sc_mode,kg_k=kg2_k,eig_k=eig2_k,occ_k=occ2_k) !, cg_k=cg2_k, 5392 5393 if (ANY( ABS(kg1_k - kg2_k) > zero)) then 5394 ierr = ierr + 2 5395 ABI_WARNING("Difference in kg_k") 5396 !write(std_out,*)"kg1_k",kg1_k 5397 !write(std_out,*)"kg2_k",kg2_k 5398 end if 5399 5400 if (ANY( ABS(eig1_k - eig2_k) > zero)) then 5401 ierr = ierr + 4 5402 ABI_WARNING("Difference in eig_k") 5403 !write(std_out,*)"eig1_k",eig1_k 5404 !write(std_out,*)"eig2_k",eig2_k 5405 end if 5406 5407 if (formeig==0) then 5408 if (ANY( ABS(occ1_k - occ2_k) > zero)) then 5409 ierr = ierr + 8 5410 ABI_WARNING("occ_k") 5411 write(std_out,*)"occ1_k",occ1_k 5412 write(std_out,*)"occ2_k",occ2_k 5413 end if 5414 end if 5415 5416 if (ANY( ABS(cg1_k - cg2_k) > zero)) then 5417 ierr = ierr + 1 5418 ABI_WARNING("Difference in cg_k") 5419 end if 5420 5421 write(msg,"(a,3(i0,2x))")" (ik_ibz, spin, ierr) ",ik_ibz,spin,ierr 5422 if (ierr/=0) then 5423 ABI_WARNING(TRIM(msg)//": FAILED") 5424 else 5425 call wrtout(std_out,TRIM(msg)//": OK") 5426 end if 5427 5428 ABI_FREE(eig1_k) 5429 ABI_FREE(occ1_k) 5430 ABI_FREE(kg1_k) 5431 ABI_FREE(cg1_k) 5432 5433 ABI_FREE(eig2_k) 5434 ABI_FREE(occ2_k) 5435 ABI_FREE(kg2_k) 5436 ABI_FREE(cg2_k) 5437 end do !ik_ibz 5438 end do !spin 5439 5440 call wfk1%close() 5441 call wfk2%close() 5442 5443 call Hdr1%free() 5444 call Hdr2%free() 5445 5446 end subroutine wfk_diff
m_wfk/wfk_findk [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_findk
FUNCTION
Find the index of the k-point in the WKF file. umklapp vectors are not allowed. Return -1 if not found.
INPUTS
wfk<class(wfk_t)> = WFK handler initialized and set in read mode kpt(3)=k-point in reduced coordinates. [ktol]=Optional tolerance for k-point comparison. For each reduced direction the absolute difference between the coordinates must be less that ktol
SOURCE
872 integer pure function wfk_findk(wfk, kpt, ktol) result(ikpt) 873 874 !Arguments ------------------------------------ 875 !scalars 876 real(dp),optional,intent(in) :: ktol 877 class(wfk_t),intent(in) :: wfk 878 !arrays 879 real(dp),intent(in) :: kpt(3) 880 881 !Local variables------------------------------- 882 !scalars 883 integer :: ik 884 real(dp) :: my_ktol 885 886 ! ************************************************************************* 887 888 my_ktol = 0.0001_dp; if (present(ktol)) my_ktol = ktol 889 890 !TODO: replace with krank type and routines, probably save mapping on init of the wfk object 891 ikpt = -1 892 do ik=1,wfk%hdr%nkpt 893 if (all(abs(wfk%hdr%kptns(:, ik) - kpt) < my_ktol)) then 894 ikpt = ik; exit 895 end if 896 end do 897 898 end function wfk_findk
m_wfk/wfk_klist2mesh [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_klist2mesh
FUNCTION
This routine receives a WFK file with u_k(G) given on a subset of k-points belonging to a k-mesh and generates a new WFK file with the complete list of k-points in the IBZ by filling the missing k-points with npw_k =1 and u(G=0) = zero. This routine is mainly used to prepare the computation of electron mobilities whose convergence with the k-point sampling is notoriously slow. Since only the electron/hole states close to the band edges contribute (say ~0.# eV), one can reduce significantly the computational cost of the NSCF run by computing a WFK file with kptopt == 0 and the explicit list of k-points located inside the pockets instead of computing all the k-points of the dense IBZ. Unfortunately, the EPH code expects a WFK on a k-mesh so we need to "convert" the initial WFK with the list of k-points to a new WFK file with k-points on the dense kmesh.
INPUTS
in_wfkpath = Input WFK file with k-point list. kerange_path = path to KERANGE.nc file. dtset <dataset_type>=all input variables for this dataset out_wfkpath = Output WFK file. comm = MPI communicator.
OUTPUT
Output is written to file out_wfkpath.
NOTES
Only GS WFK files are supported (formeig==0)
SOURCE
5484 subroutine wfk_klist2mesh(in_wfkpath, kerange_path, dtset, comm) 5485 5486 !Arguments ------------------------------------ 5487 !scalars 5488 character(len=*),intent(in) :: in_wfkpath, kerange_path 5489 type(dataset_type),intent(in) :: dtset 5490 integer,intent(in) :: comm 5491 !arrays 5492 5493 !Local variables------------------------------- 5494 !scalars 5495 integer,parameter :: formeig0 = 0, master = 0 5496 integer :: spin, ikf, ikin, nband_k, mpw, mband, nspinor, ierr, fine_mband 5497 integer :: nsppol, iomode, npw_k, ii, my_rank, ncid, fform, fform_kerange 5498 real(dp) :: cpu, wall, gflops, mae_meV, merr 5499 character(len=500) :: msg 5500 character(len=fnlen) :: my_inpath, out_wfkpath 5501 type(wfk_t),target :: iwfk 5502 type(wfk_t) :: owfk 5503 type(crystal_t) :: cryst 5504 type(hdr_type) :: fine_hdr 5505 type(hdr_type),pointer :: ihdr 5506 type(ebands_t) :: iwfk_ebands, fine_ebands 5507 !arrays 5508 integer,allocatable :: kf2kin(:), kg_k(:,:) !, kshe_mask(:,:,:) 5509 real(dp),allocatable :: cg_k(:,:), eig_k(:), occ_k(:), fine_eigen(:,:,:) 5510 5511 ! ************************************************************************* 5512 5513 call cwtime(cpu, wall, gflops, "start") 5514 5515 ! IO section are executed by master only, all other procs wait for the new WFK before returning. 5516 my_rank = xmpi_comm_rank(comm); if (my_rank /= master) goto 100 5517 5518 ! Read interpolated ebands and kshe_mask from KERANGE file and build fine_ebands object. 5519 ! NOTE: KERANGE is written by sigtk_kpts_in_erange in m_sigtk module. 5520 NCF_CHECK(nctk_open_read(ncid, kerange_path, xmpi_comm_self)) 5521 ! Read header associated to the fine k-mesh 5522 call hdr_ncread(fine_hdr, ncid, fform) 5523 fform_kerange = fform_from_ext("KERANGE.nc") 5524 ABI_CHECK(fform == fform_kerange, sjoin("Wrong fform. Got: ", itoa(fform), ", Expecting: ", itoa(fform_kerange))) 5525 ! Read eigenvalues and kmask 5526 fine_mband = maxval(fine_hdr%nband) 5527 ABI_MALLOC(fine_eigen, (fine_mband, fine_hdr%nkpt, fine_hdr%nsppol)) 5528 NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eigenvalues"), fine_eigen)) 5529 !NCF_CHECK(nctk_get_dim(ncid, "nkpt_inerange", nkpt_inerage)) 5530 !ABI_MALLOC(kshe_mask, (fine_hdr%nkpt, fine_hdr%nsppol, 2)) 5531 !NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "kshe_mask"), kshe_mask)) 5532 !ABI_MALLOC(krange2ibz, (nkpt_inerange)) 5533 !NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "krange2ibz"), krange2ibz)) 5534 !ABI_FREE(krange2ibz) 5535 NCF_CHECK(nf90_close(ncid)) 5536 ! Build fine_ebands 5537 fine_ebands = ebands_from_hdr(fine_hdr, fine_mband, fine_eigen) 5538 !call ebands_print(fine_ebands, header="SKW interpolated energies", prtvol=dtset%prtvol) 5539 ABI_FREE(fine_eigen) 5540 5541 if (my_rank == master) then 5542 write(std_out, "(2a)")ch10, repeat("=", 92) 5543 !call wrtout([std_out, ab_out], msg) 5544 write(std_out, "(a)")" Generating new WKF file with dense k-mesh:" 5545 write(std_out, "(2a)")" Taking ab-initio wavefunctions with k-point list from WFK file: ", trim(in_wfkpath) 5546 write(std_out, "(a)")" When the routine returns, this file will be replaced by a new one with the dense k-mesh" 5547 write(std_out, "(2a)")" Taking eigenvalues and k-point tables from KERANGE file: ", trim(kerange_path) 5548 write(std_out, "(a, 9(i0, 1x))")" fine_kptrlatt: ", fine_hdr%kptrlatt 5549 do ii=1,fine_hdr%nshiftk 5550 write(std_out, "(a, 3(f5.2, 1x))")" fine_shiftk: ", fine_hdr%shiftk(:, ii) 5551 end do 5552 write(std_out, "(2a)")repeat("=", 92), ch10 5553 end if 5554 5555 ! Open WFK file with k-point list, extract dimensions and allocate workspace arrays. 5556 my_inpath = in_wfkpath 5557 if (nctk_try_fort_or_ncfile(my_inpath, msg) /= 0) then 5558 ABI_ERROR(msg) 5559 end if 5560 iwfk_ebands = wfk_read_ebands(my_inpath, xmpi_comm_self) 5561 !call ebands_print(iwfk_ebands, header="iwfk_ebands", unit=std_out, prtvol=dtset%prtvol) 5562 5563 iomode = iomode_from_fname(my_inpath) 5564 call wfk_open_read(iwfk, my_inpath, formeig0, iomode, get_unit(), xmpi_comm_self) 5565 5566 if (my_rank == master .and. dtset%prtvol > 0) then 5567 fform = 0 5568 call iwfk%hdr%echo(fform, 3, unit=std_out, header="Header of iwfk file") 5569 call fine_hdr%echo(fform, 3, unit=std_out, header="Header of fine_hdr") 5570 end if 5571 5572 ihdr => iwfk%hdr 5573 mband = iwfk%mband; nsppol = iwfk%nsppol; nspinor = iwfk%nspinor 5574 5575 cryst = iwfk%hdr%get_crystal() 5576 5577 ! Find correspondence fine kmesh --> input WFK and handle possible mismatch 5578 !TODO: Write specialized routine wrapping listkk to find mapping without O(N2) scaling. 5579 ABI_MALLOC(kf2kin, (fine_ebands%nkpt)) 5580 kf2kin = -1 5581 !call kpts_map(iwfk_ebands%nkpt, iwfk_ebands%kptns, fine_ebands%nkpt, fine_ebands%kptns, kf2kin, xmpi_comm_self) 5582 do ikf=1,fine_ebands%nkpt 5583 do ii=1,iwfk_ebands%nkpt 5584 if (all(abs(fine_ebands%kptns(:, ikf) - iwfk_ebands%kptns(:, ii)) < tol12)) then 5585 kf2kin(ikf) = ii; exit 5586 end if 5587 end do 5588 end do 5589 5590 if (count(kf2kin /= -1) /= iwfk_ebands%nkpt) then 5591 write(msg, "(2a, 2(a,i0))")"Something wrong in the computation of fine_mesh --> input_mesh table.",ch10, & 5592 "Expecting: ", iwfk_ebands%nkpt, " matches, found: ", count(kf2kin /= -1) 5593 ABI_ERROR(msg) 5594 end if 5595 5596 ! Check weights (the list of k-points should be a subset of the kmesh specified by sigma_ngkpt). 5597 ierr = 0 5598 do ikf=1,fine_ebands%nkpt 5599 ikin = kf2kin(ikf) 5600 if (ikin == -1) cycle 5601 if (abs(ihdr%wtk(ikin) - fine_ebands%wtk(ikf)) > tol12) then 5602 ierr = ierr + 1 5603 if (ierr <= 10) write(std_out, *) "ihdr%wtk:", ihdr%wtk(ikin), "fine_ebands%wtk", fine_ebands%wtk(ikf) 5604 end if 5605 end do 5606 if (ierr /= 0) then 5607 write(msg, "(3a)") & 5608 "Mismatch between input k-weights and weigths associated to the fine mesh. ", ch10, & 5609 "Possible inconsistency between k-mesh defined by sigma_nshiftk and the list of k-points found in file." 5610 ABI_ERROR(msg) 5611 end if 5612 5613 ! TODO 5614 !fine_hdr%fermie ! EVOLVING variable 5615 !fine_hdr%residm ! EVOLVING variable 5616 5617 ! Build new header for output WFK. This is the most delicate part since all the arrays in fine_hdr 5618 ! that depend on k-points must be consistent with the fine k-mesh. 5619 mae_meV = zero 5620 do ikf=1,fine_ebands%nkpt 5621 ikin = kf2kin(ikf) 5622 5623 if (ikin == -1) then 5624 ! Set npwarr to 1 if k-point is not in input set to reduce file size. 5625 fine_ebands%npwarr(ikf) = 1 5626 fine_hdr%npwarr(ikf) = 1 5627 else 5628 fine_ebands%npwarr(ikf) = iwfk_ebands%npwarr(ikin) 5629 fine_hdr%npwarr(ikf) = iwfk_ebands%npwarr(ikin) 5630 5631 ! Insert ab-initio eigenvalues in the SKW-interpolated fine k-mesh. 5632 do spin=1,nsppol 5633 nband_k = iwfk_ebands%nband(ikin + (spin - 1) * iwfk_ebands%nkpt) 5634 merr = Ha_meV * maxval(abs(fine_ebands%eig(1:nband_k, ikf, spin) - iwfk_ebands%eig(1:nband_k, ikin, spin))) 5635 write(std_out, "(a, es12.4,a)")" MERR: ", merr, " (meV)" 5636 !if merr > 5637 !write(std_out, *)fine_ebands%eig(1:nband_k, ikf, spin) * Ha_eV 5638 !write(std_out, *)iwfk_ebands%eig(1:nband_k, ikin, spin) * Ha_eV 5639 !write(std_out, *) Ha_meV * (fine_ebands%eig(1:nband_k, ikf, spin) - iwfk_ebands%eig(1:nband_k, ikin, spin)) 5640 !end if 5641 mae_meV = max(mae_meV, merr) 5642 fine_ebands%eig(1:nband_k, ikf, spin) = iwfk_ebands%eig(1:nband_k, ikin, spin) 5643 fine_ebands%occ(1:nband_k, ikf, spin) = iwfk_ebands%occ(1:nband_k, ikin, spin) 5644 end do 5645 end if 5646 5647 end do 5648 5649 write(std_out, "(a, es12.4,a)") & 5650 " Max error between SKW interpolated energies and ab-initio quantities:", mae_meV, " (meV)" 5651 5652 !if (mae_meV > ten) then 5653 ! write(msg,"(2a,2(a,es12.4),a)") & 5654 ! "Large error in SKW interpolation!",ch10," MARE: ",mare, ", MAE: ", mae_meV, " (meV)" 5655 ! call wrtout(ab_out, msg) 5656 ! ABI_WARNING(msg) 5657 !end if 5658 5659 call ebands_update_occ(fine_ebands, dtset%spinmagntarget, prtvol=dtset%prtvol) 5660 !call pack_eneocc(nkpt, nsppol, mband, nband, bantot, array3d, vect) 5661 !fine_hdr%occ = reshape(fine_ebands%occ, fine_ebands%mband (1:nband_k, ikin, spin) 5662 call ebands_print(fine_ebands, header="fine_ebands", unit=std_out, prtvol=dtset%prtvol) 5663 5664 out_wfkpath = strcat(in_wfkpath, ".tmp") 5665 if (iomode == IO_MODE_ETSF) out_wfkpath = strcat(out_wfkpath, ".nc") 5666 call owfk%open_write(fine_hdr, out_wfkpath, iwfk%formeig, iomode, get_unit(), xmpi_comm_self) 5667 5668 if (iomode == IO_MODE_ETSF) then 5669 ! Add crystal structure and ebands if netcdf output. 5670 NCF_CHECK(cryst%ncwrite(owfk%fh)) 5671 NCF_CHECK(ebands_ncwrite(fine_ebands, owfk%fh)) 5672 end if 5673 5674 call fine_hdr%free() 5675 5676 ! Allocate workspace arrays for wavefunction block. 5677 mpw = maxval(fine_ebands%npwarr) 5678 ABI_MALLOC(kg_k, (3, mpw)) 5679 ABI_MALLOC(cg_k, (2, mpw * nspinor * mband)) 5680 ABI_MALLOC(eig_k, ((2*mband)**iwfk%formeig * mband) ) 5681 ABI_MALLOC(occ_k, (mband)) 5682 5683 do spin=1,nsppol 5684 do ikf=1,fine_ebands%nkpt 5685 ikin = kf2kin(ikf) 5686 nband_k = owfk%nband(ikf, spin) 5687 npw_k = owfk%hdr%npwarr(ikf) 5688 5689 !cg_k = zero 5690 if (ikin /= -1) then 5691 5692 ! Consistency check 5693 if (nband_k /= iwfk%nband(ikin, spin)) then 5694 ABI_ERROR(sjoin("Mismatch in nband_k", itoa(nband_k), "/=", itoa(iwfk%nband(ikin, spin)))) 5695 end if 5696 if (npw_k /= iwfk%hdr%npwarr(ikin)) then 5697 ABI_ERROR(sjoin("Mismatch in npw_k", itoa(npw_k), "/=", itoa(iwfk%hdr%npwarr(ikin)))) 5698 end if 5699 if (owfk%hdr%istwfk(ikf) /= iwfk%hdr%istwfk(ikin)) then 5700 ABI_ERROR(sjoin("Mismatch in istwfk_k", itoa(owfk%hdr%istwfk(ikf)), "/=", itoa(iwfk%hdr%istwfk(ikin)))) 5701 end if 5702 5703 ! Read wavefunctions from input WFK file. 5704 call iwfk%read_band_block([1, nband_k], ikin, spin, xmpio_single, kg_k=kg_k, cg_k=cg_k) !, eig_k=eig_k, occ_k=occ_k) 5705 else 5706 ! Fill wavefunctions with fake data (npw_k == 1) 5707 kg_k = 0 5708 cg_k = zero 5709 end if 5710 5711 ! Write (kpt, spin) block 5712 eig_k(1:nband_k) = fine_ebands%eig(1:nband_k, ikf, spin) 5713 occ_k(1:nband_k) = fine_ebands%occ(1:nband_k, ikf, spin) 5714 5715 call owfk%write_band_block([1, nband_k], ikf, spin, xmpio_single, kg_k=kg_k, cg_k=cg_k, eig_k=eig_k, occ_k=occ_k) 5716 end do 5717 end do 5718 5719 ! Free memory 5720 ABI_FREE(kg_k) 5721 ABI_FREE(cg_k) 5722 ABI_FREE(eig_k) 5723 ABI_FREE(occ_k) 5724 ABI_FREE(kf2kin) 5725 !ABI_FREE(kshe_mask) 5726 5727 call cryst%free() 5728 call ebands_free(iwfk_ebands) 5729 call ebands_free(fine_ebands) 5730 call iwfk%close() 5731 call owfk%close() 5732 5733 ! Rename files, keep backup copy of input WFK file. 5734 call delete_file(my_inpath, ierr) 5735 ABI_CHECK(ierr == 0, sjoin("Cannot remove OLD file:", my_inpath)) 5736 !ABI_CHECK(clib_rename(my_inpath, strcat(my_inpath, ".bkp")) == 0, "Failed to rename input WFK file.") 5737 ABI_CHECK(clib_rename(out_wfkpath, my_inpath) == 0, "Failed to rename output WFK file.") 5738 5739 call cwtime_report(" WFK with fine k-mesh written to file.", cpu, wall, gflops) 5740 5741 ! All procs wait here. 5742 100 call xmpi_barrier(comm) 5743 5744 end subroutine wfk_klist2mesh
m_wfk/wfk_nc2fort [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_nc2fort
FUNCTION
Convert a netcdf WFK file (nc_path) to a Fortran WFK file (fort_path).
NOTES
- This routine should be called by a single processor. - Only GS WFK files are supported (formeig==0)
SOURCE
4685 subroutine wfk_nc2fort(nc_path, fort_path) 4686 4687 !Arguments ------------------------------------ 4688 !scalars 4689 character(len=*),intent(in) :: nc_path,fort_path 4690 4691 !Local variables------------------------------- 4692 !scalars 4693 integer :: ik,spin,mband,mpw,nband_k 4694 type(wfk_t) :: iwfk,owfk 4695 !arrays 4696 integer,parameter :: formeig0=0 4697 integer,allocatable :: kg_k(:,:) 4698 real(dp),allocatable :: cg_k(:,:),eig_k(:),occ_k(:) 4699 4700 ! ************************************************************************* 4701 4702 call wrtout(std_out, sjoin("Converting:", nc_path, "to", fort_path)) 4703 4704 ! Open input file, extract dimensions and allocate workspace arrays. 4705 call wfk_open_read(iwfk,nc_path,formeig0,IO_MODE_ETSF,get_unit(),xmpi_comm_self) 4706 4707 mpw = maxval(iwfk%hdr%npwarr); mband = iwfk%mband 4708 ABI_MALLOC(kg_k, (3, mpw)) 4709 ABI_MALLOC(cg_k, (2, mpw*iwfk%nspinor*mband)) 4710 ABI_MALLOC(eig_k, ((2*mband)**iwfk%formeig*mband) ) 4711 ABI_MALLOC(occ_k, (mband)) 4712 4713 ! Open output file. 4714 call owfk%open_write(iwfk%hdr,fort_path,formeig0,IO_MODE_FORTRAN,get_unit(),xmpi_comm_self) 4715 4716 do spin=1,iwfk%nsppol 4717 do ik=1,iwfk%nkpt 4718 nband_k = iwfk%nband(ik,spin) 4719 4720 call iwfk%read_band_block([1,nband_k],ik,spin,xmpio_single,& 4721 kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 4722 4723 call owfk%write_band_block([1,nband_k],ik,spin,xmpio_single,& 4724 kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 4725 end do 4726 end do 4727 4728 ABI_FREE(kg_k) 4729 ABI_FREE(cg_k) 4730 ABI_FREE(eig_k) 4731 ABI_FREE(occ_k) 4732 4733 call iwfk%close() 4734 call owfk%close() 4735 4736 end subroutine wfk_nc2fort
m_wfk/wfk_ncdef_dims_vars [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_ncdef_dims_vars
FUNCTION
Write the heeder, fform as well as the etsf-io dimensions and variables.
INPUTS
ncid=netcdf file handler. hdr<hdr_tyep>=Abinit header fform=File type format of the header [write_hdr]=True if the header should be written (default) [iskss]=True if this is a KSS file (activate kdependent=No)
SOURCE
919 subroutine wfk_ncdef_dims_vars(ncid, hdr, fform, write_hdr, iskss) 920 921 !Arguments ------------------------------------ 922 !scalars 923 integer,intent(in) :: ncid,fform 924 type(hdr_type),intent(in) :: hdr 925 logical,optional,intent(in) :: write_hdr,iskss 926 927 !Local variables------------------------------- 928 !scalars 929 character(len=500) :: title,history 930 logical :: do_write_hdr,my_iskss 931 integer :: ivar,mpw,ncerr 932 ! ************************************************************************* 933 934 do_write_hdr = .True.; if (present(write_hdr)) do_write_hdr = write_hdr 935 my_iskss = .False.; if (present(iskss)) my_iskss = iskss 936 if (do_write_hdr) then 937 NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.)) 938 end if 939 940 ! Add the etsf header. 941 title = sjoin("WFK file generated by Abinit, version: ", abinit_version) 942 if (my_iskss) title = sjoin("KSS file generated by Abinit, version: ", abinit_version) 943 history = sjoin("Generated on: ", asctime()) 944 NCF_CHECK(nctk_add_etsf_header(ncid, title=title, history=history)) 945 946 mpw = MAXVAL(hdr%npwarr) 947 ncerr = nctk_def_dims(ncid, [& 948 nctkdim_t("real_or_complex_coefficients", 2), nctkdim_t("max_number_of_coefficients", mpw)]) 949 NCF_CHECK(ncerr) 950 951 ! Define kg_k 952 ncerr = nctk_def_arrays(ncid, [& 953 nctkarr_t("reduced_coordinates_of_plane_waves", "int", & 954 "number_of_reduced_dimensions, max_number_of_coefficients, number_of_kpoints")& 955 ]) 956 NCF_CHECK(ncerr) 957 958 NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", ivar)) 959 if (my_iskss) then 960 NCF_CHECK(nf90_put_att(ncid, ivar, "k_dependent", "no")) 961 else 962 NCF_CHECK(nf90_put_att(ncid, ivar, "k_dependent", "yes")) 963 end if 964 965 ncerr = nctk_def_arrays(ncid, [& 966 nctkarr_t("eigenvalues", "dp", "max_number_of_states, number_of_kpoints, number_of_spins") & 967 ]) 968 NCF_CHECK(ncerr) 969 NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues")) 970 971 ncerr = nctk_def_arrays(ncid, [& 972 nctkarr_t("h1_matrix_elements", "dp", "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins") & 973 ]) 974 NCF_CHECK(ncerr) 975 NCF_CHECK(nctk_set_atomic_units(ncid, "h1_matrix_elements")) 976 977 ncerr = nctk_def_arrays(ncid, nctkarr_t("coefficients_of_wavefunctions", "dp", & 978 "real_or_complex_coefficients, max_number_of_coefficients, number_of_spinor_components, & 979 &max_number_of_states, number_of_kpoints, number_of_spins")) 980 NCF_CHECK(ncerr) 981 982 !NF90_DEF_VAR_FILL(INTEGER NCID, INTEGER VARID, INTEGER NO_FILL, FILL_VALUE) 983 !NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", ivar)) 984 !NCF_CHECK(nf90_def_var_fill(ncid, ivar, 0, -one)) 985 986 end subroutine wfk_ncdef_dims_vars
m_wfk/wfk_open_read [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_open_read
FUNCTION
Open the WFK file in read mode.
INPUTS
fname = Name of the file formeig = 0 for GS wavefunctions, 1 for RF wavefunctions. iomode = access mode funt = Fortran unit numer. Only used if iomode == IO_MODE_FORTRAN comm = MPI communicator (used for collective parallel IO)
OUTPUT
Wfk<class(wfk_t)> = WFK handler initialized and set in read mode [Hdr_out]=Copy of the abinit header NOTES TODO it would be better if formeig and iomode could be determined from the file itself! e.g iomode from the file extension, and formeig from whether it is WFK or 1WF
SOURCE
330 subroutine wfk_open_read(Wfk, fname, formeig, iomode, funt, comm, Hdr_out) 331 332 !Arguments ------------------------------------ 333 !scalars 334 class(wfk_t),intent(inout) :: Wfk 335 integer,intent(in) :: iomode,comm,formeig,funt 336 character(len=*),intent(in) :: fname 337 type(hdr_type),optional,intent(inout) :: Hdr_out ! should be intent(out), but psc miscompiles the call! 338 339 !Local variables------------------------------- 340 !scalars 341 integer :: ierr,mpierr 342 character(len=500) :: msg 343 #ifdef HAVE_MPI_IO 344 integer :: fform, nfrec !,ncerr 345 integer(XMPI_OFFSET_KIND) :: offset 346 integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecords(:) 347 #endif 348 349 !************************************************************************ 350 351 DBG_ENTER("COLL") 352 353 Wfk%comm = comm 354 Wfk%master = 0 355 Wfk%my_rank = xmpi_comm_rank(comm) 356 Wfk%nproc = xmpi_comm_size(comm) 357 358 !Initialize the mandatory data of the Wfk datastructure 359 !@wfk_t 360 Wfk%rw_mode = WFK_READMODE 361 Wfk%chunk_bsize = WFK_CHUNK_BSIZE 362 Wfk%fname = fname 363 364 ! Master checks the existence of data file 365 if (wfk%my_rank == wfk%master) then 366 if (.not. file_exists(fname)) then 367 ! Trick needed to run Abinit test suite in netcdf mode. 368 if (file_exists(nctk_ncify(fname))) then 369 write(std_out,"(3a)")"- File: ",trim(fname)," does not exist but found netcdf file with similar name." 370 Wfk%fname = nctk_ncify(fname) 371 end if 372 if (.not. file_exists(Wfk%fname)) then 373 ABI_ERROR('Missing data file: '//TRIM(Wfk%fname)) 374 end if 375 end if 376 end if 377 call xmpi_bcast(wfk%fname, wfk%master, comm, ierr) 378 379 !TODO: owfk%get_mem_mb() 380 381 Wfk%formeig = formeig 382 Wfk%iomode = iomode 383 if (endswith(fname, ".nc")) wfk%iomode = IO_MODE_ETSF 384 ! This to test the different versions. 385 !wfk%iomode = IO_MODE_MPI 386 !if (.not. endswith(fname, ".nc") .and. xmpi_comm_size == 1) wfk%iomode == IO_MODE_FORTRAN 387 388 ! Reads fform and the Header. 389 call hdr_read_from_fname(Wfk%Hdr,fname,Wfk%fform,comm) 390 ABI_CHECK(Wfk%fform /= 0, "fform == 0") 391 392 if (Wfk%debug) call Wfk%Hdr%echo(Wfk%fform, 4, unit=std_out) 393 394 ! Copy the header if required. 395 if (present(Hdr_out)) call hdr_copy(Wfk%Hdr,Hdr_out) 396 397 ! Useful dimensions 398 Wfk%mband = MAXVAL(Wfk%Hdr%nband) 399 Wfk%nkpt = Wfk%Hdr%nkpt 400 Wfk%nsppol = Wfk%Hdr%nsppol 401 Wfk%nspinor = Wfk%Hdr%nspinor 402 403 ABI_MALLOC(Wfk%nband, (Wfk%nkpt,Wfk%nsppol)) 404 Wfk%nband = RESHAPE(Wfk%Hdr%nband, (/Wfk%nkpt,Wfk%nsppol/)) 405 406 ierr=0 407 select case (wfk%iomode) 408 case (IO_MODE_FORTRAN) 409 ! All processors see a local Fortran binary file. 410 ! Each node opens the file, skip the header and set f90_fptr. 411 Wfk%fh = funt 412 if (open_file(Wfk%fname,msg,unit=Wfk%fh,form="unformatted", status="old", action="read") /= 0) then 413 ABI_ERROR(msg) 414 end if 415 416 ! Precompute number of records for Fortran IO. 417 call wfk_compute_offsets(Wfk) 418 419 call hdr_skip(Wfk%fh,ierr) 420 ABI_CHECK(ierr==0, "hdr_skip returned ierr! /= 0") 421 Wfk%f90_fptr = [1,1,REC_NPW] 422 423 #ifdef HAVE_MPI_IO 424 case (IO_MODE_MPI) 425 call MPI_FILE_OPEN(Wfk%comm, Wfk%fname, MPI_MODE_RDONLY, xmpio_info, Wfk%fh, mpierr) 426 ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN") 427 !call MPI_FILE_SET_VIEW(Wfk%fh,origin,MPI_BYTE,MPI_BYTE,'native',xmpio_info,mpierr) 428 429 call hdr_mpio_skip(Wfk%fh,fform,Wfk%hdr_offset) 430 ! Precompute offsets for MPI-IO access 431 if (Wfk%hdr_offset > 0) then 432 call wfk_compute_offsets(Wfk) 433 else 434 ABI_ERROR("hdr_offset <=0") 435 end if 436 if (Wfk%debug) then 437 !print *, 'checking offsets upon open_read : ', trim(Wfk%fname) 438 offset = Wfk%hdr_offset 439 call hdr_bsize_frecords(Wfk%Hdr,Wfk%formeig,nfrec,bsize_frecords) 440 call xmpio_check_frmarkers(Wfk%fh,offset,xmpio_collective,nfrec,bsize_frecords,ierr) 441 end if 442 #endif 443 444 case (IO_MODE_ETSF) 445 NCF_CHECK(nctk_open_read(wfk%fh, wfk%fname, wfk%comm)) 446 447 case default 448 ABI_ERROR(sjoin('Wrong or unsupported iomode:', itoa(wfk%iomode))) 449 end select 450 451 DBG_EXIT("COLL") 452 453 end subroutine wfk_open_read
m_wfk/wfk_open_write [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_open_write
FUNCTION
Open the WFK file in write mode.
INPUTS
fname = Name of the file formeig = 0 for GS wavefunctions, 1 for RF wavefunctions. iomode = access mode funt = Fortran unit numer for Only used if iomode == IO_MODE_FORTRAN comm = MPI communicator (used for MPI-IO) [write_hdr]=True if the header should be written (default) [write_frm]=True if the fortran record markers should be written (default). Only if Fortran binary file.
OUTPUT
Wfk<class(wfk_t)> = WFK handler initialized and set in read mode
SOURCE
479 subroutine wfk_open_write(Wfk, Hdr, fname, formeig, iomode, funt, comm, write_hdr, write_frm) 480 481 !Arguments ------------------------------------ 482 !scalars 483 class(wfk_t),intent(out) :: Wfk 484 integer,intent(in) :: iomode,comm,formeig,funt 485 character(len=*),intent(in) :: fname 486 logical,optional,intent(in) :: write_hdr,write_frm 487 type(hdr_type),intent(in) :: Hdr 488 489 !Local variables------------------------------- 490 !scalars 491 integer :: mpierr,ierr, hdroffset(1) 492 real(dp) :: cpu,wall,gflops 493 logical :: do_write_frm,do_write_hdr 494 character(len=500) :: msg 495 #ifdef HAVE_MPI_IO 496 integer :: fform,nfrec,sc_mode 497 integer(XMPI_OFFSET_KIND) :: offset 498 integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecords(:) 499 #endif 500 integer :: ncerr 501 502 !************************************************************************ 503 504 DBG_ENTER("COLL") 505 506 do_write_hdr = .TRUE.; if (present(write_hdr)) do_write_hdr = write_hdr 507 do_write_frm = .TRUE.; if (present(write_frm)) do_write_frm = write_frm 508 509 !Initialize mandatory data of the Wfk datastructure 510 Wfk%rw_mode = WFK_WRITEMODE 511 Wfk%chunk_bsize = WFK_CHUNK_BSIZE 512 513 Wfk%fname = fname 514 Wfk%formeig = formeig 515 Wfk%iomode = iomode; if (endswith(fname, ".nc")) wfk%iomode = IO_MODE_ETSF 516 ! This to test the different versions. 517 !wfk%iomode = IO_MODE_MPI; if (.not. endswith(fname, ".nc") .and. xmpi_comm_size == 1) wfk%iomode == IO_MODE_FORTRAN 518 519 Wfk%comm = comm 520 Wfk%master = 0 521 Wfk%my_rank = xmpi_comm_rank(comm) 522 Wfk%nproc = xmpi_comm_size(comm) 523 Wfk%fform = 2 524 525 ! Copy the header 526 call hdr_copy(Hdr,Wfk%Hdr) 527 528 ! Master writes fform and the Header (write it afterwards if IO_MODE_ETSF) 529 if (Wfk%my_rank==Wfk%master .and. do_write_hdr .and. iomode /= IO_MODE_ETSF) then 530 call Wfk%Hdr%write_to_fname(Wfk%fname, Wfk%fform) 531 if (Wfk%debug) call Wfk%Hdr%echo(Wfk%fform, 4, unit=std_out) 532 end if 533 call xmpi_barrier(Wfk%comm) 534 535 ! Useful dimensions 536 Wfk%mband = MAXVAL(Wfk%Hdr%nband) 537 Wfk%nkpt = Wfk%Hdr%nkpt 538 Wfk%nsppol = Wfk%Hdr%nsppol 539 Wfk%nspinor = Wfk%Hdr%nspinor 540 541 ABI_MALLOC(Wfk%nband, (Wfk%nkpt,Wfk%nsppol)) 542 Wfk%nband = RESHAPE(Wfk%Hdr%nband, [Wfk%nkpt, Wfk%nsppol]) 543 544 ierr = 0 545 546 select case (wfk%iomode) 547 case (IO_MODE_FORTRAN) 548 ABI_CHECK(wfk%nproc == 1, "Cannot use Fortran-IO to write WFK file with nprocs > 1") 549 Wfk%fh = funt 550 if (open_file(Wfk%fname,msg,unit=Wfk%fh,form="unformatted", status="unknown", action="readwrite") /= 0) then 551 ABI_ERROR(msg) 552 end if 553 554 ! Precompute number of records for Fortran IO. 555 call wfk_compute_offsets(Wfk) 556 557 call hdr_skip(Wfk%fh,ierr) 558 Wfk%f90_fptr = [1, 1, REC_NPW] 559 560 #ifdef HAVE_MPI_IO 561 case (IO_MODE_MPI) 562 call cwtime(cpu, wall, gflops, "start") 563 564 ! FIXME: mode flags should be rationalized 565 !call MPI_FILE_OPEN(Wfk%comm, Wfk%fname, MPI_MODE_CREATE + MPI_MODE_WRONLY, xmpio_info, Wfk%fh, mpierr) 566 !call MPI_FILE_OPEN(Wfk%comm, Wfk%fname, MPI_MODE_CREATE + MPI_MODE_RDWR, xmpio_info, Wfk%fh, mpierr) 567 call MPI_FILE_OPEN(Wfk%comm, Wfk%fname, MPI_MODE_RDWR, xmpio_info, Wfk%fh, mpierr) 568 ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN") 569 570 !call MPI_FILE_SET_VIEW(Wfk%fh,origin,MPI_BYTE,MPI_BYTE,'native',xmpio_info,mpierr) 571 ! TODO 572 !%% call MPI_File_set_size(Wfk%fh, MPI_Offset size, mpierr) 573 !ABI_CHECK_MPI(mpierr, "MPI_FILE_SET_SIZE") 574 575 hdroffset = -1 576 if (Wfk%my_rank==Wfk%master) then 577 call hdr_mpio_skip(Wfk%fh,fform,Wfk%hdr_offset) 578 ABI_CHECK(fform == Wfk%fform,"fform != Wfk%fform") 579 !call wfk%Hdr%echo(wfk%fform, 4, unit=std_out) 580 hdroffset = Wfk%hdr_offset 581 end if 582 call xmpi_bcast(hdroffset, Wfk%master, Wfk%comm, ierr) 583 Wfk%hdr_offset = hdroffset(1) 584 585 ! Precompute offsets for MPI-IO access 586 if (Wfk%hdr_offset > 0) then 587 call wfk_compute_offsets(Wfk) 588 else 589 ABI_ERROR("hdr_offset <=0") 590 end if 591 call cwtime_report(" FILE_OPEN", cpu, wall, gflops) 592 593 ! Write Fortran record markers. 594 if (do_write_frm) then 595 call cwtime(cpu, wall, gflops, "start") 596 call hdr_bsize_frecords(Wfk%Hdr,Wfk%formeig,nfrec,bsize_frecords) 597 598 sc_mode = xmpio_collective 599 offset = Wfk%hdr_offset 600 601 if (sc_mode == xmpio_collective) then 602 call xmpio_write_frmarkers(Wfk%fh,offset,sc_mode,nfrec,bsize_frecords,ierr) 603 else 604 ierr = 0 605 if (Wfk%my_rank == Wfk%master) then 606 call xmpio_write_frmarkers(Wfk%fh,offset,xmpio_single,nfrec,bsize_frecords,ierr) 607 end if 608 end if 609 ABI_CHECK(ierr == 0, "xmpio_write_frmarkers returned ierr!=0") 610 611 !call MPI_FILE_SYNC(Wfk%fh,mpierr) 612 !ABI_CHECK_MPI(mpierr, "FILE_SYNC") 613 614 if (Wfk%debug) then 615 call xmpio_check_frmarkers(Wfk%fh,offset,sc_mode,nfrec,bsize_frecords,ierr) 616 ABI_CHECK(ierr == 0, "xmpio_check_frmarkers returned ierr!=0") 617 end if 618 619 ABI_FREE(bsize_frecords) 620 call cwtime_report(" write_frmarkers", cpu, wall, gflops) 621 end if 622 #endif 623 624 CASE (IO_MODE_ETSF) 625 !NCF_CHECK(nctk_open_modify(wfk%fh, wfk%fname, wfk%comm)) 626 627 if (nctk_has_mpiio) then 628 #ifdef HAVE_NETCDF_MPI 629 ncerr = nf90_create(wfk%fname, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), & 630 comm=wfk%comm, info=xmpio_info, ncid=wfk%fh) 631 NCF_CHECK_MSG(ncerr, sjoin("nf90_create: ", wfk%fname)) 632 #else 633 ABI_ERROR("You should not be here") 634 #endif 635 else 636 if (wfk%nproc > 1) then 637 ABI_ERROR("Your netcdf library does not support MPI-IO. Cannot write WFK file with nprocs > 1") 638 end if 639 640 ncerr = nf90_create(wfk%fname, nf90_write, wfk%fh) 641 NCF_CHECK_MSG(ncerr, sjoin("nf90_create: ", wfk%fname)) 642 end if 643 644 call wfk_ncdef_dims_vars(wfk%fh, hdr, wfk%fform, write_hdr=.True.) 645 NCF_CHECK(nctk_def_basedims(wfk%fh)) 646 647 ! Switch to data mode. 648 NCF_CHECK(nctk_set_datamode(wfk%fh)) 649 650 case default 651 ABI_ERROR(sjoin('Wrong/unsupported iomode: ', itoa(wfk%iomode))) 652 end select 653 654 DBG_EXIT("COLL") 655 656 end subroutine wfk_open_write
m_wfk/wfk_print [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_print
FUNCTION
Print information on the object.
INPUTS
wfk<class(wfk_t)> = WFK handler [header]=String to be printed as header for additional info. [unit]=Unit number for output. Defaults to std_out [prtvol]=Verbosity level
SOURCE
758 subroutine wfk_print(wfk,unit,header,prtvol) 759 760 !Arguments ------------------------------------ 761 !scalars 762 class(wfk_t),intent(inout) :: wfk 763 integer,optional,intent(in) :: unit,prtvol 764 character(len=*),optional,intent(in) :: header 765 766 !Local variables------------------------------- 767 integer,parameter :: rdwr4=4 768 integer :: my_unt,my_prtvol 769 character(len=500) :: msg 770 771 ! ************************************************************************* 772 773 my_unt = std_out; if (present(unit)) my_unt = unit 774 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 775 776 msg=' ==== Info on the wfk_t object ==== ' 777 if (present(header)) msg=' ==== '//trim(adjustl(header))//' ==== ' 778 call wrtout(my_unt,msg) 779 call wrtout(my_unt, sjoin(" iomode = ",itoa(wfk%iomode))) 780 781 call wfk%hdr%echo(wfk%fform, rdwr4 ,unit=my_unt) 782 783 end subroutine wfk_print
m_wfk/wfk_prof [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_prof
FUNCTION
Profiling tool for IO routines
INPUTS
wfk_fname=Filename formeig=0 for GS file, 1 for DFPT file nband=Number of bands to read. comm=MPI communicator
SOURCE
4756 subroutine wfk_prof(wfk_fname, formeig, nband, comm) 4757 4758 !Arguments ------------------------------------ 4759 integer,intent(in) :: nband,formeig,comm 4760 character(len=*),intent(in) :: wfk_fname 4761 4762 !Local variables------------------------------- 4763 !scalars 4764 integer,parameter :: rdwr1=1,master=0,optkg1=1,option1=1,tim_rwwf0=0,icg0=0,headform0=0 4765 integer :: iomode,wfk_unt,ik_ibz,spin,ierr,ii,option,mband 4766 integer :: npw_disk,nband_disk,mcg,fform,nband_read,sc_mode,my_rank,nproc 4767 real(dp) :: cpu,wall,gflops 4768 character(len=500) :: msg 4769 type(hdr_type) :: Hdr 4770 type(wfk_t) :: Wfk 4771 type(wffile_type) :: wff 4772 type(MPI_type) :: MPI_enreg_seq 4773 !arrays 4774 !integer,parameter :: io_modes(2) = (/IO_MODE_FORTRAN, IO_MODE_MPI/) 4775 integer,parameter :: io_modes(1) = (/IO_MODE_MPI/) 4776 integer :: ngfft(18) 4777 logical,allocatable :: my_bmask(:) 4778 integer,allocatable :: kg_k(:,:) 4779 real(dp),allocatable :: eig_k(:),cg_k(:,:),occ_k(:) 4780 4781 ! ************************************************************************* 4782 4783 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 4784 sc_mode = xmpio_collective 4785 4786 call hdr_read_from_fname(hdr,wfk_fname,fform,comm) 4787 4788 ! nband_read is the max number of bands we can read from this file. 4789 nband_read = nband 4790 if (nband_read <= 0) then 4791 nband_read = minval(hdr%nband) 4792 call wrtout(std_out, sjoin("nband == 0 --> Setting nband_read to:",itoa(nband_read))) 4793 end if 4794 if (nband_read > minval(hdr%nband)) then 4795 nband_read = minval(hdr%nband) 4796 call wrtout(std_out, sjoin("nband > hdr%nband --> Setting nband_read to:",itoa(nband_read))) 4797 end if 4798 4799 wfk_unt = get_unit() 4800 4801 do ii=1,SIZE(io_modes) 4802 iomode = io_modes(ii) 4803 !do option=1,3 4804 do option=1,3,2 4805 write(std_out,*)"iomode, option",iomode,option 4806 call cwtime(cpu,wall,gflops,"start") 4807 4808 select case (option) 4809 case (1) 4810 call wfk_open_read(Wfk,wfk_fname,formeig,iomode,wfk_unt,comm) 4811 4812 do spin=1,Hdr%nsppol 4813 do ik_ibz=1,Hdr%nkpt 4814 npw_disk = Hdr%npwarr(ik_ibz) 4815 nband_disk = Wfk%nband(ik_ibz,spin) 4816 4817 mcg = npw_disk*Hdr%nspinor*nband_read 4818 4819 ABI_MALLOC(eig_k,((2*Wfk%mband)**formeig*Wfk%mband)) 4820 ABI_MALLOC(occ_k,(Wfk%mband)) 4821 4822 ABI_MALLOC(kg_k,(3,npw_disk)) 4823 ABI_MALLOC_OR_DIE(cg_k,(2,mcg), ierr) 4824 4825 ! Read the block of bands for this (k,s). 4826 call wfk%read_band_block([1,nband_read],ik_ibz,spin,xmpio_collective,& 4827 kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 4828 4829 ABI_FREE(eig_k) 4830 ABI_FREE(occ_k) 4831 ABI_FREE(kg_k) 4832 ABI_FREE(cg_k) 4833 end do !ik_ibz 4834 end do !spin 4835 4836 call wfk%close() 4837 4838 case (2) 4839 call wfk_open_read(Wfk,wfk_fname,formeig,iomode,wfk_unt,comm) 4840 4841 do spin=1,Hdr%nsppol 4842 do ik_ibz=1,Hdr%nkpt 4843 npw_disk = Hdr%npwarr(ik_ibz) 4844 nband_disk = Hdr%nband(ik_ibz+(spin-1)*Hdr%nkpt) 4845 4846 ABI_MALLOC(my_bmask,(MAXVAL(Hdr%nband))) 4847 my_bmask=.False.; my_bmask(1:nband_read) = .True. 4848 4849 ABI_MALLOC(eig_k,((2*nband_disk)**formeig*nband_disk)) 4850 ABI_MALLOC(kg_k,(3,npw_disk)) 4851 ABI_MALLOC(occ_k,(nband_disk)) 4852 4853 mcg = npw_disk*Hdr%nspinor*COUNT(my_bmask) 4854 ABI_MALLOC_OR_DIE(cg_k,(2,mcg), ierr) 4855 4856 call wfk%read_bmask(my_bmask,ik_ibz,spin,sc_mode,kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 4857 !call wfk%read_band_block((/1,nband_read/),ik_ibz,spin,sc_mode,kg_k=kg_k,cg_k=cg_k,eig_k=eig_k,occ_k=occ_k) 4858 4859 ABI_FREE(my_bmask) 4860 ABI_FREE(eig_k) 4861 ABI_FREE(occ_k) 4862 ABI_FREE(kg_k) 4863 ABI_FREE(cg_k) 4864 end do !ik_ibz 4865 end do !spin 4866 4867 call wfk%close() 4868 4869 case (3) 4870 !Fake MPI_type for the sequential part. 4871 ngfft(1:6) = (/12,12,12,13,13,13/) 4872 call initmpi_seq(MPI_enreg_seq) 4873 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 4874 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft(2),ngfft(3),'all') 4875 4876 call WffOpen(iomode,comm,wfk_fname,ierr,wff,master,my_rank,wfk_unt) !,spaceComm_mpiio) ! optional argument 4877 ABI_CHECK(ierr==0,"ierr!=0") 4878 4879 call Hdr%free() 4880 call hdr_io(fform,Hdr,1,wff) 4881 call WffKg(wff,optkg1) 4882 4883 do spin=1,Hdr%nsppol 4884 do ik_ibz=1,Hdr%nkpt 4885 4886 npw_disk = Hdr%npwarr(ik_ibz) 4887 nband_disk = Hdr%nband(ik_ibz+(spin-1)*Hdr%nkpt) 4888 4889 mband = MAXVAL(Hdr%nband) 4890 mcg = npw_disk*Hdr%nspinor*nband_read 4891 4892 ABI_MALLOC(eig_k,((2*mband)**formeig*mband)) 4893 ABI_MALLOC(occ_k,(mband)) 4894 4895 ABI_MALLOC(kg_k,(3,optkg1*npw_disk)) 4896 ABI_MALLOC_OR_DIE(cg_k,(2,mcg), ierr) 4897 ! 4898 ! Read the block of bands for this (k,s). 4899 call rwwf(cg_k,eig_k,formeig,headform0,icg0,ik_ibz,spin,kg_k,mband,mcg,MPI_enreg_seq,nband_read,& 4900 nband_disk,npw_disk,Hdr%nspinor,occ_k,option1,optkg1,tim_rwwf0,Wff) 4901 4902 ABI_FREE(eig_k) 4903 ABI_FREE(occ_k) 4904 ABI_FREE(kg_k) 4905 ABI_FREE(cg_k) 4906 4907 end do !ik_ibz 4908 end do !spin 4909 4910 call WffClose(wff,ierr) 4911 call destroy_mpi_enreg(MPI_enreg_seq) 4912 4913 case default 4914 ABI_ERROR("Wrong method") 4915 end select 4916 4917 call cwtime(cpu,wall,gflops,"stop") 4918 write(msg,'(3(a,i2),2(a,f8.2))')& 4919 " iomode: ",iomode,", nproc: ",nproc,", option: ",option,", cpu: ",cpu,", wall:",wall 4920 call wrtout(std_out, msg) 4921 !call cwtime_report(" FULL_WFK written to file. ", cpu, wall, gflops) 4922 end do 4923 end do 4924 4925 call Hdr%free() 4926 4927 end subroutine wfk_prof
m_wfk/wfk_read_band_block [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_band_block
FUNCTION
Read a block of contiguous bands at a given (k-point, spin)
INPUTS
Wfk<class(wfk_t)>= WFK file handler object. band_block(2)=Initial and final band index. ik_ibz=Index of the k-point in the IBZ. spin=Spin index sc_mode= MPI-IO option xmpio_single ==> for reading by current proc. xmpio_collective ==> for collective reading in wfk%comm (use it wisely!) OUTPUTS [kg_k=(:,:)] = G-vectors [eig_k(:)] = Eigenvectors [cg_k(:,:)] = Fourier coefficients
NOTES
The output arrays eig_k and occ_k contain the *full* set of eigenvalues and occupation factors stored in the file and are dimensioned with wfk%mband.
SOURCE
1103 subroutine wfk_read_band_block(Wfk, band_block, ik_ibz, spin, sc_mode, kg_k, cg_k, eig_k, occ_k) 1104 1105 !Arguments ------------------------------------ 1106 !scalars 1107 integer,intent(in) :: ik_ibz,spin,sc_mode 1108 class(wfk_t),intent(inout) :: Wfk 1109 !arrays 1110 integer,intent(in) :: band_block(2) 1111 integer,intent(out), DEV_CONTARRD optional :: kg_k(:,:) ! (3,npw_k) 1112 real(dp),intent(out), DEV_CONTARRD optional :: cg_k(:,:) ! (2,npw_k*nspinor*nband) 1113 real(dp),intent(inout),optional :: eig_k((2*Wfk%mband)**Wfk%formeig*Wfk%mband) 1114 real(dp),intent(out),optional :: occ_k(Wfk%mband) 1115 1116 !Local variables------------------------------- 1117 !scalars 1118 integer :: ierr,npw_disk,nspinor_disk,nband_disk,band 1119 integer :: nband_disk_keep 1120 integer :: ipw,my_bcount,npwso,npw_tot_disk,nb_block,base 1121 integer :: npw_read,nspinor_read,nband_read 1122 character(len=500) :: msg,errmsg 1123 !arrays 1124 real(dp),ABI_CONTIGUOUS pointer :: tmp_eigk(:),tmp_occk(:) 1125 #ifdef HAVE_MPI_IO 1126 integer :: mpierr,bufsz,gkk_type,cgblock_type 1127 integer(XMPI_OFFSET_KIND) :: my_offset,my_offpad 1128 integer :: sizes(2),subsizes(2),starts(2),types(2) 1129 #endif 1130 integer :: kg_varid,eig_varid,occ_varid,cg_varid,ncerr,h1_varid,idx,ib1,ib2 1131 real(dp),allocatable :: h1mat(:,:,:) 1132 1133 !************************************************************************ 1134 1135 DBG_ENTER("COLL") 1136 1137 ABI_CHECK_IEQ(Wfk%rw_mode, WFK_READMODE, "Wfk must be in READMODE") 1138 1139 if (wfk_validate_ks(wfk, ik_ibz, spin) /= 0) then 1140 ABI_ERROR("Wrong (ik_ibz, spin) args, Aborting now") 1141 end if 1142 1143 ! Look before you leap. 1144 npw_disk = Wfk%Hdr%npwarr(ik_ibz) 1145 nspinor_disk = Wfk%nspinor 1146 1147 nband_disk = Wfk%nband(ik_ibz,spin) 1148 ! there are several cases here, reading in fewer than mband bands, 1149 ! or possibly more than you have allocated, and truncating 1150 ! nband_disk_keep could be used to distinguish these cases 1151 nband_disk_keep = nband_disk 1152 if (present(occ_k)) then 1153 nband_disk_keep = min(nband_disk, size(occ_k)) 1154 end if 1155 if (present(eig_k)) then 1156 if (Wfk%formeig == 0) then 1157 nband_disk_keep = min( nband_disk_keep, size(eig_k) ) 1158 else if (Wfk%formeig == 1) then 1159 nband_disk_keep = min( nband_disk_keep, int(sqrt(size(eig_k)/two)) ) 1160 end if 1161 end if 1162 1163 nb_block = (band_block(2) - band_block(1) + 1) 1164 ABI_CHECK(nb_block >0, "nband <=0") 1165 npw_tot_disk = npw_disk * nspinor_disk * nb_block 1166 1167 if (present(kg_k)) then 1168 ABI_CHECK(SIZE(kg_k,DIM=2) >= npw_disk,"kg_k too small") 1169 kg_k = zero 1170 end if 1171 if (present(cg_k)) then 1172 ABI_CHECK(SIZE(cg_k, DIM=2) >= npw_tot_disk,"cg_k too small") 1173 cg_k = zero 1174 end if 1175 1176 if (present(eig_k)) then 1177 if (Wfk%formeig==0) then 1178 ABI_CHECK(SIZE(eig_k) >= nband_disk, "GS eig_k too small") 1179 else if (Wfk%formeig==1) then 1180 ABI_CHECK(SIZE(eig_k) >= 2*nband_disk**2, "DFPT eig_k too small") 1181 else 1182 ABI_ERROR("formeig != [0,1]") 1183 end if 1184 end if 1185 1186 if (present(occ_k)) then 1187 !ABI_CHECK(SIZE(occ_k) <= nband_disk, "GS occ_k too large, not enough data on disk") 1188 if (Wfk%formeig==1) then 1189 ABI_ERROR("occ_k cannot be used when formeig ==1") 1190 end if 1191 end if 1192 1193 select case (Wfk%iomode) 1194 case (IO_MODE_FORTRAN) 1195 1196 ! Rewind the file to have the correct (k,s) block (if needed) 1197 call wfk_seek(Wfk,ik_ibz,spin) 1198 ! 1199 ! Read the first record: npw, nspinor, nband_disk 1200 read(Wfk%fh, err=10, iomsg=errmsg) npw_read, nspinor_read, nband_read 1201 1202 if (any( [npw_read, nspinor_read, nband_read] /= [npw_disk, nspinor_disk, nband_disk])) then 1203 write(msg,"(a,6(i0,2x))")"Mismatch between (npw, nspinor, nband) read from WFK and those found in HDR ",& 1204 npw_read, nspinor_read, nband_read, npw_disk, nspinor_disk, nband_disk 1205 ABI_ERROR(msg) 1206 end if 1207 1208 ! The second record: (k+G) vectors 1209 if (present(kg_k)) then 1210 read(Wfk%fh, err=10, iomsg=errmsg) kg_k(1:3,1:npw_disk) 1211 else 1212 read(Wfk%fh, err=10, iomsg=errmsg) ! kg_k(1:3,1:npw_disk) 1213 end if 1214 1215 select case (Wfk%formeig) 1216 case (0) 1217 ! The third record: eigenvalues and occupation factors. 1218 ! write(unitwf) (eigen(iband),iband=1,nband_disk),(occ(iband),iband=1,nband_disk) 1219 if (present(eig_k) .or. present(occ_k)) then 1220 ABI_MALLOC(tmp_eigk, (nband_disk)) 1221 ABI_MALLOC(tmp_occk, (nband_disk)) 1222 1223 read(Wfk%fh, err=10, iomsg=errmsg) tmp_eigk, tmp_occk 1224 1225 if (present(eig_k)) then 1226 eig_k = zero 1227 eig_k(1:nband_disk_keep) = tmp_eigk(1:nband_disk_keep) 1228 end if 1229 if (present(occ_k)) then 1230 occ_k = zero 1231 occ_k(1:nband_disk_keep) = tmp_occk(1:nband_disk_keep) 1232 end if 1233 1234 ABI_FREE(tmp_eigk) 1235 ABI_FREE(tmp_occk) 1236 1237 else 1238 read(Wfk%fh, err=10, iomsg=errmsg) ! eig_k(1:nband_disk), occ_k(1:nband_k) 1239 end if 1240 1241 ! Fourth record with the wave-functions. 1242 if (present(cg_k)) then 1243 npwso = npw_disk*nspinor_disk 1244 my_bcount = 0 1245 do band=1,nband_disk 1246 if (band >= band_block(1) .and. band <= band_block(2)) then 1247 ipw = my_bcount * npwso 1248 my_bcount = my_bcount + 1 1249 read(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 1250 else 1251 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 1252 end if 1253 end do 1254 1255 else 1256 do band=1,nband_disk 1257 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 1258 end do 1259 end if 1260 1261 case (1) 1262 ! formeig 1 for DFPT WF file 1263 npwso = npw_disk*nspinor_disk 1264 my_bcount = 0 1265 do band=1,nband_disk 1266 1267 if (present(eig_k)) then 1268 ! Read column matrix of size (2*nband_k) 1269 base = 2*(band-1)*nband_disk 1270 read(Wfk%fh, err=10, iomsg=errmsg) eig_k(base+1:base+2*nband_disk) 1271 else 1272 read(Wfk%fh, err=10, iomsg=errmsg ) ! eig_k(2*nband_disk) 1273 end if 1274 1275 if (present(cg_k) .and. (band >= band_block(1) .and. band <= band_block(2)) ) then 1276 ipw = my_bcount * npwso 1277 my_bcount = my_bcount + 1 1278 read(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 1279 else 1280 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 1281 end if 1282 end do 1283 1284 case default 1285 ABI_ERROR("formeig != [0,1]") 1286 end select 1287 1288 ! Reached the end of the (k,s) block. Update f90_fptr 1289 if (ik_ibz < Wfk%nkpt) then 1290 Wfk%f90_fptr = [ik_ibz+1,spin,REC_NPW] 1291 else 1292 ABI_CHECK(ik_ibz == wfk%nkpt, "ik_ibz != nkpt") 1293 if (spin==Wfk%nsppol) then 1294 Wfk%f90_fptr = FPTR_EOF ! EOF condition 1295 else 1296 Wfk%f90_fptr = [1,spin+1,REC_NPW] 1297 end if 1298 end if 1299 1300 #ifdef HAVE_MPI_IO 1301 case (IO_MODE_MPI) 1302 if (present(kg_k)) then 1303 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_KG) + xmpio_bsize_frm 1304 1305 call mpio_read_kg_k(Wfk%fh,my_offset,npw_disk,sc_mode,kg_k,mpierr) 1306 ABI_CHECK_MPI(mpierr, "reading kg") 1307 end if 1308 1309 ! formeig=0 => Read both eig and occ in tmp_eigk 1310 ! formeig=1 => Read (nband_k,nband_k) matrix of complex numbers. 1311 select case (Wfk%formeig) 1312 case (0) 1313 if (present(eig_k) .or. present(occ_k)) then 1314 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) + xmpio_bsize_frm 1315 1316 call mpio_read_eigocc_k(Wfk%fh,my_offset,nband_disk,Wfk%formeig,sc_mode,tmp_eigk,mpierr) 1317 ABI_CHECK_MPI(mpierr, "reading eigocc") 1318 1319 if (present(eig_k)) then 1320 eig_k(1:nband_disk_keep) = tmp_eigk(1:nband_disk_keep) 1321 end if 1322 if (present(occ_k)) then 1323 occ_k(1:nband_disk_keep) = tmp_eigk(nband_disk+1:nband_disk+nband_disk_keep) 1324 end if 1325 1326 ABI_FREE(tmp_eigk) 1327 end if 1328 1329 if (present(cg_k)) then 1330 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) 1331 sizes = [npw_disk*nspinor_disk, nband_disk] 1332 subsizes = [npw_disk*nspinor_disk, band_block(2)-band_block(1)+1] 1333 bufsz = 2 * npw_disk * nspinor_disk * nb_block 1334 starts = [1, band_block(1)] 1335 1336 call mpiotk_read_fsuba_dp2D(Wfk%fh,my_offset,sizes,subsizes,starts,bufsz,cg_k,& 1337 wfk%chunk_bsize,sc_mode,Wfk%comm,ierr) 1338 ABI_CHECK(ierr==0,"Fortran record too big") 1339 end if 1340 1341 case (1) 1342 if (present(eig_k)) then 1343 sizes = [nband_disk, npw_disk*nspinor_disk] 1344 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 1345 1346 call xmpio_create_fstripes(nband_disk,sizes,types,gkk_type,my_offpad,mpierr) 1347 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 1348 1349 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) + my_offpad 1350 1351 call MPI_FILE_SET_VIEW(Wfk%fh,my_offset,MPI_BYTE,gkk_type,'native',xmpio_info,mpierr) 1352 ABI_CHECK_MPI(mpierr, "SET_VIEW") 1353 call MPI_TYPE_FREE(gkk_type,mpierr) 1354 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 1355 1356 bufsz = nband_disk**2 1357 1358 if (sc_mode==xmpio_collective) then 1359 call MPI_FILE_READ_ALL(Wfk%fh,eig_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1360 else if (sc_mode==xmpio_single) then 1361 call MPI_FILE_READ(Wfk%fh,eig_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1362 else 1363 ABI_ERROR("Wrong sc_mode") 1364 end if 1365 ABI_CHECK_MPI(mpierr, "FILE_READ") 1366 end if 1367 1368 if (present(cg_k)) then 1369 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 1370 sizes = [npw_disk*nspinor_disk, nband_disk] 1371 1372 call xmpio_create_fstripes(nb_block,sizes,types,cgblock_type,my_offpad,mpierr) 1373 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 1374 1375 ! Increment my_offset to account for previous eigen and cg records if band_block(1) != 1 1376 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) + my_offpad 1377 my_offset = my_offset + (band_block(1) - 1) * ( & 1378 (2*npw_disk*wfk%nspinor*xmpi_bsize_dp + 2*xmpio_bsize_frm) + & 1379 (2*nband_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) ) 1380 1381 call MPI_FILE_SET_VIEW(Wfk%fh,my_offset,MPI_BYTE,cgblock_type,'native',xmpio_info,mpierr) 1382 ABI_CHECK_MPI(mpierr, "SET_VIEW") 1383 1384 call MPI_TYPE_FREE(cgblock_type,mpierr) 1385 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 1386 1387 bufsz = npw_disk * nspinor_disk * nb_block 1388 1389 if (sc_mode==xmpio_collective) then 1390 call MPI_FILE_READ_ALL(Wfk%fh,cg_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1391 else if (sc_mode==xmpio_single) then 1392 call MPI_FILE_READ(Wfk%fh,cg_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1393 else 1394 ABI_ERROR("Wrong sc_mode") 1395 end if 1396 ABI_CHECK_MPI(mpierr, "FILE_READ") 1397 end if 1398 1399 case default 1400 ABI_ERROR("formeig != [0,1]") 1401 end select 1402 #endif 1403 1404 case (IO_MODE_ETSF) 1405 if (present(kg_k)) then 1406 ! Read the reduced_coordinates_of_plane_waves for this k point. 1407 NCF_CHECK(nf90_inq_varid(wfk%fh, "reduced_coordinates_of_plane_waves", kg_varid)) 1408 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1409 NCF_CHECK(nctk_set_collective(wfk%fh, kg_varid)) 1410 end if 1411 ncerr = nf90_get_var(wfk%fh, kg_varid, kg_k, start=[1,1,ik_ibz], count=[3,npw_disk,1]) 1412 NCF_CHECK(ncerr) 1413 end if 1414 1415 if (Wfk%formeig==0) then 1416 ! Read eigenvalues and occupations. 1417 if (present(eig_k)) then 1418 NCF_CHECK(nf90_inq_varid(wfk%fh, "eigenvalues", eig_varid)) 1419 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1420 NCF_CHECK(nctk_set_collective(wfk%fh, eig_varid)) 1421 end if 1422 ncerr = nf90_get_var(wfk%fh, eig_varid, eig_k, start=[1,ik_ibz,spin], count=[nband_disk_keep,1,1]) 1423 NCF_CHECK(ncerr) 1424 end if 1425 1426 if (present(occ_k)) then 1427 NCF_CHECK(nf90_inq_varid(wfk%fh, "occupations", occ_varid)) 1428 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1429 NCF_CHECK(nctk_set_collective(wfk%fh, occ_varid)) 1430 end if 1431 ncerr = nf90_get_var(wfk%fh, occ_varid, occ_k, start=[1,ik_ibz,spin], count=[nband_disk,1,1]) 1432 NCF_CHECK_MSG(ncerr, "getting occ_k") 1433 end if 1434 1435 else ! formeig == 1 1436 1437 if (present(eig_k)) then 1438 ! Read h1 matrix elements. The netcdf array has shape: 1439 ! [complex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins] 1440 ABI_MALLOC(h1mat, (2, wfk%mband, wfk%mband)) 1441 NCF_CHECK(nf90_inq_varid(wfk%fh, "h1_matrix_elements", h1_varid)) 1442 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1443 NCF_CHECK(nctk_set_collective(wfk%fh, h1_varid)) 1444 end if 1445 ncerr = nf90_get_var(wfk%fh, h1_varid, h1mat, start=[1,1,1,ik_ibz,spin], & 1446 count=[2, wfk%mband, wfk%mband, 1, 1]) 1447 NCF_CHECK_MSG(ncerr, "getting h1mat_k") 1448 1449 ! For legacy reasons, I have to pack nband_k**2 elements in the first positions in eig_k 1450 ! This is important only if nband(:) depends on k i.e. mband != nband_k 1451 idx=1 1452 do ib2=1,nband_disk 1453 do ib1=1,nband_disk 1454 eig_k(idx:idx+1) = h1mat(:2,ib1,ib2) 1455 idx=idx+2 1456 end do 1457 end do 1458 ABI_FREE(h1mat) 1459 end if 1460 end if 1461 1462 if (present(cg_k)) then 1463 ! Read the nb_block bands starting from band_block(1) 1464 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 1465 NCF_CHECK(nf90_inq_varid(wfk%fh, "coefficients_of_wavefunctions", cg_varid)) 1466 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1467 NCF_CHECK(nctk_set_collective(wfk%fh, cg_varid)) 1468 end if 1469 1470 ncerr = nf90_get_var(wfk%fh, cg_varid, cg_k, start=[1,1,1,band_block(1),ik_ibz,spin], & 1471 count=[2,npw_disk,wfk%nspinor,nb_block,1,1]) 1472 NCF_CHECK_MSG(ncerr, "getting cg_k") 1473 end if 1474 1475 case default 1476 ABI_ERROR(sjoin('Wrong/unsupported iomode: ', itoa(Wfk%iomode))) 1477 end select 1478 1479 DBG_EXIT("COLL") 1480 1481 return 1482 1483 ! Handle Fortran IO error 1484 10 continue 1485 ABI_ERROR(errmsg) 1486 1487 end subroutine wfk_read_band_block
m_wfk/wfk_read_bks [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_bks
FUNCTION
Read the wavefunction and, optionally, the *DFPT* matrix elements for a given (band, k-point, spin).
INPUTS
Wfk<class(wfk_t)>=WFK file handler. band=Band index ik_ibz=Index of the k-point in the IBZ. spin=Spin index sc_mode= MPI-IO option xmpio_single ==> for reading by current proc. xmpio_collective ==> for collective reading. OUTPUTS cg_bks(2,npw_k*nspinor) = Fourier coefficients of the wavefunction [eig1_bks(2*wfk%mband)] = Matrix elements of the DFPT H1 Hamiltonian at the specified (k, spin).
SOURCE
1515 subroutine wfk_read_bks(wfk, band, ik_ibz, spin, sc_mode, cg_bks, eig1_bks) 1516 1517 !Arguments ------------------------------------ 1518 !scalars 1519 integer,intent(in) :: band,ik_ibz,spin,sc_mode 1520 class(wfk_t),intent(inout) :: wfk 1521 !arrays 1522 real(dp),DEV_CONTARRD intent(out) :: cg_bks(:,:) 1523 real(dp),optional,intent(inout) :: eig1_bks(2*wfk%mband) 1524 1525 !Local variables------------------------------- 1526 !scalars 1527 integer :: start,ib1,nspinor_disk,npw_disk,nband_disk !ierr, 1528 #ifdef HAVE_MPI_IO 1529 integer :: mpierr,bufsz,gkk_type !,cg_type 1530 integer(XMPI_OFFSET_KIND) :: my_offset,my_offpad 1531 integer :: sizes(2),types(2) 1532 #endif 1533 integer :: h1_varid,cg_varid,ncerr 1534 character(len=500) :: errmsg 1535 !arrays 1536 real(dp),allocatable :: all_eigk(:) 1537 1538 !************************************************************************ 1539 1540 if (wfk_validate_ks(wfk, ik_ibz, spin, band=band) /= 0) then 1541 ABI_ERROR("Wrong (ik_ibz, spin, band) args, Aborting now") 1542 end if 1543 npw_disk = wfk%Hdr%npwarr(ik_ibz) 1544 nband_disk = wfk%nband(ik_ibz, spin) 1545 nspinor_disk = wfk%nspinor 1546 !cg_bks = one; if (present(eig1_bks)) eig1_bks = zero ; return 1547 1548 if (.not. present(eig1_bks)) then 1549 call wfk%read_band_block([band, band], ik_ibz, spin, sc_mode, cg_k=cg_bks) 1550 return 1551 1552 else 1553 ABI_CHECK(wfk%formeig==1, "formeig must be 1 if eig1_bks is present") 1554 ABI_CHECK(size(cg_bks, dim=2) >= npw_disk*wfk%nspinor,"cg_bks too small") 1555 ABI_CHECK(size(eig1_bks) >= 2*nband_disk, "eig1_bks too small") 1556 1557 if (.False.) then 1558 !if (.True.) then 1559 ! Due to the API of wfk_read_band_block, we have to read the full set of eigenvalues 1560 ! and then extract the relevant band 1561 ! TODO: Should write another routine to avoid doing that. 1562 ABI_MALLOC(all_eigk, (2*wfk%mband**2)) 1563 call wfk%read_band_block([band, band], ik_ibz, spin, sc_mode, cg_k=cg_bks, eig_k=all_eigk) 1564 1565 ! Find the index of the slice. 1566 ! Remember that data in all_eigk does not have a constant stride if nband_disk != mband 1567 start = (band-1)*2*nband_disk 1568 1569 !write(std_out,*)size(eig1_bks), nband_disk 1570 eig1_bks(1:2*nband_disk) = all_eigk(start+1:start+2*nband_disk) 1571 ABI_FREE(all_eigk) 1572 1573 else 1574 ! Improved version 1575 select case (wfk%iomode) 1576 case (IO_MODE_FORTRAN) 1577 ! This code is not optimal because I'm rewinding the (k, s) block 1578 ! at each call but I'm not gonna spend time on this because I should 1579 ! refactor a lot of stuff. Use MPI-IO or HDF5! 1580 call wfk_seek(wfk,ik_ibz,spin) 1581 1582 ! Read the first record: npw, nspinor, nband_disk 1583 read(Wfk%fh, err=10, iomsg=errmsg) !npw_read, nspinor_read, nband_read 1584 ! The second record: (k+G) vectors 1585 read(wfk%fh, err=10, iomsg=errmsg) ! kg_k(1:3,1:npw_disk) 1586 1587 do ib1=1,nband_disk 1588 if (ib1 /= band) then 1589 read(Wfk%fh, err=10, iomsg=errmsg) ! eig_k(base+1:base+2*nband_disk) 1590 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 1591 else 1592 if (present(eig1_bks)) then 1593 read(wfk%fh, err=10, iomsg=errmsg) eig1_bks(1:2*nband_disk) 1594 else 1595 read(wfk%fh, err=10, iomsg=errmsg) 1596 end if 1597 read(wfk%fh, err=10, iomsg=errmsg) cg_bks(:,:npw_disk*wfk%nspinor) 1598 end if 1599 end do 1600 1601 ! Reached the end of the (k,s) block. Update f90_fptr 1602 call wfk_update_f90ptr(wfk, ik_ibz, spin) 1603 1604 #ifdef HAVE_MPI_IO 1605 case (IO_MODE_MPI) 1606 ! MPI-IO operations are done with file views even for contiguous data of the same type. 1607 ! See NOTES at the beginning of this module. 1608 1609 sizes = [nband_disk, npw_disk*nspinor_disk] 1610 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 1611 1612 call xmpio_create_fstripes(1,sizes,types,gkk_type,my_offpad,mpierr) 1613 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 1614 1615 !call MPI_TYPE_CONTIGUOUS(nband_disk,MPI_DOUBLE_COMPLEX,gkk_type,mpierr) 1616 !ABI_CHECK_MPI(mpierr, "type_contigous") 1617 !call MPI_TYPE_COMMIT(gkk_type,mpierr) 1618 !ABI_CHECK_MPI(mpierr, "mpi_commit") 1619 !my_offpad = 0 1620 1621 ! Increment my_offset to account for previous (iband -1) bands. 1622 my_offset = wfk%offset_ks(ik_ibz,spin,REC_EIG) + my_offpad 1623 my_offset = my_offset + (band - 1) * ( & 1624 (2*npw_disk*wfk%nspinor*xmpi_bsize_dp + 2*xmpio_bsize_frm) + & 1625 (2*nband_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) ) 1626 1627 #if 1 1628 bufsz = nband_disk 1629 if (sc_mode==xmpio_collective) then 1630 call MPI_FILE_READ_AT_ALL(wfk%fh,my_offset,eig1_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1631 else 1632 call MPI_FILE_READ_AT(wfk%fh,my_offset,eig1_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1633 end if 1634 1635 ! Read the cg_ks(G) 1636 ! Increment my_offset to account for the previous eigen and cg records of (iband-1) bands. 1637 my_offset = wfk%offset_ks(ik_ibz,spin,REC_CG) + xmpio_bsize_frm 1638 my_offset = my_offset + (band - 1) * ( & 1639 (2*npw_disk*wfk%nspinor*xmpi_bsize_dp + 2*xmpio_bsize_frm) + & 1640 (2*nband_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) ) 1641 1642 bufsz = npw_disk * nspinor_disk 1643 if (sc_mode==xmpio_collective) then 1644 call MPI_FILE_READ_AT_ALL(wfk%fh,my_offset,cg_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1645 else 1646 call MPI_FILE_READ_AT(wfk%fh,my_offset,cg_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1647 end if 1648 1649 #else 1650 call MPI_FILE_SET_VIEW(wfk%fh,my_offset,MPI_BYTE,gkk_type,'native',xmpio_info,mpierr) 1651 ABI_CHECK_MPI(mpierr, "SET_VIEW") 1652 call MPI_TYPE_FREE(gkk_type, mpierr) 1653 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 1654 1655 bufsz = nband_disk 1656 if (sc_mode==xmpio_collective) then 1657 call MPI_FILE_READ_ALL(wfk%fh,eig1_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1658 else if (sc_mode==xmpio_single) then 1659 call MPI_FILE_READ(wfk%fh,eig1_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1660 else 1661 ABI_ERROR("Wrong sc_mode") 1662 end if 1663 ABI_CHECK_MPI(mpierr, "FILE_READ") 1664 1665 ! Read the cg_ks(G) 1666 ! Increment my_offset to account for the previous eigen and cg records of (iband-1) bands. 1667 my_offset = wfk%offset_ks(ik_ibz,spin,REC_CG) + xmpio_bsize_frm 1668 my_offset = my_offset + (band - 1) * ( & 1669 (2*npw_disk*wfk%nspinor*xmpi_bsize_dp + 2*xmpio_bsize_frm) + & 1670 (2*nband_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) ) 1671 1672 call MPI_TYPE_CONTIGUOUS(npw_disk*nspinor_disk,MPI_DOUBLE_COMPLEX,cg_type,mpierr) 1673 ABI_CHECK_MPI(mpierr, "type_contigous") 1674 call MPI_TYPE_COMMIT(cg_type,mpierr) 1675 ABI_CHECK_MPI(mpierr, "mpi_commit") 1676 1677 call MPI_FILE_SET_VIEW(wfk%fh,my_offset,MPI_BYTE,cg_type,'native',xmpio_info,mpierr) 1678 ABI_CHECK_MPI(mpierr, "SET_VIEW") 1679 call MPI_TYPE_FREE(cg_type, mpierr) 1680 ABI_CHECK_MPI(mpierr, "MPI_TYPE_FREE") 1681 1682 bufsz = npw_disk * nspinor_disk 1683 if (sc_mode==xmpio_collective) then 1684 call MPI_FILE_READ_ALL(wfk%fh,cg_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1685 else if (sc_mode==xmpio_single) then 1686 call MPI_FILE_READ(wfk%fh,cg_bks,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 1687 else 1688 ABI_ERROR("Wrong sc_mode") 1689 end if 1690 ABI_CHECK_MPI(mpierr, "FILE_READ") 1691 #endif 1692 #endif 1693 1694 case (IO_MODE_ETSF) 1695 ! Read h1 matrix elements. The netcdf array has shape: 1696 ! [complex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins] 1697 NCF_CHECK(nf90_inq_varid(wfk%fh, "h1_matrix_elements", h1_varid)) 1698 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1699 NCF_CHECK(nctk_set_collective(wfk%fh, h1_varid)) 1700 end if 1701 ncerr = nf90_get_var(wfk%fh, h1_varid, eig1_bks, start=[1,1,band,ik_ibz,spin], count=[2, nband_disk, 1, 1, 1]) 1702 NCF_CHECK_MSG(ncerr, "getting h1mat_k") 1703 1704 ! Read the wavefunction. The coefficients_of_wavefunctions on file have shape: 1705 ! [cplex, mpw, nspinor, mband, nkpt, nsppol] 1706 NCF_CHECK(nf90_inq_varid(wfk%fh, "coefficients_of_wavefunctions", cg_varid)) 1707 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 1708 NCF_CHECK(nctk_set_collective(wfk%fh, cg_varid)) 1709 end if 1710 1711 ncerr = nf90_get_var(wfk%fh, cg_varid, cg_bks, start=[1,1,1,band,ik_ibz,spin], & 1712 count=[2,npw_disk,wfk%nspinor,1,1,1]) 1713 NCF_CHECK_MSG(ncerr, "getting cg_k") 1714 1715 case default 1716 ABI_ERROR(sjoin('Wrong value for iomode:', itoa(Wfk%iomode))) 1717 end select 1718 end if 1719 1720 end if 1721 1722 return 1723 1724 ! Handle Fortran IO error 1725 10 continue 1726 ABI_ERROR(errmsg) 1727 1728 end subroutine wfk_read_bks
m_wfk/wfk_read_bmask [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_bmask
FUNCTION
Read a set of bands at a given k-point, spin. The bands to be read are specified by the logical mask `bmask`.
INPUTS
Wfk<class(wfk_t)>= ik_ibz=Index of the k-point in the IBZ. spin=Spin index sc_mode= MPI-IO option xmpio_single ==> for reading by current proc. xmpio_collective ==> for collective reading. OUTPUTS [kg_k=(:,:)] = G-vectors [cg_k(:,:)] = Fourier coefficients [eig_k(:)] = Eigenvectors [occ_k(:)] = Occupation
NOTES
The output arrays eig_k and occ_k contain the *full* set of eigenvalues and occupation factors stored in the file and are dimensioned with wfk%mband.
SOURCE
2194 subroutine wfk_read_bmask(Wfk, bmask, ik_ibz, spin, sc_mode, kg_k, cg_k, eig_k, occ_k) 2195 2196 !Arguments ------------------------------------ 2197 !scalars 2198 integer,intent(in) :: ik_ibz,spin,sc_mode 2199 class(wfk_t),intent(inout) :: Wfk 2200 !arrays 2201 logical,intent(in) :: bmask(Wfk%mband) 2202 integer,intent(out), DEV_CONTARRD optional :: kg_k(:,:) !(3,npw_k) 2203 real(dp),intent(out), DEV_CONTARRD optional :: cg_k(:,:) !(2,npw_k*nspinor*nband) 2204 real(dp),intent(out),optional :: eig_k((2*Wfk%mband)**Wfk%formeig*Wfk%mband) 2205 real(dp),intent(out),optional :: occ_k(Wfk%mband) 2206 2207 !Local variables------------------------------- 2208 !scalars 2209 integer :: npw_disk,nspinor_disk,nband_disk,ipw,my_bcount,cnt,npwso,npw_tot,pt1,pt2,band 2210 integer :: npw_read,nspinor_read,nband_read,nb_tot,ncount,my_bcnt,my_maxb,base,nb, ierr 2211 character(len=500) :: msg,errmsg 2212 !arrays 2213 real(dp),ABI_CONTIGUOUS pointer :: tmp_eigk(:),tmp_occk(:) 2214 integer :: mpierr,cgscatter_type,cg_type,method,block,nblocks,nbxblock 2215 integer :: bstart,bstop,bufsz,ugsz,brest,max_nband 2216 integer(XMPI_OFFSET_KIND) :: my_offset,base_ofs,my_offpad 2217 integer :: band_block(2),sizes(2),subsizes(2),starts(2),types(2) 2218 integer,allocatable :: block_length(:),block_type(:) 2219 integer(XMPI_ADDRESS_KIND),allocatable :: block_displ(:) 2220 real(dp),allocatable :: buffer(:,:) 2221 integer :: kg_varid,eig_varid,occ_varid,cg_varid,ncerr 2222 integer,allocatable :: blocks(:,:) 2223 2224 !************************************************************************ 2225 2226 DBG_ENTER("COLL") 2227 2228 ABI_CHECK_IEQ(Wfk%rw_mode, WFK_READMODE, "Wfk must be in READMODE") 2229 2230 !do band=1,wfk%mband 2231 ! if (.not. bmask(band)) continue 2232 ! if (wfk_validate_ks(wfk, ik_ibz, spin, band=band) /= 0) then 2233 ! ABI_ERROR("Wrong (ik_ibz, spin, band) args, Aborting now") 2234 ! end if 2235 !end if 2236 2237 ! Look before you leap. 2238 npw_disk = Wfk%Hdr%npwarr(ik_ibz) 2239 nspinor_disk = Wfk%nspinor 2240 nband_disk = Wfk%nband(ik_ibz,spin) 2241 nb_tot = COUNT(bmask) 2242 npw_tot = npw_disk * nspinor_disk * nb_tot 2243 2244 if (present(kg_k)) then 2245 ABI_CHECK((SIZE(kg_k,DIM=2) >= npw_disk),"kg_k too small") 2246 end if 2247 2248 if (present(cg_k)) then 2249 ABI_CHECK(SIZE(cg_k, DIM=2) >= npw_tot, "Too small cg_k") 2250 end if 2251 2252 if (present(eig_k)) then 2253 if (Wfk%formeig==0) then 2254 ABI_CHECK(SIZE(eig_k) >= nband_disk, "GS eig_k too small") 2255 else if (Wfk%formeig==1) then 2256 ABI_CHECK(SIZE(eig_k) >= 2*nband_disk**2, "DFPT eig_k too small") 2257 else 2258 ABI_ERROR("formeig != [0,1]") 2259 end if 2260 end if 2261 2262 if (present(occ_k)) then 2263 ABI_CHECK(Wfk%formeig==0,"occ_k with formeig != 0") 2264 ABI_CHECK(SIZE(occ_k) >= nband_disk, "GS eig_k too small") 2265 end if 2266 2267 select case (Wfk%iomode) 2268 case (IO_MODE_FORTRAN) 2269 2270 ! Rewind the file to have the correct (k,s) block (if needed) 2271 call wfk_seek(Wfk, ik_ibz, spin) 2272 2273 ! Read the first record: npw, nspinor, nband_disk 2274 read(Wfk%fh, err=10, iomsg=errmsg) npw_read, nspinor_read, nband_read 2275 2276 if (any([npw_read, nspinor_read, nband_read] /= [npw_disk, nspinor_disk, nband_disk])) then 2277 write(msg,"(a,6(i0,2x))")"Mismatch between (npw, nspinor, nband) read from WFK and those found in HDR ",& 2278 npw_read, nspinor_read, nband_read, npw_disk, nspinor_disk, nband_disk 2279 ABI_ERROR(msg) 2280 end if 2281 2282 ! The second record: (k+G) vectors 2283 if (present(kg_k)) then 2284 read(Wfk%fh, err=10, iomsg=errmsg) kg_k(1:3,1:npw_disk) 2285 else 2286 read(Wfk%fh, err=10, iomsg=errmsg) ! kg_k(1:3,1:npw_disk) 2287 end if 2288 2289 ! The third record: eigenvalues and occupation factors. 2290 if (Wfk%formeig == 0) then 2291 2292 if (present(eig_k) .or. present(occ_k)) then 2293 ABI_MALLOC(tmp_eigk, (nband_disk)) 2294 ABI_MALLOC(tmp_occk, (nband_disk)) 2295 2296 read(Wfk%fh, err=10, iomsg=errmsg) tmp_eigk, tmp_occk 2297 2298 if (present(eig_k)) eig_k = tmp_eigk 2299 if (present(occ_k)) occ_k = tmp_occk 2300 ABI_FREE(tmp_eigk) 2301 ABI_FREE(tmp_occk) 2302 2303 else 2304 read(Wfk%fh, err=10, iomsg=errmsg) ! eig_k(1:nband_disk) 2305 end if 2306 2307 ! The wave-functions. 2308 if (present(cg_k)) then 2309 npwso = npw_disk*nspinor_disk 2310 my_bcount = 0 2311 do band=1,nband_disk 2312 if (bmask(band)) then 2313 ipw = my_bcount * npwso 2314 my_bcount = my_bcount + 1 2315 read(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 2316 else 2317 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 2318 end if 2319 end do 2320 2321 else 2322 do band=1,nband_disk 2323 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 2324 end do 2325 end if 2326 2327 else if (Wfk%formeig == 1) then 2328 ! Read matrix of size (2*nband_k**2) 2329 npwso = npw_disk*nspinor_disk 2330 my_bcount = 0 2331 2332 do band=1,nband_disk 2333 base = 2*(band-1)*nband_disk 2334 if (present(eig_k)) then 2335 read(Wfk%fh, err=10, iomsg=errmsg) eig_k(base+1:base+2*nband_disk) 2336 else 2337 read(Wfk%fh, err=10, iomsg=errmsg) ! eig_k(base+1:base+2*nband_disk) 2338 end if 2339 2340 if (bmask(band) .and. present(cg_k)) then 2341 ipw = my_bcount * npwso 2342 my_bcount = my_bcount + 1 2343 read(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 2344 else 2345 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 2346 end if 2347 end do 2348 2349 else 2350 ABI_ERROR("formeig != [0,1]") 2351 end if 2352 2353 ! Reached the end of the (k,s) block. Update f90_fptr 2354 call wfk_update_f90ptr(wfk, ik_ibz, spin) 2355 2356 #ifdef HAVE_MPI_IO 2357 case (IO_MODE_MPI) 2358 if (present(kg_k)) then 2359 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_KG) + xmpio_bsize_frm 2360 2361 call mpio_read_kg_k(Wfk%fh,my_offset,npw_disk,sc_mode,kg_k,mpierr) 2362 ABI_CHECK_MPI(mpierr, "mpio_read_kg_k") 2363 end if 2364 2365 ! The third record: eigenvalues and occupation factors. 2366 if (present(eig_k) .or. present(occ_k)) then 2367 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) + xmpio_bsize_frm 2368 ! 2369 ! formeig=0 => Read both eig and occ in tmp_eigk. 2370 ! formeig=1 => Read (nband_k,nband_k) matrix of complex numbers. 2371 ! 2372 call mpio_read_eigocc_k(Wfk%fh,my_offset,nband_disk,Wfk%formeig,sc_mode,tmp_eigk,mpierr) 2373 ABI_CHECK_MPI(mpierr, "mpio_read_eigocc") 2374 2375 if (Wfk%formeig == 0) then 2376 if (present(eig_k)) eig_k(1:nband_disk) = tmp_eigk(1:nband_disk) 2377 if (present(occ_k)) occ_k(1:nband_disk) = tmp_eigk(nband_disk+1:) 2378 else if (Wfk%formeig == 1) then 2379 if (present(eig_k)) eig_k(1:2*nband_disk**2) = tmp_eigk(1:2*nband_disk**2) 2380 else 2381 ABI_ERROR("formeig not in [0,1]") 2382 end if 2383 2384 ABI_FREE(tmp_eigk) 2385 end if 2386 2387 if (present(cg_k)) then 2388 method = 0 2389 2390 select case (method) 2391 case (0) 2392 ! DATA SIEVING: 2393 ! read max_nband states in chuncks of nbxblock, then extract my states according to bmask. 2394 ! 2395 ! MAX number of bands read by the procs in the communicator 2396 my_maxb = nband_disk 2397 do band=nband_disk,1,-1 2398 if (bmask(band)) then 2399 my_maxb = band 2400 EXIT 2401 end if 2402 end do 2403 call xmpi_max(my_maxb, max_nband, Wfk%comm, mpierr) 2404 !max_nband = nband_disk 2405 ! 2406 ! MPI-IO crashes if we try to read a large number of bands in a single call. 2407 nbxblock = max_nband 2408 if ((two * npw_disk *nspinor_disk * nbxblock * xmpi_bsize_dp) > Wfk%chunk_bsize) then 2409 nbxblock = Wfk%chunk_bsize / (2*npw_disk*nspinor_disk*xmpi_bsize_dp) 2410 if (nbxblock == 0) nbxblock = 50 2411 end if 2412 !nbxblock = 2 2413 2414 nblocks = max_nband / nbxblock 2415 brest = MOD(max_nband, nbxblock) 2416 if (brest /= 0) nblocks = nblocks + 1 2417 2418 !write(std_out, *) "full_size:", 2 * npw_disk *nspinor_disk *nbxblock * xmpi_bsize_dp, Wfk%chunk_bsize 2419 !write(std_out,*)"in buffered bmask with nblocks:", nblocks, ", nbxblock: ", nbxblock 2420 2421 base_ofs = Wfk%offset_ks(ik_ibz, spin, REC_CG) 2422 sizes = [npw_disk * nspinor_disk, nband_disk] 2423 2424 my_bcnt = 0 ! index of my band in cg_k 2425 do block=1,nblocks 2426 bstart = 1 + (block-1) * nbxblock 2427 bstop = bstart + nbxblock - 1 2428 if (bstop > max_nband) bstop = max_nband 2429 nb = bstop - bstart + 1 2430 2431 ! Allocate and read the buffer 2432 ! Note that in the API calls we mix real and complex. 2433 ! bufsz is the size in terms of complex numbers. 2434 band_block = [bstart, bstop] 2435 ugsz = npw_disk*nspinor_disk 2436 bufsz = ugsz * (bstop - bstart + 1) 2437 !write(std_out,*)" bstart, bstop:", band_block 2438 ABI_MALLOC_OR_DIE(buffer, (2, bufsz), ierr) 2439 2440 ! Read the cg_ks(G). different versions depending on formeig 2441 if (wfk%formeig == 0) then 2442 subsizes = [npw_disk*nspinor_disk, band_block(2)-band_block(1)+1] 2443 starts = [1, bstart] 2444 2445 call mpiotk_read_fsuba_dp2D(Wfk%fh,base_ofs,sizes,subsizes,starts,& 2446 2 * bufsz,buffer,Wfk%chunk_bsize,sc_mode,Wfk%comm,mpierr) 2447 ABI_CHECK(mpierr == 0, "Fortran record too big") 2448 2449 ! New version based: master reads and broadcasts the buffer. 2450 !if (wfk%my_rank == 0) then 2451 ! call mpiotk_read_fsuba_dp2D(Wfk%fh,base_ofs,sizes,subsizes,starts,& 2452 ! 2 * bufsz,buffer,Wfk%chunk_bsize,xmpio_single,wfk%comm,mpierr) 2453 ! ABI_CHECK(mpierr == 0, "Fortran record too big") 2454 !end if 2455 !call xmpi_bcast(buffer, 0, wfk%comm, ierr) 2456 2457 else if (wfk%formeig == 1) then 2458 2459 ! Increment my_offset to account for the previous eigen and cg records of (iband-1) bands. 2460 my_offset = wfk%offset_ks(ik_ibz,spin,REC_CG) + xmpio_bsize_frm 2461 my_offset = my_offset + (bstart - 1) * ( & 2462 (2*npw_disk*wfk%nspinor*xmpi_bsize_dp + 2*xmpio_bsize_frm) + & 2463 (2*nband_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) ) 2464 2465 sizes = [npw_disk*nspinor_disk, nband_disk] 2466 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 2467 2468 call xmpio_create_fstripes(nb,sizes,types,cg_type,my_offpad,mpierr) 2469 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 2470 2471 call MPI_FILE_SET_VIEW(wfk%fh,my_offset,MPI_BYTE,cg_type,'native',xmpio_info,mpierr) 2472 ABI_CHECK_MPI(mpierr, "SET_VIEW") 2473 call MPI_TYPE_FREE(cg_type,mpierr) 2474 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 2475 2476 if (sc_mode == xmpio_collective) then 2477 call MPI_FILE_READ_ALL(wfk%fh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2478 else if (sc_mode == xmpio_single) then 2479 call MPI_FILE_READ(wfk%fh,buffer,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2480 else 2481 ABI_ERROR("Wrong sc_mode") 2482 end if 2483 ABI_CHECK_MPI(mpierr, "FILE_READ") 2484 end if 2485 2486 ! Extract my bands from buffer. 2487 do band=bstart,bstop 2488 if (bmask(band)) then 2489 my_bcnt = my_bcnt + 1 2490 pt1 = 1 + (my_bcnt - 1) * ugsz 2491 pt2 = 1 + (band - bstart) * ugsz 2492 cg_k(:,pt1:pt1+ugsz-1) = buffer(:,pt2:pt2+ugsz-1) 2493 end if 2494 end do 2495 2496 ABI_FREE(buffer) 2497 end do 2498 2499 case (1, 2) 2500 ABI_CHECK(wfk%formeig == 0, "formeig == 1 not coded") 2501 call MPI_TYPE_CONTIGUOUS(npw_disk*nspinor_disk,MPI_DOUBLE_COMPLEX,cg_type,mpierr) 2502 ABI_CHECK_MPI(mpierr, "type_contigous") 2503 2504 if (method == 1) then 2505 ncount = nb_tot 2506 ABI_MALLOC(block_length, (ncount+2)) 2507 ABI_MALLOC(block_type, (ncount+2)) 2508 ABI_MALLOC(block_displ, (ncount+2)) 2509 2510 block_length(1)=1 2511 block_displ (1)=0 2512 block_type (1)=MPI_LB 2513 2514 my_bcount = 1 2515 do band=1,Wfk%mband 2516 if (bmask(band)) then 2517 my_bcount = my_bcount + 1 2518 block_length(my_bcount) = 1 2519 block_type(my_bcount) = cg_type 2520 block_displ(my_bcount) = xmpio_bsize_frm + & 2521 (band-1) * (2*npw_disk*nspinor_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) 2522 end if 2523 end do 2524 2525 block_length(ncount+2) = 1 2526 block_displ (ncount+2) = block_displ(my_bcount) 2527 block_type (ncount+2) = MPI_UB 2528 2529 else if (method == 2) then 2530 ! this file view is not efficient but it's similar to the 2531 ! one used in wff_readwrite. Let's see if MPI-IO likes it! 2532 ncount = nb_tot* nspinor_disk * npw_disk 2533 2534 ABI_MALLOC(block_length, (ncount+2)) 2535 ABI_MALLOC(block_type, (ncount+2)) 2536 ABI_MALLOC(block_displ, (ncount+2)) 2537 2538 block_length(1)=1 2539 block_displ (1)=0 2540 block_type (1)=MPI_LB 2541 ! 2542 ! The view starts at REC_CG 2543 cnt = 1 2544 do band=1,Wfk%mband 2545 if (bmask(band)) then 2546 base_ofs = xmpio_bsize_frm + & 2547 (band-1) * (2*npw_disk*nspinor_disk*xmpi_bsize_dp + 2*xmpio_bsize_frm) 2548 do ipw=1,npw_disk*nspinor_disk 2549 cnt = cnt + 1 2550 block_length(cnt) = 1 2551 block_type(cnt) = MPI_DOUBLE_COMPLEX 2552 block_displ(cnt) = base_ofs + 2*(ipw-1)*xmpi_bsize_dp 2553 end do 2554 end if 2555 end do 2556 2557 block_length(ncount+2) = 1 2558 block_displ (ncount+2) = block_displ(cnt) 2559 block_type (ncount+2) = MPI_UB 2560 end if 2561 2562 call xmpio_type_struct(ncount+2,block_length,block_displ,block_type,cgscatter_type,mpierr) 2563 ABI_CHECK_MPI(mpierr, "type_struct") 2564 2565 ABI_FREE(block_length) 2566 ABI_FREE(block_type) 2567 ABI_FREE(block_displ) 2568 2569 call MPI_TYPE_FREE(cg_type, mpierr) 2570 ABI_CHECK_MPI(mpierr, "MPI_TYPE_FREE") 2571 2572 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) 2573 2574 call MPI_FILE_SET_VIEW(Wfk%fh, my_offset, MPI_BYTE, cgscatter_type, 'native', xmpio_info, mpierr) 2575 ABI_CHECK_MPI(mpierr, "SET_VIEW") 2576 2577 call MPI_TYPE_FREE(cgscatter_type, mpierr) 2578 ABI_CHECK_MPI(mpierr, "MPI_TYPE_FREE") 2579 2580 call MPI_FILE_READ_ALL(Wfk%fh, cg_k, npw_tot, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr) 2581 ABI_CHECK_MPI(mpierr, "FILE_READ_ALL") 2582 2583 case default 2584 ABI_ERROR("Wrong method") 2585 end select 2586 end if 2587 #endif 2588 2589 case (IO_MODE_ETSF) 2590 ABI_CHECK(wfk%formeig == 0, "formeig != 0 not coded") 2591 !write(std_out,*)"bmask: ",bmask 2592 2593 ! TODO: extract routines (see other similar calls) 2594 if (present(kg_k)) then 2595 ! Read the reduced_coordinates_of_plane_waves for this k point. 2596 NCF_CHECK(nf90_inq_varid(wfk%fh, "reduced_coordinates_of_plane_waves", kg_varid)) 2597 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2598 NCF_CHECK(nctk_set_collective(wfk%fh, kg_varid)) 2599 end if 2600 2601 ncerr = nf90_get_var(wfk%fh, kg_varid, kg_k, start=[1,1,ik_ibz], count=[3,npw_disk,1]) 2602 NCF_CHECK(ncerr) 2603 end if 2604 2605 ! Read eigenvalues and occupations. 2606 if (Wfk%formeig==0) then 2607 if (present(eig_k)) then 2608 NCF_CHECK(nf90_inq_varid(wfk%fh, "eigenvalues", eig_varid)) 2609 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2610 NCF_CHECK(nctk_set_collective(wfk%fh, eig_varid)) 2611 end if 2612 2613 ncerr = nf90_get_var(wfk%fh, eig_varid, eig_k, start=[1,ik_ibz,spin], count=[nband_disk,1,1]) 2614 NCF_CHECK(ncerr) 2615 end if 2616 2617 if (present(occ_k)) then 2618 NCF_CHECK(nf90_inq_varid(wfk%fh, "occupations", occ_varid)) 2619 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2620 NCF_CHECK(nctk_set_collective(wfk%fh, occ_varid)) 2621 end if 2622 2623 ncerr = nf90_get_var(wfk%fh, occ_varid, occ_k, start=[1,ik_ibz,spin], count=[nband_disk,1,1]) 2624 NCF_CHECK_MSG(ncerr, "getting occ_k") 2625 end if 2626 else 2627 ABI_ERROR("formeig !=0 not compatible with ETSF-IO") 2628 end if 2629 2630 if (present(cg_k)) then 2631 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 2632 NCF_CHECK(nf90_inq_varid(wfk%fh, "coefficients_of_wavefunctions", cg_varid)) 2633 ! TODO: Collective 2634 !if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2635 ! NCF_CHECK(nctk_set_collective(wfk%fh, cg_varid)) 2636 !end if 2637 #if 0 2638 ! Simple and very inefficient version for debugging. 2639 ipw = 1 2640 do band=1,wfk%mband 2641 if (.not. bmask(band)) cycle 2642 ncerr = nf90_get_var(wfk%fh, cg_varid, cg_k(:,ipw:), start=[1,1,1,band,ik_ibz,spin], & 2643 count=[2,npw_disk,wfk%nspinor,1,1,1]) 2644 NCF_CHECK_MSG(ncerr, "getting cg_k block") 2645 ipw = ipw + wfk%nspinor * npw_disk 2646 end do 2647 #else 2648 ! Read bands in blocks defined by bmask. 2649 ! be careful when in collective mode because processors may call the routine with nblocks==0 2650 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 2651 call mask2blocks(bmask, nblocks, blocks) 2652 !ABI_CHECK(nblocks /= 0, "nblocks==0") 2653 2654 ipw = 1 2655 do block=1,nblocks 2656 band_block = blocks(:,block) 2657 nb = band_block(2) - band_block(1) + 1 2658 ncerr = nf90_get_var(wfk%fh, cg_varid, cg_k(:,ipw:), start=[1,1,1,band_block(1),ik_ibz,spin], & 2659 count=[2,npw_disk,wfk%nspinor,nb,1,1]) 2660 NCF_CHECK_MSG(ncerr, "getting cg_k block") 2661 ipw = ipw + wfk%nspinor * npw_disk * nb 2662 end do 2663 ABI_FREE(blocks) 2664 #endif 2665 2666 ! Prototype for collective version. 2667 !min_band = lfind(bmask); if min_ 2668 !max_band = lfind(bmask, back=.True.) ! n+1 2669 !! TODO: xmpi_min_max 2670 !call xmpi_max(my_min_band, min_band, comm_cell, ierr) 2671 !call xmpi_min(my_min_band, max_band, comm_cell, ierr) 2672 !nb = max_band - min_band + 1 2673 !ncalls = nb / 2674 2675 !NCF_CHECK(nf90_var_par_access(ncid, cg_varid, nf90_collective)) 2676 !do block=1,nblocks 2677 ! ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 2678 ! ncerr = nf90_get_var(wfk%fh, cg_varid, cg_k, start=[1,1,1,band_block(1),ik_ibz,spin], & 2679 ! count=[2,npw_disk,wfk%nspinor,band_block(2)-band_block(1)+1,1,1]) 2680 ! NCF_CHECK_MSG(ncerr, "getting cg_k block") 2681 ! do band=band_start,band_end 2682 ! if (.not. bmask(band)) cycle 2683 ! cg_k(:,:) = 2684 ! end do 2685 !end do 2686 end if 2687 2688 case default 2689 ABI_ERROR(sjoin('Wrong/unsupported value of iomode: ', itoa(wfk%iomode))) 2690 end select 2691 2692 DBG_EXIT("COLL") 2693 2694 return 2695 2696 ! Handle Fortran IO error 2697 10 continue 2698 ABI_ERROR(errmsg) 2699 2700 end subroutine wfk_read_bmask
m_wfk/wfk_read_ebands [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_ebands
FUNCTION
Read the GS eigenvalues and return ebands_t object.
INPUTS
path=WFK file name comm=MPI communicator OUTPUTS ebands<ebands_t>=GS band-structure. [out_hdr]=Abinit header.
SOURCE
2722 type(ebands_t) function wfk_read_ebands(path, comm, out_hdr) result(ebands) 2723 2724 !Arguments ------------------------------------ 2725 !scalars 2726 character(len=*),intent(in) :: path 2727 integer,intent(in) :: comm 2728 type(hdr_type),optional,intent(inout) :: out_hdr ! ifort and others are buggy for optional intent(out) structured types 2729 2730 !Local variables------------------------------- 2731 !scalars 2732 type(hdr_type) :: hdr 2733 !arrays 2734 real(dp),pointer :: eigen(:,:,:) 2735 2736 !************************************************************************ 2737 2738 call wfk_read_eigenvalues(path, eigen, hdr, comm) 2739 ebands = ebands_from_hdr(hdr, maxval(hdr%nband), eigen) 2740 if (present(out_hdr)) call hdr_copy(hdr, out_hdr) 2741 2742 ABI_FREE(eigen) 2743 call hdr%free() 2744 2745 end function wfk_read_ebands
m_wfk/wfk_read_eigenvalues [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_eigenvalues
FUNCTION
Read all the GS eigenvalues stored in the WFK file fname.
INPUTS
fname=Name of the file comm=MPI communicator. OUTPUTS eigen = In input: nullified pointer In output: eigen(mband,nkpt,nsppol) contains the GS eigevalues. Hdr_out<hdr_type>=The header of the file
SOURCE
2822 subroutine wfk_read_eigenvalues(fname, eigen, Hdr_out, comm, occ) 2823 2824 !Arguments ------------------------------------ 2825 !scalars 2826 integer,intent(in) :: comm 2827 character(len=*),intent(in) :: fname 2828 type(hdr_type),intent(out) :: Hdr_out 2829 !arrays 2830 !TODO: Replace pointers with allocatable. 2831 real(dp),pointer :: eigen(:,:,:) 2832 real(dp),pointer,optional :: occ(:,:,:) 2833 2834 !Local variables------------------------------- 2835 !scalars 2836 integer,parameter :: master = 0, formeig0 = 0 2837 integer :: ik_ibz,spin,my_rank,ierr,iomode,funt,sc_mode,mband 2838 real(dp) :: cpu, wall, gflops, cpu_io, wall_io, gflops_io 2839 type(wfk_t) :: Wfk 2840 2841 !************************************************************************ 2842 2843 call cwtime(cpu, wall, gflops, "start") 2844 my_rank = xmpi_comm_rank(comm) 2845 iomode = iomode_from_fname(fname) 2846 2847 ! Should not do this but it seems that read_eigk with IO_MODE_MPI is very inefficient when 2848 ! we read files with lots of k-points. 2849 if (iomode == IO_MODE_MPI) iomode = IO_MODE_FORTRAN 2850 2851 call wrtout(std_out, sjoin(" Reading eigenvalues from:", fname, ", with iomode:", iomode2str(iomode))) 2852 2853 if (my_rank == master) then 2854 ! Master reads and broadcasts 2855 call cwtime(cpu_io, wall_io, gflops_io, "start") 2856 2857 ! Open the file. 2858 sc_mode = xmpio_single 2859 funt = get_unit() 2860 call wfk_open_read(Wfk, fname, formeig0, iomode, funt, xmpi_comm_self, Hdr_out=Hdr_out) 2861 2862 ! Read the eigenvalues and optionally the occupation factors. 2863 ABI_MALLOC(eigen, (Wfk%mband, Wfk%nkpt, Wfk%nsppol)) 2864 eigen = HUGE(zero) 2865 if (present(occ)) then 2866 ABI_MALLOC(occ, (Wfk%mband, Wfk%nkpt, Wfk%nsppol)) 2867 occ = HUGE(zero) 2868 end if 2869 2870 do spin=1,Wfk%nsppol 2871 do ik_ibz=1,Wfk%nkpt 2872 if (present(occ)) then 2873 call wfk%read_eigk(ik_ibz, spin, sc_mode, eigen(:,ik_ibz,spin), occ_k=occ(:,ik_ibz,spin)) 2874 else 2875 call wfk%read_eigk(ik_ibz, spin, sc_mode, eigen(:,ik_ibz,spin)) 2876 end if 2877 end do 2878 end do 2879 2880 ! Close the file. 2881 call wfk%close() 2882 call cwtime_report(" wfk_read_eigenvalues_io", cpu_io, wall_io, gflops_io) 2883 end if 2884 2885 ! Broadcast data 2886 if (xmpi_comm_size(comm) > 1) then 2887 call Hdr_out%bcast(master, my_rank, comm) 2888 mband = MAXVAL(Hdr_out%nband) 2889 if (my_rank /= master) then 2890 ABI_MALLOC(eigen, (mband, Hdr_out%nkpt, Hdr_out%nsppol)) 2891 if (present(occ)) then 2892 ABI_MALLOC(occ, (mband, Hdr_out%nkpt, Hdr_out%nsppol)) 2893 end if 2894 end if 2895 call xmpi_bcast(eigen, master, comm, ierr) 2896 if (present(occ)) call xmpi_bcast(occ, master, comm, ierr) 2897 end if 2898 2899 call cwtime_report(" wfk_read_eigenvalues", cpu, wall, gflops) 2900 2901 end subroutine wfk_read_eigenvalues
m_wfk/wfk_read_eigk [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_eigk
FUNCTION
Helper function to read all the eigenvalues for a given (k-point,spin)
INPUTS
Wfk<class(wfk_t)>= WFK file handler ik_ibz=Index of the k-point in the IBZ. spin=spin index sc_mode= MPI-IO option xmpio_single ==> for reading by current proc. xmpio_collective ==> for collective reading. OUTPUTS eig_k(1:nband_k) = GS Eigenvalues for the given (k,s) occ_k(1:nband_k) = Occupation factors for the given (k,s)
NOTES
The buffers eig_k and occ_k are dimensions with wfk%mband. The routine will fill the first nband_k positions with data read from file where nband_k is the number of bands on file i.e. wfk%nband(ik_ibz,spin)
SOURCE
2776 subroutine wfk_read_eigk(Wfk,ik_ibz,spin,sc_mode,eig_k,occ_k) 2777 2778 !Arguments ------------------------------------ 2779 !scalars 2780 integer,intent(in) :: ik_ibz,spin,sc_mode 2781 class(wfk_t),intent(inout) :: Wfk 2782 !arrays 2783 real(dp),intent(out) :: eig_k((2*Wfk%mband)**Wfk%formeig*Wfk%mband) 2784 real(dp),optional,intent(out) :: occ_k(Wfk%mband) 2785 2786 !Local variables------------------------------- 2787 !scalars 2788 integer,parameter :: band_block00(2) = [0, 0] 2789 2790 !************************************************************************ 2791 2792 if (present(occ_k)) then 2793 ABI_CHECK(Wfk%formeig == 0, "formeig !=0") 2794 call wfk%read_band_block(band_block00,ik_ibz,spin,sc_mode,eig_k=eig_k,occ_k=occ_k) 2795 else 2796 call wfk%read_band_block(band_block00,ik_ibz,spin,sc_mode,eig_k=eig_k) 2797 end if 2798 2799 end subroutine wfk_read_eigk
m_wfk/wfk_read_h1mat [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_h1mat
FUNCTION
Read all H1 matrix elements in the WFK file fname inside the MPI communicator comm.
INPUTS
path=File name comm=MPI communicator. OUTPUTS eigen(2*hdr_out%mband**2*hdr_out%nkpt*hdr_out%nsppol)=Array with the matrix elements of H1 packed in the first positions. The array is allocated by the procedure. Hdr_out<hdr_type>=The header of the file
SOURCE
3581 subroutine wfk_read_h1mat(fname, eigen, hdr_out, comm) 3582 3583 !Arguments ------------------------------------ 3584 !scalars 3585 character(len=*),intent(in) :: fname 3586 integer,intent(in) :: comm 3587 type(hdr_type),intent(out) :: Hdr_out 3588 !arrays 3589 real(dp),allocatable,intent(out) :: eigen(:) 3590 3591 !Local variables------------------------------- 3592 !scalars 3593 integer,parameter :: master=0,formeig1=1 3594 integer :: spin,ik_ibz,nband_k,ptr,ierr,iomode,mband,my_rank 3595 type(wfk_t) :: wfk 3596 !arrays 3597 integer,parameter :: band_block00(2)=[0,0] 3598 3599 !************************************************************************ 3600 3601 my_rank = xmpi_comm_rank(comm) 3602 3603 if (my_rank==master) then 3604 ! Open the file. 3605 iomode = iomode_from_fname(fname) 3606 call wfk_open_read(wfk, fname, formeig1, iomode, get_unit(), xmpi_comm_self, hdr_out=hdr_out) 3607 3608 ! Read h1 mat and pack them in the first positions. 3609 ABI_MALLOC(eigen, (2*wfk%mband**2*wfk%nkpt*wfk%nsppol)) 3610 3611 ptr=1 3612 do spin=1,wfk%nsppol 3613 do ik_ibz=1,wfk%nkpt 3614 nband_k = wfk%nband(ik_ibz,spin) 3615 call wfk_read_band_block(wfk, band_block00, ik_ibz, spin, xmpio_single, eig_k=eigen(ptr:)) 3616 ptr = ptr + 2*nband_k**2 3617 end do 3618 end do 3619 3620 call wfk%close() 3621 end if 3622 3623 ! Broadcast data 3624 if (xmpi_comm_size(comm) > 1) then 3625 call hdr_out%bcast(master, my_rank, comm) 3626 3627 mband = maxval(Hdr_out%nband) 3628 if (my_rank/=master) then 3629 ABI_MALLOC(eigen, (2*mband**2*hdr_out%nkpt*hdr_out%nsppol)) 3630 end if 3631 call xmpi_bcast(eigen,master,comm,ierr) 3632 end if 3633 3634 end subroutine wfk_read_h1mat
m_wfk/wfk_read_my_kptbands [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_read_my_kptbands
FUNCTION
Fill a cg (kg, eigen, occ) array with wavefunctions in a given BZ based on a distribution of k, b, s attributed to present processor
INPUTS
inpath_ = file name distrb_flags = logical mask for band, k, spins on this processor comm = mpi communicator formeig = flag for GS or response function format of eigenvalues istwfk_in = reciprocal space storage (reduced PW sphere or not) kptns_in = requested k points, to be extracted from file or completed nkpt_in = number of requested k npwarr = array of number of plane waves at each k istwfk_in = storage flag for plane waves mcg = max size of cg array mband_in = max number of bands over all k mband_mem_in = max number of bands stored on each processor nkpt_in = total number of k-points nspinor_in = number of spinor components 1 or 2 nsppol_in = number of spin polarization channels usepaw_in = enable PAW or not? (1/0)
OUTPUT
cg = plane wave coefficients kg = plane wave coordinates eigen = eigenvectors at all bands and my k occ = occupations of all bands at my k pawrhoij = PAW matrix elements in projectors
SOURCE
2941 subroutine wfk_read_my_kptbands(inpath_, distrb_flags, comm, ecut_eff_in, & 2942 formeig, istwfk_in, kptns_in, mcg, mband_in, mband_mem_in, mkmem_in, mpw_in, & 2943 natom_in, nkpt_in, npwarr, nspinor_in, nsppol_in, usepaw_in, & 2944 cg, kg, eigen, occ, pawrhoij, ask_accurate_) 2945 2946 !Arguments ------------------------------------ 2947 !scalars 2948 integer, intent(in) :: comm, nkpt_in, formeig 2949 integer, intent(in) :: mcg, mpw_in, mkmem_in 2950 integer, intent(in) :: mband_in, mband_mem_in, natom_in, nspinor_in, nsppol_in, usepaw_in 2951 real(dp), intent(in) :: ecut_eff_in 2952 !=dtset%ecut*(dtset%dilatmx)**2 ! ecut * dilatmx**2 2953 !arrays 2954 integer, intent(in) :: istwfk_in(nkpt_in) 2955 integer, intent(in) :: npwarr(nkpt_in) 2956 character(len=fnlen), intent(in) :: inpath_ 2957 logical, intent(in) :: distrb_flags(nkpt_in,mband_in,nsppol_in) 2958 real(dp), intent(in),target :: kptns_in(3,nkpt_in) 2959 real(dp), intent(out) :: cg(2,mcg) 2960 integer, intent(out), optional :: kg(3,mpw_in*mkmem_in) 2961 real(dp), intent(out), optional :: eigen(mband_in*(2*mband_in)**formeig*nkpt_in*nsppol_in) 2962 real(dp), intent(out), optional :: occ(mband_in*nkpt_in*nsppol_in) 2963 type(pawrhoij_type),intent(inout),optional,target :: pawrhoij(natom_in) 2964 integer, intent(in), optional :: ask_accurate_ 2965 2966 !Local variables------------------------------- 2967 !scalars 2968 integer,parameter :: formeig0 = 0, master = 0 2969 integer :: spin,ikf,ik_disk,nband_k,mpw_disk,mband,nspinor 2970 integer :: iomode,nsppol,isym,itimrev 2971 integer :: npw_disk,npw_kf,istwf_disk,istwf_kf 2972 integer :: ikpt,ii,jj,kk,ll,iqst,nqst 2973 integer :: ibdoff, ierr, my_rank 2974 integer :: wfk_unt, iband, nband_me, nband_me_disk 2975 integer :: nband_me_saved, iband_saved 2976 integer :: spin_saved, spin_sym 2977 integer :: mpierr 2978 integer :: ask_accurate, sppoldbl 2979 real(dp) :: cpu, wall, gflops 2980 real(dp) :: ecut_eff_disk 2981 real(dp) :: dksqmax 2982 character(len=fnlen) :: inpath 2983 logical :: isirred_kf 2984 logical :: needthisk 2985 logical :: convnsppol1to2 2986 type(wfk_t),target :: wfk_disk 2987 type(crystal_t) :: cryst 2988 !arrays 2989 integer :: g0(3),work_ngfft(18),gmax_disk(3),gmax_kf(3),gmax(3) 2990 integer,allocatable :: kg_kf(:,:), icg(:,:), ikg(:), ibdeig(:,:), ibdocc(:,:) 2991 integer,allocatable :: symrelT(:,:,:) 2992 integer,allocatable :: rbz2disk(:,:),kg_disk(:,:),iperm(:),rbz2disk_sort(:) 2993 real(dp) :: kf(3),k_disk(3), ksym(3) 2994 real(dp),allocatable :: cg_disk(:,:),eig_disk(:),occ_disk(:),work(:,:,:,:) 2995 2996 ! ************************************************************************* 2997 2998 call cwtime(cpu, wall, gflops, "start") 2999 3000 my_rank = xmpi_comm_rank(comm) 3001 3002 ! Master checks the existence of data file 3003 if (my_rank == master) then 3004 inpath = inpath_ 3005 if (.not. file_exists(inpath)) then 3006 ! Trick needed to run Abinit test suite in netcdf mode. 3007 if (file_exists(nctk_ncify(inpath))) then 3008 write(std_out,"(3a)")"- File: ",trim(inpath)," does not exist but found netcdf file with similar name." 3009 inpath = nctk_ncify(inpath) 3010 end if 3011 if (.not. file_exists(inpath)) then 3012 ABI_ERROR('Missing data file: '//TRIM(inpath)) 3013 end if 3014 end if 3015 end if 3016 3017 call xmpi_bcast(inpath, master, comm, ierr) 3018 call wrtout(std_out, sjoin(" About to read wavefunctions from:", inpath)) 3019 3020 ! now attack the cg reading 3021 iomode = iomode_from_fname(inpath) 3022 wfk_unt = get_unit() 3023 3024 ! TODO: this still does not read in parallel properly: 3025 ! if I use xmpi_comm_self only the mother thread gets eigen and cg 3026 ! if I use comm and MPIO_stuff then it hangs on this call 3027 ! if I impose FORTRAN_IO and xmpio_single it complains the file is already opened by another proc 3028 ABI_UNUSED(comm) 3029 call wfk_open_read(wfk_disk,inpath,formeig,iomode,wfk_unt,xmpi_comm_self) 3030 3031 if(present(eigen)) eigen = zero 3032 if(present(occ)) occ = zero 3033 if(present(kg)) kg = 0 3034 3035 ! this initialization is needed in case we read a file with fewer bands and only fill part of cg 3036 cg = zero 3037 3038 ! ABI_CHECK(wfk_disk%mband >= mband_in, "input mband too large for this file") 3039 mband = wfk_disk%mband; 3040 ABI_CHECK(wfk_disk%nspinor == nspinor_in, "input nspinor does not agree with file") 3041 nspinor = wfk_disk%nspinor 3042 !checks: impose each individual nband conserved wrt disk? 3043 3044 ABI_CHECK(wfk_disk%nsppol <= nsppol_in, "nsppol can not decrease when reading from disk") 3045 !ABI_CHECK(wfk_disk%nsppol == nsppol_in, "nsppol does not agree with file") 3046 nsppol = nsppol_in; 3047 convnsppol1to2=.false. 3048 if (wfk_disk%nsppol < nsppol_in) convnsppol1to2 = .true. 3049 3050 ! NB: npw can differ as can istwfk 3051 mpw_disk = maxval(wfk_disk%Hdr%npwarr) 3052 ecut_eff_disk = wfk_disk%hdr%ecut_eff ! ecut * dilatmx**2 3053 3054 ABI_MALLOC(kg_disk, (3, mpw_disk)) 3055 ABI_MALLOC(cg_disk, (2, mpw_disk*nspinor*mband_mem_in)) 3056 ABI_MALLOC(eig_disk, ((2*mband)**wfk_disk%formeig*mband) ) 3057 ABI_MALLOC(occ_disk, (mband)) 3058 3059 itimrev = kpts_timrev_from_kptopt(wfk_disk%hdr%kptopt) 3060 cryst = wfk_disk%hdr%get_crystal(itimrev + 1) 3061 3062 sppoldbl = 1 3063 ABI_MALLOC (rbz2disk, (sppoldbl*nkpt_in, 6)) 3064 3065 ABI_MALLOC (symrelT, (3,3,cryst%nsym)) 3066 ! TODO: from Matteo, this should be symrel straight, not transposed. Perhaps the logic in mapkptsets is transposed? 3067 do isym=1,cryst%nsym 3068 symrelT(:,:,isym) = transpose(cryst%symrel(:,:,isym)) 3069 end do 3070 3071 ask_accurate=1 3072 if (present(ask_accurate_)) ask_accurate=ask_accurate_ 3073 3074 ! Use listkk instead of rank-based routines since in DFPT we may receive a k+q mesh 3075 ! with q along a path --> max_linear_density in krank becomes large e.g. 1440 3076 ! and the computation of the rank overflows. 3077 ! Note also that ctgk_rotate assumes use_symrec=False and symrel in input. 3078 call listkk(dksqmax, cryst%gmet, rbz2disk, wfk_disk%hdr%kptns, kptns_in, wfk_disk%hdr%nkpt, nkpt_in, cryst%nsym, & 3079 sppoldbl, cryst%symafm, cryst%symrel, cryst%timrev-1, xmpi_comm_self, use_symrec=.False.) 3080 3081 !call xmpi_barrier(comm) 3082 if (ask_accurate == 1) then 3083 ABI_CHECK(dksqmax < tol8, sjoin("WFK file read but k-points too far from requested set, dksqmax:", ftoa(dksqmax))) 3084 end if 3085 3086 ! More efficienct algorithm based on random access IO: 3087 ! For each point in the irred disk set: 3088 ! - Read wavefunctions from wfk_disk 3089 ! - For each k-point in the star of kpt_disk: 3090 ! - Rotate wavefunctions in G-space to get the k-point in the requested BZ. 3091 ! - save kbz data. 3092 3093 ! Construct sorted mapping RBZ --> irred kdisk set, to speedup qbz search below. 3094 ABI_MALLOC(iperm, (nkpt_in)) 3095 ABI_MALLOC(rbz2disk_sort, (nkpt_in)) 3096 iperm = [(ii, ii=1,nkpt_in)] 3097 rbz2disk_sort = rbz2disk(:,1) 3098 call sort_int(nkpt_in, rbz2disk_sort, iperm) 3099 3100 ! prepare offsets for k-points, which could arrive in a random order from the irred k 3101 ! these are valid in the output arrays, not in the disk file 3102 !TODO: if nband_me is not constant over the k-points, this becomes a huge pain to predict... 3103 ABI_MALLOC(icg, (nkpt_in,nsppol)) 3104 ABI_MALLOC(ikg, (nkpt_in)) 3105 ABI_MALLOC(ibdeig, (nkpt_in,nsppol)) 3106 ABI_MALLOC(ibdocc, (nkpt_in,nsppol)) 3107 icg = 0 3108 ikg = 0 3109 ibdeig = 0 3110 ibdocc = 0 3111 ii = 0 3112 kk = 0 3113 ll = 0 3114 do spin=1,nsppol 3115 jj = 0 3116 do ikpt=1,nkpt_in 3117 ik_disk = rbz2disk(ikpt,1) 3118 3119 ! conversion of single spin AFM wfk file to full 2 component one in memory 3120 spin_sym=spin 3121 if (convnsppol1to2) spin_sym=1 3122 ! this allows for reading fewer bands from disk than the disk version of nband 3123 nband_k = min(wfk_disk%nband(ik_disk,spin_sym), mband_in) 3124 ibdeig(ikpt,spin) = kk 3125 ibdocc(ikpt,spin) = ll 3126 kk = kk+nband_k*(2*nband_k)**formeig 3127 ll = ll+nband_k 3128 3129 if (.not. any(distrb_flags(ikpt,:,spin))) cycle 3130 ! TODO: this does not take into account variable nband(ik) 3131 icg(ikpt,spin) = ii 3132 ikg(ikpt) = jj 3133 ! this allows for variable nband_k < mband_mem 3134 ii = ii+min(nband_k,mband_mem_in)*npwarr(ikpt)*nspinor_in 3135 jj = jj+npwarr(ikpt) 3136 end do 3137 end do 3138 3139 ! main loop reading in wfk and spinning them out to all kptns_in which need them 3140 3141 ! MG TODO: I believe this is not the most efficient way to implement the IO algorithm 3142 ! One might have only the master proc reading all the (ik_ibz, spin, mband) states 3143 ! perhaps blocking on the band dimension to reduce memory and then broadcast the block of bands. 3144 ! At this point, each proc rotates the wavefunctions and store it in memory if these states are needed. 3145 3146 do spin=1,nsppol 3147 ! for nsppol=1 input and nsppol=2 run, no need to continue the spin loop 3148 if (convnsppol1to2 .and. spin > 1) exit 3149 3150 do ik_disk=1,wfk_disk%hdr%nkpt 3151 k_disk = wfk_disk%hdr%kptns(:, ik_disk) 3152 3153 ! this allows for reading fewer bands from disk than the maximum 3154 nband_k = min(wfk_disk%nband(ik_disk,spin), mband_in) 3155 istwf_disk = wfk_disk%hdr%istwfk(ik_disk) 3156 npw_disk = wfk_disk%hdr%npwarr(ik_disk) 3157 3158 ! Find number of symmetric k-points associated to ik_disk 3159 nqst = 0 3160 needthisk=.false. 3161 iqst = 0 3162 ! scan to the first point which uses this kdisk 3163 do iqst = 1, nkpt_in 3164 if (rbz2disk_sort(iqst) == ik_disk) exit 3165 end do 3166 ! how many equivalent k? Could be 0, and we will not necessarily use them all if their bands are on other cpus 3167 do ii=iqst, nkpt_in 3168 if (rbz2disk_sort(ii) /= ik_disk) exit 3169 nqst = nqst + 1 3170 if (any(distrb_flags(iperm(ii),:,spin))) needthisk=.true. 3171 if (convnsppol1to2 .and. any(distrb_flags(iperm(ii),:,nsppol+1-spin))) needthisk=.true. 3172 end do ! loop over equivalent k 3173 3174 ! do we need the present kdisk, or one of its images? 3175 ! TODO: check if the eigenvalues are correct all the same 3176 if (.not. needthisk) cycle 3177 3178 ABI_CHECK(nqst > 0 .and. rbz2disk_sort(iqst) == ik_disk, "Wrong iqst") 3179 3180 ! loop over equivalent images found in rbz set for current k_disk point 3181 iband_saved = -1 3182 nband_me_saved = -1 3183 do jj=0,nqst-1 3184 ikf = iperm(iqst+jj) 3185 ABI_CHECK(ik_disk == rbz2disk(ikf,1), "ik_disk !/ ind qq(1)") 3186 3187 kf = kptns_in(:,ikf) 3188 istwf_kf = istwfk_in(ikf) 3189 npw_kf = npwarr(ikf) 3190 3191 do spin_sym = 1, nsppol 3192 if (.not. convnsppol1to2 .and. spin_sym /= spin) cycle 3193 3194 ! how many bands in memory for this cpu_ 3195 nband_me = count(distrb_flags(ikf,:,spin_sym)) 3196 ! no need to put wfk at this k for this processor into memory 3197 if (nband_me == 0) cycle 3198 3199 ! find starting band index 3200 do iband = 1, nband_k 3201 if (distrb_flags(ikf,iband,spin_sym)) exit 3202 end do 3203 ! check bands are contiguous in distrb_flags for this ikf and find first band needed, iband 3204 if (.not. distrb_flags(ikf,iband+nband_me-1,spin_sym)) then 3205 ABI_ERROR("wfk_read_my_kptbands: bands not contiguous in distrb_flags") 3206 end if 3207 3208 ! if nband_me goes beyond the end of the bands on disk, just read those we have 3209 nband_me_disk = min(nband_k,nband_me) 3210 3211 ! In parallel, iband+nband_me-1 could be larger than mband_disk 3212 ! we want to limit nband_me_disk in that case too, just for the last band procs 3213 if (iband+nband_me-1 > nband_k) then 3214 nband_me_disk = nband_k+1-iband 3215 end if 3216 3217 ! may need to re-read if for a different equivalent k if I need other bands 3218 if (iband /= iband_saved .or. nband_me_disk /= nband_me_saved .or. spin /= spin_saved) then 3219 if (formeig > 0) then 3220 call wfk_disk%read_band_block([iband,iband+nband_me_disk-1],ik_disk,spin,xmpio_single,& 3221 kg_k=kg_disk,cg_k=cg_disk,eig_k=eig_disk) 3222 else 3223 call wfk_disk%read_band_block([iband,iband+nband_me_disk-1],ik_disk,spin,xmpio_single,& 3224 kg_k=kg_disk,cg_k=cg_disk,eig_k=eig_disk,occ_k=occ_disk) 3225 end if 3226 ! in nsppol=1 nspden=2 case the occupations are doubled 3227 if (convnsppol1to2) then 3228 occ_disk = half * occ_disk 3229 end if 3230 iband_saved = iband 3231 nband_me_saved = nband_me_disk 3232 spin_saved = spin 3233 end if 3234 3235 ! reset isym for each spin_sym 3236 isym = rbz2disk(ikf,2) 3237 ! there is a first time reversal possible from the irred set found above to the kptns in input. 3238 ! a second possible time reversal if the irred k is not explicitly in the disk file, but only it's time reversed image 3239 itimrev = rbz2disk(ikf,6) 3240 g0 = rbz2disk(ikf,3:5) ! IS(k_disk) + g0 = k_bz 3241 3242 ! complete the spin down wfk with an AFM symop 3243 if (spin_sym /= spin) then 3244 ! try next symop to find afm operation to get the spin component we want 3245 do isym = 1, cryst%nsym 3246 if (cryst%symafm(isym) == 1) cycle 3247 ksym = matmul(symrelT(:,:,isym), k_disk) 3248 if (sum(abs(ksym-kf)) < tol8) exit 3249 end do 3250 ABI_CHECK(isym <= cryst%nsym, "did not find the AFM symop I need to get isppol=2 wave functions from disk") 3251 end if 3252 3253 isirred_kf = (isym == 1 .and. itimrev == 0 .and. all(g0 == 0) .and. cryst%symafm(isym) == 1) 3254 if (present(eigen)) then 3255 ibdoff = ibdeig(ikf,spin_sym)+(iband-1)*(2*nband_k)**formeig 3256 eigen(ibdoff+1:ibdoff+nband_me_disk*(2*nband_k)**formeig) = & 3257 eig_disk((iband-1)*(2*nband_k)**formeig+1:(iband-1+nband_me_disk)*(2*nband_k)**formeig) 3258 end if 3259 if (present(occ)) then 3260 ibdoff = ibdocc(ikf,spin_sym)+(iband-1) 3261 occ(ibdoff+1:ibdoff+nband_me_disk) = occ_disk(iband:iband-1+nband_me_disk) 3262 end if 3263 3264 ! The test on npwarr is needed because we may change istwfk e.g. gamma. 3265 if (isirred_kf .and. wfk_disk%hdr%npwarr(ik_disk) == npwarr(ikf)) then 3266 if (present(kg)) then 3267 kg(:,ikg(ikf)+1:ikg(ikf)+npw_kf) = kg_disk (:,1:npw_kf) 3268 end if 3269 cg(:,icg(ikf,spin_sym)+1:icg(ikf,spin_sym)+npw_kf*nband_me_disk*nspinor_in) = & 3270 cg_disk(:,1:npw_kf*nband_me_disk*nspinor_in) 3271 else 3272 ! Compute G-sphere centered on kf 3273 call get_kg(kf,istwf_kf,ecut_eff_in,cryst%gmet,npw_kf,kg_kf) 3274 ! npw found for the present sphere must be equal to size of array for output 3275 ABI_CHECK(npw_kf == npwarr(ikf), "Wrong npw_kf") 3276 3277 if (present(kg)) then 3278 kg(:,ikg(ikf)+1:ikg(ikf)+npw_kf) = kg_kf (:,1:npw_kf) 3279 end if 3280 3281 ! FFT box must enclose the two spheres centered on kdisk and kf 3282 gmax_disk = maxval(abs(kg_disk(:,1:npw_disk)), dim=2) 3283 gmax_kf = maxval(abs(kg_kf), dim=2) 3284 do ii=1,3 3285 gmax(ii) = max(gmax_disk(ii), gmax_kf(ii)) 3286 end do 3287 gmax = 2*gmax + 1 3288 call ngfft_seq(work_ngfft, gmax) 3289 ABI_CALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6))) 3290 3291 ! Rotate nband_k wavefunctions (output in cg) 3292 call cgtk_rotate(cryst,k_disk,isym,itimrev,g0,nspinor,nband_me_disk,& 3293 npw_disk,kg_disk,npw_kf,kg_kf,istwf_disk,istwf_kf,cg_disk,& 3294 cg(:,icg(ikf,spin_sym)+1:icg(ikf,spin_sym)+npw_kf*nband_me_disk*nspinor_in),& 3295 work_ngfft,work) 3296 3297 ABI_FREE(work) 3298 ABI_FREE(kg_kf) 3299 end if 3300 end do ! spin_sym 3301 end do ! equiv kpt jj 3302 end do ! kpt disk 3303 end do ! sppol 3304 3305 ! this sums over the whole kpt communicator, so also the band procs. 3306 ! need to 0 out bands which are not mine 3307 if(present(eigen)) call xmpi_sum(eigen,comm,mpierr) 3308 if(present(occ)) call xmpi_sum(occ,comm,mpierr) 3309 if(present(kg)) call xmpi_sum(kg,comm,mpierr) 3310 3311 if(present(pawrhoij) .and. usepaw_in==1) call pawrhoij_copy(wfk_disk%hdr%pawrhoij,pawrhoij) 3312 3313 ABI_FREE(icg) 3314 ABI_FREE(ikg) 3315 ABI_FREE(ibdeig) 3316 ABI_FREE(ibdocc) 3317 ABI_FREE(symrelT) 3318 ABI_FREE(iperm) 3319 ABI_FREE(rbz2disk_sort) 3320 ABI_FREE(rbz2disk) 3321 ABI_FREE(kg_disk) 3322 ABI_FREE(cg_disk) 3323 ABI_FREE(eig_disk) 3324 ABI_FREE(occ_disk) 3325 3326 call cryst%free() 3327 call wfk_disk%close() 3328 3329 call cwtime_report(" wfk_read_my_kptbands:", cpu, wall, gflops) 3330 3331 end subroutine wfk_read_my_kptbands
m_wfk/wfk_rewind [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_rewind
FUNCTION
Rewind the file, skip the header and modifies Wfk%f90_fptr $ Mainly used for debugging purposes when IO_MODE_FORTRAN is used.
SOURCE
3649 subroutine wfk_rewind(wfk) 3650 3651 !Arguments ------------------------------------ 3652 class(wfk_t),intent(inout) :: wfk 3653 3654 !Local variables------------------------------- 3655 integer :: ierr 3656 3657 ! ************************************************************************* 3658 3659 select case (wfk%iomode) 3660 case (IO_MODE_FORTRAN) 3661 rewind(wfk%fh) 3662 call hdr_skip(wfk%fh,ierr) 3663 ABI_CHECK(ierr==0, "hdr_skip returned ierr! /= 0") 3664 wfk%f90_fptr = [1,1,REC_NPW] 3665 3666 case default 3667 ABI_ERROR("should not be called when wfk%iomode /= IO_MODE_FORTRAN") 3668 end select 3669 3670 end subroutine wfk_rewind
m_wfk/wfk_seek [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_seek
FUNCTION
Move the internal file pointer so that it points to the block (ik_ibz, spin). Needed only if iomode==IO_MODE_FORTRAN
INPUTS
ik_ibz,spin = (k-point,spin) indices
SIDE EFFECTS
Wfk<class(wfk_t)> : modifies Wfk%f90_fptr and the internal F90 file pointer.
SOURCE
3691 subroutine wfk_seek(Wfk,ik_ibz,spin) 3692 3693 !Arguments ------------------------------------ 3694 integer,intent(in) :: ik_ibz,spin 3695 class(wfk_t),intent(inout) :: Wfk 3696 3697 !Local variables------------------------------- 3698 integer :: ierr,ik_fpt,spin_fpt,recn_wanted,recn_fpt,rec_type 3699 character(len=500) :: msg 3700 3701 ! ************************************************************************* 3702 3703 select case (Wfk%iomode) 3704 case (IO_MODE_FORTRAN) 3705 ! 3706 ! Find the position inside the file. 3707 if (ALL(Wfk%f90_fptr==FPTR_EOF)) then ! handle the EOF condition 3708 if (Wfk%debug) call wrtout(std_out,"EOF condition","PERS") 3709 recn_fpt = Wfk%recn_eof 3710 else 3711 ik_fpt = Wfk%f90_fptr(1) 3712 spin_fpt = Wfk%f90_fptr(2) 3713 rec_type = Wfk%f90_fptr(3) 3714 recn_fpt = Wfk%recn_ks(ik_fpt,spin_fpt, rec_type) 3715 end if 3716 recn_wanted = Wfk%recn_ks(ik_ibz,spin, REC_NPW) 3717 3718 if (Wfk%debug) then 3719 write(msg,'(a,3(i0,2x))')"seeking ik_ibz, spin, recn_wanted-recn_fpt: ",ik_ibz,spin,recn_wanted - recn_fpt 3720 call wrtout(std_out,msg,"PERS") 3721 end if 3722 3723 call mvrecord(Wfk%fh, (recn_wanted - recn_fpt) ,ierr) 3724 ABI_CHECK(ierr == 0, "error in mvrecord") 3725 3726 Wfk%f90_fptr = [ik_ibz, spin, REC_NPW] 3727 3728 case default 3729 ABI_ERROR("should not be called when Wfk%iomode /= IO_MODE_FORTRAN") 3730 end select 3731 3732 end subroutine wfk_seek
m_wfk/wfk_t [ Types ]
NAME
wfk_t
FUNCTION
File handler for the WFK file.
SOURCE
111 type,public :: wfk_t 112 113 integer :: fh 114 ! unit number if IO_MODE_FORTRAN 115 ! MPI file handler if IO_MODE_MPI 116 ! Netcdf file handler if IO_MODE_ETSF 117 118 integer :: iomode 119 ! Method used to access the WFK file: 120 ! IO_MODE_FORTRAN for usual Fortran IO routines 121 ! IO_MODE_MPI if MPI/IO routines. 122 ! IO_MODE_ETSF, NetCDF format read/written via etsf-io. 123 124 integer :: mband 125 ! Max number of bands stored on file (MAX(Hdr%nband)) 126 127 integer :: nkpt 128 ! Number of k-points. 129 130 integer :: nsppol 131 ! Number of spins 132 133 integer :: nspinor 134 ! Number of spinor components. 135 136 integer :: formeig 137 ! format of the eigenvalues 138 ! 0 => vector of eigenvalues (GS case) 139 ! 1 => Hermitian matrix of eigenvalues (DFPT case) 140 ! TODO: this should be reported somewhere in the WFK file, at present is passed to wfk_open 141 142 integer :: fform 143 ! File type format of the header 144 145 integer :: rw_mode = WFK_NOMODE 146 ! (Read|Write) mode 147 148 character(len=fnlen) :: fname = ABI_NOFILE 149 ! File name 150 151 integer :: master 152 ! master node of the IO procedure 153 154 integer :: my_rank 155 ! index of my processor in the MPI communicator comm 156 157 integer :: nproc 158 ! number of processors in comm 159 160 integer :: comm 161 ! MPI communicator 162 163 integer :: recn_eof 164 ! EOF record number (used for Fortran IO) 165 166 integer(XMPI_OFFSET_KIND) :: offset_eof 167 ! EOF offset (used for MPI-IO access) 168 169 logical :: debug = .FALSE. 170 !logical :: debug=.TRUE. 171 172 type(hdr_type) :: Hdr 173 ! Abinit header. 174 175 integer,allocatable :: nband(:,:) 176 ! nband(nkpt,nsppol) = Number of bands at each (k,s) 177 178 integer :: f90_fptr(3) = [0,0,0] 179 ! The position of the file pointer used for sequential access with Fortran-IO. 180 ! f90_fprt(1) = Index of the k-point associated to the block. 181 ! f90_fprt(2) = the spin associated to the block. 182 ! f90_fprt(3) = Record Type (see REC_* variables). 183 ! [0,0,0] corresponds to the beginning of the file. 184 ! FPTR_EOF signals the end of file 185 186 integer,allocatable :: recn_ks(:,:,:) 187 ! recn_ks(k,s,1) : record number of (npw, nspinor, nband_disk) 188 ! recn_ks(k,s,2) : record number of the (k+G) vectors. 189 ! recn_ks(k,s,3) : record number of the eigenvalues. 190 ! recn_ks(k,s,4) : record number of the first wavefunction in the wf coefficients block. 191 192 integer(XMPI_OFFSET_KIND),allocatable :: offset_ks(:,:,:) 193 ! offset_ks(k,s,1) : offset of the record: npw, nspinor, nband_disk. 194 ! offset_ks(k,s,2) : offset of the Second record: (k+G) vectors. 195 ! offset_ks(k,s,3) : offset of the third record eigenvalues. 196 ! offset_ks(k,s,4) : offset of the fourth record (wavefunction coefficients). 197 ! 198 ! ********************************************************************** 199 ! NB: The offset point to the Fortran record marker and not to the data 200 ! ********************************************************************** 201 202 integer(XMPI_OFFSET_KIND) :: hdr_offset 203 ! offset of the header 204 ! TODO this should be the output of a hdr method! 205 206 integer(XMPI_OFFSET_KIND) :: chunk_bsize 207 ! IO is performed in chunks of max size chunk_bsize [bytes] 208 209 contains 210 211 procedure :: open_write => wfk_open_write 212 ! Open the WFK file in write mode. 213 214 procedure :: close => wfk_close 215 ! Close the WFK file and release the memory allocated in wfk_t. 216 217 procedure :: print => wfk_print 218 ! Print info on the wfk_t object 219 220 procedure :: findk => wfk_findk 221 ! Returns the index of the k-point in the WFK file. 222 223 procedure :: compare => wfk_compare 224 ! Test two wfk_t objects for consistency. 225 226 procedure :: read_band_block => wfk_read_band_block 227 ! Read a contiguous block of bands for a given (kpoint, spin) 228 229 procedure :: read_bks => wfk_read_bks 230 ! Read the wavefunction and the eigenvalues for a given (band, k-point, spin) 231 232 procedure :: write_band_block => wfk_write_band_block 233 ! Write a contiguous block of bands for a given (kpoint, spin) 234 235 procedure :: read_bmask => wfk_read_bmask 236 ! Read a scattered set of bands for a given (kpoint, spin). 237 238 procedure :: read_eigk => wfk_read_eigk 239 ! Read the eigenvalues at a given (kpoint,spin). 240 241 procedure :: write_h1mat => wfk_write_h1mat 242 ! Write all the H1 matrix elements. 243 end type wfk_t 244 245 246 public :: wfk_open_read ! Open the WFK file in read mode. 247 public :: wfk_tofullbz ! Generate a new WFK file with wavefunctions in the full BZ and istwfk==1 248 ! Mainly used to interface ABINIT with other codes that 249 ! cannot handle symmetries e.g. lobster 250 public :: wfk_nc2fort ! Convert a netcdf WFK file to a Fortran WFK file. 251 public :: wfk_ncdef_dims_vars ! Define basic dimensions for netcdf file format. 252 public :: wfk_read_ebands ! Read the GS eigenvalues and return ebands_t object. 253 public :: wfk_read_eigenvalues ! Read all the GS eigenvalues stored in the WFK file. 254 public :: wfk_read_h1mat ! Read all the H1 matrix elements. 255 public :: wfk_read_my_kptbands ! Read in all of my bands and k, depending on a distribution flag array 256 public :: wfk_write_my_kptbands ! Write all of my bands and k to a file, depending on a distribution flag array 257 public :: wfk_klist2mesh ! Generate a full WFK file with k in the IBZ from a file with a subset of k-points 258 ! Mainly used in the transport part when the kerange trick is employed. 259 260 ! Profiling tools 261 public :: wfk_prof ! Profiling tool. 262 263 ! Unit tests 264 public :: wfk_diff ! Compare two WFK file for binary equality. 265 public :: wfk_create_wfkfile ! Create a FAKE WFK file. 266 public :: wfk_check_wfkfile ! Read a FAKE WFK file and perform basic tests. 267 public :: wfk_check_symtab
m_wfk/wfk_tofullbz [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_tofullbz
FUNCTION
Generate a new WFK file with wavefunctions in the full BZ and istwfk==1 Mainly used to interface ABINIT with other codes that cannot handle symmetries e.g. lobster
INPUTS
in_path = Input WFK file dtset <dataset_type>=all input variables for this dataset psps <pseudopotential_type>=all the information about psps pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data out_path = Output WFK file.
OUTPUT
Output is written to file out_path hdr_kfull: header of the WFK file.
NOTES
- This routine should be called by a single processor. - Only GS WFK files are supported (formeig==0)
SOURCE
4391 subroutine wfk_tofullbz(in_path, dtset, psps, pawtab, out_path, hdr_kfull) 4392 4393 !Arguments ------------------------------------ 4394 !scalars 4395 character(len=*),intent(in) :: in_path,out_path 4396 type(pseudopotential_type),intent(in) :: psps 4397 type(dataset_type),intent(in) :: dtset 4398 type(hdr_type),intent(out) :: hdr_kfull 4399 !arrays 4400 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 4401 4402 !Local variables------------------------------- 4403 !scalars 4404 integer,parameter :: formeig0=0,kptopt3=3 4405 integer :: spin,ikf,ik_ibz,nband_k,mpw_ki,mpw_kf,mband,nspinor,nkfull 4406 integer :: in_iomode,nsppol,nkibz,out_iomode,isym,itimrev 4407 integer :: npw_ki,npw_kf,istwf_ki,istwf_kf,ii,jj,iqst,nqst 4408 real(dp) :: ecut_eff,dksqmax,cpu,wall,gflops 4409 character(len=500) :: msg 4410 character(len=fnlen) :: my_inpath 4411 logical :: isirred_kf 4412 logical,parameter :: force_istwfk1=.True. 4413 type(wfk_t),target :: iwfk 4414 type(wfk_t) :: owfk 4415 type(crystal_t) :: cryst 4416 4417 type(hdr_type),pointer :: ihdr 4418 type(ebands_t) :: ebands_ibz 4419 type(ebands_t),target :: ebands_full 4420 type(wvl_internal_type) :: dummy_wvl 4421 !arrays 4422 integer :: g0(3),work_ngfft(18),gmax_ki(3),gmax_kf(3),gmax(3) 4423 integer,allocatable :: bz2ibz(:,:),kg_ki(:,:),kg_kf(:,:),iperm(:),bz2ibz_sort(:) 4424 real(dp) :: kf(3),kibz(3) 4425 real(dp),allocatable :: cg_ki(:,:),cg_kf(:,:),eig_ki(:),occ_ki(:),work(:,:,:,:) 4426 real(dp), ABI_CONTIGUOUS pointer :: kfull(:,:) 4427 4428 ! ************************************************************************* 4429 4430 if (all(dtset%kptrlatt == 0)) then 4431 write(msg,"(5a)")& 4432 "Cannot produce full WFK file because kptrlatt == 0",ch10,& 4433 "Please use nkgpt and shiftk to define a homogeneous k-mesh.",ch10,& 4434 "Returning to caller" 4435 ABI_WARNING(msg) 4436 return 4437 end if 4438 4439 call cwtime(cpu, wall, gflops, "start") 4440 my_inpath = in_path 4441 4442 if (nctk_try_fort_or_ncfile(my_inpath, msg) /= 0) then 4443 ABI_ERROR(msg) 4444 end if 4445 call wrtout(std_out, sjoin(" Converting:", my_inpath, "to", out_path)) 4446 4447 in_iomode = iomode_from_fname(my_inpath) 4448 4449 ebands_ibz = wfk_read_ebands(my_inpath, xmpi_comm_self) 4450 4451 ! Open input file, extract dimensions and allocate workspace arrays. 4452 call wfk_open_read(iwfk,my_inpath,formeig0,in_iomode,get_unit(),xmpi_comm_self) 4453 ihdr => iwfk%hdr 4454 4455 mband = iwfk%mband; mpw_ki = maxval(iwfk%Hdr%npwarr); nkibz = iwfk%nkpt 4456 nsppol = iwfk%nsppol; nspinor = iwfk%nspinor 4457 ecut_eff = iwfk%hdr%ecut_eff ! ecut * dilatmx**2 4458 4459 ABI_MALLOC(kg_ki, (3, mpw_ki)) 4460 ABI_MALLOC(cg_ki, (2, mpw_ki*nspinor*mband)) 4461 ABI_MALLOC(eig_ki, ((2*mband)**iwfk%formeig*mband) ) 4462 ABI_MALLOC(occ_ki, (mband)) 4463 4464 cryst = iwfk%hdr%get_crystal() 4465 4466 ! Build new header for owfk. This is the most delicate part since all the arrays in hdr_full 4467 ! that depend on k-points must be consistent with kfull and nkfull. 4468 call ebands_expandk(ebands_ibz, cryst, ecut_eff, force_istwfk1, dksqmax, bz2ibz, ebands_full) 4469 4470 if (dksqmax > tol12) then 4471 write(msg, '(3a,es16.6,4a)' )& 4472 'At least one of the k points could not be generated from a symmetrical one.',ch10,& 4473 'dksqmax=',dksqmax,ch10,& 4474 'Action: check your WFK file and k-point input variables',ch10,& 4475 ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.' 4476 ABI_ERROR(msg) 4477 end if 4478 4479 nkfull = ebands_full%nkpt 4480 kfull => ebands_full%kptns 4481 4482 ! Build new header and update pawrhoij. 4483 call hdr_init_lowlvl(hdr_kfull,ebands_full,psps,pawtab,dummy_wvl,abinit_version,& 4484 ihdr%pertcase,ihdr%natom,ihdr%nsym,ihdr%nspden,ihdr%ecut,dtset%pawecutdg,ihdr%ecutsm,dtset%dilatmx,& 4485 ihdr%intxc,ihdr%ixc,ihdr%stmbias,ihdr%usewvl,dtset%pawcpxocc,dtset%pawspnorb,dtset%ngfft,dtset%ngfftdg,ihdr%so_psp,& 4486 ihdr%qptn,cryst%rprimd,cryst%xred,ihdr%symrel,ihdr%tnons,ihdr%symafm,ihdr%typat,ihdr%amu,ihdr%icoulomb,& 4487 kptopt3,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,dtset%cellcharge(1),& 4488 dtset%kptrlatt_orig,dtset%kptrlatt,& 4489 dtset%nshiftk_orig,dtset%nshiftk,dtset%shiftk_orig,dtset%shiftk) 4490 4491 if (psps%usepaw == 1) call pawrhoij_copy(iwfk%hdr%pawrhoij, hdr_kfull%pawrhoij) 4492 4493 out_iomode = iomode_from_fname(out_path) 4494 call owfk%open_write(hdr_kfull, out_path, iwfk%formeig, out_iomode, get_unit(), xmpi_comm_self) 4495 4496 ! workspace array for BZ wavefunction block. 4497 mpw_kf = maxval(ebands_full%npwarr) 4498 ABI_MALLOC(cg_kf, (2,mpw_kf*nspinor*mband)) 4499 4500 if (out_iomode == IO_MODE_FORTRAN) then 4501 call wrtout(std_out,"Using (slow) Fortran IO version to generate full WFK file", do_flush=.True.) 4502 4503 ! Fortran IO does not support random access hence the external loop is on the k-points in the full BZ. 4504 ! 4505 ! For each point in the BZ: 4506 ! - Find symmetric k-point in the IBZ and read IBZ wavefunctions from iwfk 4507 ! - Rotate wavefunctions in G-space and write kbz data 4508 ! 4509 ! Inefficient since we are reading the same IBZ block several times. 4510 do spin=1,nsppol 4511 do ikf=1,nkfull 4512 ik_ibz = bz2ibz(ikf,1); isym = bz2ibz(ikf,2); itimrev = bz2ibz(ikf,6); g0 = bz2ibz(ikf,3:5) ! IS(k_ibz) + g0 = k_bz 4513 isirred_kf = (isym == 1 .and. itimrev == 0 .and. all(g0 == 0)) 4514 4515 nband_k = iwfk%nband(ik_ibz,spin) 4516 kf = kfull(:,ikf) 4517 kibz = ebands_ibz%kptns(:,ik_ibz) 4518 4519 istwf_ki = iwfk%hdr%istwfk(ik_ibz) 4520 istwf_kf = owfk%hdr%istwfk(ikf) 4521 npw_ki = iwfk%hdr%npwarr(ik_ibz) 4522 4523 ! Read IBZ data. 4524 call iwfk%read_band_block([1,nband_k],ik_ibz,spin,xmpio_single,& 4525 kg_k=kg_ki,cg_k=cg_ki,eig_k=eig_ki,occ_k=occ_ki) 4526 4527 ! The test on npwarr is needed because we may change istwfk e.g. gamma. 4528 if (isirred_kf .and. iwfk%hdr%npwarr(ik_ibz) == owfk%hdr%npwarr(ikf)) then 4529 4530 call owfk%write_band_block([1,nband_k],ikf,spin,xmpio_single,& 4531 kg_k=kg_ki,cg_k=cg_ki,eig_k=eig_ki,occ_k=occ_ki) 4532 4533 else 4534 ! Compute G-sphere centered on kf 4535 call get_kg(kf,istwf_kf,ecut_eff,cryst%gmet,npw_kf,kg_kf) 4536 ABI_CHECK(npw_kf == owfk%hdr%npwarr(ikf), "Wrong npw_kf") 4537 4538 ! FFT box must enclose the two spheres centered on ki and kf 4539 gmax_ki = maxval(abs(kg_ki(:,1:npw_ki)), dim=2) 4540 gmax_kf = maxval(abs(kg_kf), dim=2) 4541 do ii=1,3 4542 gmax(ii) = max(gmax_ki(ii), gmax_kf(ii)) 4543 end do 4544 gmax = 2*gmax + 1 4545 call ngfft_seq(work_ngfft, gmax) 4546 ABI_CALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6))) 4547 4548 ! Rotate nband_k wavefunctions (output in cg_kf) 4549 call cgtk_rotate(cryst,kibz,isym,itimrev,g0,nspinor,nband_k,& 4550 npw_ki,kg_ki,npw_kf,kg_kf,istwf_ki,istwf_kf,cg_ki,cg_kf,work_ngfft,work) 4551 4552 ABI_FREE(work) 4553 4554 ! Write data 4555 call owfk%write_band_block([1,nband_k],ikf,spin,xmpio_single,& 4556 kg_k=kg_kf,cg_k=cg_kf,eig_k=eig_ki,occ_k=occ_ki) 4557 4558 ABI_FREE(kg_kf) 4559 end if 4560 end do 4561 end do 4562 4563 else 4564 ! 4565 ! More efficienct algorithm based on random access IO: 4566 ! For each point in the IBZ: 4567 ! - Read wavefunctions from iwfk 4568 ! - For each k-point in the star of kpt_ibz: 4569 ! - Rotate wavefunctions in G-space to get the k-point in the full BZ. 4570 ! - Write kbz data to file. 4571 if (out_iomode == IO_MODE_MPI) call wrtout(std_out," Using MPI-IO to generate full WFK file", do_flush=.True.) 4572 if (out_iomode == IO_MODE_ETSF) call wrtout(std_out, "Using Netcdf-IO to generate full WFK file", do_flush=.True.) 4573 4574 ! Construct sorted mapping BZ --> IBZ to speedup qbz search below. 4575 ABI_MALLOC(iperm, (nkfull)) 4576 ABI_MALLOC(bz2ibz_sort, (nkfull)) 4577 iperm = [(ii, ii=1,nkfull)] 4578 bz2ibz_sort = bz2ibz(:,1) 4579 call sort_int(nkfull, bz2ibz_sort, iperm) 4580 4581 do spin=1,nsppol 4582 iqst = 0 4583 do ik_ibz=1,iwfk%nkpt 4584 nband_k = iwfk%nband(ik_ibz,spin) 4585 kibz = ebands_ibz%kptns(:,ik_ibz) 4586 istwf_ki = iwfk%hdr%istwfk(ik_ibz) 4587 npw_ki = iwfk%hdr%npwarr(ik_ibz) 4588 4589 call iwfk%read_band_block([1,nband_k],ik_ibz,spin,xmpio_single,& 4590 kg_k=kg_ki,cg_k=cg_ki,eig_k=eig_ki,occ_k=occ_ki) 4591 4592 ! Find number of symmetric q-points associated to ik_ibz 4593 nqst = 0 4594 do ii=iqst+1,nkfull 4595 if (bz2ibz_sort(ii) /= ik_ibz) exit 4596 nqst = nqst + 1 4597 end do 4598 ABI_CHECK(nqst > 0 .and. bz2ibz_sort(iqst+1) == ik_ibz, "Wrong iqst") 4599 4600 do jj=1,nqst 4601 iqst = iqst + 1 4602 ikf = iperm(iqst) 4603 ABI_CHECK(ik_ibz == bz2ibz(ikf,1), "ik_ibz !/ ind qq(1)") 4604 4605 isym = bz2ibz(ikf,2); itimrev = bz2ibz(ikf,6); g0 = bz2ibz(ikf,3:5) ! IS(k_ibz) + g0 = k_bz 4606 isirred_kf = (isym == 1 .and. itimrev == 0 .and. all(g0 == 0)) 4607 4608 kf = kfull(:,ikf) 4609 istwf_kf = owfk%hdr%istwfk(ikf) 4610 4611 ! The test on npwarr is needed because we may change istwfk e.g. gamma. 4612 if (isirred_kf .and. iwfk%hdr%npwarr(ik_ibz) == owfk%hdr%npwarr(ikf)) then 4613 4614 call owfk%write_band_block([1,nband_k],ikf,spin,xmpio_single,& 4615 kg_k=kg_ki,cg_k=cg_ki,eig_k=eig_ki,occ_k=occ_ki) 4616 4617 else 4618 ! Compute G-sphere centered on kf 4619 call get_kg(kf,istwf_kf,ecut_eff,cryst%gmet,npw_kf,kg_kf) 4620 ABI_CHECK(npw_kf == owfk%hdr%npwarr(ikf), "Wrong npw_kf") 4621 4622 ! FFT box must enclose the two spheres centered on ki and kf 4623 gmax_ki = maxval(abs(kg_ki(:,1:npw_ki)), dim=2) 4624 gmax_kf = maxval(abs(kg_kf), dim=2) 4625 do ii=1,3 4626 gmax(ii) = max(gmax_ki(ii), gmax_kf(ii)) 4627 end do 4628 gmax = 2*gmax + 1 4629 call ngfft_seq(work_ngfft, gmax) 4630 ABI_CALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6))) 4631 4632 ! Rotate nband_k wavefunctions (output in cg_kf) 4633 call cgtk_rotate(cryst,kibz,isym,itimrev,g0,nspinor,nband_k,& 4634 npw_ki,kg_ki,npw_kf,kg_kf,istwf_ki,istwf_kf,cg_ki,cg_kf,work_ngfft,work) 4635 4636 ABI_FREE(work) 4637 4638 ! Write data 4639 call owfk%write_band_block([1,nband_k],ikf,spin,xmpio_single,& 4640 kg_k=kg_kf,cg_k=cg_kf,eig_k=eig_ki,occ_k=occ_ki) 4641 4642 ABI_FREE(kg_kf) 4643 end if 4644 end do 4645 end do 4646 end do 4647 4648 ABI_FREE(iperm) 4649 ABI_FREE(bz2ibz_sort) 4650 end if 4651 4652 call cwtime_report(sjoin(" FULL_WFK written to: ", out_path), cpu, wall, gflops) 4653 4654 ABI_FREE(kg_ki) 4655 ABI_FREE(cg_ki) 4656 ABI_FREE(eig_ki) 4657 ABI_FREE(occ_ki) 4658 ABI_FREE(bz2ibz) 4659 ABI_FREE(cg_kf) 4660 4661 call cryst%free() 4662 call ebands_free(ebands_ibz) 4663 call ebands_free(ebands_full) 4664 call iwfk%close() 4665 call owfk%close() 4666 4667 end subroutine wfk_tofullbz
m_wfk/wfk_update_f90ptr [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_update_f90ptr
FUNCTION
Update wfk%f90_ptr. Used if wfk%iomode == IO_MODE_FORTRAN.
INPUTS
ik_ibz=K-point index, spin=Spin index.
SOURCE
3750 subroutine wfk_update_f90ptr(wfk, ik_ibz, spin) 3751 3752 !Arguments ------------------------------------ 3753 class(wfk_t),intent(inout) :: wfk 3754 integer,intent(in) :: ik_ibz,spin 3755 3756 ! ************************************************************************* 3757 3758 if (ik_ibz < wfk%nkpt) then 3759 wfk%f90_fptr = [ik_ibz+1,spin,REC_NPW] 3760 else 3761 ABI_CHECK(ik_ibz == wfk%nkpt, "ik_ibz != nkpt") 3762 if (spin == wfk%nsppol) then 3763 wfk%f90_fptr = FPTR_EOF ! EOF condition 3764 else 3765 wfk%f90_fptr = [1,spin+1,REC_NPW] 3766 end if 3767 end if 3768 3769 end subroutine wfk_update_f90ptr
m_wfk/wfk_validate_ks [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_validate_ks
FUNCTION
Validate the k-point, the spin index and, optionally, the band index. Return non-zero value if error.
INPUTS
wfk<class(wfk_t)> = WFK handler ik_ibz=k-point index. spin=Spin index. [band]=Band index.
SOURCE
804 integer function wfk_validate_ks(wfk, ik_ibz, spin, band) result(ierr) 805 806 !Arguments ------------------------------------ 807 !scalars 808 integer,intent(in) :: ik_ibz, spin 809 integer,optional,intent(in) :: band 810 class(wfk_t),intent(in) :: wfk 811 812 !Local variables------------------------------- 813 !scalars 814 character(len=500) :: msg 815 816 ! ************************************************************************* 817 ierr = 0 818 819 if (ik_ibz <= 0 .or. ik_ibz > wfk%nkpt) then 820 ierr = ierr + 1 821 write(msg, '(2(a,i0))')'ik_ibz = ',ik_ibz,' whereas it should be between 1 and ',wfk%nkpt 822 ABI_WARNING(msg) 823 end if 824 825 if (spin <= 0 .or. spin > wfk%nsppol) then 826 ierr = ierr + 1 827 write(msg, '(2(a,i0))')'spin = ',spin,' whereas it should be between 1 and ',wfk%nsppol 828 ABI_WARNING(msg) 829 end if 830 831 if (present(band)) then 832 if (band <=0) then 833 ierr = ierr + 1 834 ABI_WARNING(sjoin('Negative band index: band = ',itoa(band))) 835 end if 836 837 ! Don't touch nband array if wrong indices. 838 if (spin > 0 .and. spin <= wfk%nsppol .and. ik_ibz > 0 .and. ik_ibz <= wfk%nkpt) then 839 if (band > wfk%nband(ik_ibz, spin)) then 840 ierr = ierr + 1 841 write(msg, '(2(a,i0))')'band = ',band,' whereas it should be between 1 and ',wfk%nband(ik_ibz,spin) 842 ABI_WARNING(msg) 843 end if 844 end if 845 end if 846 847 !if (ierr /= 0) then 848 ! ABI_ERROR("Wrong (ik_ibz, spin) args, Aborting now") 849 !end if 850 851 end function wfk_validate_ks
m_wfk/wfk_write_band_block [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_write_band_block
FUNCTION
Write a block of contigous bands.
INPUTS
Wfk<class(wfk_t)>= band_block(2)=Initial and final band index. ik_ibz=Index of the k-point in the IBZ. spin=Spin index sc_mode= MPI-IO option xmpio_single ==> for writing by current proc. xmpio_collective ==> for collective writing. [kg_k=(:,:)] = G-vectors [cg_k(:,:)] = Fourier coefficients [eig_k(:)] = Eigenvalues (dimensioned with wfk%mband, see below) [occ_k(:)] = Occupancies ((dimensioned with wfk%mband, see below)
SOURCE
1753 subroutine wfk_write_band_block(Wfk, band_block, ik_ibz, spin, sc_mode, kg_k, cg_k, eig_k, occ_k) 1754 1755 !Arguments ------------------------------------ 1756 !scalars 1757 class(wfk_t),intent(inout) :: Wfk 1758 integer,intent(in) :: ik_ibz,spin,sc_mode 1759 !arrays 1760 integer,intent(in) :: band_block(2) 1761 integer,intent(in),optional :: kg_k(:,:) ! (3, npw_k) 1762 real(dp),intent(in),optional :: cg_k(:,:) ! (2, npw_k*nspinor*nband) 1763 real(dp),intent(in),optional :: eig_k((2*Wfk%mband)**Wfk%formeig*Wfk%mband) 1764 real(dp),intent(in),optional :: occ_k(Wfk%mband) 1765 1766 !Local variables------------------------------- 1767 !scalars 1768 integer :: npw_disk,nspinor_disk,nband_disk,band 1769 integer :: ipw,my_bcount,npwso,npw_tot,nb_block,base 1770 character(len=500) :: errmsg !msg, 1771 real(dp) :: cpu, wall, gflops 1772 !arrays 1773 real(dp),ABI_CONTIGUOUS pointer :: tmp_eigk(:) 1774 !real(dp), allocatable :: eig_buffer(:), cg_buffer(:,:) 1775 #ifdef HAVE_MPI_IO 1776 integer :: mpierr,bufsz,recnpw_type,gkk_type,cgblock_type 1777 integer(XMPI_OFFSET_KIND) :: my_offset,my_offpad 1778 integer :: sizes(2),subsizes(2),starts(2),dims(3),types(2) 1779 !integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecords(:) 1780 #endif 1781 integer :: kg_varid,eig_varid,occ_varid,cg_varid,ncerr,h1_varid 1782 1783 !************************************************************************ 1784 1785 DBG_ENTER("COLL") 1786 ABI_CHECK_IEQ(Wfk%rw_mode, WFK_WRITEMODE, "Wfk must be in WRITEMODE") 1787 1788 call cwtime(cpu, wall, gflops, "start") 1789 1790 ! Look before you leap. 1791 npw_disk = Wfk%Hdr%npwarr(ik_ibz) 1792 nspinor_disk = Wfk%nspinor 1793 nband_disk = Wfk%nband(ik_ibz,spin) 1794 nb_block = (band_block(2) - band_block(1) + 1) 1795 npw_tot = npw_disk * nspinor_disk * nb_block 1796 1797 ! MG: We don't need to allocate memory if we want to skip records. 1798 ! Plain read without variable will do the job. 1799 !ABI_MALLOC (eig_buffer, (2*nband_disk)) 1800 !ABI_MALLOC (cg_buffer, (2,npw_disk*nspinor_disk)) 1801 1802 if (PRESENT(kg_k)) then 1803 ABI_CHECK_IGEQ(SIZE(kg_k, DIM=2), npw_disk, "kg_k too small") 1804 end if 1805 1806 if (PRESENT(cg_k)) then 1807 ABI_CHECK_IGEQ(SIZE(cg_k, DIM=2), npw_tot, "cg_k too small") 1808 end if 1809 1810 if (PRESENT(eig_k)) then 1811 if (Wfk%formeig == 0) then 1812 ABI_CHECK_IGEQ(SIZE(eig_k), nband_disk, "GS eig_k too small") 1813 ABI_CHECK(PRESENT(occ_k), "both eig_k and occ_k must be present") 1814 else if (Wfk%formeig == 1) then 1815 ABI_CHECK_IGEQ(SIZE(eig_k), 2*nband_disk**2, "DFPT eig_k too small") 1816 else 1817 ABI_ERROR("formeig != [0,1]") 1818 end if 1819 end if 1820 1821 if (PRESENT(occ_k)) then 1822 ABI_CHECK_IGEQ(SIZE(occ_k), nband_disk, "GS eig_k too small") 1823 ABI_CHECK(PRESENT(eig_k), "both eig_k and occ_k must be present") 1824 ABI_CHECK(Wfk%formeig == 0, "formeig /=0 with occ_k in input!") 1825 end if 1826 1827 select case (Wfk%iomode) 1828 case (IO_MODE_FORTRAN) 1829 1830 ! Rewind the file to have the correct (k,s) block (if needed) 1831 call wfk_seek(Wfk,ik_ibz,spin) 1832 1833 ! Write the first record: npw, nspinor, nband_disk 1834 write(Wfk%fh, err=10, iomsg=errmsg) npw_disk, nspinor_disk, nband_disk 1835 1836 ! The second record: (k+G) vectors 1837 if (PRESENT(kg_k)) then 1838 write(Wfk%fh, err=10, iomsg=errmsg) kg_k(1:3,1:npw_disk) 1839 else 1840 read(Wfk%fh, err=10, iomsg=errmsg) ! kg_k(1:3,1:npw_disk) 1841 end if 1842 1843 ! The third record: eigenvalues occupation factors and wavefunctions. 1844 select case (Wfk%formeig) 1845 case (0) 1846 !write(unitwf) (eigen(iband),iband=1,nband_disk),(occ(iband),iband=1,nband_disk) 1847 1848 if (present(eig_k) .and. present(occ_k)) then 1849 write(Wfk%fh, err=10, iomsg=errmsg) eig_k(1:nband_disk), occ_k(1:nband_disk) 1850 else 1851 ABI_ERROR("Not coded") 1852 write(Wfk%fh, err=10, iomsg=errmsg) ! eig_k(1:nband_disk), occ_k(1:nband_disk) 1853 end if 1854 1855 ! The wave-functions. 1856 if (present(cg_k)) then 1857 npwso = npw_disk*nspinor_disk 1858 ! fast forward the bands which are not mine 1859 ! could do in a single read, but need to check if cg_k is big enough as a buffer 1860 ! e.g. for band_block(1)=100 and band_block(2)=105 1861 do band=1,band_block(1)-1 1862 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_buffer(1:2,1:npwso) 1863 end do 1864 1865 ! MJV 2021/02: I think my coding is correct - the previous one would only accept band_block(:) = 1,nband_disk 1866 my_bcount = 0 1867 do band=band_block(1), band_block(2) 1868 ipw = my_bcount * npwso 1869 my_bcount = my_bcount + 1 1870 write(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 1871 end do 1872 1873 else 1874 ABI_ERROR("Not coded") 1875 do band=1,nband_disk 1876 write(Wfk%fh, err=10, iomsg=errmsg) ! cg_k(1:2,ipw+1:ipw+npwso) 1877 end do 1878 end if 1879 1880 case (1) 1881 ! Write column of matrix of total size (2*nband_k**2) 1882 ! And the wave-functions. 1883 npwso = npw_disk*nspinor_disk 1884 1885 ! fast forward the bands which are not mine 1886 ! could do in a single read, but need to check if cg_k is big enough as a buffer 1887 ! e.g. for band_block(1)=100 and band_block(2)=105 1888 do band=1, band_block(1)-1 1889 read(Wfk%fh, err=10, iomsg=errmsg) ! eig_buffer(1:2*nband_disk) 1890 read(Wfk%fh, err=10, iomsg=errmsg) ! cg_buffer(1:2,1:npwso) 1891 end do 1892 1893 my_bcount = 0 1894 do band=band_block(1),band_block(2) 1895 base = 2*(band-1)*nband_disk 1896 !NB: interleaves the arrays eig_k and cg_k in the RF case with formeig 1 1897 write(Wfk%fh, err=10, iomsg=errmsg) eig_k(base+1:base+2*nband_disk) 1898 ! MJV 2021/02: I think my coding is correct - the previous one would only accept band_block(:) = 1,nband_disk 1899 ipw = my_bcount * npwso 1900 my_bcount = my_bcount + 1 1901 write(Wfk%fh, err=10, iomsg=errmsg) cg_k(1:2,ipw+1:ipw+npwso) 1902 end do 1903 1904 case default 1905 ABI_ERROR("formeig != [0,1]") 1906 end select 1907 1908 ! Reached the end of the (k,s) block. Update f90_fptr 1909 call wfk_update_f90ptr(wfk, ik_ibz, spin) 1910 1911 #ifdef HAVE_MPI_IO 1912 case (IO_MODE_MPI) 1913 ! record 1 npw, nspinor, nband of length 3 1914 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_NPW) 1915 1916 ! bsize_rec(1) = 3 * xmpi_bsize_int 1917 ! call xmpio_write_frmarkers(Wfk%fh,my_offset,sc_mode,1,bsize_rec,mpierr) 1918 ! ABI_CHECK(mpierr==0,"mpierr!=0") 1919 1920 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_NPW) + xmpio_bsize_frm 1921 1922 call MPI_TYPE_CONTIGUOUS(3, MPI_INTEGER, recnpw_type, mpierr) 1923 ABI_CHECK_MPI(mpierr, "writing REC_NPW") 1924 1925 call MPI_TYPE_COMMIT(recnpw_type, mpierr) 1926 ABI_CHECK_MPI(mpierr, "writing REC_NPW") 1927 1928 ! NB: This is a collection operation so all proch in wfk%comm must call the routine. 1929 call MPI_FILE_SET_VIEW(Wfk%fh, my_offset, MPI_BYTE, recnpw_type, 'native', xmpio_info, mpierr) 1930 ABI_CHECK_MPI(mpierr, "writing REC_NPW") 1931 1932 call MPI_TYPE_FREE(recnpw_type, mpierr) 1933 ABI_CHECK_MPI(mpierr, "writing REC_NPW") 1934 1935 dims = [npw_disk, nspinor_disk, nband_disk] 1936 1937 if (sc_mode == xmpio_collective) then 1938 call MPI_FILE_WRITE_ALL(Wfk%fh, dims, SIZE(dims), MPI_INTEGER, MPI_STATUS_IGNORE, mpierr) 1939 else if (sc_mode == xmpio_single) then 1940 call MPI_FILE_WRITE(Wfk%fh, dims, SIZE(dims), MPI_INTEGER, MPI_STATUS_IGNORE, mpierr) 1941 else 1942 ABI_ERROR("Wrong sc_mode") 1943 end if 1944 ABI_CHECK_MPI(mpierr, "writing REC_NPW") 1945 1946 !---------------------------------------------------------------------------- 1947 ! record 2 kg 1948 if (present(kg_k)) then 1949 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_KG) 1950 1951 ! bsize_rec(1) = 3 * npw_disk * xmpi_bsize_int 1952 ! call xmpio_write_frmarkers(Wfk%fh,my_offset,sc_mode,1,bsize_rec,mpierr) 1953 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_KG) + xmpio_bsize_frm 1954 1955 call mpio_write_kg_k(Wfk%fh, my_offset, npw_disk, sc_mode, kg_k, mpierr) 1956 ABI_CHECK_MPI(mpierr, "mpio_write_kg_k") 1957 end if 1958 1959 if (Wfk%formeig==0) then 1960 !---------------------------------------------------------------------------- 1961 ! record 3 eigk occk 1962 if (present(eig_k) .and. present(occ_k)) then 1963 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) 1964 1965 ! bsize_rec(1) = 2 * nband_disk * xmpi_bsize_dp 1966 ! call xmpio_write_frmarkers(Wfk%fh,my_offset,sc_mode,1,bsize_rec,mpierr) 1967 1968 !TODO: check if we need 2*bsize_frm here 1969 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) + xmpio_bsize_frm 1970 ! 1971 ! Write both eig and occ in tmp_eigk 1972 bufsz = 2*nband_disk 1973 ABI_MALLOC(tmp_eigk, (bufsz)) 1974 1975 tmp_eigk(1:nband_disk) = eig_k(1:nband_disk) 1976 tmp_eigk(nband_disk+1:) = occ_k(1:nband_disk) 1977 1978 call mpio_write_eigocc_k(Wfk%fh, my_offset, nband_disk, Wfk%formeig, sc_mode, tmp_eigk, mpierr) 1979 ABI_CHECK_MPI(mpierr, "mpio_write_eigocc_k") 1980 1981 ABI_FREE(tmp_eigk) 1982 end if 1983 1984 !---------------------------------------------------------------------------- 1985 ! record 4 cg 1986 if (present(cg_k)) then 1987 !TODO: in principle these markers are written when the file is opened, no need here. 1988 !ABI_MALLOC(bsize_frecords, (nb_block)) 1989 !bsize_frecords = 2 * npw_disk * nspinor_disk * xmpi_bsize_dp 1990 !! TODO: why 2*frm size here? Each band cg is a single record! 1991 !my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) + (band_block(1)-1) * (bsize_frecords(1) + 2*xmpio_bsize_frm) 1992 !call xmpio_write_frmarkers(Wfk%fh,my_offset,sc_mode,nb_block,bsize_frecords,mpierr) 1993 !ABI_CHECK(mpierr==0,"mpierr!=0") 1994 !ABI_FREE(bsize_frecords) 1995 !print *, "Writing cg" 1996 1997 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) 1998 sizes = [npw_disk*nspinor_disk, nband_disk] 1999 subsizes = [npw_disk*nspinor_disk, band_block(2)-band_block(1)+1] 2000 bufsz = 2 * npw_disk * nspinor_disk * nb_block 2001 starts = [1, band_block(1)] 2002 2003 call mpiotk_write_fsuba_dp2D(Wfk%fh,my_offset,sizes,subsizes,starts,bufsz,cg_k,Wfk%chunk_bsize,sc_mode,Wfk%comm,mpierr) 2004 ABI_CHECK(mpierr == 0, "mpierr != 0") 2005 end if 2006 2007 else if (Wfk%formeig == 1) then 2008 2009 !---------------------------------------------------------------------------- 2010 ! record 3 eigk occk 2011 if (present(eig_k)) then 2012 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 2013 sizes = [nband_disk, npw_disk*nspinor_disk] 2014 2015 call xmpio_create_fstripes(nband_disk,sizes,types,gkk_type,my_offpad,mpierr) 2016 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 2017 2018 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_EIG) + my_offpad 2019 2020 call MPI_FILE_SET_VIEW(Wfk%fh,my_offset,MPI_BYTE,gkk_type,'native',xmpio_info,mpierr) 2021 ABI_CHECK_MPI(mpierr, "SET_VIEW") 2022 2023 call MPI_TYPE_FREE(gkk_type,mpierr) 2024 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 2025 2026 ! NB: bufsz is not 2*nband**2 because we use COMPLEX below 2027 bufsz = nband_disk**2 2028 2029 if (sc_mode == xmpio_collective) then 2030 call MPI_FILE_WRITE_ALL(Wfk%fh,eig_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2031 else if (sc_mode == xmpio_single) then 2032 call MPI_FILE_WRITE(Wfk%fh,eig_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2033 else 2034 ABI_ERROR("Wrong sc_mode") 2035 end if 2036 2037 ABI_CHECK_MPI(mpierr, "FILE_WRITE") 2038 end if 2039 2040 !---------------------------------------------------------------------------- 2041 ! record 4 cg 2042 if (present(cg_k)) then 2043 !ABI_CHECK(band_block(1)==1,"band_block(1) !=1 not coded") 2044 2045 types = [MPI_DOUBLE_COMPLEX, MPI_DOUBLE_COMPLEX] 2046 !sizes = [npw_disk*nspinor_disk, band_block(2)-band_block(1)] 2047 sizes = [npw_disk*nspinor_disk, nband_disk] 2048 2049 call xmpio_create_fstripes(nb_block,sizes,types,cgblock_type,my_offpad,mpierr) 2050 ABI_CHECK_MPI(mpierr, "xmpio_create_fstripes") 2051 2052 ! TODO: check that the following offset is correct 2053 ! check that the 4 * xmpio_bsize_frm is correct: 1 record marker for eigen and 1 for cg in principle! 2054 ! even if the cg is followed by 2 frm, and eig 1, then it should be 3, not 4 2055 my_offset = Wfk%offset_ks(ik_ibz,spin,REC_CG) + my_offpad & 2056 + (band_block(1)-1) * (2 * nband_disk * xmpi_bsize_dp & 2057 + 2 * npw_disk * nspinor_disk * xmpi_bsize_dp & 2058 + 4 * xmpio_bsize_frm) 2059 2060 call MPI_FILE_SET_VIEW(Wfk%fh,my_offset,MPI_BYTE,cgblock_type,'native',xmpio_info,mpierr) 2061 ABI_CHECK_MPI(mpierr, "SET_VIEW") 2062 2063 call MPI_TYPE_FREE(cgblock_type,mpierr) 2064 ABI_CHECK_MPI(mpierr, "TYPE_FREE") 2065 2066 bufsz = npw_disk * nspinor_disk * nb_block 2067 if (sc_mode == xmpio_collective) then 2068 call MPI_FILE_WRITE_ALL(Wfk%fh,cg_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2069 else if (sc_mode == xmpio_single) then 2070 call MPI_FILE_WRITE (Wfk%fh,cg_k,bufsz,MPI_DOUBLE_COMPLEX,MPI_STATUS_IGNORE,mpierr) 2071 else 2072 ABI_ERROR("Wrong sc_mode") 2073 end if 2074 ABI_CHECK_MPI(mpierr, "FILE_WRITE") 2075 end if 2076 2077 else 2078 ABI_ERROR("formeig not in [0,1]") 2079 end if 2080 #endif 2081 2082 case (IO_MODE_ETSF) 2083 if (present(kg_k)) then 2084 ! Write the reduced_coordinates_of_plane_waves for this k point. 2085 NCF_CHECK(nf90_inq_varid(wfk%fh, "reduced_coordinates_of_plane_waves", kg_varid)) 2086 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2087 NCF_CHECK(nctk_set_collective(wfk%fh, kg_varid)) 2088 end if 2089 ncerr = nf90_put_var(wfk%fh, kg_varid, kg_k, start=[1,1,ik_ibz], count=[3,npw_disk,1]) 2090 NCF_CHECK_MSG(ncerr, "putting kg_k") 2091 end if 2092 2093 ! Write eigenvalues and occupation factors. 2094 if (Wfk%formeig == 0) then 2095 2096 if (present(eig_k)) then 2097 NCF_CHECK(nf90_inq_varid(wfk%fh, "eigenvalues", eig_varid)) 2098 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2099 NCF_CHECK(nctk_set_collective(wfk%fh, eig_varid)) 2100 end if 2101 ncerr = nf90_put_var(wfk%fh, eig_varid, eig_k, start=[1,ik_ibz,spin], count=[nband_disk,1,1]) 2102 NCF_CHECK_MSG(ncerr, "putting eig_k") 2103 end if 2104 2105 if (present(occ_k)) then 2106 NCF_CHECK(nf90_inq_varid(wfk%fh, "occupations", occ_varid)) 2107 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2108 NCF_CHECK(nctk_set_collective(wfk%fh, occ_varid)) 2109 end if 2110 ncerr = nf90_put_var(wfk%fh, occ_varid, occ_k, start=[1,ik_ibz,spin], count=[nband_disk,1,1]) 2111 NCF_CHECK_MSG(ncerr, "putting occ_k") 2112 end if 2113 2114 else if (Wfk%formeig == 1) then 2115 if (present(occ_k)) then 2116 ABI_ERROR("Don't pass occ_k when formeig==1 and ETSF-IO") 2117 end if 2118 if (present(eig_k)) then 2119 !ABI_WARNING("Don't pass eig_k when formeig==1 and ETSF-IO") 2120 NCF_CHECK(nf90_inq_varid(wfk%fh, "h1_matrix_elements", h1_varid)) 2121 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2122 NCF_CHECK(nctk_set_collective(wfk%fh, h1_varid)) 2123 end if 2124 ncerr = nf90_put_var(wfk%fh, h1_varid, eig_k, start=[1,1,1,ik_ibz,spin], count=[2, nband_disk, nband_disk, 1, 1]) 2125 NCF_CHECK_MSG(ncerr, "puting h1mat_k") 2126 end if 2127 2128 else 2129 ABI_ERROR("formeig != [0,1]") 2130 end if 2131 2132 if (present(cg_k)) then 2133 ! Write the nb_block bands starting from band_block(1) 2134 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 2135 NCF_CHECK(nf90_inq_varid(wfk%fh, "coefficients_of_wavefunctions", cg_varid)) 2136 if (sc_mode == xmpio_collective .and. wfk%nproc > 1) then 2137 NCF_CHECK(nctk_set_collective(wfk%fh, cg_varid)) 2138 end if 2139 2140 ncerr = nf90_put_var(wfk%fh, cg_varid, cg_k, start=[1,1,1,band_block(1),ik_ibz,spin], & 2141 count=[2, npw_disk, wfk%nspinor, nb_block, 1, 1]) 2142 NCF_CHECK_MSG(ncerr, "putting cg_k") 2143 end if 2144 2145 case default 2146 ABI_ERROR(sjoin('Wrong value of iomode:', itoa(Wfk%iomode))) 2147 end select 2148 2149 !ABI_FREE(eig_buffer) 2150 !ABI_FREE(cg_buffer) 2151 2152 call cwtime_report(" wfk_write_band_block", cpu, wall, gflops) 2153 DBG_EXIT("COLL") 2154 2155 return 2156 2157 ! Handle Fortran IO error 2158 10 continue 2159 ABI_ERROR(errmsg) 2160 2161 end subroutine wfk_write_band_block
m_wfk/wfk_write_h1mat [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_write_h1mat
FUNCTION
Write all H1 matrix elements in the WFK file fname.
INPUTS
OUTPUTS
SOURCE
3531 subroutine wfk_write_h1mat(Wfk,sc_mode,eigen) 3532 3533 !Arguments ------------------------------------ 3534 !scalars 3535 integer,intent(in) :: sc_mode 3536 class(wfk_t),intent(inout) :: Wfk 3537 !arrays 3538 real(dp),intent(in) :: eigen(2*Wfk%mband**2*Wfk%nkpt*Wfk%nsppol) 3539 3540 !Local variables------------------------------- 3541 !scalars 3542 integer :: spin,ik_ibz,nband_k,ptr 3543 !arrays 3544 integer,parameter :: band_block00(2)=[0,0] 3545 3546 !************************************************************************ 3547 3548 ptr=1 3549 do spin=1,Wfk%nsppol 3550 do ik_ibz=1,Wfk%nkpt 3551 nband_k = Wfk%nband(ik_ibz,spin) 3552 call wfk%write_band_block(band_block00,ik_ibz,spin,sc_mode,eig_k=eigen(ptr:)) 3553 ptr = ptr + 2*nband_k**2 3554 end do 3555 end do 3556 3557 end subroutine wfk_write_h1mat
m_wfk/wfk_write_my_kptbands [ Functions ]
[ Top ] [ m_wfk ] [ Functions ]
NAME
wfk_write_my_kptbands
FUNCTION
From a cg (kg, eigen, occ) array write the corresponding file distributed bands on all procs, not just k-points
INPUTS
outpath_ = file name distrb_flags = logical mask for band, k, spins on this processor comm = mpi communicator formeig = flag for GS or response function format of eigenvalues kptns_in = requested k points, to be extracted from file or completed nkpt_in = number of requested k npwarr = array of number of plane waves at each k cg = plane wave coefficients kg = plane wave coordinates eigen = eigenvectors at all bands and my k occ = occupations of all bands at my k
OUTPUT
writes to file
SOURCE
3362 subroutine wfk_write_my_kptbands(outpath_, distrb_flags, comm, formeig, hdr,& 3363 iomode_, mband_in, mband_mem_in, mkmem_in, mpw_in, nkpt_in, nspinor_in, nsppol_in, & 3364 cg_in, kg_in, eigen, occ) 3365 3366 !Arguments ------------------------------------ 3367 !scalars 3368 integer, intent(in) :: comm, nkpt_in, formeig, iomode_ 3369 integer, intent(in) :: mband_in,mband_mem_in,mpw_in, nspinor_in, nsppol_in 3370 integer, intent(in) :: mkmem_in 3371 type(hdr_type),intent(in) :: hdr 3372 !arrays 3373 character(len=fnlen), intent(in) :: outpath_ 3374 logical, intent(in) :: distrb_flags(nkpt_in,mband_in,nsppol_in) 3375 3376 real(dp), intent(in), target :: cg_in(2,mpw_in*nspinor_in*mband_mem_in*mkmem_in*nsppol_in) 3377 integer, intent(in), target :: kg_in(3,mpw_in*mkmem_in) 3378 real(dp), intent(in) :: eigen((mband_in*(2*mband_in)**formeig)*nkpt_in*nsppol_in) 3379 real(dp), intent(in),optional :: occ(mband_in*nkpt_in*nsppol_in) 3380 3381 !Local variables------------------------------- 3382 !scalars 3383 integer :: spin,ik_rbz,nband_k 3384 integer :: npw_k 3385 integer :: iomode 3386 integer :: wfk_unt, iband, nband_me 3387 integer :: ii,jj,kk,ll 3388 integer, allocatable :: icg(:,:) 3389 integer, allocatable :: ikg(:) 3390 integer, allocatable :: ibdeig(:,:) 3391 integer, allocatable :: ibdocc(:,:) 3392 character(len=fnlen) :: outpath 3393 real(dp) :: cpu,wall,gflops 3394 type(wfk_t),target :: wfk_disk 3395 real(dp), pointer :: cg(:,:) 3396 integer, pointer :: kg(:,:) 3397 3398 ! ************************************************************************* 3399 3400 call cwtime(cpu, wall, gflops, "start") 3401 3402 ! if iomode ncdf check that outpath has the correct termination 3403 outpath = outpath_ 3404 iomode = iomode_ 3405 if (iomode==IO_MODE_ETSF .and. .not. endswith(outpath, ".nc")) then 3406 outpath = nctk_ncify(outpath) 3407 else 3408 ! adjust for mpiio if needed 3409 iomode = iomode_from_fname(outpath) 3410 end if 3411 3412 wfk_unt = get_unit() 3413 wfk_disk%debug = .true. 3414 call wfk_disk%open_write(hdr,outpath,formeig,iomode,wfk_unt,comm) !xmpi_comm_self) 3415 3416 ! no kpt on this proc, make local dummies for cg and kg 3417 if (mkmem_in == 0) then 3418 ABI_MALLOC(cg, (2,mpw_in)) 3419 ABI_MALLOC(kg, (3,mpw_in)) 3420 else 3421 cg => cg_in 3422 kg => kg_in 3423 end if 3424 3425 ABI_MALLOC(icg, (nkpt_in,nsppol_in)) 3426 ABI_MALLOC(ikg, (nkpt_in)) 3427 ABI_MALLOC(ibdeig, (nkpt_in,nsppol_in)) 3428 ABI_MALLOC(ibdocc, (nkpt_in,nsppol_in)) 3429 icg = 0 3430 ikg = 0 3431 ibdeig = 0 3432 ibdocc = 0 3433 ii = 0 3434 kk = 0 3435 ll = 0 3436 do spin=1,nsppol_in 3437 jj = 0 3438 do ik_rbz=1,nkpt_in 3439 ! this allows for reading fewer bands from disk than the disk version of nband 3440 nband_k = hdr%nband(ik_rbz+(spin-1)*hdr%nkpt) 3441 ibdeig(ik_rbz,spin) = kk 3442 ibdocc(ik_rbz,spin) = ll 3443 kk = kk+nband_k*(2*nband_k)**formeig 3444 ll = ll+nband_k 3445 3446 if (.not. any(distrb_flags(ik_rbz,:,spin))) cycle 3447 ! TODO: this does not take into account variable nband(ik) 3448 icg(ik_rbz,spin) = ii 3449 ikg(ik_rbz) = jj 3450 ! this allows for variable nband_k < mband_mem 3451 ii = ii+min(nband_k,mband_mem_in)*hdr%npwarr(ik_rbz)*nspinor_in 3452 jj = jj+hdr%npwarr(ik_rbz) 3453 end do 3454 end do 3455 3456 do spin=1,nsppol_in 3457 do ik_rbz=1,nkpt_in 3458 3459 nband_k = hdr%nband(ik_rbz+(spin-1)*hdr%nkpt) 3460 npw_k = hdr%npwarr(ik_rbz) 3461 3462 ! even if I do not have any bands to run, go through the mpio calls to avoid deadlocks 3463 ! if (.not. any(distrb_flags(ik_rbz,:,spin))) then 3464 ! ibdeig = ibdeig + nband_k*(2*nband_k)**formeig 3465 ! ibdocc = ibdocc + nband_k 3466 ! cycle 3467 ! end if 3468 3469 ! in case nband is not constant with k this creates chaos in the file writing 3470 ! as distrb_flags is allocated for mband, and true 3471 nband_me = min(count(distrb_flags(ik_rbz,:,spin)), nband_k) 3472 if (nband_me == 0) then 3473 iband = 1 ! does write_band_block accept the range [1,0]? 3474 else 3475 do iband = 1, nband_k 3476 if (distrb_flags(ik_rbz,iband,spin)) exit 3477 end do 3478 !TODO: check all nband_me entries in distrib_flags - the distribution could be random but with iband+nband_me-1 .true. 3479 if (.not. distrb_flags(ik_rbz,iband+nband_me-1,spin)) then 3480 ABI_ERROR("wfk_write_my_kptbands: bands not contiguous in distrb_flags") 3481 end if 3482 end if 3483 3484 if (present(occ)) then 3485 call wfk_disk%write_band_block([iband,iband+nband_me-1],ik_rbz,spin,xmpio_collective,& 3486 kg_k=kg(:,ikg(ik_rbz)+1:ikg(ik_rbz)+npw_k), & 3487 cg_k=cg(:,icg(ik_rbz,spin)+1:icg(ik_rbz,spin)+npw_k*nband_me*nspinor_in),& 3488 eig_k=eigen(ibdeig(ik_rbz,spin)+1:ibdeig(ik_rbz,spin)+nband_k*(2*nband_k)**formeig), & 3489 occ_k=occ(ibdocc(ik_rbz,spin)+1:ibdocc(ik_rbz,spin)+nband_k)) 3490 else 3491 call wfk_disk%write_band_block([iband,iband+nband_me-1],ik_rbz,spin,xmpio_collective,& 3492 kg_k=kg(:,ikg(ik_rbz)+1:ikg(ik_rbz)+npw_k), & 3493 cg_k=cg(:,icg(ik_rbz,spin)+1:icg(ik_rbz,spin)+npw_k*nband_me*nspinor_in),& 3494 eig_k=eigen(ibdeig(ik_rbz,spin)+1:ibdeig(ik_rbz,spin)+nband_k*(2*nband_k)**formeig)) 3495 end if 3496 3497 end do ! kpt 3498 end do ! sppol 3499 3500 call wfk_disk%close() 3501 3502 call cwtime_report(" wfk_write_my_kptbands. ", cpu, wall, gflops) 3503 3504 ABI_FREE(icg) 3505 ABI_FREE(ikg) 3506 ABI_FREE(ibdeig) 3507 ABI_FREE(ibdocc) 3508 if (mkmem_in == 0) then 3509 ABI_FREE(cg) 3510 ABI_FREE(kg) 3511 end if 3512 3513 end subroutine wfk_write_my_kptbands
m_wfkfile/wfk_compute_offsets [ Functions ]
NAME
wfk_compute_offsets
FUNCTION
Compute the offsets corresponding to the different sections of the file (G-vectors, eigenvalues, u(G). Needed only for Fortran-IO or MPI-IO.
SOURCE
3784 subroutine wfk_compute_offsets(Wfk) 3785 3786 !Arguments ------------------------------------ 3787 class(wfk_t),intent(inout) :: Wfk 3788 3789 !Local variables------------------------------- 3790 !scalars 3791 integer :: spin,ik_ibz,npw_k,nband_k,bsize_frm,mpi_type_frm,base !,band 3792 integer(XMPI_OFFSET_KIND) :: offset 3793 ! this variable is needed to force arithmetic in the right kind 3794 ! and avoid integer overflows with large nband npw. 3795 ! TODO: check if same is needed elsewhere for offsets 3796 integer(XMPI_OFFSET_KIND) :: increment 3797 3798 ! ************************************************************************* 3799 3800 select case (Wfk%iomode) 3801 case (IO_MODE_FORTRAN) 3802 ! Compute record number for Fortran IO 3803 ABI_MALLOC(Wfk%recn_ks, (Wfk%nkpt,Wfk%nsppol,REC_NUM)) 3804 3805 ! We start to count the number of Fortran records from the end of the Header 3806 ! Hence recn gives the relative position from the header, it's not an absolute position. 3807 base = 0 3808 do spin=1,Wfk%nsppol 3809 do ik_ibz=1,Wfk%nkpt 3810 nband_k = Wfk%nband(ik_ibz,spin) 3811 Wfk%recn_ks(ik_ibz,spin, REC_NPW) = base + 1 3812 Wfk%recn_ks(ik_ibz,spin, REC_KG) = base + 2 3813 Wfk%recn_ks(ik_ibz,spin, REC_EIG) = base + 3 3814 Wfk%recn_ks(ik_ibz,spin, REC_CG) = base + 4 3815 base = Wfk%recn_ks(ik_ibz,spin,REC_CG) 3816 if (Wfk%formeig==0) then 3817 ! add records for each cg (iband), and account for offset of 1 added for REC_NPW 3818 !TODO check if variable nband(k) works here 3819 base = base + nband_k - 1 3820 else if (Wfk%formeig==1) then 3821 ! add records for each eig1(:,iband) and cg (iband), and account for offset of 1 added for REC_NPW 3822 base = base + 2*(nband_k-1) 3823 else 3824 ABI_ERROR("formeig != [0,1]") 3825 end if 3826 end do 3827 end do 3828 3829 ! Save EOF position 3830 Wfk%recn_eof = base + 1 3831 3832 case (IO_MODE_MPI) 3833 ! Compute offsets for MPI-IO. 3834 ABI_MALLOC(Wfk%offset_ks, (Wfk%nkpt,Wfk%nsppol,REC_NUM)) 3835 3836 bsize_frm = xmpio_bsize_frm ! Byte length of the Fortran record marker. 3837 mpi_type_frm = xmpio_mpi_type_frm ! MPI type of the record marker. 3838 3839 ! The offset of the Header. TODO 3840 ! hdr_offset(Hdr) 3841 offset = Wfk%hdr_offset 3842 3843 do spin=1,Wfk%nsppol 3844 do ik_ibz=1,Wfk%nkpt 3845 npw_k = Wfk%Hdr%npwarr(ik_ibz) 3846 nband_k = Wfk%nband(ik_ibz,spin) 3847 !--------------------------------------------------------------------------- 3848 ! First record: npw, nspinor, nband_disk 3849 !--------------------------------------------------------------------------- 3850 Wfk%offset_ks(ik_ibz,spin,REC_NPW) = offset 3851 3852 if (Wfk%Hdr%headform>=40) then 3853 ! npw, nspinor, nband_disk 3854 offset = offset + 3*xmpi_bsize_int + 2*bsize_frm 3855 else 3856 ABI_ERROR("Old headforms < 40 are not supported") 3857 end if 3858 Wfk%offset_ks(ik_ibz,spin,REC_KG) = offset 3859 3860 !--------------------------------------------------------------------------- 3861 ! Second record: (k+G) vectors 3862 ! kg_k(1:3,1:npw_k) 3863 !--------------------------------------------------------------------------- 3864 offset = offset + 3*npw_k*xmpi_bsize_int + 2*bsize_frm 3865 Wfk%offset_ks(ik_ibz,spin,REC_EIG) = offset 3866 ! 3867 !--------------------------------------------------------------------------- 3868 ! Third record: eigenvalues 3869 !--------------------------------------------------------------------------- 3870 if (Wfk%formeig==0) then 3871 ! eigen(1:nband_k), occ(1:nband_k) 3872 !offset = offset + 2*Wfk%mband*xmpi_bsize_dp + 2*bsize_frm 3873 offset = offset + 2*nband_k*xmpi_bsize_dp + 2*bsize_frm 3874 Wfk%offset_ks(ik_ibz,spin,REC_CG) = offset 3875 3876 !--------------------------------------------------------------------------- 3877 ! Fourth record: Wavefunction coefficients 3878 !--------------------------------------------------------------------------- 3879 ! do band=1,nband_k; write(unitwf) cg_k(1:2,npw_k*nspinor); end do 3880 !offset = offset + Wfk%mband * (2*npw_k*Wfk%nspinor*xmpi_bsize_dp + 2*bsize_frm) 3881 increment = 2*npw_k*Wfk%nspinor*xmpi_bsize_dp + 2*bsize_frm 3882 increment = nband_k * increment 3883 offset = offset + increment 3884 !offset = offset + nband_k * (2*npw_k*Wfk%nspinor*xmpi_bsize_dp + 2*bsize_frm) 3885 3886 else if (Wfk%formeig==1) then 3887 ! read(unitwf) eigen(2*nband_k) 3888 !Wfk%offset_ks(ik_ibz,spin,REC_CG) = offset + 2*Wfk%mband*xmpi_bsize_dp + 2*bsize_frm 3889 Wfk%offset_ks(ik_ibz,spin,REC_CG) = offset + 2*nband_k*xmpi_bsize_dp + 2*bsize_frm 3890 3891 !--------------------------------------------------------------------------- 3892 ! Fourth record: Wavefunction coefficients 3893 !--------------------------------------------------------------------------- 3894 increment = (2*npw_k*Wfk%nspinor*xmpi_bsize_dp + 2*bsize_frm) + & 3895 (2*nband_k*xmpi_bsize_dp + 2*bsize_frm) 3896 !Wfk%mband * (2*npw_k*Wfk%nspinor*xmpi_bsize_dp + 2*bsize_frm) + & 3897 !Wfk%mband * (2*Wfk%mband*xmpi_bsize_dp + 2*bsize_frm) 3898 increment = nband_k * increment 3899 offset = offset + increment 3900 3901 else 3902 ABI_ERROR("Wrong formeig") 3903 end if 3904 3905 end do ! ik_ibz 3906 end do ! spin 3907 3908 ! Save EOF offset 3909 Wfk%offset_eof = offset 3910 3911 ! Check for possible wraparound errors. 3912 if (ANY(Wfk%offset_ks <= 0) .or. Wfk%offset_eof < 0) then 3913 ABI_ERROR("Found negative offset. File too large for MPI-IO!!!") 3914 end if 3915 end select 3916 3917 if (Wfk%debug) call wfk_show_offsets(Wfk) 3918 3919 end subroutine wfk_compute_offsets
m_wfkfile/wfk_show_offsets [ Functions ]
NAME
wfk_show_offsets
FUNCTION
Print the offsets.
SOURCE
3933 subroutine wfk_show_offsets(Wfk) 3934 3935 !Arguments ------------------------------------ 3936 class(wfk_t),intent(inout) :: Wfk 3937 3938 !Local variables------------------------------- 3939 !scalars 3940 integer :: spin,ik_ibz 3941 3942 ! ************************************************************************* 3943 3944 select case (Wfk%iomode) 3945 3946 case (IO_MODE_FORTRAN) 3947 write(std_out,*)"Record number relative to the header." 3948 do spin=1,Wfk%nsppol 3949 do ik_ibz=1,Wfk%nkpt 3950 write(std_out,"(a,2(i0,2x),a,4(a,i0,a))") & 3951 "(ik_ibz, spin) ",ik_ibz,spin,ch10, & 3952 " recn(REC_NPW): ",Wfk%recn_ks(ik_ibz,spin,REC_NPW),ch10,& 3953 " recn(REC_KG) : ",Wfk%recn_ks(ik_ibz,spin,REC_KG), ch10,& 3954 " recn(REC_EIG): ",Wfk%recn_ks(ik_ibz,spin,REC_EIG),ch10,& 3955 " recn(REC_CG) : ",Wfk%recn_ks(ik_ibz,spin,REC_CG),ch10 3956 end do 3957 end do 3958 3959 write(std_out,"(a,i0)")"EOS position: ",Wfk%recn_eof 3960 3961 case (IO_MODE_MPI) 3962 write(std_out,"(a,i0)")"hdr_offset ",Wfk%hdr_offset 3963 3964 do spin=1,Wfk%nsppol 3965 do ik_ibz=1,Wfk%nkpt 3966 write(std_out,"(a,2(i0,2x),a,4(a,i0,a))") & 3967 "(ik_ibz, spin) ",ik_ibz,spin,ch10, & 3968 " offset(REC_NPW): ",Wfk%offset_ks(ik_ibz,spin,REC_NPW),ch10,& 3969 " offset(REC_KG) : ",Wfk%offset_ks(ik_ibz,spin,REC_KG), ch10,& 3970 " offset(REC_EIG): ",Wfk%offset_ks(ik_ibz,spin,REC_EIG),ch10,& 3971 " offset(REC_CG) : ",Wfk%offset_ks(ik_ibz,spin,REC_CG),ch10 3972 end do ! ik_ibz 3973 end do ! spin 3974 ! 3975 ! Write EOF position 3976 write(std_out,"(a,i0)")"offset_eof ",Wfk%offset_eof 3977 end select 3978 3979 end subroutine wfk_show_offsets