Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_eig66_hdf
       8             : use m_juDFT
       9             :    !*****************************************************************
      10             :    ! DESC:Module for hdf-io of eig-file
      11             :    !      To be compatible with f90 interface of HDF, use kind for vars
      12             :    !
      13             :    !      !ATTENTION before calling openeig and after calling closeeig!
      14             :    !      !the hdf library has to be initialized or finalized, respectively
      15             :    !
      16             :    !      CONTAINS the following subroutines:
      17             :    !      openeig        opens file
      18             :    !      closeeig       closes file
      19             :    !      read_keb       reads kpt, enpara and basis data
      20             :    !      read_neig      read no of eigenvalues (and eigenvalues itself)
      21             :    !      read_eig       reads eigenvectors
      22             :    !      writeeig       saves all data for kpt
      23             :    !      writesingleeig saves data for one kpt and energy
      24             :    !
      25             :    !
      26             :    !                          Daniel Wortmann
      27             :    !*****************************************************************
      28             :    USE m_eig66_data
      29             :    USE m_types_mat
      30             : #ifdef CPP_HDF
      31             :    USE hdf5
      32             :    USE m_hdf_tools
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             :    INTEGER, PARAMETER :: one = 1, two = 2, three = 3, zero = 0
      37             :    !to have the correct
      38             :    !type for array constructors
      39             : 
      40             : #endif
      41             :    PUBLIC open_eig, close_eig
      42             :    PUBLIC read_eig
      43             :    PUBLIC write_eig!,writesingleeig,writeeigc,writebas
      44             : 
      45             : CONTAINS
      46           0 :    SUBROUTINE priv_find_data(id, d)
      47             :       INTEGER, INTENT(IN)::id
      48             :       TYPE(t_data_hdf), POINTER:: d
      49             : 
      50             :       CLASS(t_data), POINTER   ::dp
      51           0 :       CALL eig66_find_data(dp, id)
      52             :       SELECT TYPE (dp)
      53             :       TYPE is (t_data_hdf)
      54           0 :          d => dp
      55             :       CLASS default
      56           0 :          CALL judft_error("BUG: wrong datatype in eig66_hdf")
      57             :       END SELECT
      58           0 :    END SUBROUTINE priv_find_data
      59             :    !----------------------------------------------------------------------
      60           0 :    SUBROUTINE open_eig(id, fmpi_comm, nmat, neig, nkpts, jspins, create, l_real, l_soc, readonly, l_olap, filename)
      61             : 
      62             :       !*****************************************************************
      63             :       !     opens hdf-file for eigenvectors+values
      64             :       !*****************************************************************
      65             : #ifdef CPP_HDFMPI
      66             :       USE mpi
      67             : #endif
      68             :       IMPLICIT NONE
      69             : 
      70             :       INTEGER, INTENT(IN) :: id, fmpi_comm
      71             :       INTEGER, INTENT(IN) :: nmat, neig, nkpts, jspins
      72             :       LOGICAL, INTENT(IN) :: create, readonly, l_real, l_soc, l_olap
      73             :       CHARACTER(LEN=*), OPTIONAL :: filename
      74             : 
      75             : #ifdef CPP_HDF
      76             : 
      77             :       INTEGER         :: hdferr, access_mode
      78             :       INTEGER(HID_T)  :: creation_prp, access_prp, spaceid
      79             :       LOGICAL         :: l_exist
      80             :       INTEGER(HSIZE_T):: dims(7)
      81             :       TYPE(t_data_HDF), POINTER::d
      82             :       !Set creation and access properties
      83             : 
      84             : #ifdef CPP_HDFMPI
      85           0 :       IF (readonly) THEN
      86           0 :          access_prp = H5P_DEFAULT_f
      87           0 :          creation_prp = H5P_DEFAULT_f
      88             :       ELSE
      89           0 :          CALL h5pcreate_f(H5P_FILE_ACCESS_F, access_prp, hdferr)
      90             :          !      CALL h5pset_fapl_mpiposix_f(access_prp,fmpi_comm,
      91             :          !     +.false.,hdferr)
      92           0 :          CALL h5pset_fapl_mpio_f(access_prp, fmpi_comm, MPI_INFO_NULL, hdferr)
      93           0 :          creation_prp = H5P_DEFAULT_f !no special creation property
      94             :       ENDIF
      95             : #else
      96             :       access_prp = H5P_DEFAULT_f
      97             :       creation_prp = H5P_DEFAULT_f
      98             : #endif
      99             : 
     100           0 :       if(l_olap) call juDFT_error("olap not implemented for hdf5")
     101             : 
     102           0 :       CALL priv_find_data(id, d)
     103           0 :       IF (PRESENT(filename)) d%fname = filename
     104           0 :       CALL eig66_data_storedefault(d, jspins, nkpts, nmat, neig, l_real, l_soc)
     105             :       !set access_flags according
     106           0 :       IF (readonly) THEN
     107           0 :          access_mode = H5F_ACC_RDONLY_F
     108             :       ELSE
     109           0 :          access_mode = H5F_ACC_RDWR_F
     110             :       ENDIF
     111             :       !     OPEN FILE and get D%FID's
     112           0 :       IF (create) THEN
     113           0 :          INQUIRE (FILE=TRIM(d%fname)//'.hdf', EXIST=l_exist)
     114           0 :          access_mode = H5F_ACC_TRUNC_F
     115             :          !         IF (l_exist) WRITE (*,*)'Warning: eig.hdf was overwritten'
     116           0 :          CALL h5fcreate_f(TRIM(d%fname)//'.hdf', access_Mode, d%fid, hdferr, creation_prp, access_prp)
     117             :          ! create dataspaces and datasets
     118             :          !   scalars
     119           0 :          dims(:2) = (/nkpts, jspins/)
     120           0 :          CALL h5screate_simple_f(2, dims(:2), spaceid, hdferr)
     121           0 :          CALL h5dcreate_f(d%fid, "neig", H5T_NATIVE_INTEGER, spaceid, d%neigsetid, hdferr)
     122           0 :          CALL h5sclose_f(spaceid, hdferr)
     123             :          !     ew
     124           0 :          dims(:3) = (/neig, nkpts, jspins/)
     125           0 :          CALL h5screate_simple_f(3, dims(:3), spaceid, hdferr)
     126           0 :          CALL h5dcreate_f(d%fid, "energy", H5T_NATIVE_DOUBLE, spaceid, d%energysetid, hdferr)
     127           0 :          CALL h5sclose_f(spaceid, hdferr)
     128             :          !     ev
     129           0 :          if (l_real .and. .not. l_soc) THEN
     130           0 :             dims(:5) = (/one, nmat, neig, nkpts, jspins/)
     131             :          else
     132           0 :             dims(:5) = (/two, nmat, neig, nkpts, jspins/)
     133             :          endif
     134           0 :          CALL h5screate_simple_f(5, dims(:5), spaceid, hdferr)
     135           0 :          CALL h5dcreate_f(d%fid, "ev", H5T_NATIVE_DOUBLE, spaceid, d%evsetid, hdferr)
     136           0 :          CALL h5sclose_f(spaceid, hdferr)
     137             :       ELSE
     138           0 :          CALL h5fopen_f(TRIM(d%fname)//'.hdf', access_Mode, d%fid, hdferr, access_prp)
     139             :          !get dataset-ids
     140           0 :          CALL h5dopen_f(d%fid, 'energy', d%energysetid, hdferr)
     141           0 :          CALL h5dopen_f(d%fid, 'neig', d%neigsetid, hdferr)
     142           0 :          CALL h5dopen_f(d%fid, 'ev', d%evsetid, hdferr)
     143             :       endif
     144           0 :       IF (.NOT. access_prp == H5P_DEFAULT_f) CALL H5Pclose_f(access_prp&
     145           0 :               &     , hdferr)
     146             : #else
     147             :       CALL juDFT_error("Could not use HDF5 for IO, please recompile")
     148             : #endif
     149           0 :    END SUBROUTINE open_eig
     150             :    !----------------------------------------------------------------------
     151           0 :    SUBROUTINE close_eig(id, filename)
     152             :       !*****************************************************************
     153             :       !     closes hdf-file for eigenvectors+values
     154             :       !*****************************************************************
     155             :       IMPLICIT NONE
     156             :       INTEGER, INTENT(IN)                   :: id
     157             :       CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: filename
     158             : 
     159             :       INTEGER::hdferr
     160             :       TYPE(t_data_HDF), POINTER::d
     161             : 
     162             :       !close datasets
     163             : #ifdef CPP_HDF
     164           0 :       CALL priv_find_data(id, d)
     165             : 
     166           0 :       CALL h5dclose_f(d%energysetid, hdferr)
     167           0 :       CALL h5dclose_f(d%wikssetid, hdferr)
     168           0 :       CALL h5dclose_f(d%neigsetid, hdferr)
     169           0 :       CALL h5dclose_f(d%evsetid, hdferr)
     170             :       !close file
     171           0 :       CALL h5fclose_f(d%fid, hdferr)
     172             :       !If a filename was given and the name is not the current filename
     173           0 :       IF (PRESENT(filename)) THEN
     174           0 :          IF (filename .NE. d%fname) THEN
     175           0 :             CALL system("mv "//TRIM(d%fname)//".hdf "//TRIM(filename)//".hdf")
     176             :          ENDIF
     177             :       ENDIF
     178           0 :       d%fname = "eig"
     179           0 :       CALL eig66_remove_data(id)
     180             : 
     181             : #endif
     182           0 :    END SUBROUTINE close_eig
     183             : #ifdef CPP_HDF
     184             :    !----------------------------------------------------------------------
     185           0 :    SUBROUTINE priv_r_vec(d, nk, jspin, list, z)
     186             : 
     187             :       USE m_hdf_tools
     188             :       IMPLICIT NONE
     189             :       TYPE(t_data_HDF), INTENT(IN)::d
     190             :       INTEGER, INTENT(IN)  :: nk, jspin
     191             :       INTEGER, OPTIONAL, INTENT(IN)  :: list(:)
     192             :       REAL, INTENT(OUT) :: z(:, :)
     193             : 
     194             :       INTEGER :: nmat
     195             :       INTEGER i
     196             : 
     197           0 :       nmat = SIZE(z, 1)
     198             :       !read eigenvectors
     199           0 :       IF (.NOT. PRESENT(list)) THEN
     200             :          ! read all eigenvectors
     201             :          CALL io_read_real2(d%evsetid, (/1, 1, 1, nk, jspin/), &
     202           0 :                             (/1, nmat, SIZE(z, 2), 1, 1/), "evec", z(:nmat, :))
     203             :       ELSE
     204           0 :          DO i = 1, SIZE(list)
     205             :             CALL io_read_real1(d%evsetid, (/1, 1, list(i), nk, jspin/),&
     206           0 :                  &                      (/1, nmat, 1, 1, 1/), "evec", z(:nmat, i))
     207             :          ENDDO
     208             :       END IF
     209             : 
     210           0 :    END SUBROUTINE priv_r_vec
     211             : 
     212             : #endif
     213             : 
     214           0 :    SUBROUTINE write_eig(id, nk, jspin, neig, neig_total, eig, n_size, n_rank, zmat, smat)
     215             : 
     216             :       !*****************************************************************
     217             :       !     writes all eignevecs for the nk-th kpoint
     218             :       !*****************************************************************
     219             :       IMPLICIT NONE
     220             : 
     221             :       INTEGER, INTENT(IN)          :: id, nk, jspin
     222             :       INTEGER, INTENT(IN), OPTIONAL :: n_size, n_rank
     223             :       INTEGER, INTENT(IN), OPTIONAL :: neig, neig_total
     224             :       REAL, INTENT(IN), OPTIONAL :: eig(:)
     225             :       TYPE(t_mat), INTENT(IN), OPTIONAL :: zmat, smat
     226             : 
     227             :       INTEGER i, j, k, nv_local, n1, n2, ne
     228             :       TYPE(t_data_HDF), POINTER::d
     229           0 :       call timestart("write_eig: HDF")
     230           0 :       CALL priv_find_data(id, d)
     231             : 
     232           0 :       if(present(smat)) call juDFT_error("writing smat in HDF not supported yet")
     233             : #ifdef CPP_HDF
     234             :       !
     235             :       !write enparas
     236             :       !
     237           0 :       nv_local = HUGE(1)
     238             : 
     239             :       !
     240             :       !write eigenvalues
     241             :       !
     242             : 
     243           0 :       call timestart("write neig_total")
     244           0 :       IF (PRESENT(neig_total)) THEN
     245           0 :          CALL io_write_integer0(d%neigsetid, (/nk, jspin/), (/1, 1/), "neig_total", neig_total)
     246             :       ENDIF
     247           0 :       call timestop("write neig_total")
     248             : 
     249           0 :       call timestart("write eig")
     250             :       IF (PRESENT(n_rank) .AND. PRESENT(n_size) .AND.&
     251           0 :            &        PRESENT(eig) .AND. PRESENT(neig)) THEN
     252             :          CALL io_write_real1s(&
     253             :               &                     d%energysetid, (/n_rank + 1, nk, jspin/),        &
     254           0 :               &                     (/neig, 1, 1/), eig(:neig), (/n_size, 1, 1/))
     255             :          !write eigenvectors
     256             :          !
     257           0 :       ELSEIF (PRESENT(eig) .AND. PRESENT(neig)) THEN
     258             :          CALL io_write_real1s(&
     259             :               &                     d%energysetid, (/1, nk, jspin/),&
     260           0 :               &                     (/neig, 1, 1/), eig(:neig), (/1, 1, 1/))
     261             :       ELSE
     262           0 :          IF (PRESENT(eig)) CALL juDFT_error("BUG in calling write_eig")
     263             :       ENDIF
     264           0 :       call timestop("write eig")
     265             : 
     266           0 :       call timestart("write zmat")
     267           0 :       IF (PRESENT(zmat) .AND. .NOT. PRESENT(neig))&
     268           0 :            &    CALL juDFT_error("BUG in calling write_eig with eigenvector")
     269             : 
     270           0 :       n1 = 1; n2 = 0
     271           0 :       IF (PRESENT(n_size)) n1 = n_size
     272           0 :       IF (PRESENT(n_rank)) n2 = n_rank
     273           0 :       IF (PRESENT(zmat)) THEN
     274           0 :          IF (zmat%l_real) THEN
     275             :             CALL io_write_real2s(&
     276             :                  &                     d%evsetid, (/1, 1, n2 + 1, nk, jspin/),&
     277           0 :                  &           (/1, SIZE(zmat%data_r, 1), neig, 1, 1/), REAL(zmat%data_r(:, :neig)), (/1, 1, n1, 1, 1/))
     278             :          ELSE
     279             :             CALL io_write_real2s(&
     280             :                  &                     d%evsetid, (/1, 1, n2 + 1, nk, jspin/),&
     281           0 :                  &           (/1, SIZE(zmat%data_c, 1), neig, 1, 1/), REAL(zmat%data_c(:, :neig)), (/1, 1, n1, 1, 1/))
     282             :             CALL io_write_real2s(&
     283             :                  &                     d%evsetid, (/2, 1, n2 + 1, nk, jspin/),&
     284             :                  &           (/1, SIZE(zmat%data_c, 1), neig, 1, 1/), AIMAG(zmat%data_c(:, :neig)),&
     285           0 :                  &           (/1, 1, n1, 1, 1/))
     286             :          ENDIF
     287             :       ENDIF
     288           0 :       call timestop("write zmat")
     289             : #endif
     290           0 :       call timestop("write_eig: HDF")
     291           0 :    END SUBROUTINE write_eig
     292             : 
     293             : #ifdef CPP_HDF
     294             : 
     295             :    !----------------------------------------------------------------------
     296           0 :    SUBROUTINE priv_r_vecc(&
     297           0 :         &                     d, nk, jspin, list, z)
     298             : 
     299             :       USE m_hdf_tools
     300             :       IMPLICIT NONE
     301             :       TYPE(t_data_HDF), INTENT(IN)::d
     302             :       INTEGER, INTENT(IN)  :: nk, jspin
     303             :       INTEGER, OPTIONAL, INTENT(IN)  :: list(:)
     304             :       COMPLEX, INTENT(OUT) :: z(:, :)
     305             : 
     306           0 :       REAL, ALLOCATABLE :: z1(:, :, :)
     307             :       INTEGER i, j
     308             :       INTEGER :: nmat
     309             : 
     310           0 :       nmat = SIZE(z, 1)
     311             : 
     312           0 :       IF (.NOT. PRESENT(list)) THEN
     313             :          ! read all eigenvectors
     314           0 :          ALLOCATE (z1(2, nmat, SIZE(z, 2)))
     315             :          CALL io_read_real3(d%evsetid, (/1, 1, 1, nk, jspin/),&
     316           0 :               &                      (/2, nmat, SIZE(z, 2), 1, 1/), "z1", z1)
     317           0 :          DO i = 1, SIZE(z, 2)
     318           0 :             z(:, i) = CMPLX(z1(1, :, i), z1(2, :, i))
     319             :          ENDDO
     320             :       ELSE
     321           0 :          ALLOCATE (z1(2, nmat, 1))
     322           0 :          DO i = 1, SIZE(list)
     323             :             CALL io_read_real3(d%evsetid, (/1, 1, list(i), nk, jspin/),&
     324           0 :              &                      (/2, nmat, 1, 1, 1/), "z1", z1)
     325           0 :             z(:, i) = CMPLX(z1(1, :, 1), z1(2, :, 1))
     326             :          ENDDO
     327             :       END IF
     328           0 :    END SUBROUTINE priv_r_vecc
     329             :    !-----------------------------------------------------------------------
     330             : 
     331             : #endif
     332             : 
     333           0 :    SUBROUTINE read_eig(id, nk, jspin, neig, eig, list, zMat, smat)
     334             :       IMPLICIT NONE
     335             :       INTEGER, INTENT(IN)            :: id, nk, jspin
     336             :       INTEGER, INTENT(OUT), OPTIONAL  :: neig
     337             :       REAL, INTENT(OUT), OPTIONAL  :: eig(:)
     338             :       INTEGER, INTENT(IN), OPTIONAL   :: list(:)
     339             :       TYPE(t_mat), OPTIONAL  :: zmat, smat
     340             : 
     341             : #ifdef CPP_HDF
     342             :       INTEGER:: n1, n, k
     343             :       TYPE(t_data_HDF), POINTER::d
     344           0 :       if(present(smat)) call juDFT_error("reading smat not supported for HDF")
     345           0 :       CALL priv_find_data(id, d)
     346             : 
     347           0 :       IF (PRESENT(neig)) THEN
     348           0 :          CALL io_read_integer0(d%neigsetid, (/nk, jspin/), (/1, 1/), neig)
     349             : 
     350           0 :          IF (PRESENT(eig)) THEN                           ! read eigenv
     351           0 :             IF (neig > SIZE(eig)) THEN
     352           0 :                WRITE (*, *) neig, SIZE(eig)
     353           0 :                CALL juDFT_error("eig66_hdf$readeig", calledby="eig66_hdf")
     354             :             ENDIF
     355             :             CALL io_read_real1(d%energysetid, (/1, nk, jspin/), (/neig, 1, 1/),&
     356           0 :                  &                      "eig", eig(:neig))
     357             :          ENDIF
     358             :       ENDIF
     359             : 
     360           0 :       IF (PRESENT(zMat)) THEN
     361           0 :          IF (zmat%l_real) THEN
     362           0 :             CALL priv_r_vec(d, nk, jspin, list, zmat%data_r)
     363             :          ELSE
     364           0 :             CALL priv_r_vecc(d, nk, jspin, list, zmat%data_c)
     365             :          ENDIF
     366             :       ENDIF
     367             : #endif
     368           0 :    END SUBROUTINE read_eig
     369             : 
     370           0 : END MODULE