Merge branch 'develop' of fleur-git:fleur into develop

io/CMakeLists.txt

@@ -15,6 +15,7 @@ io/nocoInputCheck.F90
 io/inpnoco.F90
 io/loddop.f90
 io/cdnpot_io_hdf.F90
+io/cdnpot_io_common.F90
 io/cdn_io.F90
 io/pot_io.F90
 io/rw_inp.f90

io/cdn_io.F90

@@ -19,6 +19,7 @@ MODULE m_cdn_io
    USE m_loddop
    USE m_wrtdop
    USE m_cdnpot_io_hdf
+   USE m_cdnpot_io_common
 #ifdef CPP_HDF
    USE hdf5
 #endif
@@ -129,12 +130,13 @@ MODULE m_cdn_io
    END SUBROUTINE printDensityFileInfo
 
 
-   SUBROUTINE readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
+   SUBROUTINE readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
                           relCdnIndex,fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
       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
@@ -155,7 +157,7 @@ MODULE m_cdn_io
 
       ! local variables
       INTEGER            :: mode, datend, k, i, iVac, j, iUnit
-      LOGICAL            :: l_exist, l_rhomatFile
+      LOGICAL            :: l_exist, l_rhomatFile, l_DimChange
       CHARACTER(LEN=30)  :: filename
 
 #ifdef CPP_HDF
@@ -216,9 +218,14 @@ MODULE m_cdn_io
                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
 
             CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, oneD, archiveName, densityType,&
-                                fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
+                                fermiEnergy,l_qfix,l_DimChange,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
             CALL closeCDNPOT_HDF(fileID)
+
+            IF(l_DimChange) THEN
+               CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
+                           1,-1.0,fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
+            END IF
             RETURN
          ELSE
             WRITE(*,*) 'cdn.hdf file or relevant density entry not found.'
@@ -346,7 +353,7 @@ MODULE m_cdn_io
 
       INTEGER           :: mode, iterTemp, k, i, iVac, j, iUnit
       INTEGER           :: d1, d10, asciioffset, iUnitTemp
-      LOGICAL           :: l_exist, l_storeIndices, l_writeNew
+      LOGICAL           :: l_exist, l_storeIndices, l_writeNew, l_same
       LOGICAL           :: l_writeAll
       CHARACTER(len=30) :: filename
       CHARACTER(len=5)  :: cdnfile
@@ -379,92 +386,10 @@ MODULE m_cdn_io
          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
                           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)
-         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
+         CALL checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, oneD, stars, sphhar, sym,&
+                                       currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
+                                       currentStepfunctionIndex,l_storeIndices)
+
          previousDensityIndex = readDensityIndex
          writeDensityIndex = readDensityIndex
          IF(relCdnIndex.NE.0) THEN
@@ -919,8 +844,10 @@ MODULE m_cdn_io
       LOGICAL, INTENT(IN)         :: l_xcExtended,l_ExtData
       LOGICAL, INTENT(OUT)        :: l_error
 
+
+      TYPE(t_stars)               :: starsTemp
       INTEGER                     :: mode, ioStatus, ngz,izmin,izmax
-      LOGICAL                     :: l_exist
+      LOGICAL                     :: l_exist, l_same
 
       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
@@ -947,7 +874,14 @@ MODULE m_cdn_io
             IF (currentStarsIndex.LT.1) THEN
                mode = CDN_DIRECT_MODE ! (no stars entry found in cdn.hdf file)
             ELSE
-               CALL readStarsHDF(fileID, currentStarsIndex, stars)
+               CALL readStarsHDF(fileID, currentStarsIndex, starsTemp)
+               CALL compareStars(stars, starsTemp, l_same)
+               WRITE(*,*) 'l_same', l_same
+               IF(l_same) THEN
+                  CALL readStarsHDF(fileID, currentStarsIndex, stars)
+               ELSE
+                 mode = CDN_DIRECT_MODE ! (no adequate stars entry found in cdn.hdf file)
+               END IF
             END IF
             CALL closeCDNPOT_HDF(fileID)
 #endif

io/cdnpot_io_common.F90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2017 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 contains common subroutines required for density IO
+!!! as well as for potential IO
+!!!
+!!!                             GM'17
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+MODULE m_cdnpot_io_common
+
+   USE m_types
+   USE m_juDFT
+   USE m_cdnpot_io_hdf
+#ifdef CPP_HDF
+   USE hdf5
+#endif
+
+   CONTAINS
+
+   SUBROUTINE compareStars(stars, refStars, l_same)
+
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_stars),INTENT(IN)  :: refStars
+
+      LOGICAL,      INTENT(OUT) :: l_same
+
+      l_same = .TRUE.
+
+      IF(ABS(stars%gmaxInit-refStars%gmaxInit).GT.1e-10) l_same = .FALSE.
+      IF(stars%ng3.NE.refStars%ng3) l_same = .FALSE.
+      IF(stars%ng2.NE.refStars%ng2) l_same = .FALSE.
+      IF(stars%mx1.NE.refStars%mx1) l_same = .FALSE.
+      IF(stars%mx2.NE.refStars%mx2) l_same = .FALSE.
+      IF(stars%mx3.NE.refStars%mx3) l_same = .FALSE.
+
+   END SUBROUTINE compareStars
+
+   SUBROUTINE compareStepfunctions(stars, refStars, l_same)
+
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_stars),INTENT(IN)  :: refStars
+
+      LOGICAL,      INTENT(OUT) :: l_same
+
+      l_same = .TRUE.
+
+      IF(stars%ng3.NE.refStars%ng3) l_same = .FALSE.
+      IF(stars%mx1.NE.refStars%mx1) l_same = .FALSE.
+      IF(stars%mx2.NE.refStars%mx2) l_same = .FALSE.
+      IF(stars%mx3.NE.refStars%mx3) l_same = .FALSE.
+
+   END SUBROUTINE compareStepfunctions
+
+   SUBROUTINE compareStructure(atoms, vacuum, cell, refAtoms, refVacuum,&
+                               refCell, l_same)
+
+      TYPE(t_atoms),INTENT(IN)  :: atoms, refAtoms
+      TYPE(t_vacuum),INTENT(IN) :: vacuum, refVacuum
+      TYPE(t_cell),INTENT(IN)   :: cell, refCell
+
+      LOGICAL,      INTENT(OUT) :: l_same
+
+      l_same = .TRUE.
+
+      IF(atoms%ntype.NE.refAtoms%ntype) l_same = .FALSE.
+      IF(atoms%nat.NE.refAtoms%nat) l_same = .FALSE.
+      IF(atoms%lmaxd.NE.refAtoms%lmaxd) l_same = .FALSE.
+      IF(atoms%jmtd.NE.refAtoms%jmtd) l_same = .FALSE.
+      IF(vacuum%dvac.NE.refVacuum%dvac) l_same = .FALSE.
+      IF(ANY(ABS(cell%amat(:,:)-refCell%amat(:,:)).GT.1e-10)) l_same = .FALSE.
+      IF(l_same) THEN
+         IF(ANY(atoms%nz(:).NE.refAtoms%nz(:))) l_same = .FALSE.
+         IF(ANY(atoms%lmax(:).NE.refAtoms%lmax(:))) l_same = .FALSE.
+      END IF
+      IF(l_same) THEN
+         DO i = 1, atoms%nat
+            IF(ANY(ABS(atoms%pos(:,i)-refAtoms%pos(:,i)).GT.1e-10)) l_same = .FALSE.
+         END DO
+      END IF
+
+      ! NOTE: This subroutine certainly is not yet complete. Especially symmetry should
+      !       also be stored and compared for structure considerations.
+
+   END SUBROUTINE compareStructure
+
+   SUBROUTINE compareLatharms(latharms, refLatharms, l_same)
+
+      TYPE(t_sphhar)       :: latharms, refLatharms
+
+      LOGICAL,      INTENT(OUT) :: l_same
+
+      l_same = .TRUE.
+
+      IF(latharms%ntypsd.NE.refLatharms%ntypsd) l_same = .FALSE.
+      IF(latharms%memd.NE.refLatharms%memd) l_same = .FALSE.
+      IF(latharms%nlhd.NE.refLatharms%nlhd) l_same = .FALSE.
+
+   END SUBROUTINE compareLatharms
+
+#ifdef CPP_HDF
+   SUBROUTINE checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, oneD, stars, latharms, sym,&
+                                       currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
+                                       currentStepfunctionIndex,l_storeIndices)
+
+      TYPE(t_input),INTENT(IN)  :: input
+      TYPE(t_atoms),INTENT(IN)  :: atoms
+      TYPE(t_cell), INTENT(IN)  :: cell
+      TYPE(t_vacuum),INTENT(IN) :: vacuum
+      TYPE(t_oneD),INTENT(IN)   :: oneD
+      TYPE(t_stars),INTENT(IN)  :: stars
+      TYPE(t_sphhar),INTENT(IN) :: latharms
+      TYPE(t_sym),INTENT(IN)    :: sym
+
+      INTEGER(HID_T), INTENT(IN) :: fileID
+      INTEGER, INTENT(INOUT)     :: currentStarsIndex,currentLatharmsIndex
+      INTEGER, INTENT(INOUT)     :: currentStructureIndex,currentStepfunctionIndex
+      LOGICAL, INTENT(OUT)       :: l_storeIndices
+
+      TYPE(t_stars)        :: starsTemp
+      TYPE(t_vacuum)       :: vacuumTemp
+      TYPE(t_atoms)        :: atomsTemp
+      TYPE(t_sphhar)       :: latharmsTemp
+      TYPE(t_input)        :: inputTemp
+      TYPE(t_sym)          :: symTemp
+      TYPE(t_cell)         :: cellTemp
+      TYPE(t_oneD)         :: oneDTemp
+
+      LOGICAL                    :: l_same, l_writeAll
+
+      l_storeIndices = .FALSE.
+      l_writeAll = .FALSE.
+
+      IF(currentStructureIndex.EQ.0) THEN
+         currentStructureIndex = 1
+         l_storeIndices = .TRUE.
+         CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, currentStructureIndex)
+      ELSE
+         CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, currentStructureIndex)
+         CALL compareStructure(atoms, vacuum, cell, atomsTemp, vacuumTemp, cellTemp, l_same)
+
+         IF(.NOT.l_same) 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)
+         CALL compareStars(stars, starsTemp, l_same)
+         IF((.NOT.l_same).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, latharms)
+      ELSE
+         CALL readLatharmsHDF(fileID, currentLatharmsIndex, latharmsTemp)
+         CALL compareLatharms(latharms, latharmsTemp, l_same)
+         IF((.NOT.l_same).OR.l_writeAll) THEN
+            currentLatharmsIndex = currentLatharmsIndex + 1
+            l_storeIndices = .TRUE.
+            CALL writeLatharmsHDF(fileID, currentLatharmsIndex, latharms)
+         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)
+         CALL compareStepfunctions(stars, starsTemp, l_same)
+         IF((.NOT.l_same).OR.l_writeAll) THEN
+            currentStepfunctionIndex = currentStepfunctionIndex + 1
+            l_storeIndices = .TRUE.
+            CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, stars)
+         END IF
+      END IF
+
+   END SUBROUTINE checkAndWriteMetadataHDF
+#endif
+
+
+END MODULE m_cdnpot_io_common

io/pot_io.F90

 !--------------------------------------------------------------------------------
-! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! Copyright (c) 2017 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.
 !--------------------------------------------------------------------------------
@@ -19,6 +19,7 @@ MODULE m_pot_io
    USE m_loddop
    USE m_wrtdop
    USE m_cdnpot_io_hdf
+   USE m_cdnpot_io_common
 #ifdef CPP_HDF
    USE hdf5
 #endif
@@ -200,27 +201,9 @@ MODULE m_pot_io
          CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,&
                           currentStructureIndex,currentStepfunctionIndex)
 
-         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
-         IF(currentStepfunctionIndex.EQ.0) THEN
-            currentStepfunctionIndex = 1
-            l_storeIndices = .TRUE.
-            CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, stars)
-         END IF
+         CALL checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, oneD, stars, sphhar, sym,&
+                                       currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
+                                       currentStepfunctionIndex,l_storeIndices)
 
          archiveName = 'illegalPotentialArchive'
          IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN

main/cdngen.f90

@@ -106,7 +106,7 @@
       IF(noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
       IF((.NOT.noco%l_noco).AND.mpi%irank.EQ.0) THEN
          ALLOCATE(cdom(1),cdomvz(1,1),cdomvxy(1,1,1))
-         CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
+         CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
                           CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
          DEALLOCATE(cdom,cdomvz,cdomvxy)
       END IF

main/mix.F90

@@ -182,7 +182,7 @@ CONTAINS
     END IF
 
     !---> reload densities of current iteration
-    CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
+    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                      CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
 
     !
@@ -194,7 +194,7 @@ CONTAINS
          intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp(-3,-3,1,1,1), nmap,nmaph,mapmt,mapvac,mapvac2,sm) 
 
     !     load output charge density
-    CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
+    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                      CDN_OUTPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
 
     !

main/totale.f90

@@ -146,7 +146,7 @@ CONTAINS
     archiveType = CDN_ARCHIVE_TYPE_CDN1_const
     IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_CDN_const
 
-    CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
+    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                      CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
 
     !+for

main/vgen.F90

@@ -155,10 +155,10 @@ CONTAINS
           IF (noco%l_noco) THEN
              CALL juDFT_error("vgen:1",calledby ="vgen")
           ENDIF
-          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,CDN_OUTPUT_DEN_const,&
+          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,CDN_OUTPUT_DEN_const,&
                            0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
        ELSE
-          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
+          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
                            0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
        END IF
 
@@ -405,7 +405,7 @@ CONTAINS
        !     ----> reload the density for calculating vxc (for spin-pol. case)
        !
        IF (input%jspins.EQ.2) THEN
-          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
+          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
                            0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
           vr(:,0:,:,2) = vr(:,0:,:,1)
           vpw(:,2) = vpw(:,1)
@@ -678,7 +678,7 @@ CONTAINS
        IF (input%total) THEN
 
           IF (noco%l_noco) THEN ! load qpw,rht,rhtxy from 'cdn'-file
-             CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
+             CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
                               0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
           ENDIF
           !

optional/bmt.f90

@@ -50,7 +50,7 @@ contains
     archiveType = CDN_ARCHIVE_TYPE_CDN1_const
     IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
 
-    CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
+    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                      CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,it,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy)
 
     IF ( typmag < atoms%ntype ) THEN 

optional/cdnsp.f90

@@ -59,7 +59,7 @@
 
           input%jspins=1
           CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
-          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
+          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,&
                            CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
           input%jspins=2
 

optional/flipcdn.f90

@@ -60,7 +60,7 @@
           ENDIF
 
           !     ---> read the charge density 
-          CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,archiveType,&
+          CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                      CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
 
           !     ---> flip cdn for each atom with nflip=-1

vgen/rhodirgen.f90

@@ -116,7 +116,7 @@ CONTAINS
     ENDIF
 
     !---> reload the density matrix from file rhomat_inp
-    CALL readDensity(stars,vacuum,atoms,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_NOCO_const,CDN_INPUT_DEN_const,&
+    CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_NOCO_const,CDN_INPUT_DEN_const,&
                      0,fermiEnergyTemp,l_qfix,iter,rho,qpw,rht,rhtxy,cdom,cdomvz,cdomvxy)
 
     CALL qfix(&