pot_io.F90 9.78 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 69 70 71 72
#ifdef CPP_HDF
      INTEGER(HID_T)    :: fileID
#endif
      INTEGER           :: currentStarsIndex,currentLatharmsIndex
      INTEGER           :: currentStructureIndex
      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
            CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex)
80

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
            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
            CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex)

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

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

154
         CLOSE(iUnit)
155 156 157 158
      END IF

   END SUBROUTINE readPotential

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

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

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

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

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

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

194 195 196
      CALL getMode(mode)

      IF(mode.EQ.POT_HDF5_MODE) THEN
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
#ifdef CPP_HDF
         CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex)

         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

         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

230 231 232 233 234 235 236 237 238 239
         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
240 241
         CALL writePotentialHDF(input, fileID, archiveName, potentialType,&
                                currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
242
                                iter,fr,fpw,fzTemp,fzxyTemp)
243 244 245 246 247 248 249 250

         IF(l_storeIndices) THEN
            CALL writePOTHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
                                    currentStructureIndex)
         END IF

         CALL closeCDNPOT_HDF(fileID)
#endif
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 279
      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
280 281 282 283 284 285 286 287 288
      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
289 290 291
   END SUBROUTINE getMode

END MODULE m_pot_io