eig66_hdf.F90 12.8 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8
MODULE m_eig66_hdf
#include "juDFT_env.h"
Matthias Redies's avatar
Matthias Redies committed
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
   !*****************************************************************
   ! DESC:Module for hdf-io of eig-file
   !      To be compatible with f90 interface of HDF, use kind for vars
   !
   !      !ATTENTION before calling openeig and after calling closeeig!
   !      !the hdf library has to be initialized or finalized, respectively
   !
   !      CONTAINS the following subroutines:
   !      openeig        opens file
   !      closeeig       closes file
   !      read_keb       reads kpt, enpara and basis data
   !      read_neig      read no of eigenvalues (and eigenvalues itself)
   !      read_eig       reads eigenvectors
   !      writeeig       saves all data for kpt
   !      writesingleeig saves data for one kpt and energy
   !
   !
   !                          Daniel Wortmann
   !*****************************************************************
   USE m_eig66_data
   USE m_types
30
#ifdef CPP_HDF
Matthias Redies's avatar
Matthias Redies committed
31 32 33
   USE hdf5
   USE m_hdf_tools
   IMPLICIT NONE
34

Matthias Redies's avatar
Matthias Redies committed
35 36 37 38
   PRIVATE
   INTEGER, PARAMETER :: one = 1, two = 2, three = 3, zero = 0
   !to have the correct
   !type for array constructors
39 40

#endif
Matthias Redies's avatar
Matthias Redies committed
41 42 43
   PUBLIC open_eig, close_eig
   PUBLIC read_eig
   PUBLIC write_eig!,writesingleeig,writeeigc,writebas
44 45

CONTAINS
Matthias Redies's avatar
Matthias Redies committed
46 47 48 49 50 51 52 53 54 55 56 57 58 59
   SUBROUTINE priv_find_data(id, d)
      INTEGER, INTENT(IN)::id
      TYPE(t_data_hdf), POINTER:: d

      CLASS(t_data), POINTER   ::dp
      CALL eig66_find_data(dp, id)
      SELECT TYPE (dp)
      TYPE is (t_data_hdf)
         d => dp
      CLASS default
         CALL judft_error("BUG: wrong datatype in eig66_hdf")
      END SELECT
   END SUBROUTINE priv_find_data
   !----------------------------------------------------------------------
Matthias Redies's avatar
Matthias Redies committed
60
   SUBROUTINE open_eig(id, mpi_comm, nmat, neig, nkpts, jspins, create, l_real, l_soc, readonly, l_olap, filename)
Matthias Redies's avatar
Matthias Redies committed
61 62 63 64 65 66 67 68

      !*****************************************************************
      !     opens hdf-file for eigenvectors+values
      !*****************************************************************
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: id, mpi_comm
      INTEGER, INTENT(IN) :: nmat, neig, nkpts, jspins
Matthias Redies's avatar
Matthias Redies committed
69
      LOGICAL, INTENT(IN) :: create, readonly, l_real, l_soc, l_olap
Matthias Redies's avatar
Matthias Redies committed
70
      CHARACTER(LEN=*), OPTIONAL :: filename
71 72

#ifdef CPP_HDF
73

Matthias Redies's avatar
Matthias Redies committed
74 75 76 77 78 79
      INTEGER         :: hdferr, access_mode
      INTEGER(HID_T)  :: creation_prp, access_prp, spaceid
      LOGICAL         :: l_exist
      INTEGER(HSIZE_T):: dims(7)
      TYPE(t_data_HDF), POINTER::d
      !Set creation and access properties
Matthias Redies's avatar
Matthias Redies committed
80

81
#ifdef CPP_HDFMPI
Matthias Redies's avatar
Matthias Redies committed
82 83 84 85 86 87 88 89 90 91 92
      INCLUDE 'mpif.h'
      IF (readonly) THEN
         access_prp = H5P_DEFAULT_f
         creation_prp = H5P_DEFAULT_f
      ELSE
         CALL h5pcreate_f(H5P_FILE_ACCESS_F, access_prp, hdferr)
         !      CALL h5pset_fapl_mpiposix_f(access_prp,MPI_COMM,
         !     +.false.,hdferr)
         CALL h5pset_fapl_mpio_f(access_prp, MPI_COMM, MPI_INFO_NULL, hdferr)
         creation_prp = H5P_DEFAULT_f !no special creation property
      ENDIF
93
#else
Matthias Redies's avatar
Matthias Redies committed
94 95 96
      access_prp = H5P_DEFAULT_f
      creation_prp = H5P_DEFAULT_f
#endif
Matthias Redies's avatar
Matthias Redies committed
97 98 99
      
      if(l_olap) call juDFT_error("olap not implemented for hdf5")
      
Matthias Redies's avatar
Matthias Redies committed
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
      CALL priv_find_data(id, d)
      IF (PRESENT(filename)) d%fname = filename
      CALL eig66_data_storedefault(d, jspins, nkpts, nmat, neig, l_real, l_soc)
      !set access_flags according
      IF (readonly) THEN
         access_mode = H5F_ACC_RDONLY_F
      ELSE
         access_mode = H5F_ACC_RDWR_F
      ENDIF
      !     OPEN FILE and get D%FID's
      IF (create) THEN
         INQUIRE (FILE=TRIM(d%fname)//'.hdf', EXIST=l_exist)
         access_mode = H5F_ACC_TRUNC_F
         !         IF (l_exist) WRITE (*,*)'Warning: eig.hdf was overwritten'
         CALL h5fcreate_f(TRIM(d%fname)//'.hdf', access_Mode, d%fid, hdferr, creation_prp, access_prp)
         ! create dataspaces and datasets
         !   scalars
         dims(:2) = (/nkpts, jspins/)
         CALL h5screate_simple_f(2, dims(:2), spaceid, hdferr)
         CALL h5dcreate_f(d%fid, "neig", H5T_NATIVE_INTEGER, spaceid, d%neigsetid, hdferr)
         CALL h5sclose_f(spaceid, hdferr)
         !     ew
         dims(:3) = (/neig, nkpts, jspins/)
         CALL h5screate_simple_f(3, dims(:3), spaceid, hdferr)
         CALL h5dcreate_f(d%fid, "energy", H5T_NATIVE_DOUBLE, spaceid, d%energysetid, hdferr)
         CALL h5sclose_f(spaceid, hdferr)
         !     ev
         if (l_real .and. .not. l_soc) THEN
            dims(:5) = (/one, nmat, neig, nkpts, jspins/)
         else
            dims(:5) = (/two, nmat, neig, nkpts, jspins/)
         endif
         CALL h5screate_simple_f(5, dims(:5), spaceid, hdferr)
         CALL h5dcreate_f(d%fid, "ev", H5T_NATIVE_DOUBLE, spaceid, d%evsetid, hdferr)
         CALL h5sclose_f(spaceid, hdferr)
      ELSE
         CALL h5fopen_f(TRIM(d%fname)//'.hdf', access_Mode, d%fid, hdferr, access_prp)
         !get dataset-ids
         CALL h5dopen_f(d%fid, 'energy', d%energysetid, hdferr)
         CALL h5dopen_f(d%fid, 'neig', d%neigsetid, hdferr)
         CALL h5dopen_f(d%fid, 'ev', d%evsetid, hdferr)
      endif
      IF (.NOT. access_prp == H5P_DEFAULT_f) CALL H5Pclose_f(access_prp&
              &     , hdferr)
144
#else
Matthias Redies's avatar
Matthias Redies committed
145
      CALL juDFT_error("Could not use HDF5 for IO, please recompile")
146
#endif
Matthias Redies's avatar
Matthias Redies committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160
   END SUBROUTINE open_eig
   !----------------------------------------------------------------------
   SUBROUTINE close_eig(id, filename)
      !*****************************************************************
      !     closes hdf-file for eigenvectors+values
      !*****************************************************************
      IMPLICIT NONE
      INTEGER, INTENT(IN)                   :: id
      CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: filename

      INTEGER::hdferr
      TYPE(t_data_HDF), POINTER::d

      !close datasets
161
#ifdef CPP_HDF
Matthias Redies's avatar
Matthias Redies committed
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
      CALL priv_find_data(id, d)

      CALL h5dclose_f(d%energysetid, hdferr)
      CALL h5dclose_f(d%wikssetid, hdferr)
      CALL h5dclose_f(d%neigsetid, hdferr)
      CALL h5dclose_f(d%evsetid, hdferr)
      !close file
      CALL h5fclose_f(d%fid, hdferr)
      !If a filename was given and the name is not the current filename
      IF (PRESENT(filename)) THEN
         IF (filename .NE. d%fname) THEN
            CALL system("mv "//TRIM(d%fname)//".hdf "//TRIM(filename)//".hdf")
         ENDIF
      ENDIF
      d%fname = "eig"
      CALL eig66_remove_data(id)
178 179

#endif
Matthias Redies's avatar
Matthias Redies committed
180
   END SUBROUTINE close_eig
181
#ifdef CPP_HDF
Matthias Redies's avatar
Matthias Redies committed
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
   !----------------------------------------------------------------------
   SUBROUTINE priv_r_vec(d, nk, jspin, list, z)

      USE m_hdf_tools
      IMPLICIT NONE
      TYPE(t_data_HDF), INTENT(IN)::d
      INTEGER, INTENT(IN)  :: nk, jspin
      INTEGER, OPTIONAL, INTENT(IN)  :: list(:)
      REAL, INTENT(OUT) :: z(:, :)

      INTEGER :: nmat
      INTEGER i

      nmat = SIZE(z, 1)
      !read eigenvectors
      IF (.NOT. PRESENT(list)) THEN
         ! read all eigenvectors
         CALL io_read_real2(d%evsetid, (/1, 1, 1, nk, jspin/), &
                            (/1, nmat, SIZE(z, 2), 1, 1/), z(:nmat, :))
      ELSE
         DO i = 1, SIZE(list)
            CALL io_read_real1(d%evsetid, (/1, 1, list(i), nk, jspin/),&
                 &                      (/1, nmat, 1, 1, 1/), z(:nmat, i))
         ENDDO
      END IF

   END SUBROUTINE priv_r_vec
209 210 211

#endif

Gregor Michalicek's avatar
Gregor Michalicek committed
212
   SUBROUTINE write_eig(id, nk, jspin, neig, neig_total, eig, n_size, n_rank, zmat, smat)
213

Matthias Redies's avatar
Matthias Redies committed
214 215 216 217
      !*****************************************************************
      !     writes all eignevecs for the nk-th kpoint
      !*****************************************************************
      IMPLICIT NONE
218

Matthias Redies's avatar
Matthias Redies committed
219 220 221
      INTEGER, INTENT(IN)          :: id, nk, jspin
      INTEGER, INTENT(IN), OPTIONAL :: n_size, n_rank
      INTEGER, INTENT(IN), OPTIONAL :: neig, neig_total
Gregor Michalicek's avatar
Gregor Michalicek committed
222
      REAL, INTENT(IN), OPTIONAL :: eig(:)
Matthias Redies's avatar
Matthias Redies committed
223
      TYPE(t_mat), INTENT(IN), OPTIONAL :: zmat, smat
224

Matthias Redies's avatar
Matthias Redies committed
225 226 227
      INTEGER i, j, k, nv_local, n1, n2, ne
      TYPE(t_data_HDF), POINTER::d
      CALL priv_find_data(id, d)
228

Matthias Redies's avatar
Matthias Redies committed
229
      if(present(smat)) call juDFT_error("writing smat in HDF not supported yet")
230
#ifdef CPP_HDF
Matthias Redies's avatar
Matthias Redies committed
231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
      !
      !write enparas
      !
      nv_local = HUGE(1)

      !
      !write eigenvalues
      !

      IF (PRESENT(neig_total)) THEN
         CALL io_write_integer0(d%neigsetid, (/nk, jspin/), (/1, 1/), neig_total)
      ENDIF

      IF (PRESENT(n_rank) .AND. PRESENT(n_size) .AND.&
           &        PRESENT(eig) .AND. PRESENT(neig)) THEN
         CALL io_write_real1s(&
              &                     d%energysetid, (/n_rank + 1, nk, jspin/),        &
              &                     (/neig, 1, 1/), eig(:neig), (/n_size, 1, 1/))
         !write eigenvectors
         !
      ELSEIF (PRESENT(eig) .AND. PRESENT(neig)) THEN
         CALL io_write_real1s(&
              &                     d%energysetid, (/1, nk, jspin/),&
              &                     (/neig, 1, 1/), eig(:neig), (/1, 1, 1/))
      ELSE
         IF (PRESENT(eig)) CALL juDFT_error("BUG in calling write_eig")
      ENDIF
      IF (PRESENT(zmat) .AND. .NOT. PRESENT(neig))&
           &    CALL juDFT_error("BUG in calling write_eig with eigenvector")

      n1 = 1; n2 = 0
      IF (PRESENT(n_size)) n1 = n_size
      IF (PRESENT(n_rank)) n2 = n_rank
      IF (PRESENT(zmat)) THEN
         IF (zmat%l_real) THEN
            CALL io_write_real2s(&
                 &                     d%evsetid, (/1, 1, n2 + 1, nk, jspin/),&
                 &           (/1, SIZE(zmat%data_r, 1), neig, 1, 1/), REAL(zmat%data_r(:, :neig)), (/1, 1, n1, 1, 1/))
         ELSE
            CALL io_write_real2s(&
                 &                     d%evsetid, (/1, 1, n2 + 1, nk, jspin/),&
                 &           (/1, SIZE(zmat%data_c, 1), neig, 1, 1/), REAL(zmat%data_c(:, :neig)), (/1, 1, n1, 1, 1/))
            CALL io_write_real2s(&
                 &                     d%evsetid, (/2, 1, n2 + 1, nk, jspin/),&
                 &           (/1, SIZE(zmat%data_c, 1), neig, 1, 1/), AIMAG(zmat%data_c(:, :neig)),&
                 &           (/1, 1, n1, 1, 1/))
         ENDIF
      ENDIF
279 280

#endif
Matthias Redies's avatar
Matthias Redies committed
281
   END SUBROUTINE write_eig
282 283 284

#ifdef CPP_HDF

Matthias Redies's avatar
Matthias Redies committed
285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
   !----------------------------------------------------------------------
   SUBROUTINE priv_r_vecc(&
        &                     d, nk, jspin, list, z)

      USE m_hdf_tools
      IMPLICIT NONE
      TYPE(t_data_HDF), INTENT(IN)::d
      INTEGER, INTENT(IN)  :: nk, jspin
      INTEGER, OPTIONAL, INTENT(IN)  :: list(:)
      COMPLEX, INTENT(OUT) :: z(:, :)

      REAL, ALLOCATABLE :: z1(:, :, :)
      INTEGER i, j
      INTEGER :: nmat

      nmat = SIZE(z, 1)

      IF (.NOT. PRESENT(list)) THEN
         ! read all eigenvectors
         ALLOCATE (z1(2, nmat, SIZE(z, 2)))
         CALL io_read_real3(d%evsetid, (/1, 1, 1, nk, jspin/),&
              &                      (/2, nmat, SIZE(z, 2), 1, 1/), z1)
         DO i = 1, SIZE(z, 2)
            z(:, i) = CMPLX(z1(1, :, i), z1(2, :, i))
         ENDDO
      ELSE
         ALLOCATE (z1(2, nmat, 1))
         DO i = 1, SIZE(list)
            CALL io_read_real3(d%evsetid, (/1, 1, list(i), nk, jspin/),&
             &                      (/2, nmat, 1, 1, 1/), z1)
            z(:, i) = CMPLX(z1(1, :, 1), z1(2, :, 1))
         ENDDO
      END IF
   END SUBROUTINE priv_r_vecc
   !-----------------------------------------------------------------------
320 321 322

#endif

Gregor Michalicek's avatar
Gregor Michalicek committed
323
   SUBROUTINE read_eig(id, nk, jspin, neig, eig, list, zMat, smat)
Matthias Redies's avatar
Matthias Redies committed
324 325 326
      IMPLICIT NONE
      INTEGER, INTENT(IN)            :: id, nk, jspin
      INTEGER, INTENT(OUT), OPTIONAL  :: neig
Gregor Michalicek's avatar
Gregor Michalicek committed
327
      REAL, INTENT(OUT), OPTIONAL  :: eig(:)
Matthias Redies's avatar
Matthias Redies committed
328
      INTEGER, INTENT(IN), OPTIONAL   :: list(:)
Matthias Redies's avatar
Matthias Redies committed
329
      TYPE(t_mat), OPTIONAL  :: zmat, smat
330 331

#ifdef CPP_HDF
Matthias Redies's avatar
Matthias Redies committed
332 333
      INTEGER:: n1, n, k
      TYPE(t_data_HDF), POINTER::d
Matthias Redies's avatar
Matthias Redies committed
334
      if(present(smat)) call juDFT_error("reading smat not supported for HDF")
Matthias Redies's avatar
Matthias Redies committed
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
      CALL priv_find_data(id, d)

      IF (PRESENT(neig)) THEN
         CALL io_read_integer0(d%neigsetid, (/nk, jspin/), (/1, 1/), neig)

         IF (PRESENT(eig)) THEN                           ! read eigenv
            IF (neig > SIZE(eig)) THEN
               WRITE (*, *) neig, SIZE(eig)
               CALL juDFT_error("eig66_hdf$readeig", calledby="eig66_hdf")
            ENDIF
            CALL io_read_real1(d%energysetid, (/1, nk, jspin/), (/neig, 1, 1/),&
                 &                      eig(:neig))
         ENDIF
      ENDIF

      IF (PRESENT(zMat)) THEN
         IF (zmat%l_real) THEN
            CALL priv_r_vec(d, nk, jspin, list, zmat%data_r)
         ELSE
            CALL priv_r_vecc(d, nk, jspin, list, zmat%data_c)
         ENDIF
      ENDIF
357
#endif
Matthias Redies's avatar
Matthias Redies committed
358
   END SUBROUTINE read_eig
359

Matthias Redies's avatar
Matthias Redies committed
360
END MODULE
361