Initial commit for the introduction of the pot.hdf file.

So far this is only tested for simple systems.

Like the cdn.hdf file the usage of pot.hdf is also enables with the -hdf_cdn
command line switch.

io/CMakeLists.txt

@@ -16,7 +16,7 @@ io/inpnoco.F90
 io/loddop.f90
 io/cdnpot_io_hdf.F90
 io/cdn_io.F90
-io/pot_io.f90
+io/pot_io.F90
 io/rw_inp.f90
 io/rw_noco.f90
 io/r_inpXML.F90

io/cdn_io.F90

@@ -125,7 +125,7 @@ MODULE m_cdn_io
                   CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="readDensity")
             END SELECT
             l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType)
-            CALL closeCDN_HDF(fileID)
+            CALL closeCDNPOT_HDF(fileID)
          END IF
 
          IF (l_exist) THEN
@@ -135,7 +135,7 @@ MODULE m_cdn_io
             CALL readDensityHDF(fileID, archiveName, densityType,&
                                 fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
-            CALL closeCDN_HDF(fileID)
+            CALL closeCDNPOT_HDF(fileID)
             RETURN
          ELSE
             WRITE(*,*) 'cdn.hdf file or relevant density entry not found.'
@@ -350,11 +350,11 @@ MODULE m_cdn_io
                               fermiEnergy,l_qfix,iter+relCdnIndex,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
          IF(l_storeIndices) THEN
-            CALL writeHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
-                                 currentStructureIndex,readDensityIndex,lastDensityIndex)
+            CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
+                                    currentStructureIndex,readDensityIndex,lastDensityIndex)
          END IF
 
-         CALL closeCDN_HDF(fileID)
+         CALL closeCDNPOT_HDF(fileID)
 #endif
       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
          ! Write density to cdn.str file
@@ -536,7 +536,7 @@ MODULE m_cdn_io
          ELSE
             l_error = .TRUE.
          END IF
-         CALL closeCDN_HDF(fileID)
+         CALL closeCDNPOT_HDF(fileID)
 #endif
       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
          STOP 'cdn.str not yet implemented!'
@@ -576,7 +576,7 @@ MODULE m_cdn_io
             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
                              readDensityIndex,lastDensityIndex)
             CALL readCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints)
-            CALL closeCDN_HDF(fileID)
+            CALL closeCDNPOT_HDF(fileID)
             RETURN
          ELSE
             WRITE(*,*) 'No core density is available in HDF5 format.'
@@ -645,7 +645,7 @@ MODULE m_cdn_io
          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
                           readDensityIndex,lastDensityIndex)
          CALL writeCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints)
-         CALL closeCDN_HDF(fileID)
+         CALL closeCDNPOT_HDF(fileID)
 #endif
       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
          ! Write core density to cdn.str file
@@ -713,9 +713,9 @@ MODULE m_cdn_io
             WRITE(*,*) 'archiveName: ', TRIM(ADJUSTL(archiveName))
             CALL juDFT_error("For selected starting density index no in-density is present.",calledby ="setStartingDensity")
          END IF
-         CALL writeHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
-                              currentStructureIndex,sdIndex,lastDensityIndex)
-         CALL closeCDN_HDF(fileID)
+         CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
+                                 currentStructureIndex,sdIndex,lastDensityIndex)
+         CALL closeCDNPOT_HDF(fileID)
 #endif
       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
          STOP 'CDN_STREAM_MODE not yet implemented!'

io/pot_io.f90 → io/pot_io.F90

@@ -18,6 +18,10 @@ MODULE m_pot_io
    USE m_juDFT
    USE m_loddop
    USE m_wrtdop
+   USE m_cdnpot_io_hdf
+#ifdef CPP_HDF
+   USE hdf5
+#endif
    IMPLICIT NONE
 
    PRIVATE
@@ -58,19 +62,57 @@ MODULE m_pot_io
       LOGICAL           :: l_exist
       CHARACTER(len=30) :: filename
 
+#ifdef CPP_HDF
+      INTEGER(HID_T)    :: fileID
+#endif
+      INTEGER           :: currentStarsIndex,currentLatharmsIndex
+      INTEGER           :: currentStructureIndex
+      INTEGER           :: potentialType
+      CHARACTER(LEN=30) :: archiveName
+
       CALL getMode(mode)
 
       IF(mode.EQ.POT_HDF5_MODE) THEN
+#ifdef CPP_HDF
          INQUIRE(FILE='pot.hdf',EXIST=l_exist)
          IF (l_exist) THEN
-            !load density from pot.hdf and exit subroutine
+            CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex)
 
-            RETURN
+            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)
          ELSE
-            WRITE(*,*) 'pot.hdf file not found.'
+            WRITE(*,*) 'Potential entry or pot.hdf file not found.'
             WRITE(*,*) 'Falling back to stream access file pot.str.'
             mode = POT_STREAM_MODE
          END IF
+#else
+         WRITE(*,*) 'Not compiled for pot.hdf file usage.'
+         WRITE(*,*) 'Falling back to stream access file pot.str.'
+         mode = POT_STREAM_MODE
+#endif
       END IF
 
       IF(mode.EQ.POT_STREAM_MODE) THEN
@@ -114,15 +156,17 @@ MODULE m_pot_io
 
    END SUBROUTINE readPotential
 
-   SUBROUTINE writePotential(stars,vacuum,atoms,sphhar,input,sym,archiveType,&
+   SUBROUTINE writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,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_cell), INTENT(IN)  :: cell
       TYPE(t_sphhar),INTENT(IN) :: sphhar
       TYPE(t_input),INTENT(IN)  :: input
       TYPE(t_sym),INTENT(IN)    :: sym
+      TYPE(t_oneD),INTENT(IN)   :: oneD
 
       INTEGER, INTENT (IN)      :: iter
       INTEGER, INTENT (IN)      :: archiveType
@@ -133,14 +177,64 @@ MODULE m_pot_io
 
       ! local variables
       INTEGER           :: mode, iUnit
-      LOGICAL           :: l_exist
+      LOGICAL           :: l_exist, l_storeIndices
       CHARACTER(len=30) :: filename
 
+#ifdef CPP_HDF
+      INTEGER(HID_T)    :: fileID
+#endif
+      INTEGER           :: currentStarsIndex,currentLatharmsIndex
+      INTEGER           :: currentStructureIndex
+      INTEGER           :: potentialType
+      CHARACTER(LEN=30) :: archiveName
+
       CALL getMode(mode)
 
       IF(mode.EQ.POT_HDF5_MODE) THEN
-         ! Write potential to pot.hdf file
-         STOP 'POT_HDF5_MODE not yet implemented!'
+#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
+
+         CALL writePotentialHDF(input, fileID, archiveName, potentialType,&
+                                currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
+                                iter,fr,fpw,fz,fzxy)
+
+         IF(l_storeIndices) THEN
+            CALL writePOTHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,&
+                                    currentStructureIndex)
+         END IF
+
+         CALL closeCDNPOT_HDF(fileID)
+#endif
       ELSE IF(mode.EQ.POT_STREAM_MODE) THEN
          ! Write potential to pot.str file
          STOP 'POT_STREAM_MODE not yet implemented!'
@@ -170,8 +264,15 @@ MODULE m_pot_io
       INTEGER, INTENT(OUT) :: mode
 
       mode = POT_DIRECT_MODE
-      IF (juDFT_was_argument("-stream_pot")) mode=POT_STREAM_MODE
-      IF (juDFT_was_argument("-hdf_pot")) mode=POT_HDF5_MODE
+      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
    END SUBROUTINE getMode
 
 END MODULE m_pot_io

main/vgen.F90

@@ -417,7 +417,7 @@ CONTAINS
              vpw_w(1:stars%ng3,js)=vpw_w(1:stars%ng3,js)/stars%nstr(1:stars%ng3)     ! the PW-coulomb part is not
              ! used otherwise anyway.
           ENDDO
-          CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_COUL_const,&
+          CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_COUL_const,&
                               iter,vr,vpw_w,vz,vxy)
           DO js = 1,input%jspins
              DO i = 1,stars%ng3
@@ -836,7 +836,7 @@ CONTAINS
                    vpw_w(i,js)=vpw_w(i,js)/stars%nstr(i)
                 ENDDO
              ENDDO
-             CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_TOT_const,&
+             CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_TOT_const,&
                                  iter,vr,vpw_w,vz,vxy)
 
              DO js=1,input%jspins
@@ -845,7 +845,7 @@ CONTAINS
                 ENDDO
              ENDDO
 
-             CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_X_const,&
+             CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_X_const,&
                                  iter,vxr,vxpw_w,vz,vxy)
           END IF