First steps towards changing structure, stars, stepfunction, lattice harmonics in cdn.hdf

There are still bugs. Probably many bugs.

io/cdn_io.F90

@@ -215,7 +215,7 @@ MODULE m_cdn_io
             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
 
-            CALL readDensityHDF(fileID, archiveName, densityType,&
+            CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, oneD, archiveName, densityType,&
                                 fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
             CALL closeCDNPOT_HDF(fileID)
@@ -338,13 +338,16 @@ MODULE m_cdn_io
       TYPE(t_sphhar)       :: sphharTemp
       TYPE(t_input)        :: inputTemp
       TYPE(t_sym)          :: symTemp
+      TYPE(t_cell)         :: cellTemp
+      TYPE(t_oneD)         :: oneDTemp
 
       COMPLEX, ALLOCATABLE :: fpwTemp(:,:), fzxyTemp(:,:,:,:)
       REAL, ALLOCATABLE    :: frTemp(:,:,:,:), fzTemp(:,:,:)
 
       INTEGER           :: mode, iterTemp, k, i, iVac, j, iUnit
       INTEGER           :: d1, d10, asciioffset, iUnitTemp
-      LOGICAL           :: l_exist, l_storeIndices
+      LOGICAL           :: l_exist, l_storeIndices, l_writeNew
+      LOGICAL           :: l_writeAll
       CHARACTER(len=30) :: filename
       CHARACTER(len=5)  :: cdnfile
 
@@ -377,25 +380,90 @@ MODULE m_cdn_io
                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
 
          l_storeIndices = .FALSE.
+         l_writeAll = .FALSE.
+         IF(currentStructureIndex.EQ.0) THEN
+            currentStructureIndex = 1
+            l_storeIndices = .TRUE.
+            l_writeAll = .TRUE.
+            CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, currentStructureIndex)
+         ELSE
+            CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, currentStructureIndex)
+            l_writeNew = .FALSE.
+            IF(atoms%ntype.NE.atomsTemp%ntype) l_writeNew = .TRUE.
+            IF(atoms%nat.NE.atomsTemp%nat) l_writeNew = .TRUE.
+            IF(atoms%lmaxd.NE.atomsTemp%lmaxd) l_writeNew = .TRUE.
+            IF(atoms%jmtd.NE.atomsTemp%jmtd) l_writeNew = .TRUE.
+            IF(vacuum%dvac.NE.vacuumTemp%dvac) l_writeNew = .TRUE.
+            IF(ANY(ABS(cell%amat(:,:)-cellTemp%amat(:,:)).GT.1e-10)) l_writeNew = .TRUE.
+            IF(.NOT.l_writeNew) THEN
+               IF(ANY(atoms%nz(:).NE.atomsTemp%nz(:))) l_writeNew = .TRUE.
+               IF(ANY(atoms%lmax(:).NE.atomsTemp%lmax(:))) l_writeNew = .TRUE.
+            END IF
+            IF(.NOT.l_writeNew) THEN
+               DO i = 1, atoms%nat
+                  IF(ANY(ABS(atoms%pos(:,i)-atomsTemp%pos(:,i)).GT.1e-10)) l_writeNew = .TRUE.
+               END DO
+            END IF
+            IF(l_writeNew) THEN
+               currentStructureIndex = currentStructureIndex + 1
+               l_storeIndices = .TRUE.
+               l_writeAll = .TRUE.
+               CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, currentStructureIndex)
+            END IF
+         END IF
          IF (currentStarsIndex.EQ.0) THEN
             currentStarsIndex = 1
             l_storeIndices = .TRUE.
             CALL writeStarsHDF(fileID, currentStarsIndex, stars)
+         ELSE
+            CALL readStarsHDF(fileID, currentStarsIndex, starsTemp)
+            l_writeNew = .FALSE.
+            IF(ABS(stars%gmax-starsTemp%gmax).GT.1e-10) l_writeNew = .TRUE.
+            IF(stars%ng3.NE.starsTemp%ng3) l_writeNew = .TRUE.
+            IF(stars%ng2.NE.starsTemp%ng2) l_writeNew = .TRUE.
+            IF(stars%mx1.NE.starsTemp%mx1) l_writeNew = .TRUE.
+            IF(stars%mx2.NE.starsTemp%mx2) l_writeNew = .TRUE.
+            IF(stars%mx3.NE.starsTemp%mx3) l_writeNew = .TRUE.
+            IF(stars%kimax.NE.starsTemp%kimax) l_writeNew = .TRUE.
+            IF(stars%kimax2.NE.starsTemp%kimax2) l_writeNew = .TRUE.
+            IF(l_writeNew.OR.l_writeAll) THEN
+               currentStarsIndex = currentStarsIndex + 1
+               l_storeIndices = .TRUE.
+               CALL writeStarsHDF(fileID, currentStarsIndex, stars)
+            END IF
          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)
+         ELSE
+            CALL readLatharmsHDF(fileID, currentLatharmsIndex, sphharTemp)
+            l_writeNew = .FALSE.
+            IF(sphhar%ntypsd.NE.sphharTemp%ntypsd) l_writeNew = .TRUE.
+            IF(sphhar%memd.NE.sphharTemp%memd) l_writeNew = .TRUE.
+            IF(sphhar%nlhd.NE.sphharTemp%nlhd) l_writeNew = .TRUE.
+            IF(l_writeNew.OR.l_writeAll) THEN
+               currentLatharmsIndex = currentLatharmsIndex + 1
+               l_storeIndices = .TRUE.
+               CALL writeLatharmsHDF(fileID, currentLatharmsIndex, sphhar)
+            END IF
          END IF
          IF(currentStepfunctionIndex.EQ.0) THEN
             currentStepfunctionIndex = 1
             l_storeIndices = .TRUE.
             CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, stars)
+         ELSE
+            CALL readStepfunctionHDF(fileID, currentStepfunctionIndex, starsTemp)
+            l_writeNew = .FALSE.
+            IF(stars%ng3.NE.starsTemp%ng3) l_writeNew = .TRUE.
+            IF(stars%mx1.NE.starsTemp%mx1) l_writeNew = .TRUE.
+            IF(stars%mx2.NE.starsTemp%mx2) l_writeNew = .TRUE.
+            IF(stars%mx3.NE.starsTemp%mx3) l_writeNew = .TRUE.
+            IF(l_writeNew.OR.l_writeAll) THEN
+               currentStepfunctionIndex = currentStepfunctionIndex + 1
+               l_storeIndices = .TRUE.
+               CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, stars)
+            END IF
          END IF
          previousDensityIndex = readDensityIndex
          writeDensityIndex = readDensityIndex

io/cdnpot_io_hdf.F90

@@ -1196,7 +1196,7 @@ MODULE m_cdnpot_io_hdf
       CALL h5dclose_f(neqSetID, hdfError)
 
       dimsInt(:1)=(/atoms%ntype/)
-      CALL h5dopen_f(groupID, 'jri', nzSetID, hdfError)
+      CALL h5dopen_f(groupID, 'jri', jriSetID, hdfError)
       CALL io_read_integer1(jriSetID,(/1/),dimsInt(:1),atoms%jri)
       CALL h5dclose_f(jriSetID, hdfError)
 
@@ -1801,9 +1801,17 @@ MODULE m_cdnpot_io_hdf
 
    END SUBROUTINE writePotentialHDF
 
-   SUBROUTINE readDensityHDF(fileID, archiveName, densityType,&
+   SUBROUTINE readDensityHDF(fileID, input, stars, latharms, atoms, vacuum, oneD,&
+                             archiveName, densityType,&
                              fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
+      TYPE(t_input),INTENT(IN)  :: input
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_sphhar),INTENT(IN) :: latharms
+      TYPE(t_atoms),INTENT(IN)  :: atoms
+      TYPE(t_vacuum),INTENT(IN) :: vacuum
+      TYPE(t_oneD),INTENT(IN)   :: oneD
+
       INTEGER(HID_T), INTENT(IN)   :: fileID
       INTEGER, INTENT(IN)          :: densityType
       CHARACTER(LEN=*), INTENT(IN) :: archiveName
@@ -1930,23 +1938,27 @@ MODULE m_cdnpot_io_hdf
       CALL io_read_attreal0(groupID,'fermiEnergy',fermiEnergy)
       CALL io_read_attlog0(groupID,'l_qfix',l_qfix)
 
-      dimsInt(:4)=(/jmtd,nlhd+1,ntype,jspins/)
+      fr = 0.0
+      dimsInt(:4)=(/MIN(jmtd,atoms%jmtd),MIN(nlhd,latharms%nlhd)+1,MIN(ntype,atoms%ntype),MIN(jspins,input%jspins)/)
       CALL h5dopen_f(groupID, 'fr', frSetID, hdfError)
       CALL io_read_real4(frSetID,(/1,1,1,1/),dimsInt(:4),fr)
       CALL h5dclose_f(frSetID, hdfError)
 
-      dimsInt(:3)=(/2,ng3,jspins/)
+      fpw = CMPLX(0.0,0.0)
+      dimsInt(:3)=(/2,MIN(ng3,stars%ng3),MIN(jspins,input%jspins)/)
       CALL h5dopen_f(groupID, 'fpw', fpwSetID, hdfError)
       CALL io_read_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),fpw)
       CALL h5dclose_f(fpwSetID, hdfError)
 
       IF (l_film) THEN
-         dimsInt(:3)=(/nmzd,2,jspins/)
+         fz = 0.0
+         dimsInt(:3)=(/MIN(nmzd,vacuum%nmzd),2,MIN(jspins,input%jspins)/)
          CALL h5dopen_f(groupID, 'fz', fzSetID, hdfError)
          CALL io_read_real3(fzSetID,(/1,1,1/),dimsInt(:3),fz)
          CALL h5dclose_f(fzSetID, hdfError)
 
-         dimsInt(:5)=(/2,nmzxyd,ng2-1,2,jspins/)
+         fzxy = CMPLX(0.0,0.0)
+         dimsInt(:5)=(/2,MIN(nmzxyd,vacuum%nmzxyd),MIN(ng2,stars%ng2)-1,2,MIN(jspins,input%jspins)/)
          CALL h5dopen_f(groupID, 'fzxy', fzxySetID, hdfError)
          CALL io_read_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),fzxy)
          CALL h5dclose_f(fzxySetID, hdfError)
@@ -1955,18 +1967,21 @@ MODULE m_cdnpot_io_hdf
       IF((localDensityType.EQ.DENSITY_TYPE_NOCO_IN_const).OR.&
          (localDensityType.EQ.DENSITY_TYPE_NOCO_OUT_const)) THEN
 
-         dimsInt(:2)=(/2,ng3/)
+         cdom = CMPLX(0.0,0.0)
+         dimsInt(:2)=(/2,MIN(ng3,stars%ng3)/)
          CALL h5dopen_f(groupID, 'cdom', cdomSetID, hdfError)
          CALL io_read_complex1(cdomSetID,(/-1,1/),dimsInt(:2),cdom)
          CALL h5dclose_f(cdomSetID, hdfError)
 
          IF (l_film) THEN
-            dimsInt(:3)=(/2,nmz,nvac/)
+            cdomvz = CMPLX(0.0,0.0)
+            dimsInt(:3)=(/2,MIN(nmz,vacuum%nmz),MIN(nvac,vacuum%nvac)/)
             CALL h5dopen_f(groupID, 'cdomvz', cdomvzSetID, hdfError)
             CALL io_read_complex2(cdomvzSetID,(/-1,1,1/),dimsInt(:3),cdomvz)
             CALL h5dclose_f(cdomvzSetID, hdfError)
 
-            dimsInt(:4)=(/2,nmzxy,od_nq2-1,nvac/)
+            cdomvxy = CMPLX(0.0,0.0)
+            dimsInt(:4)=(/2,MIN(nmzxy,vacuum%nmzxy),MIN(od_nq2,oneD%odi%nq2)-1,MIN(nvac,vacuum%nvac)/)
             CALL h5dopen_f(groupID, 'cdomvxy', cdomvxySetID, hdfError)
             CALL io_read_complex3(cdomvxySetID,(/-1,1,1,1/),dimsInt(:4),cdomvxy)
             CALL h5dclose_f(cdomvxySetID, hdfError)