pot_io.F90 10.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! This module is a wrapper for the potential I/O
!!!
!!!                             GM'17
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

MODULE m_pot_io

   USE m_types
   USE m_juDFT
   USE m_loddop
   USE m_wrtdop
21 22 23 24
   USE m_cdnpot_io_hdf
#ifdef CPP_HDF
   USE hdf5
#endif
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
   IMPLICIT NONE

   PRIVATE
   PUBLIC readPotential, writePotential
   PUBLIC POT_ARCHIVE_TYPE_TOT_const, POT_ARCHIVE_TYPE_COUL_const
   PUBLIC POT_ARCHIVE_TYPE_X_const

   INTEGER,          PARAMETER :: POT_ARCHIVE_TYPE_TOT_const = 1
   INTEGER,          PARAMETER :: POT_ARCHIVE_TYPE_COUL_const = 2
   INTEGER,          PARAMETER :: POT_ARCHIVE_TYPE_X_const = 3

   INTEGER,          PARAMETER :: POT_DIRECT_MODE = 1
   INTEGER,          PARAMETER :: POT_STREAM_MODE = 2
   INTEGER,          PARAMETER :: POT_HDF5_MODE   = 3

   CONTAINS

   SUBROUTINE readPotential(stars,vacuum,atoms,sphhar,input,sym,archiveType,&
                            iter,fr,fpw,fz,fzxy)

      TYPE(t_stars),INTENT(IN)  :: stars
      TYPE(t_vacuum),INTENT(IN) :: vacuum
      TYPE(t_atoms),INTENT(IN)  :: atoms
      TYPE(t_sphhar),INTENT(IN) :: sphhar
      TYPE(t_input),INTENT(IN)  :: input
      TYPE(t_sym),INTENT(IN)    :: sym

      INTEGER, INTENT (OUT)     :: iter
      INTEGER, INTENT (IN)      :: archiveType

      !     ..
      !     .. Array Arguments ..
57
      COMPLEX, INTENT (OUT) :: fpw(stars%ng3,input%jspins), fzxy(vacuum%nmzxyd,stars%ng2-1,2,input%jspins)
58 59 60
      REAL,    INTENT (OUT) :: fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins), fz(vacuum%nmzd,2,input%jspins)

      ! local variables
61
      INTEGER           :: mode, iUnit
62 63 64
      LOGICAL           :: l_exist
      CHARACTER(len=30) :: filename

65 66 67 68
#ifdef CPP_HDF
      INTEGER(HID_T)    :: fileID
#endif
      INTEGER           :: currentStarsIndex,currentLatharmsIndex
69
      INTEGER           :: currentStructureIndex,currentStepfunctionIndex
70 71 72
      INTEGER           :: potentialType
      CHARACTER(LEN=30) :: archiveName

73 74 75
      CALL getMode(mode)

      IF(mode.EQ.POT_HDF5_MODE) THEN
76
#ifdef CPP_HDF
77 78
         INQUIRE(FILE='pot.hdf',EXIST=l_exist)
         IF (l_exist) THEN
79 80
            CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,&
                             currentStructureIndex,currentStepfunctionIndex)
81

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
            archiveName = 'illegalPotentialArchive'
            IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN
               archiveName = 'pottot'
            END IF
            IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN
               archiveName = 'potcoul'
            END IF
            IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN
               archiveName = 'potx'
            END IF

            potentialType = POTENTIAL_TYPE_IN_const

            l_exist = isPotentialEntryPresentHDF(fileID,archiveName,potentialType)

            CALL closeCDNPOT_HDF(fileID)
         END IF

         IF(l_exist) THEN
101 102
            CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,&
                             currentStructureIndex,currentStepfunctionIndex)
103 104 105 106 107

            CALL readPotentialHDF(fileID, archiveName, potentialType,&
                                  iter,fr,fpw,fz,fzxy)

            CALL closeCDNPOT_HDF(fileID)
108
         ELSE
109
            WRITE(*,*) 'Potential entry or pot.hdf file not found.'
110 111 112
            WRITE(*,*) 'Falling back to stream access file pot.str.'
            mode = POT_STREAM_MODE
         END IF
113 114 115 116 117
#else
         WRITE(*,*) 'Not compiled for pot.hdf file usage.'
         WRITE(*,*) 'Falling back to stream access file pot.str.'
         mode = POT_STREAM_MODE
#endif
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 144 145 146 147 148 149 150 151 152 153 154 155
      END IF

      IF(mode.EQ.POT_STREAM_MODE) THEN
         INQUIRE(FILE='pot.str',EXIST=l_exist)
         IF (l_exist) THEN
            !load density from cdn.str and exit subroutine

            RETURN
         ELSE
            WRITE(*,*) 'pot.str file not found.'
            WRITE(*,*) 'Falling back to direct access file.'
            mode = POT_DIRECT_MODE
         END IF
      END IF

      IF (mode.EQ.POT_DIRECT_MODE) THEN
         filename = 'illegalPotentialArchive'
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN
            filename = 'pottot'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN
            filename = 'potcoul'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN
            filename = 'potx'
         END IF

         INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
         IF(.NOT.l_exist) THEN
            CALL juDFT_error("potential file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readPotential")
         END IF

         iUnit = 11
         OPEN (iUnit,file=TRIM(ADJUSTL(filename)),form='unformatted',status='unknown')

         CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
                     iUnit,iter,fr,fpw,fz,fzxy)

156
         CLOSE(iUnit)
157 158 159 160
      END IF

   END SUBROUTINE readPotential

161
   SUBROUTINE writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
162 163 164 165 166
                             iter,fr,fpw,fz,fzxy)

      TYPE(t_stars),INTENT(IN)  :: stars
      TYPE(t_vacuum),INTENT(IN) :: vacuum
      TYPE(t_atoms),INTENT(IN)  :: atoms
167
      TYPE(t_cell), INTENT(IN)  :: cell
168 169 170
      TYPE(t_sphhar),INTENT(IN) :: sphhar
      TYPE(t_input),INTENT(IN)  :: input
      TYPE(t_sym),INTENT(IN)    :: sym
171
      TYPE(t_oneD),INTENT(IN)   :: oneD
172 173 174 175 176

      INTEGER, INTENT (IN)      :: iter
      INTEGER, INTENT (IN)      :: archiveType
      !     ..
      !     .. Array Arguments ..
177
      COMPLEX, INTENT (IN) :: fpw(stars%ng3,input%jspins), fzxy(vacuum%nmzxyd,stars%ng2-1,2,input%jspins)
178 179 180 181
      REAL,    INTENT (IN) :: fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins), fz(vacuum%nmzd,2,input%jspins)

      ! local variables
      INTEGER           :: mode, iUnit
182
      LOGICAL           :: l_exist, l_storeIndices
183 184
      CHARACTER(len=30) :: filename

185 186 187 188
#ifdef CPP_HDF
      INTEGER(HID_T)    :: fileID
#endif
      INTEGER           :: currentStarsIndex,currentLatharmsIndex
189
      INTEGER           :: currentStructureIndex,currentStepfunctionIndex
190 191 192
      INTEGER           :: potentialType
      CHARACTER(LEN=30) :: archiveName

193 194 195
      REAL              :: fzTemp(vacuum%nmzd,2,input%jspins)
      COMPLEX           :: fzxyTemp(vacuum%nmzxyd,stars%ng2-1,2,input%jspins)

196 197 198
      CALL getMode(mode)

      IF(mode.EQ.POT_HDF5_MODE) THEN
199
#ifdef CPP_HDF
200 201
         CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,&
                          currentStructureIndex,currentStepfunctionIndex)
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

         l_storeIndices = .FALSE.
         IF (currentStarsIndex.EQ.0) THEN
            currentStarsIndex = 1
            l_storeIndices = .TRUE.
            CALL writeStarsHDF(fileID, currentStarsIndex, stars)
         END IF
         IF (currentLatharmsIndex.EQ.0) THEN
            currentLatharmsIndex = 1
            l_storeIndices = .TRUE.
            CALL writeLatharmsHDF(fileID, currentLatharmsIndex, sphhar)
         END IF
         IF(currentStructureIndex.EQ.0) THEN
            currentStructureIndex = 1
            l_storeIndices = .TRUE.
            CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, currentStructureIndex)
         END IF
219 220 221 222 223
         IF(currentStepfunctionIndex.EQ.0) THEN
            currentStepfunctionIndex = 1
            l_storeIndices = .TRUE.
            CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, stars)
         END IF
224 225 226 227 228 229 230 231 232 233 234 235 236 237

         archiveName = 'illegalPotentialArchive'
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN
            archiveName = 'pottot'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN
            archiveName = 'potcoul'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN
            archiveName = 'potx'
         END IF

         potentialType = POTENTIAL_TYPE_IN_const

238 239 240 241 242 243 244 245 246 247
         fzTemp(:,:,:) = fz(:,:,:)
         fzxyTemp(:,:,:,:) = fzxy(:,:,:,:)
         IF(vacuum%nvac.EQ.1) THEN
            fzTemp(:,2,:)=fzTemp(:,1,:)
            IF (sym%invs) THEN
               fzxyTemp(:,:,2,:) = CONJG(fzxyTemp(:,:,1,:))
            ELSE
               fzxyTemp(:,:,2,:) = fzxyTemp(:,:,1,:)
            END IF
         END IF
248 249
         CALL writePotentialHDF(input, fileID, archiveName, potentialType,&
                                currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
250
                                currentStepfunctionIndex,iter,fr,fpw,fzTemp,fzxyTemp)
251 252 253

         IF(l_storeIndices) THEN
            CALL writePOTHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
254
                                    currentStructureIndex,currentStepfunctionIndex)
255 256 257 258
         END IF

         CALL closeCDNPOT_HDF(fileID)
#endif
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
      ELSE IF(mode.EQ.POT_STREAM_MODE) THEN
         ! Write potential to pot.str file
         STOP 'POT_STREAM_MODE not yet implemented!'
      ELSE
         ! Direct mode
         filename = 'illegalPotentialArchive'
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN
            filename = 'pottot'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN
            filename = 'potcoul'
         END IF
         IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN
            filename = 'potx'
         END IF

         iUnit = 11
         OPEN (iUnit,file=TRIM(ADJUSTL(filename)),form='unformatted',status='unknown')
         CALL wrtdop(stars,vacuum,atoms,sphhar,input,sym,&
                     iUnit,iter,fr,fpw,fz,fzxy)
         CLOSE(iUnit)
      END IF

   END SUBROUTINE writePotential

   SUBROUTINE getMode(mode)
      INTEGER, INTENT(OUT) :: mode

      mode = POT_DIRECT_MODE
288 289 290 291 292 293 294 295 296
      IF (juDFT_was_argument("-stream_cdn")) mode=POT_STREAM_MODE
      IF (juDFT_was_argument("-hdf_cdn")) THEN
#ifdef CPP_HDF
         mode=POT_HDF5_MODE
#else
         WRITE(*,*) 'Code not compiled with HDF5 support.'
         WRITE(*,*) 'Falling back to direct access.'
#endif
      END IF
297 298 299
   END SUBROUTINE getMode

END MODULE m_pot_io