1!>========================================================================= 2!! 3!! Module: 4!! 5!! epswrite_hdf5_m Originally by JRD Last Modified 12/2014 (FHJ) 6!! 7!! Routines to write header info for epsmat files in HDF5 format. 8!! 9!!========================================================================= 10 11#include "f_defs.h" 12 13module epswrite_hdf5_m 14#ifdef HDF5 15 use hdf5 16 use global_m 17 use hdf5_io_m 18 use wfn_io_hdf5_m 19 implicit none 20 21 private 22 23 public :: & 24 eps_hdf5_setup, & 25 set_subspace_neigen_q, & 26 set_qpt_done, & 27 is_qpt_done 28 29contains 30 31 32subroutine eps_hdf5_setup(kp, gvec, syms, crys, pol, qgrid, nq, qpts, nmtx, & 33 nmtx_max, nband, name, restart) 34 type(kpoints), intent(in) :: kp 35 type(gspace), intent(in) :: gvec 36 type(symmetry), intent(in) :: syms 37 type(crystal), intent(in) :: crys 38 type(polarizability), intent(in) :: pol 39 integer, intent(in) :: qgrid(3) 40 integer, intent(in) :: nq 41 real(DP), intent(in) :: qpts(:,:) !< (3,nq) 42 integer, intent(in) :: nmtx(:) !< (nq) 43 integer, intent(in) :: nmtx_max 44 integer, intent(in) :: nband 45 character(len=*), intent(in) :: name 46 logical, intent(inout), optional :: restart 47 48 integer(HID_T) :: file_id 49 integer(HID_T) :: dspace 50 integer(HID_T) :: dset_id 51 integer(HSIZE_T) :: dims(6) 52 53 integer :: error, ii, intdata(1) 54 integer :: neigen_of_q(nq) 55 logical :: qpt_done(nq) 56 real(DP), allocatable :: tmp(:,:) 57 real(DP) :: freqs_tmp(2,pol%nfreq) 58 real(DP) :: realdata(1) 59 logical :: restart_, file_exists, file_ok 60 character(len=3) :: sheader='WFN' 61 62 PUSH_SUB(eps_hdf5_setup) 63 64 restart_=.false. 65 if (present(restart)) restart_ = restart 66 67 ! FHJ: try to restart the calculation, if possible and desired. 68 ! We ignore the restart flags if the file doesn`t exist. However, the code 69 ! dies if the file exists and is incompatible or looks corrupted. 70 if (restart_) then 71 call try_restart() 72 if (file_ok) return 73 endif 74 75 ! FHJ: Set up file: write MF header and create groups 76 write(6,'(1x,2a)') "Initializing ", trim(name) 77 call setup_hdf5_mf_file(trim(name)) 78 call write_hdf5_header_type(trim(name), sheader, SCALARSIZE, kp, gvec, syms, crys) 79 call write_hdf5_gvectors(trim(name), gvec%ng, gvec%components) 80 call h5fopen_f(trim(name), H5F_ACC_RDWR_F, file_id, error) 81 call hdf5_create_group(file_id, 'eps_header', error) 82 call hdf5_create_group(file_id, 'eps_header/params', error) 83 call hdf5_create_group(file_id, 'eps_header/qpoints', error) 84 call hdf5_create_group(file_id, 'eps_header/freqs', error) 85 call hdf5_create_group(file_id, 'eps_header/gspace', error) 86 call hdf5_create_group(file_id, 'mats', error) 87 88 if( pol%subspace ) then 89 call hdf5_create_group(file_id, 'eps_header/subspace', error) 90 endif 91 92 call hdf5_write_int(file_id, 'eps_header/versionnumber', VER_EPS_HDF5, error) 93 call hdf5_write_int(file_id, 'eps_header/flavor', SCALARSIZE, error) 94 95 ! FHJ: General datasets 96 call hdf5_write_int(file_id, 'eps_header/params/matrix_type', pol%matrix_type, error) 97 call hdf5_write_logical(file_id, 'eps_header/params/has_advanced', pol%has_advanced, error) 98 call hdf5_write_int(file_id, 'eps_header/params/nmatrix', pol%nmatrix, error) 99 call hdf5_write_int(file_id, 'eps_header/params/matrix_flavor', pol%matrix_flavor, error) 100 call hdf5_write_int(file_id, 'eps_header/params/icutv', pol%icutv, error) 101 call hdf5_write_double(file_id, 'eps_header/params/ecuts', pol%ecuts, error) 102 call hdf5_write_int(file_id, 'eps_header/params/nband', nband, error) 103 call hdf5_write_double(file_id, 'eps_header/params/efermi', pol%efermi/ryd, error) 104 call hdf5_write_int(file_id, 'eps_header/params/intraband_flag', pol%intraband_flag, error) 105 call hdf5_write_double(file_id, 'eps_header/params/intraband_overlap_min', pol%intraband_overlap_min, error) 106 call hdf5_write_logical(file_id, 'eps_header/params/subsample', pol%subsample, error) 107 call hdf5_write_logical(file_id, 'eps_header/params/subspace', pol%subspace, error) 108 109 ! FHJ: Q-points-related datasets 110 qpt_done(:) = .false. 111 call hdf5_write_int(file_id, 'eps_header/qpoints/nq', nq, error) 112 call hdf5_write_double_array(file_id, 'eps_header/qpoints/qpts', (/3,nq/), qpts, error) 113 call hdf5_write_int_array(file_id, 'eps_header/qpoints/qgrid', (/3/), qgrid, error) 114 call hdf5_write_logical_array(file_id, 'eps_header/qpoints/qpt_done', (/nq/), qpt_done, error) 115 116 ! FHJ: Frequency-related datasets 117 call hdf5_write_int(file_id, 'eps_header/freqs/freq_dep', pol%freq_dep, error) 118 call hdf5_write_int(file_id, 'eps_header/freqs/nfreq', pol%nfreq, error) 119 call hdf5_write_int(file_id, 'eps_header/freqs/nfreq_imag', pol%nfreq_imag, error) 120 do ii=1,pol%nfreq !FHJ: TODO - have an unique complex freq grid in the future! 121 freqs_tmp(1,ii) = pol%dFreqGrid(ii) + dble(pol%dFreqBrd(ii)) 122 freqs_tmp(2,ii) = IMAG(pol%dFreqBrd(ii)) 123 enddo 124 call hdf5_write_double_array(file_id, 'eps_header/freqs/freqs', (/2, pol%nfreq/), freqs_tmp, error) 125 126 ! FHJ: G-vectors-related datasets 127 call hdf5_write_int_array(file_id, 'eps_header/gspace/nmtx', (/nq/), nmtx, error) 128 call hdf5_write_int(file_id, 'eps_header/gspace/nmtx_max', nmtx_max, error) 129 call hdf5_create_dset(file_id, 'eps_header/gspace/ekin', H5T_NATIVE_DOUBLE, (/gvec%ng,nq/), error) 130 ! FHJ: Epsmat version 3+ includes the bare Coulomb interaction as well. 131 ! We use this to get the advanced epsilon from the retarded matrix. 132 call hdf5_create_dset(file_id, 'eps_header/gspace/vcoul', H5T_NATIVE_DOUBLE, (/nmtx_max,nq/), error) 133 call hdf5_create_dset(file_id, 'eps_header/gspace/gind_eps2rho', H5T_NATIVE_INTEGER, (/gvec%ng,nq/), error) 134 call hdf5_create_dset(file_id, 'eps_header/gspace/gind_rho2eps', H5T_NATIVE_INTEGER, (/gvec%ng,nq/), error) 135 136 !DVF: subspace approximation-related datasets 137 if( pol%subspace ) then 138 call hdf5_write_logical(file_id, 'eps_header/subspace/keep_full_eps_static', pol%keep_full_eps_static, error) 139 call hdf5_write_logical(file_id, 'eps_header/subspace/matrix_in_subspace_basis', pol%matrix_in_subspace_basis, error) 140 call hdf5_write_double(file_id, 'eps_header/subspace/eps_eigenvalue_cutoff', pol%chi_eigenvalue_cutoff, error) 141 call hdf5_write_int(file_id, 'eps_header/subspace/neig_max', pol%neig_sub_input, error) 142 ! use same strategy as qpt_done, since the number of eigenvector will never exceed neig_sub_input, we 143 ! initially set all of them to neig_sub_input and we update the actual value when calculated 144 neigen_of_q(:) = pol%neig_sub_input 145 call hdf5_write_int_array(file_id, 'eps_header/subspace/neig', (/nq/), neigen_of_q, error) 146 endif 147 148 ! FHJ: Matrix-elements-related datasets 149 call hdf5_create_dset(file_id, 'mats/matrix', H5T_NATIVE_DOUBLE, & 150 (/pol%matrix_flavor,nmtx_max,nmtx_max,pol%nfreq,pol%nmatrix,nq/), error) 151 call hdf5_create_dset(file_id, 'mats/matrix-diagonal', H5T_NATIVE_DOUBLE, & 152 (/pol%matrix_flavor,nmtx_max,nq/), error) 153 154 if( pol%subspace .and. pol%matrix_in_subspace_basis ) then 155 ! 156 ! here pol%nmatrix should always be ONE, this may change if we want to write 157 ! the subspace chi, which can be spin polarized, which not sure if make any sense. 158 ! 159 call hdf5_create_dset(file_id, 'mats/matrix_fulleps0', H5T_NATIVE_DOUBLE, & 160 (/pol%matrix_flavor,nmtx_max,nmtx_max,1,pol%nmatrix,nq/), error) 161 !XXX 162 !XXX actual (max) size 163 call hdf5_create_dset(file_id, 'mats/matrix_eigenvec', H5T_NATIVE_DOUBLE, & 164 (/pol%matrix_flavor,nmtx_max,pol%neig_sub_input,1,pol%nmatrix,nq/), error) 165 !XXX 166 !XXX full matrix 167 !XXX call hdf5_create_dset(file_id, 'mats/matrix_eigenvec', H5T_NATIVE_DOUBLE, & 168 !XXX (/pol%matrix_flavor,nmtx_max,nmtx_max,1,pol%nmatrix,nq/), error) 169 !XXX 170 call hdf5_create_dset(file_id, 'mats/matrix_subspace', H5T_NATIVE_DOUBLE, & 171 (/pol%matrix_flavor,pol%neig_sub_input,pol%neig_sub_input,pol%nfreq,pol%nmatrix,nq/), error) 172 end if 173 174 call h5fclose_f(file_id, error) 175 176 POP_SUB(eps_hdf5_setup) 177 178contains 179 180 subroutine try_restart() 181 integer :: nspin_old, nq_old 182 real(DP) :: qpts_old(3,nq) 183 integer :: freq_dep_old, nfreq_old 184 integer :: ng_old, gvecs_old(3,gvec%ng), nmtx_old(nq), matrix_flavor_old, nmatrix_old 185 integer :: neig_max_old 186 logical :: subspace_exists, subspace_old, matrix_in_subspace_basis_old, keep_full_eps_static_old 187 logical :: subspace_error 188 189 subspace_error = .false. 190 write(6,'(1x,2a)') "Trying to restart file ", trim(name) 191 inquire(file=trim(name), exist=file_exists) 192 if (file_exists) then 193 call h5fopen_f(trim(name), H5F_ACC_RDONLY_F, file_id, error) 194 if (error==0) then 195 ! FHJ: Consistency check. 196 call h5lexists_f(file_id, 'eps_header/qpoints/qpt_done', file_ok, error) 197 if (file_ok) then 198 call hdf5_read_int(file_id, 'mf_header/kpoints/nspin', nspin_old, error) 199 if (error==0) call hdf5_read_int(file_id, & 200 'eps_header/qpoints/nq', nq_old, error) 201 if (error==0) call hdf5_read_double_array(file_id, & 202 'eps_header/qpoints/qpts', (/3,nq/), qpts_old, error) 203 if (error==0) call hdf5_read_int(file_id, & 204 'eps_header/freqs/freq_dep', freq_dep_old, error) 205 if (error==0) call hdf5_read_int(file_id, & 206 'eps_header/freqs/nfreq', nfreq_old, error) 207 if (error==0) call hdf5_read_int(file_id, & 208 'mf_header/gspace/ng', ng_old, error) 209 if (error==0) call hdf5_read_int_array(file_id, & 210 'mf_header/gspace/components', (/3,gvec%ng/), gvecs_old, error) 211 if (error==0) call hdf5_read_int_array(file_id, & 212 'eps_header/gspace/nmtx', (/nq/), nmtx_old, error) 213 if (error==0) call hdf5_read_int(file_id, & 214 'eps_header/params/matrix_flavor', matrix_flavor_old, error) 215 if (error==0) call hdf5_read_int(file_id, & 216 'eps_header/params/nmatrix', nmatrix_old, error) 217 file_ok = (error==0) .and. (nspin_old==kp%nspin .and. nq_old==nq .and. & 218 all(dabs(qpts_old-qpts)<TOL_ZERO) .and. freq_dep_old==pol%freq_dep .and. & 219 nfreq_old==pol%nfreq .and. ng_old==gvec%ng .and. all(gvecs_old==gvec%components) .and. & 220 all(nmtx_old==nmtx) .and. matrix_flavor_old==pol%matrix_flavor .and. & 221 nmatrix_old==pol%nmatrix) 222 ! additional check in case of subspace 223 if ( file_ok ) then 224 ! check the existence of the parameter (this ensure that the routine will 225 ! work also for old .h5 files) 226 call h5lexists_f(file_id, 'eps_header/params/subspace', subspace_exists, error) 227 if (subspace_exists .and. error == 0) then 228 call hdf5_read_logical(file_id, 'eps_header/params/subspace', subspace_old, error) 229 ! eps_header/params/subspace exists, check new and old flags 230 file_ok = (error==0) .and. (pol%subspace .eqv. subspace_old) 231 if( file_ok ) then 232 if( subspace_old .and. error == 0) then 233 ! if the flags are equivalent and we are doing a subspace calculation 234 ! check the remaining parameters 235 ! this complication is due to the fact that the subspace section is not 236 ! created in non subspace calculation... 237 call hdf5_read_int(file_id, 'eps_header/subspace/neig_max', neig_max_old, error) 238 if (error==0) call hdf5_read_logical(file_id, & 239 'eps_header/subspace/matrix_in_subspace_basis', matrix_in_subspace_basis_old, error) 240 if (error==0) call hdf5_read_logical(file_id, & 241 'eps_header/subspace/keep_full_eps_static', keep_full_eps_static_old, error) 242 file_ok = (error==0) & 243 .and. (neig_max_old == pol%neig_sub_input) & 244 .and. (matrix_in_subspace_basis_old .eqv. pol%matrix_in_subspace_basis) & 245 .and. (keep_full_eps_static_old .eqv. pol%keep_full_eps_static) 246 if( .not. file_ok ) subspace_error = .true. 247 end if 248 else 249 if(pol%subspace) then 250 write(6,'(1x,A)') "ERROR: Trying to restart a subspace calc from a non-subspace calc" 251 else 252 write(6,'(1x,A)') "ERROR: Trying to restart a non-subspace calc from a subspace calc" 253 end if 254 subspace_error = .true. 255 file_ok = .false. 256 end if 257 else 258 if (pol%subspace) then 259 ! eps_header/params/subspace don`t exist but we are trying to restart a subspace calc. 260 write(6,'(1x,A)') "ERROR: Not a subspace initialized .h5 file" 261 subspace_error = .true. 262 file_ok = .false. 263 end if 264 end if 265 end if 266 endif 267 call h5fclose_f(file_id, error) 268 if (file_ok) then 269 ! FHJ: File *seems* alright, we don`t have to initialize it 270 write(6,'(1x,2a)') "Everything looks ok: restarting file ", trim(name) 271 return 272 endif 273 write(0,*) 274 write(0,'(3a)') "ERROR: file ", trim(name), " is incompatible with the current calculation." 275 if( subspace_error ) then 276 write(0,*) "ERROR is most likely due to incompatible .h5 files within subspace approximation" 277 end if 278 write(0,*) 'Values from file vs. calculation:' 279 write(0,*) 'nspin:', nspin_old, kp%nspin 280 write(0,*) 'nq:', nq_old, nq 281 write(0,*) 'qpts (same?):', all(dabs(qpts_old-qpts)<TOL_ZERO) 282 write(0,*) 'freq_dep:', freq_dep_old, pol%freq_dep 283 write(0,*) 'nfreq:', nfreq_old, pol%nfreq 284 285 write(0,*) 'ng:', ng_old, gvec%ng 286 write(0,*) 'gvecs (same?):', all(gvecs_old==gvec%components) 287 write(0,*) 'nmtx (same?):', all(nmtx_old==nmtx) 288 write(0,*) 'matrix_flavor:', matrix_flavor_old, pol%matrix_flavor 289 write(0,*) 'nmatrix:', nmatrix_old, pol%nmatrix 290 write(0,*) 291 write(0,*) 'NOTE: you should only trust the first pair of values that disagree.' 292 write(0,*) 293 call die("file "//trim(name)//" is incompatible with the current calculation.", only_root_writes=.true.) 294 else 295 write(0,'(3a,i0,a)') "ERROR: ", trim(name), " is not a valid HDF5 file (error code: ", error," )." 296 write(0,'(a)') 'Make sure the file is not corrupted' 297 call die("file "//trim(name)//" looks corrupted", only_root_writes=.true.) 298 endif 299 endif 300 file_ok = .false. 301 write(0,'(3a)') "WARNING: file ", trim(name), " doesn`t exist. We`ll start the calculation from scratch." 302 restart_ = .false. 303 if (present(restart)) restart = .false. 304 305 end subroutine try_restart 306 307end subroutine eps_hdf5_setup 308 309 310subroutine setup_subspace_mats_hdf5(fname, pol,iq) 311 character(len=*), intent(in) :: fname 312 integer, intent(in) :: iq 313 type(polarizability), intent(in) :: pol 314 315 integer(HID_T) :: file_id 316 integer :: error 317 character(len=20) :: subspace_name 318! DVF: we create different subgroups in the mats group for the 319! subspace matrices and basis vectors at each q-point. This is 320! necessary because we don`t know the size of the subspace basis 321! at each q-point before we actually do the calculation like we do 322! for the reciprocal space basis. So, we need to allocate the matrices 323! right after we have run scalapack/ELPA, and then create the needed 324! subgroups and dsets. That is the purpose of this routine. 325 326! The above description assumes we specify the eigenvector tolerance and 327! not the number of eigenvectors. Specifying the tolerance is best practice 328! since the eigenvector tolerance is directly related to the convergence of 329! one`s calculation, and specifying the tolerance requires no prior knowledge 330! of the size of the g-space, and hence is more automatic for users. 331 332 PUSH_SUB(setup_subspace_mats_hdf5) 333 334 call h5fopen_f(trim(fname), H5F_ACC_RDWR_F, file_id, error) 335 ! DVF: subspace matrix-elements-related datasets 336 write(subspace_name,'(a14,i4)') "mats/subspace_",iq 337 call hdf5_create_group(file_id, trim(subspace_name), error) 338 call hdf5_create_dset(file_id, trim(subspace_name) // '/matrix_sub', H5T_NATIVE_DOUBLE, & 339 (/pol%matrix_flavor,pol%neig_sub,pol%neig_sub,pol%nfreq,pol%nmatrix,1/), error) 340 call hdf5_create_dset(file_id, trim(subspace_name) // '/basis_sub', H5T_NATIVE_DOUBLE, & 341 (/pol%matrix_flavor,pol%nmtx,pol%neig_sub,1,pol%nmatrix,1/), error) 342 call h5fclose_f(file_id, error) 343 344 POP_SUB(setup_subspace_mats_hdf5) 345 346end subroutine setup_subspace_mats_hdf5 347 348subroutine set_subspace_neigen_q(fname, iq, neigen) 349 character(len=*), intent(in) :: fname 350 integer, intent(in) :: iq, neigen 351 352 integer(HID_T) :: file_id 353 integer :: nq, error 354 integer, allocatable :: neigen_of_q(:) 355 356 PUSH_SUB(set_subspace_neigen_q) 357 358 call open_file(99, trim(fname), status='old') 359 call close_file(99) 360 361 call h5fopen_f(trim(fname), H5F_ACC_RDWR_F, file_id, error) 362 call hdf5_read_int(file_id, 'eps_header/qpoints/nq', nq, error) 363 SAFE_ALLOCATE(neigen_of_q, (nq)) 364 365 call hdf5_read_int_array(file_id, 'eps_header/subspace/neig', (/nq/), neigen_of_q, error) 366 neigen_of_q(iq) = neigen 367 call hdf5_write_int_array(file_id, 'eps_header/subspace/neig', (/nq/), neigen_of_q, error) 368 369 call h5fclose_f(file_id, error) 370 SAFE_DEALLOCATE(neigen_of_q) 371 372 POP_SUB(set_subspace_neigen_q) 373 374end subroutine set_subspace_neigen_q 375 376subroutine set_qpt_done(fname, iq) 377 character(len=*), intent(in) :: fname 378 integer, intent(in) :: iq 379 380 integer(HID_T) :: file_id 381 integer :: nq, error 382 logical, allocatable :: qpt_done(:) 383 384 PUSH_SUB(set_qpt_done) 385 386 call open_file(99, trim(fname), status='old') 387 call close_file(99) 388 389 call h5fopen_f(trim(fname), H5F_ACC_RDWR_F, file_id, error) 390 call hdf5_read_int(file_id, 'eps_header/qpoints/nq', nq, error) 391 SAFE_ALLOCATE(qpt_done, (nq)) 392 call hdf5_read_logical_array(file_id, 'eps_header/qpoints/qpt_done', (/nq/), qpt_done, error) 393 qpt_done(iq) = .true. 394 call hdf5_write_logical_array(file_id, 'eps_header/qpoints/qpt_done', (/nq/), qpt_done, error) 395 call h5fclose_f(file_id, error) 396 SAFE_DEALLOCATE(qpt_done) 397 398 POP_SUB(set_qpt_done) 399 400end subroutine set_qpt_done 401 402 403logical function is_qpt_done(fname, iq) 404 character(len=*), intent(in) :: fname 405 integer, intent(in) :: iq 406 407 integer(HID_T) :: file_id 408 integer :: nq, error 409 logical, allocatable :: qpt_done(:) 410 411 PUSH_SUB(is_qpt_done) 412 413 call open_file(99, trim(fname), status='old') 414 call close_file(99) 415 416 call h5fopen_f(trim(fname), H5F_ACC_RDONLY_F, file_id, error) 417 call hdf5_read_int(file_id, 'eps_header/qpoints/nq', nq, error) 418 SAFE_ALLOCATE(qpt_done, (nq)) 419 call hdf5_read_logical_array(file_id, 'eps_header/qpoints/qpt_done', (/nq/), qpt_done, error) 420 is_qpt_done = qpt_done(iq) 421 call h5fclose_f(file_id, error) 422 SAFE_DEALLOCATE(qpt_done) 423 424 POP_SUB(is_qpt_done) 425 426end function is_qpt_done 427 428#endif 429end module epswrite_hdf5_m 430