LCOV - code coverage report
Current view: top level - io - cdn_io.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 267 676 39.5 %
Date: 2019-09-08 04:53:50 Functions: 14 17 82.4 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2017 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       8             : !!!
       9             : !!! This module is a wrapper for the charge density I/O
      10             : !!!
      11             : !!!                             GM'17
      12             : !!!
      13             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      14             : 
      15             : MODULE m_cdn_io
      16             : 
      17             :    USE m_types
      18             :    USE m_juDFT
      19             :    USE m_loddop
      20             :    USE m_wrtdop
      21             :    USE m_cdnpot_io_hdf
      22             :    USE m_cdnpot_io_common
      23             :    USE m_constants
      24             : #ifdef CPP_HDF
      25             :    USE hdf5
      26             : #endif
      27             :    IMPLICIT NONE
      28             : 
      29             :    PRIVATE
      30             :    PUBLIC printDensityFileInfo
      31             :    PUBLIC readDensity, writeDensity
      32             :    PUBLIC isDensityFilePresent, isCoreDensityPresent
      33             :    PUBLIC readCoreDensity, writeCoreDensity
      34             :    PUBLIC readStars, writeStars
      35             :    PUBLIC readStepfunction, writeStepfunction
      36             :    PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
      37             :    PUBLIC storeStructureIfNew,transform_by_moving_atoms
      38             :    PUBLIC getIOMode
      39             :    PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
      40             :    PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
      41             :    PUBLIC CDN_ARCHIVE_TYPE_CDN_const
      42             : 
      43             :    INTEGER,          PARAMETER :: CDN_INPUT_DEN_const = 1
      44             :    INTEGER,          PARAMETER :: CDN_OUTPUT_DEN_const = 2
      45             : 
      46             :    INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_CDN1_const = 1
      47             :    INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_NOCO_const = 2
      48             :    INTEGER,          PARAMETER :: CDN_ARCHIVE_TYPE_CDN_const  = 3
      49             : 
      50             :    INTEGER,          PARAMETER :: CDN_DIRECT_MODE = 1
      51             :    INTEGER,          PARAMETER :: CDN_STREAM_MODE = 2
      52             :    INTEGER,          PARAMETER :: CDN_HDF5_MODE   = 3
      53             : 
      54             :    CONTAINS
      55             : 
      56           0 :    SUBROUTINE printDensityFileInfo()
      57             : 
      58             :       INTEGER            :: mode, i
      59             :       LOGICAL            :: l_exist
      60             : 
      61             : #ifdef CPP_HDF
      62             :       INTEGER(HID_T)    :: fileID
      63             : #endif
      64             :       INTEGER           :: currentStarsIndex,currentLatharmsIndex
      65             :       INTEGER           :: currentStructureIndex,currentStepfunctionIndex
      66             :       INTEGER           :: readDensityIndex, lastDensityIndex
      67             :       CHARACTER(LEN=30) :: archiveName
      68             : 
      69             :       INTEGER           :: dateTemp, timeTemp
      70             :       INTEGER           :: iterTemp, starsIndexTemp, latharmsIndexTemp 
      71             :       INTEGER           :: structureIndexTemp,stepfunctionIndexTemp
      72             :       INTEGER           :: previousDensityIndex, jspinsTemp
      73             :       REAL              :: fermiEnergyTemp, distanceTemp
      74             :       LOGICAL           :: l_qfixTemp
      75             :       CHARACTER(LEN=10) :: dateString
      76             :       CHARACTER(LEN=10) :: timeString
      77             :       CHARACTER(LEN=19) :: timeStampString
      78             :       CHARACTER(LEN=15) :: distanceString
      79             : 
      80           0 :       CALL getIOMode(mode)
      81             : 
      82           0 :       WRITE(*,*) 'Available densities info:'
      83           0 :       WRITE(*,*)
      84             : 
      85           0 :       IF(mode.EQ.CDN_HDF5_MODE) THEN
      86             : #ifdef CPP_HDF
      87           0 :          INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
      88           0 :          IF (l_exist) THEN
      89             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
      90           0 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
      91           0 :             WRITE(*,*) 'densityIndex   iteration   prevDensity   prevDistance        timeStamp'
      92           0 :             DO i = 1, lastDensityIndex
      93           0 :                archiveName = ''
      94           0 :                WRITE(archiveName,'(a,i0)') '/cdn-', i
      95             : 
      96           0 :                l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
      97           0 :                IF(.NOT.l_exist) THEN
      98             :                   CYCLE
      99             :                END IF
     100             : 
     101             :                CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_UNDEFINED_const,&
     102             :                                         iterTemp, starsIndexTemp, latharmsIndexTemp, structureIndexTemp,&
     103             :                                         stepfunctionIndexTemp,previousDensityIndex, jspinsTemp,&
     104           0 :                                         dateTemp, timeTemp, distanceTemp, fermiEnergyTemp, l_qfixTemp)
     105             : 
     106           0 :                WRITE(dateString,'(i8)') dateTemp
     107           0 :                WRITE(timeString,'(i6)') timeTemp
     108             : 
     109           0 :                distanceString = ''
     110           0 :                IF (distanceTemp.GE.-1e-10) THEN
     111           0 :                   WRITE(distanceString,'(f15.8)') distanceTemp
     112             :                END IF
     113             : 
     114             :                WRITE(timeStampString,'(a4,a1,a2,a1,a2,1x,a2,a1,a2,a1,a2)') &
     115           0 :                   dateString(1:4),'/',dateString(5:6),'/',dateString(7:8),&
     116           0 :                   timeString(1:2),':',timeString(3:4),':',timeString(5:6)
     117             : 
     118           0 :                WRITE(*,'(1x,i7,6x,i7,7x,i7,4x,a15,3x,a)') i, iterTemp, previousDensityIndex, distanceString,&
     119           0 :                                                             TRIM(ADJUSTL(timeStampString))
     120             :             END DO
     121           0 :             CALL closeCDNPOT_HDF(fileID)
     122             :          ELSE
     123           0 :             WRITE(*,'(a)') "No cdn.hdf file found. Density file info is not available."
     124             :          END IF
     125             : #else
     126             :          WRITE(*,'(a)') "Fleur is not compiled with HDF5 support. Density file info is not available."
     127             : #endif
     128             :       ELSE
     129           0 :          WRITE(*,'(a)') "Density file info is only available if '-hdf_cdn' switch is used."
     130             :       END IF
     131             : 
     132           0 :    END SUBROUTINE printDensityFileInfo
     133             : 
     134             : 
     135          35 :    SUBROUTINE readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
     136             :                           relCdnIndex,fermiEnergy,l_qfix,den,inFilename)
     137             : 
     138             :       TYPE(t_stars),INTENT(IN)     :: stars
     139             :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
     140             :       TYPE(t_atoms),INTENT(IN)     :: atoms
     141             :       TYPE(t_cell), INTENT(IN)     :: cell
     142             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
     143             :       TYPE(t_input),INTENT(IN)     :: input
     144             :       TYPE(t_sym),INTENT(IN)       :: sym
     145             :       TYPE(t_oneD),INTENT(IN)      :: oneD
     146             :       TYPE(t_potden),INTENT(INOUT) :: den
     147             : 
     148             :       INTEGER, INTENT (IN)      :: inOrOutCDN
     149             :       INTEGER, INTENT (IN)      :: relCdnIndex
     150             :       INTEGER, INTENT (IN)      :: archiveType
     151             :       REAL,    INTENT (OUT)     :: fermiEnergy
     152             :       LOGICAL, INTENT (OUT)     :: l_qfix
     153             : 
     154             :       CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: inFilename
     155             : 
     156             :       ! local variables
     157             :       INTEGER            :: mode, datend, k, i, iVac, j, iUnit, l, numLines, ioStatus, iofl
     158             :       LOGICAL            :: l_exist, l_rhomatFile, l_DimChange
     159             :       CHARACTER(LEN=30)  :: filename
     160             : 
     161             : #ifdef CPP_HDF
     162             :       INTEGER(HID_T) :: fileID
     163             : #endif
     164             :       INTEGER           :: currentStarsIndex,currentLatharmsIndex
     165             :       INTEGER           :: currentStructureIndex,currentStepfunctionIndex
     166             :       INTEGER           :: readDensityIndex, lastDensityIndex
     167             :       INTEGER           :: previousDensityIndex, densityType
     168             :       CHARACTER(LEN=30) :: archiveName
     169             :       TYPE(t_cell)      :: cellTemp
     170             : 
     171          35 :       COMPLEX, ALLOCATABLE :: cdomvz(:,:)
     172             : 
     173          35 :       fermiEnergy = 0.0
     174          35 :       l_qfix = .FALSE.
     175             : 
     176          35 :       CALL getIOMode(mode)
     177             : 
     178             : #ifndef CPP_HDF
     179             :       filename = 'cdn.hdf'
     180             :       IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
     181             :       INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     182             :       IF (l_exist) THEN
     183             :          CALL juDFT_warn('Fleur not compiled for HDF5, but '//TRIM(ADJUSTL(filename))//' present',calledby='readDensity')
     184             :       END IF
     185             : #endif
     186             : 
     187             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     188             : #ifdef CPP_HDF
     189             : 
     190          35 :          densityType = 0
     191          35 :          archiveName = ''
     192             : 
     193          35 :          filename = 'cdn.hdf'
     194          35 :          IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
     195             : 
     196          35 :          INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     197          35 :          IF (l_exist) THEN
     198             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     199          35 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
     200             : 
     201          35 :             IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
     202           1 :                archiveName = 'cdn'
     203             :             ELSE
     204          34 :                WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
     205             :             END IF
     206             : 
     207          70 :             SELECT CASE (inOrOutCDN)
     208             :                CASE (CDN_INPUT_DEN_const)
     209          35 :                   IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     210           7 :                      densityType = DENSITY_TYPE_NOCO_IN_const
     211             :                   ELSE
     212          28 :                      densityType = DENSITY_TYPE_IN_const
     213             :                   END IF
     214             :                CASE (CDN_OUTPUT_DEN_const)
     215           0 :                   IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     216           0 :                      densityType = DENSITY_TYPE_NOCO_OUT_const
     217             :                   ELSE
     218           0 :                      densityType = DENSITY_TYPE_OUT_const
     219             :                   END IF
     220             :                CASE DEFAULT
     221           0 :                   WRITE(*,*) 'inOrOutCDN = ', inOrOutCDN
     222          35 :                   CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="readDensity")
     223             :             END SELECT
     224          35 :             l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType)
     225          35 :             CALL closeCDNPOT_HDF(fileID)
     226             :          END IF
     227             : 
     228          35 :          IF (l_exist) THEN
     229             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     230          35 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
     231             : 
     232             :             CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, oneD, archiveName, densityType,&
     233          35 :                                 fermiEnergy,l_qfix,l_DimChange,den)
     234             : 
     235          35 :             CALL closeCDNPOT_HDF(fileID)
     236             : 
     237          35 :             IF(l_DimChange) THEN
     238             :                CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
     239           0 :                            1,-1.0,fermiEnergy,l_qfix,den)
     240             :             END IF
     241             :          ELSE
     242           0 :             INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     243           0 :             IF(l_exist) THEN
     244           0 :                WRITE(*,*) 'densityType is ', densityType
     245           0 :                WRITE(*,*) 'Relevant density entry is '//TRIM(ADJUSTL(archiveName))//'.'
     246           0 :                WRITE(*,*) 'Entry not found in '//TRIM(ADJUSTL(filename))//'.'
     247             :             ELSE
     248           0 :                WRITE(*,*) TRIM(ADJUSTL(filename))//' file not found.'
     249             :             END IF
     250           0 :             WRITE(*,*) 'Falling back to stream access.'
     251           0 :             mode = CDN_STREAM_MODE
     252             :          END IF
     253             : #endif
     254             :       END IF
     255             : 
     256          35 :       IF(mode.EQ.CDN_STREAM_MODE) THEN
     257           0 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
     258           0 :          IF (l_exist) THEN
     259             :             !load density from cdn.str and exit subroutine
     260             : 
     261             :          ELSE
     262           0 :             WRITE(*,*) 'cdn.str file not found.'
     263           0 :             WRITE(*,*) 'Falling back to direct access file cdn1.'
     264           0 :             mode = CDN_DIRECT_MODE
     265             :          END IF
     266             :       END IF
     267             : 
     268          35 :       IF (mode.EQ.CDN_DIRECT_MODE) THEN
     269           0 :          filename = 'cdn1'
     270           0 :          l_rhomatFile = .FALSE.
     271           0 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     272           0 :             INQUIRE(file="rhomat_inp",EXIST=l_exist)
     273           0 :             IF (l_exist) filename = 'rhomat_inp'
     274           0 :             IF (inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) THEN
     275           0 :                INQUIRE(file="rhomat_out",EXIST=l_exist)
     276           0 :                IF (l_exist) filename = 'rhomat_out'
     277             :             END IF
     278           0 :             IF(l_exist) l_rhomatFile = .TRUE.
     279             :          END IF
     280           0 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
     281           0 :             filename = 'cdn'
     282             :          END IF
     283             : 
     284           0 :          IF(PRESENT(inFilename)) filename = inFilename
     285             : 
     286           0 :          INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     287           0 :          IF(.NOT.l_exist) THEN
     288           0 :             CALL juDFT_error("charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readDensity")
     289             :          END IF
     290             : 
     291           0 :          iUnit = 93
     292           0 :          OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='old')
     293             : 
     294           0 :          IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
     295             :             ! call loddop to move the file position to the output density
     296             :             CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
     297           0 :                         iUnit,den%iter,den%mt,den%pw,den%vacz,den%vacxy)
     298             :          END IF
     299             : 
     300             :          ! read in the density
     301             :          CALL loddop(stars,vacuum,atoms,sphhar,input,sym,&
     302           0 :                      iUnit,den%iter,den%mt,den%pw,den%vacz,den%vacxy)
     303             : 
     304             :          ! read in additional data if l_noco and data is present
     305           0 :          IF ((archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const).AND.l_rhomatFile) THEN
     306           0 :             READ (iUnit,iostat=datend) (den%pw(k,3),k=1,stars%ng3)
     307           0 :             IF (datend == 0) THEN
     308           0 :                IF (input%film) THEN
     309           0 :                   ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
     310           0 :                   READ (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
     311           0 :                   DO iVac = 1, vacuum%nvac
     312           0 :                      DO i = 1, vacuum%nmz
     313           0 :                         den%vacz(i,iVac,3) = REAL(cdomvz(i,iVac))
     314           0 :                         den%vacz(i,iVac,4) = AIMAG(cdomvz(i,iVac))
     315             :                      END DO
     316             :                   END DO
     317           0 :                   DEALLOCATE(cdomvz)
     318           0 :                   READ (iUnit) (((den%vacxy(i,j-1,iVac,3),i=1,vacuum%nmzxy),j=2,stars%ng2), iVac=1,vacuum%nvac)
     319             :                END IF
     320             :             ELSE
     321             :                ! (datend < 0)  =>  no off-diagonal magnetisation stored
     322             :                !                   in "rhomat_inp"
     323           0 :                IF (datend > 0) THEN
     324           0 :                   WRITE(*,*) 'datend = ', datend
     325           0 :                   CALL juDFT_error("density file has illegal state.",calledby ="readDensity")
     326             :                END IF
     327           0 :                den%pw(:,3) = CMPLX(0.0,0.0)
     328           0 :                IF (input%film) THEN
     329           0 :                   den%vacz(:,:,3:4) = 0.0
     330           0 :                   den%vacxy(:,:,:,3) = CMPLX(0.0,0.0)
     331             :                END IF
     332             :             END IF
     333           0 :          ELSE IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     334           0 :             den%pw(:,3) = CMPLX(0.0,0.0)
     335           0 :             IF (input%film) THEN
     336           0 :                den%vacz(:,:,3:4) = 0.0
     337           0 :                den%vacxy(:,:,:,3) = CMPLX(0.0,0.0)
     338             :             END IF
     339             :          END IF
     340           0 :          CLOSE(iUnit)
     341             :       END IF
     342             : 
     343          35 :       INQUIRE(FILE='n_mmp_mat',EXIST=l_exist)
     344          35 :       IF(l_exist.AND.atoms%n_u.GT.0) THEN
     345           1 :          OPEN (69,file='n_mmp_mat',status='unknown',form='formatted')
     346           1 :          READ (69,'(7f20.13)',IOSTAT=ioStatus) den%mmpMat
     347           1 :          REWIND(69)
     348           1 :          numLines = 0
     349          56 :          DO
     350          57 :             READ (69,*,iostat=iofl)
     351          57 :             IF (iofl < 0) EXIT
     352          57 :             numLines = numLines + 1
     353             :          END DO
     354           1 :          IF (MOD(numLines,14*input%jspins).NE.0) THEN
     355           0 :             WRITE(*,*) 'The n_mmp_mat file could not be read.'
     356           0 :             WRITE(*,*) 'Was it an old style file with linear mixing parameter specification'
     357           0 :             WRITE(*,*) 'in the last line? Linear mixing for the density matrix can now be'
     358           0 :             WRITE(*,*) 'activated and specified in the inp.xml file.'
     359           0 :             CALL juDFT_error("strange n_mmp_mat-file...", calledby = "readDensity")
     360             :          END IF
     361           1 :          CLOSE(69)
     362             :       END IF
     363             : 
     364          70 :    END SUBROUTINE readDensity
     365             : 
     366         223 :    SUBROUTINE writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
     367             :                            relCdnIndex,distance,fermiEnergy,l_qfix,den,inFilename)
     368             : 
     369             :       TYPE(t_stars),INTENT(IN)     :: stars
     370             :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
     371             :       TYPE(t_atoms),INTENT(IN)     :: atoms
     372             :       TYPE(t_cell), INTENT(IN)     :: cell
     373             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
     374             :       TYPE(t_input),INTENT(IN)     :: input
     375             :       TYPE(t_sym),INTENT(IN)       :: sym
     376             :       TYPE(t_oneD),INTENT(IN)      :: oneD
     377             :       TYPE(t_potden),INTENT(INOUT) :: den
     378             : 
     379             :       INTEGER, INTENT (IN)      :: inOrOutCDN
     380             :       INTEGER, INTENT (IN)      :: relCdnIndex
     381             :       INTEGER, INTENT (IN)      :: archiveType
     382             :       REAL,    INTENT (IN)      :: fermiEnergy, distance
     383             :       LOGICAL, INTENT (IN)      :: l_qfix
     384             : 
     385             :       CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: inFilename
     386             : 
     387         223 :       TYPE(t_stars)        :: starsTemp
     388         223 :       TYPE(t_vacuum)       :: vacuumTemp
     389         223 :       TYPE(t_atoms)        :: atomsTemp
     390         223 :       TYPE(t_sphhar)       :: sphharTemp
     391             :       TYPE(t_input)        :: inputTemp
     392         223 :       TYPE(t_sym)          :: symTemp
     393             :       TYPE(t_cell)         :: cellTemp
     394             :       TYPE(t_oneD)         :: oneDTemp
     395             : 
     396         446 :       COMPLEX, ALLOCATABLE :: fpwTemp(:,:), fzxyTemp(:,:,:,:)
     397         223 :       REAL, ALLOCATABLE    :: frTemp(:,:,:,:), fzTemp(:,:,:)
     398             : 
     399             :       INTEGER           :: mode, iterTemp, k, i, iVac, j, iUnit
     400             :       INTEGER           :: d1, d10, asciioffset, iUnitTemp
     401             :       LOGICAL           :: l_exist, l_storeIndices, l_writeNew, l_same
     402             :       LOGICAL           :: l_writeAll
     403             :       CHARACTER(len=30) :: filename
     404             :       CHARACTER(len=5)  :: cdnfile
     405             : 
     406             : #ifdef CPP_HDF
     407             :       INTEGER(HID_T) :: fileID
     408             : #endif
     409             :       INTEGER           :: currentStarsIndex,currentLatharmsIndex
     410             :       INTEGER           :: currentStructureIndex,currentStepfunctionIndex
     411             :       INTEGER           :: readDensityIndex, writeDensityIndex, lastDensityIndex
     412             :       INTEGER           :: previousDensityIndex, densityType
     413             :       INTEGER           :: starsIndexTemp, latharmsIndexTemp, structureIndexTemp
     414             :       INTEGER           :: stepfunctionIndexTemp
     415             :       INTEGER           :: jspinsTemp
     416             :       INTEGER           :: date, time, dateTemp, timeTemp
     417             :       REAL              :: fermiEnergyTemp, distanceTemp
     418             :       LOGICAL           :: l_qfixTemp, l_CheckBroyd
     419             :       CHARACTER(LEN=30) :: archiveName
     420             :       CHARACTER(LEN=8)  :: dateString
     421             :       CHARACTER(LEN=10) :: timeString
     422             :       CHARACTER(LEN=10) :: zone
     423             : 
     424         223 :       COMPLEX, ALLOCATABLE :: cdomvz(:,:)
     425             : 
     426         223 :       CALL getIOMode(mode)
     427         223 :       CALL DATE_AND_TIME(dateString,timeString,zone)
     428         223 :       READ(dateString,'(i8)') date
     429         223 :       READ(timeString,'(i6)') time
     430             : 
     431         223 :       l_CheckBroyd = .TRUE.
     432         223 :       IF(PRESENT(inFilename)) l_CheckBroyd = .FALSE.
     433             : 
     434         223 :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     435             : #ifdef CPP_HDF
     436             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     437         223 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
     438             : 
     439             :          CALL checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, oneD, stars, sphhar, sym,&
     440             :                                        currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     441         223 :                                        currentStepfunctionIndex,l_storeIndices,l_CheckBroyd)
     442             : 
     443         223 :          previousDensityIndex = readDensityIndex
     444         223 :          writeDensityIndex = readDensityIndex
     445         223 :          IF(relCdnIndex.NE.0) THEN
     446         189 :             writeDensityIndex = lastDensityIndex+relCdnIndex
     447         189 :             lastDensityIndex = writeDensityIndex
     448         189 :             readDensityIndex = writeDensityIndex
     449         189 :             l_storeIndices = .TRUE.
     450             :          END IF
     451             : 
     452         223 :          archiveName = ''
     453         223 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
     454           1 :             archiveName = 'cdn'
     455             :          ELSE
     456         222 :             WRITE(archiveName,'(a,i0)') '/cdn-', writeDensityIndex
     457             :          END IF
     458             : 
     459         223 :          densityType = 0
     460         446 :          SELECT CASE (inOrOutCDN)
     461             :             CASE (CDN_INPUT_DEN_const)
     462         223 :                IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     463         131 :                   densityType = DENSITY_TYPE_NOCO_IN_const
     464             :                ELSE
     465          92 :                   densityType = DENSITY_TYPE_IN_const
     466             :                END IF
     467             :             CASE (CDN_OUTPUT_DEN_const)
     468           0 :                IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     469           0 :                   densityType = DENSITY_TYPE_NOCO_OUT_const
     470             :                ELSE
     471           0 :                   densityType = DENSITY_TYPE_OUT_const
     472             :                END IF
     473             :             CASE DEFAULT
     474           0 :                WRITE(*,*) 'inOrOutCDN = ', inOrOutCDN
     475         223 :                CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="writeDensity")
     476             :          END SELECT
     477             : 
     478         223 :          IF(relCdnIndex.EQ.0) THEN
     479          34 :             l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
     480          34 :             IF(l_exist) THEN
     481             :                CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_UNDEFINED_const,&
     482             :                                         iterTemp, starsIndexTemp, latharmsIndexTemp, structureIndexTemp,&
     483             :                                         stepfunctionIndexTemp,previousDensityIndex, jspinsTemp,&
     484          33 :                                         dateTemp, timeTemp, distanceTemp, fermiEnergyTemp, l_qfixTemp)
     485             :             END IF
     486             :          END IF
     487             : 
     488         223 :          IF(vacuum%nvac.EQ.1) THEN
     489         163 :             den%vacz(:,2,:)=den%vacz(:,1,:)
     490         163 :             IF (sym%invs) THEN
     491         163 :                den%vacxy(:,:,2,:) = CONJG(den%vacxy(:,:,1,:))
     492             :             ELSE
     493           0 :                den%vacxy(:,:,2,:) = den%vacxy(:,:,1,:)
     494             :             END IF
     495             :          END IF
     496             : 
     497             :          CALL writeDensityHDF(input, fileID, archiveName, densityType, previousDensityIndex,&
     498             :                               currentStarsIndex, currentLatharmsIndex, currentStructureIndex,&
     499             :                               currentStepfunctionIndex,date,time,distance,fermiEnergy,l_qfix,&
     500         223 :                               den%iter+relCdnIndex,den)
     501             : 
     502         223 :          IF(l_storeIndices) THEN
     503             :             CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     504         190 :                                     currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     505             :          END IF
     506             : 
     507         223 :          CALL closeCDNPOT_HDF(fileID)
     508             : #endif
     509             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
     510             :          ! Write density to cdn.str file
     511             :          STOP 'CDN_STREAM_MODE not yet implemented!'
     512             :       ELSE
     513           0 :          filename = 'cdn1'
     514           0 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     515           0 :             filename = 'rhomat_inp'
     516           0 :             IF(inOrOutCDN.EQ.CDN_OUTPUT_DEN_const) filename = 'rhomat_out'
     517             :          END IF
     518           0 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
     519           0 :             filename = 'cdn'
     520             :          END IF
     521             : 
     522           0 :          IF(PRESENT(inFilename)) filename = inFilename
     523             : 
     524           0 :          IF ((relCdnIndex.EQ.1).AND.(archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(den%iter.EQ.0)) THEN
     525           0 :             INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     526           0 :             IF(l_exist) THEN
     527           0 :                CALL juDFT_error("Trying to generate starting density while a density exists.",calledby ="writeDensity")
     528             :             END IF
     529             :          END IF
     530             : 
     531           0 :          iUnit = 93
     532           0 :          OPEN (iUnit,file=TRIM(ADJUSTL(filename)),FORM='unformatted',STATUS='unknown')
     533             : 
     534           0 :          IF ((relCdnIndex.EQ.1).AND.(archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(den%iter.GE.1)) THEN
     535           0 :             inputTemp%jspins = input%jspins
     536           0 :             vacuumTemp%nmzxyd = vacuum%nmzxyd
     537           0 :             atomsTemp%jmtd = atoms%jmtd
     538           0 :             sphharTemp%nlhd = sphhar%nlhd
     539           0 :             vacuumTemp%nmzd = vacuum%nmzd
     540           0 :             atomsTemp%ntype = atoms%ntype
     541           0 :             ALLOCATE (sphharTemp%nlh(SIZE(sphhar%nlh)))
     542           0 :             sphharTemp%nlh(:) = sphhar%nlh(:)
     543           0 :             ALLOCATE (atomsTemp%ntypsy(SIZE(atoms%ntypsy)))
     544           0 :             atomsTemp%ntypsy(:) = atoms%ntypsy(:)
     545           0 :             ALLOCATE (atomsTemp%jri(SIZE(atoms%jri)))
     546           0 :             atomsTemp%jri(:) = atoms%jri(:)
     547           0 :             ALLOCATE (atomsTemp%neq(SIZE(atoms%neq)))
     548           0 :             atomsTemp%neq(:) = atoms%neq(:)
     549           0 :             ALLOCATE (atomsTemp%zatom(SIZE(atoms%zatom)))
     550           0 :             atomsTemp%zatom(:) = atoms%zatom(:)
     551           0 :             ALLOCATE (atomsTemp%rmt(SIZE(atoms%rmt)))
     552           0 :             atomsTemp%rmt(:) = atoms%rmt(:)
     553           0 :             ALLOCATE (atomsTemp%dx(SIZE(atoms%dx)))
     554           0 :             atomsTemp%dx(:) = atoms%dx(:)
     555           0 :             starsTemp%ng3 = stars%ng3
     556           0 :             symTemp%invs = sym%invs
     557           0 :             inputTemp%film = input%film
     558           0 :             vacuumTemp%nvac = vacuum%nvac
     559           0 :             vacuumTemp%nmzxy = vacuum%nmzxy
     560           0 :             vacuumTemp%nmz = vacuum%nmz
     561           0 :             vacuumTemp%dvac = vacuum%dvac
     562           0 :             vacuumTemp%delz = vacuum%delz
     563           0 :             starsTemp%ng2 = stars%ng2
     564           0 :             symTemp%invs2 = sym%invs2
     565           0 :             ALLOCATE (fpwTemp(stars%ng3,input%jspins))
     566           0 :             ALLOCATE (fzxyTemp(vacuum%nmzxyd,stars%ng2-1,2,input%jspins))
     567           0 :             ALLOCATE (frTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins))
     568           0 :             ALLOCATE (fzTemp(vacuum%nmzd,2,input%jspins))
     569             : 
     570             :             !--->    generate name of file to hold the results of this iteration
     571           0 :             d1 = MOD(den%iter,10)
     572           0 :             d10 = MOD(INT((den%iter+0.5)/10),10)
     573           0 :             asciioffset = IACHAR('1')-1
     574             :             IF ( d10.GE.10 ) asciioffset = IACHAR('7')
     575           0 :             cdnfile = 'cdn'//ACHAR(d10+asciioffset)//ACHAR(d1+IACHAR('1')-1)
     576             : 
     577           0 :             iUnitTemp = 72
     578           0 :             OPEN (iUnitTemp,file=cdnfile,form='unformatted',status='unknown')
     579           0 :             REWIND iUnitTemp
     580             : 
     581             :             CALL loddop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
     582           0 :                         iUnit,iterTemp,frTemp,fpwTemp,fzTemp,fzxyTemp)
     583             :             CALL wrtdop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
     584           0 :                         iUnitTemp,iterTemp,frTemp,fpwTemp,fzTemp,fzxyTemp)
     585             : 
     586           0 :             CLOSE(iUnitTemp)
     587           0 :             REWIND iUnit
     588             : 
     589           0 :             DEALLOCATE (fzTemp, frTemp, fzxyTemp, fpwTemp)
     590           0 :             DEALLOCATE (atomsTemp%neq, atomsTemp%jri, atomsTemp%zatom, atomsTemp%ntypsy, sphharTemp%nlh)
     591           0 :             DEALLOCATE (atomsTemp%rmt, atomsTemp%dx)
     592             :          END IF
     593             : 
     594           0 :          IF ((inOrOutCDN.EQ.CDN_OUTPUT_DEN_const).AND.(archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const)) THEN
     595             : 
     596             :             ! Generate data in temp arrays and variables to be able to perform loddop call.
     597             :             ! loddop is called to move the file position to the output density position.
     598           0 :             inputTemp%jspins = input%jspins
     599           0 :             vacuumTemp%nmzxyd = vacuum%nmzxyd
     600           0 :             atomsTemp%jmtd = atoms%jmtd
     601           0 :             sphharTemp%nlhd = sphhar%nlhd
     602           0 :             vacuumTemp%nmzd = vacuum%nmzd
     603           0 :             atomsTemp%ntype = atoms%ntype
     604           0 :             ALLOCATE (sphharTemp%nlh(SIZE(sphhar%nlh)))
     605           0 :             sphharTemp%nlh(:) = sphhar%nlh(:)
     606           0 :             ALLOCATE (atomsTemp%ntypsy(SIZE(atoms%ntypsy)))
     607           0 :             atomsTemp%ntypsy(:) = atoms%ntypsy(:)
     608           0 :             ALLOCATE (atomsTemp%jri(SIZE(atoms%jri)))
     609           0 :             atomsTemp%jri(:) = atoms%jri(:)
     610           0 :             ALLOCATE (atomsTemp%neq(SIZE(atoms%neq)))
     611           0 :             atomsTemp%neq(:) = atoms%neq(:)
     612           0 :             starsTemp%ng3 = stars%ng3
     613           0 :             symTemp%invs = sym%invs
     614           0 :             inputTemp%film = input%film
     615           0 :             vacuumTemp%nvac = vacuum%nvac
     616           0 :             vacuumTemp%nmzxy = vacuum%nmzxy
     617           0 :             vacuumTemp%nmz = vacuum%nmz
     618           0 :             vacuumTemp%dvac = vacuum%dvac
     619           0 :             vacuumTemp%delz = vacuum%delz
     620           0 :             starsTemp%ng2 = stars%ng2
     621           0 :             symTemp%invs2 = sym%invs2
     622           0 :             ALLOCATE (fpwTemp(stars%ng3,input%jspins))
     623           0 :             ALLOCATE (fzxyTemp(vacuum%nmzxyd,stars%ng2-1,2,input%jspins))
     624           0 :             ALLOCATE (frTemp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins))
     625           0 :             ALLOCATE (fzTemp(vacuum%nmzd,2,input%jspins))
     626             : 
     627             :             CALL loddop(starsTemp,vacuumTemp,atomsTemp,sphharTemp,inputTemp,symTemp,&
     628           0 :                         iUnit,iterTemp,frTemp,fpwTemp,fzTemp,fzxyTemp)
     629             : 
     630           0 :             DEALLOCATE (fzTemp, frTemp, fzxyTemp, fpwTemp)
     631           0 :             DEALLOCATE (atomsTemp%neq, atomsTemp%jri, atomsTemp%ntypsy, sphharTemp%nlh)
     632             :          END IF
     633             : 
     634             :          ! Write the density
     635             :          CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
     636           0 :                      iUnit,den%iter+relCdnIndex,den%mt,den%pw,den%vacz,den%vacxy)
     637             : 
     638             :          ! Write additional data if l_noco
     639           0 :          IF (archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const) THEN
     640           0 :             WRITE (iUnit) (den%pw(k,3),k=1,stars%ng3)
     641           0 :             IF (input%film) THEN
     642           0 :                ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
     643           0 :                DO iVac = 1, vacuum%nvac
     644           0 :                   DO i = 1, vacuum%nmz
     645           0 :                      cdomvz(i,iVac) = CMPLX(den%vacz(i,iVac,3),den%vacz(i,iVac,4))
     646             :                   END DO
     647             :                END DO
     648           0 :                WRITE (iUnit) ((cdomvz(i,iVac),i=1,vacuum%nmz),iVac=1,vacuum%nvac)
     649           0 :                WRITE (iUnit) (((den%vacxy(i,j-1,iVac,3),i=1,vacuum%nmzxy),j=2,stars%ng2), iVac=1,vacuum%nvac)
     650           0 :                DEALLOCATE(cdomvz)
     651             :             END IF
     652             :          END IF
     653             : 
     654           0 :          CLOSE(iUnit)
     655             :       END IF
     656             : 
     657             :       !write density matrix to n_mmp_mat_out file
     658         223 :       IF((inOrOutCDN.EQ.CDN_INPUT_DEN_const).AND.(relCdnIndex.EQ.1).AND.&
     659             :          ((archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).OR.(archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const))) THEN
     660         189 :          IF(atoms%n_u.GT.0) THEN
     661          15 :             filename = 'n_mmp_mat'
     662          15 :             IF (mode.EQ.CDN_HDF5_MODE) THEN
     663          15 :                filename = 'n_mmp_mat_out'
     664             :             END IF
     665          15 :             IF(ANY(den%mmpMat(:,:,:,:).NE.0.0)) THEN
     666          14 :                IF ((mode.EQ.CDN_HDF5_MODE).AND..NOT.(input%ldauLinMix.AND.(input%ldauMixParam.EQ.0.0))) THEN
     667          14 :                   INQUIRE(file='n_mmp_mat',exist=l_exist)
     668          14 :                   IF(l_exist) THEN
     669           1 :                      CALL system('mv n_mmp_mat n_mmp_mat_old')
     670           1 :                      PRINT *,"n_mmp_mat moved to n_mmp_mat_old"
     671             :                   END IF
     672             :                END IF
     673          14 :                OPEN (69,file=TRIM(ADJUSTL(filename)),status='replace',form='formatted')
     674          14 :                WRITE (69,'(7f20.13)') den%mmpMat(:,:,:,:)
     675          14 :                CLOSE (69)
     676             :             END IF
     677             :          END IF
     678             :       END IF
     679             : 
     680         669 :    END SUBROUTINE writeDensity
     681             : 
     682          10 :    SUBROUTINE readPrevEFermi(eFermiPrev,l_error)
     683             : 
     684             :       REAL,    INTENT(OUT) :: eFermiPrev
     685             :       LOGICAL, INTENT(OUT) :: l_error
     686             : 
     687             :       INTEGER        :: mode
     688             : #ifdef CPP_HDF
     689             :       INTEGER(HID_T) :: fileID
     690             : #endif
     691             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex
     692             :       INTEGER        :: currentStructureIndex,currentStepfunctionIndex
     693             :       INTEGER        :: readDensityIndex, lastDensityIndex
     694             : 
     695             :       INTEGER           :: starsIndex, latharmsIndex, structureIndex
     696             :       INTEGER           :: stepfunctionIndex
     697             :       INTEGER           :: date, time, iter, jspins, previousDensityIndex
     698             :       REAL              :: fermiEnergy, distance
     699             :       LOGICAL           :: l_qfix, l_exist
     700             :       CHARACTER(LEN=30) :: archiveName
     701             : 
     702          10 :       CALL getIOMode(mode)
     703             : 
     704          10 :       eFermiPrev = 0.0
     705          10 :       l_error = .FALSE.
     706          10 :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     707             : #ifdef CPP_HDF
     708             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     709          10 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     710          10 :          WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
     711             :          CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_NOCO_IN_const,&
     712             :                                   iter, starsIndex, latharmsIndex, structureIndex, stepfunctionIndex,&
     713          10 :                                   previousDensityIndex, jspins, date, time, distance, fermiEnergy, l_qfix)
     714          10 :          eFermiPrev = fermiEnergy
     715          10 :          CALL closeCDNPOT_HDF(fileID)
     716             : #endif
     717             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
     718             :          STOP 'cdn.str not yet implemented!'
     719             :       ELSE
     720           0 :          l_error = .TRUE.
     721             :       END IF
     722             : 
     723          10 :    END SUBROUTINE
     724             : 
     725           2 :    SUBROUTINE readCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
     726             : 
     727             :       TYPE(t_atoms),INTENT(IN)     :: atoms
     728             :       TYPE(t_input),INTENT(IN)     :: input
     729             :       TYPE(t_dimension),INTENT(IN) :: DIMENSION
     730             : 
     731             :       REAL, INTENT(OUT) :: rhcs(:,:,:)!(atoms%jmtd,atoms%ntype,input%jspins)
     732             :       REAL, INTENT(OUT) :: tecs(:,:)!(atoms%ntype,input%jspins)
     733             :       REAL, INTENT(OUT) :: qints(:,:)!(atoms%ntype,input%jspins)
     734             : 
     735             :       INTEGER            :: mode, iUnit, iSpin, iAtom, i
     736             :       LOGICAL            :: l_exist
     737             :       CHARACTER(LEN=30)  :: filename
     738             : 
     739             : #ifdef CPP_HDF
     740             :       INTEGER(HID_T) :: fileID
     741             : #endif
     742             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex
     743             :       INTEGER        :: currentStructureIndex,currentStepfunctionIndex
     744             :       INTEGER        :: readDensityIndex, lastDensityIndex
     745             : 
     746           2 :       CALL getIOMode(mode)
     747             : 
     748             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     749             : #ifdef CPP_HDF
     750           2 :          l_exist = isCoreDensityPresentHDF()
     751           2 :          IF (l_exist) THEN
     752             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     753           2 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     754           2 :             CALL readCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints)
     755           2 :             CALL closeCDNPOT_HDF(fileID)
     756           2 :             RETURN
     757             :          ELSE
     758           0 :             WRITE(*,*) 'No core density is available in HDF5 format.'
     759           0 :             WRITE(*,*) 'Falling back to stream access file cdn.str.'
     760           0 :             mode = CDN_STREAM_MODE
     761             :          END IF
     762             : #endif
     763             :       END IF
     764             : 
     765             :       IF(mode.EQ.CDN_STREAM_MODE) THEN
     766           0 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
     767           0 :          IF (l_exist) THEN
     768             :             !load density from cdn.str and exit subroutine
     769             : 
     770             :             RETURN
     771             :          ELSE
     772           0 :             WRITE(*,*) 'cdn.str file not found.'
     773           0 :             WRITE(*,*) 'Falling back to direct access file cdnc.'
     774           0 :             mode = CDN_DIRECT_MODE
     775             :          END IF
     776             :       END IF
     777             : 
     778             :       IF (mode.EQ.CDN_DIRECT_MODE) THEN
     779           0 :          filename = 'cdnc'
     780           0 :          INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
     781           0 :          IF(.NOT.l_exist) THEN
     782           0 :             CALL juDFT_error("core charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readCoreDensity")
     783             :          END IF
     784           0 :          iUnit = 17
     785           0 :          OPEN (iUnit,file=TRIM(ADJUSTL(filename)),form='unformatted',status='unknown')
     786           0 :          DO iSpin = 1,input%jspins
     787           0 :             DO iAtom = 1,atoms%ntype
     788           0 :                READ (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
     789           0 :                READ (iUnit) tecs(iAtom,iSpin)
     790             :             END DO
     791           0 :             READ (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
     792             :          END DO
     793           0 :          CLOSE (iUnit)
     794             :       END IF
     795             : 
     796             :    END SUBROUTINE readCoreDensity
     797             : 
     798         162 :    SUBROUTINE writeCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
     799             : 
     800             :       TYPE(t_atoms),INTENT(IN)     :: atoms
     801             :       TYPE(t_input),INTENT(IN)     :: input
     802             :       TYPE(t_dimension),INTENT(IN) :: DIMENSION
     803             : 
     804             :       REAL, INTENT(IN) :: rhcs(:,:,:)!(atoms%jmtd,atoms%ntype,input%jspins)
     805             :       REAL, INTENT(IN) :: tecs(:,:)!(atoms%ntype,input%jspins)
     806             :       REAL, INTENT(IN) :: qints(:,:)!(atoms%ntype,input%jspins)
     807             : 
     808             :       INTEGER :: mode, iUnit, iSpin, iAtom, i
     809             : 
     810             : #ifdef CPP_HDF
     811             :       INTEGER(HID_T) :: fileID
     812             : #endif
     813             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex
     814             :       INTEGER        :: currentStructureIndex,currentStepfunctionIndex
     815             :       INTEGER        :: readDensityIndex, lastDensityIndex
     816             : 
     817         162 :       CALL getIOMode(mode)
     818             : 
     819             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     820             : #ifdef CPP_HDF
     821             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     822         162 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     823         162 :          CALL writeCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints)
     824         162 :          CALL closeCDNPOT_HDF(fileID)
     825             : #endif
     826             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
     827             :          ! Write core density to cdn.str file
     828             :          STOP 'CDN_STREAM_MODE not yet implemented!'
     829             :       ELSE
     830           0 :          iUnit = 17
     831           0 :          OPEN (iUnit,file='cdnc',form='unformatted',status='unknown')
     832           0 :          DO iSpin = 1,input%jspins
     833           0 :             DO iAtom = 1,atoms%ntype
     834           0 :                WRITE (iUnit) (rhcs(i,iAtom,iSpin),i=1,atoms%jri(iAtom))
     835           0 :                WRITE (iUnit) tecs(iAtom,iSpin)
     836             :             END DO
     837           0 :             WRITE (iUnit) (qints(iAtom,iSpin),iAtom=1,atoms%ntype)
     838             :          END DO
     839           0 :          CLOSE (iUnit)
     840             :       END IF
     841             : 
     842         162 :    END SUBROUTINE writeCoreDensity
     843             : 
     844          38 :    SUBROUTINE storeStructureIfNew(input,stars, atoms, cell, vacuum, oneD, sym,mpi,sphhar,noco)
     845             : 
     846             :       TYPE(t_input),INTENT(IN)   :: input
     847             :       TYPE(t_atoms), INTENT(IN)  :: atoms
     848             :       TYPE(t_cell), INTENT(IN)   :: cell
     849             :       TYPE(t_vacuum), INTENT(IN) :: vacuum
     850             :       TYPE(t_oneD),INTENT(IN)    :: oneD
     851             :       TYPE(t_sym),INTENT(IN)     :: sym
     852             :       TYPE(t_mpi),INTENT(IN)      :: mpi
     853             :       TYPE(t_sphhar),INTENT(IN)   :: sphhar
     854             :       TYPE(t_noco),INTENT(IN)     :: noco
     855             :       TYPE(t_stars),INTENT(IN)    :: stars
     856             : 
     857             :       TYPE(t_input)              :: inputTemp
     858          38 :       TYPE(t_atoms)              :: atomsTemp
     859             :       TYPE(t_cell)               :: cellTemp
     860          38 :       TYPE(t_vacuum)             :: vacuumTemp
     861             :       TYPE(t_oneD)               :: oneDTemp
     862          38 :       TYPE(t_sym)                :: symTemp
     863             : 
     864             :       INTEGER        :: mode
     865             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
     866             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
     867             :       LOGICAL        :: l_writeStructure, l_same
     868             : #ifdef CPP_HDF
     869             :       INTEGER(HID_T) :: fileID
     870             : #endif
     871             : 
     872          38 :       CALL getIOMode(mode)
     873             : 
     874             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     875             : #ifdef CPP_HDF
     876          38 :          l_writeStructure = .FALSE.
     877             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     878          38 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     879             : 
     880          38 :          IF (currentStructureIndex.EQ.0) THEN
     881          27 :             currentStructureIndex = currentStructureIndex + 1
     882          27 :             l_writeStructure = .TRUE.
     883             :          ELSE
     884          11 :             CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, currentStructureIndex)
     885          11 :             CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same)
     886          11 :             IF(.NOT.l_same) THEN
     887           0 :                currentStructureIndex = currentStructureIndex + 1
     888           0 :                l_writeStructure = .TRUE.
     889             :             END IF
     890             :          END IF
     891             : 
     892             :  
     893             :          IF (l_writeStructure) THEN
     894          27 :             CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, sym, currentStructureIndex,.TRUE.)
     895             :             CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     896          27 :                                     currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     897             :          END IF
     898             : 
     899          38 :          CALL closeCDNPOT_HDF(fileID)
     900             : #endif
     901             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
     902             :          ! Write stars to stars file
     903             :          STOP 'CDN_STREAM_MODE not yet implemented!'
     904             :       ELSE
     905             :          ! In direct access mode no structure information is written to any file.
     906             :       END IF
     907             : 
     908          76 :    END SUBROUTINE storeStructureIfNew
     909             : 
     910          38 :    SUBROUTINE transform_by_moving_atoms(mpi,stars,atoms,vacuum, cell, sym, sphhar,input,oned,noco)
     911             :      USE m_types
     912             :      USE m_qfix
     913             :      use m_fix_by_gaussian
     914             :      IMPLICIT NONE
     915             :      TYPE(t_mpi),INTENT(IN)      :: mpi
     916             :      TYPE(t_atoms),INTENT(IN)    :: atoms
     917             :      TYPE(t_sym),INTENT(IN)      :: sym
     918             :      TYPE(t_vacuum),INTENT(IN)   :: vacuum
     919             :      TYPE(t_sphhar),INTENT(IN)   :: sphhar
     920             :      TYPE(t_input),INTENT(IN)    :: input
     921             :      TYPE(t_oneD),INTENT(IN)     :: oneD
     922             :      TYPE(t_cell),INTENT(IN)     :: cell
     923             :      TYPE(t_noco),INTENT(IN)     :: noco
     924             :      TYPE(t_stars),INTENT(IN)    :: stars
     925             :      !Locals
     926             :      INTEGER :: archiveType
     927             :      LOGICAL :: l_qfix
     928             :      REAL    :: fermiEnergy,fix
     929          76 :      REAL    :: shifts(3,atoms%nat)
     930             : 
     931          38 :      TYPE(t_potden):: den
     932             : 
     933             :      TYPE(t_input)              :: inputTemp
     934          38 :      TYPE(t_atoms)              :: atomsTemp
     935             :      TYPE(t_cell)               :: cellTemp
     936          38 :      TYPE(t_vacuum)             :: vacuumTemp
     937             :      TYPE(t_oneD)               :: oneDTemp
     938          38 :      TYPE(t_sym)                :: symTemp
     939             :      
     940             : 
     941             :      INTEGER :: mode,currentStarsIndex,currentLatharmsIndex,currentStructureIndex
     942             :      INTEGER :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex,structureindex
     943             :      LOGICAL :: l_same,l_structure_by_shift
     944             : 
     945             : #ifdef CPP_HDF
     946             :       INTEGER(HID_T) :: fileID
     947             :       character(len=50) :: archivename
     948             : #endif
     949             : #ifdef CPP_MPI
     950             :       INCLUDE 'mpif.h'
     951             :       integer :: ierr
     952             : #endif
     953          38 :       l_same=.TRUE.;l_structure_by_shift=.TRUE.
     954             :           
     955          38 :       CALL getIOMode(mode)
     956             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
     957             : #ifdef CPP_HDF
     958          38 :          IF (mpi%irank==0) THEN
     959             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
     960          19 :                  currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
     961          19 :             IF (currentStructureIndex>0.AND.lastdensityindex>0) THEN
     962             :                !Determine structure of last density
     963           5 :                WRITE(archivename,"(a,i0)") "cdn-",lastdensityindex
     964           5 :                CALL peekDensityEntryHDF(fileID, archivename, DENSITY_TYPE_IN_const, structureIndex=structureIndex)
     965             :                !Read that structure
     966           5 :                CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, StructureIndex)
     967           5 :                CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same,l_structure_by_shift)
     968             :             ENDIF
     969          19 :             CALL closeCDNPOT_HDF(fileID)
     970             :          ENDIF
     971             : #ifdef CPP_MPI
     972          38 :          CALL mpi_bcast(l_same,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
     973          38 :          CALL mpi_bcast(l_structure_by_shift,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
     974             : #endif
     975          38 :          IF (l_same.OR..NOT.l_structure_by_shift) RETURN ! nothing to do
     976             :          
     977           0 :          IF (mpi%irank==0) THEN
     978           0 :             WRITE(6,*) "Atomic movement detected, trying to adjust charge density" 
     979             :             
     980             :             !Calculate shifts
     981           0 :             shifts=atomsTemp%taual-atoms%taual
     982             :             
     983             :             !Determine type of charge
     984           0 :             archiveType = MERGE(CDN_ARCHIVE_TYPE_NOCO_const,CDN_ARCHIVE_TYPE_CDN1_const,noco%l_noco)
     985             :             !read the current density
     986           0 :             CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
     987             :             CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
     988           0 :                  0,fermiEnergy,l_qfix,den)
     989             :          ENDIF
     990             :          !Now fix the density
     991           0 :          SELECT CASE(input%qfix)
     992             :          CASE (0,1) !just qfix the density
     993           0 :             IF (mpi%irank==0) WRITE(6,*) "Using qfix to adjust density"
     994           0 :             if (mpi%irank==0) CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
     995           0 :                  den,noco%l_noco,mpi%isize==1,force_fix=.TRUE.,fix=fix)
     996             :          CASE(2,3)
     997           0 :             if (mpi%irank==0) CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,&
     998           0 :                  den,noco%l_noco,mpi%isize==1,force_fix=.TRUE.,fix=fix,fix_pw_only=.true.)
     999             :          CASE(4,5)
    1000           0 :             if (mpi%irank==0) CALL fix_by_gaussian(shifts,atoms,stars,mpi,sym,vacuum,sphhar,input,oned,cell,noco,den)
    1001             :          CASE default
    1002           0 :             CALL judft_error("Wrong choice of qfix in input")
    1003             :          END SELECT
    1004             :          !Now write the density to file
    1005           0 :          IF (mpi%irank==0) CALL writedensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
    1006           0 :               0,-1.0,fermiEnergy,l_qfix,den)
    1007             :          
    1008             : #endif
    1009             :       END IF
    1010         114 :     END SUBROUTINE transform_by_moving_atoms
    1011             : 
    1012          27 :    SUBROUTINE writeStars(stars,l_xcExtended,l_ExtData)
    1013             : 
    1014             :       TYPE(t_stars),INTENT(IN)   :: stars
    1015             :       LOGICAL, INTENT(IN)        :: l_xcExtended, l_ExtData
    1016             : 
    1017             :       INTEGER        :: mode, ngz, izmin, izmax
    1018             : 
    1019             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
    1020             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
    1021             : #ifdef CPP_HDF
    1022             :       INTEGER(HID_T) :: fileID
    1023             : #endif
    1024             : 
    1025          54 :       INTEGER :: igz(stars%ng3)
    1026             : 
    1027          27 :       ngz = 0
    1028          27 :       izmin = 0
    1029          27 :       izmax = 0
    1030       92924 :       igz = 0
    1031             : 
    1032          27 :       CALL getIOMode(mode)
    1033             : 
    1034             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1035             : #ifdef CPP_HDF
    1036             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1037          27 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1038             : 
    1039          27 :          currentStarsIndex = currentStarsIndex + 1
    1040          27 :          CALL writeStarsHDF(fileID, currentStarsIndex, currentStructureIndex, stars,.TRUE.)
    1041             : 
    1042             :          CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1043          27 :                                  currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1044             : 
    1045          27 :          CALL closeCDNPOT_HDF(fileID)
    1046             : #endif
    1047             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
    1048             :          ! Write stars to stars file
    1049             :          STOP 'CDN_STREAM_MODE not yet implemented!'
    1050             :       ELSE
    1051           0 :          OPEN (51,file='stars',form='unformatted',status='unknown')
    1052           0 :          WRITE (51) stars%gmax,stars%ng3,stars%ng2,ngz,izmin,izmax,stars%mx1,stars%mx2,stars%mx3
    1053           0 :          IF(l_ExtData) THEN
    1054           0 :             IF (.NOT.l_xcExtended) THEN
    1055           0 :                WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
    1056           0 :                           stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1057           0 :                           stars%igfft2,stars%pgfft2
    1058             :             ELSE
    1059           0 :                WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
    1060           0 :                           stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1061           0 :                           stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
    1062             :             END IF
    1063             :          ELSE
    1064           0 :             IF (.NOT.l_xcExtended) THEN
    1065           0 :                WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
    1066           0 :                           stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1067           0 :                           stars%igfft2,stars%pgfft2
    1068             :             ELSE
    1069           0 :                WRITE (51) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
    1070           0 :                           stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1071           0 :                           stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
    1072             :             END IF
    1073             :          END IF
    1074           0 :          CLOSE (51)
    1075             :       END IF
    1076          27 :    END SUBROUTINE writeStars
    1077             : 
    1078          38 :    SUBROUTINE readStars(stars,l_xcExtended,l_ExtData,l_error)
    1079             : 
    1080             :       TYPE(t_stars),INTENT(INOUT) :: stars
    1081             :       LOGICAL, INTENT(IN)         :: l_xcExtended,l_ExtData
    1082             :       LOGICAL, INTENT(OUT)        :: l_error
    1083             : 
    1084             : 
    1085          38 :       TYPE(t_stars)               :: starsTemp
    1086             :       INTEGER                     :: mode, ioStatus, ngz,izmin,izmax
    1087             :       LOGICAL                     :: l_exist, l_same
    1088             : 
    1089             :       INTEGER        :: structureIndexTemp
    1090             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
    1091             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
    1092             : #ifdef CPP_HDF
    1093             :       INTEGER(HID_T) :: fileID
    1094             : #endif
    1095             : 
    1096          76 :       INTEGER :: igz(stars%ng3)
    1097             : 
    1098          38 :       l_error = .FALSE.
    1099          38 :       ngz = 0
    1100          38 :       izmin = 0
    1101          38 :       izmax = 0
    1102      106903 :       igz = 0
    1103             : 
    1104          38 :       CALL getIOMode(mode)
    1105             : 
    1106             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1107          38 :          INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
    1108          38 :          IF (l_exist) THEN
    1109             : #ifdef CPP_HDF
    1110             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1111          38 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1112          38 :             IF (currentStarsIndex.LT.1) THEN
    1113             :                mode = CDN_DIRECT_MODE ! (no stars entry found in cdn.hdf file)
    1114             :             ELSE
    1115          11 :                CALL peekStarsHDF(fileID, currentStarsIndex, structureIndexTemp)
    1116          11 :                l_same = structureIndexTemp.EQ.currentStructureIndex
    1117          11 :                IF(l_same) THEN
    1118          11 :                   CALL readStarsHDF(fileID, currentStarsIndex, starsTemp)
    1119          11 :                   CALL compareStars(stars, starsTemp, l_same)
    1120             :                END IF
    1121          11 :                IF(l_same) THEN
    1122          11 :                   CALL readStarsHDF(fileID, currentStarsIndex, stars)
    1123             :                ELSE
    1124             :                  mode = CDN_DIRECT_MODE ! (no adequate stars entry found in cdn.hdf file)
    1125             :                END IF
    1126             :             END IF
    1127          38 :             CALL closeCDNPOT_HDF(fileID)
    1128             : #endif
    1129             :          END IF
    1130          38 :          IF(.NOT.l_exist) THEN
    1131             :             mode = CDN_STREAM_MODE
    1132             :          END IF
    1133             :       END IF
    1134             : 
    1135          38 :       IF(mode.EQ.CDN_STREAM_MODE) THEN
    1136           0 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
    1137           0 :          IF (l_exist) THEN
    1138           0 :             STOP 'cdn.str code path not yet implemented!'
    1139             :          END IF
    1140             :          IF (.NOT.l_exist) THEN
    1141             :             mode = CDN_DIRECT_MODE
    1142             :          END IF
    1143             :       END IF
    1144             : 
    1145          38 :       IF (mode.EQ.CDN_DIRECT_MODE) THEN
    1146          27 :          INQUIRE(FILE='stars',EXIST=l_exist)
    1147          27 :          IF(.NOT.l_exist) THEN
    1148          27 :             l_error = .TRUE.
    1149          27 :             RETURN
    1150             :          END IF
    1151           0 :          OPEN (51,file='stars',form='unformatted',status='unknown')
    1152             : 
    1153           0 :          READ (51,IOSTAT=ioStatus) stars%gmax,stars%ng3,stars%ng2,ngz,izmin,izmax,stars%mx1,stars%mx2,stars%mx3
    1154           0 :          IF (ioStatus.NE.0) THEN
    1155           0 :             l_error = .TRUE.
    1156           0 :             RETURN
    1157             :          END IF
    1158             : 
    1159           0 :          IF (l_ExtData) THEN
    1160           0 :             IF (.NOT.l_xcExtended) THEN
    1161           0 :                READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
    1162           0 :                                          stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1163           0 :                                          stars%igfft2,stars%pgfft2
    1164           0 :                stars%ft2_gfx = 0.0
    1165           0 :                stars%ft2_gfy = 0.0
    1166             :             ELSE
    1167           0 :                READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%phi2,stars%kv3,stars%kv2,&
    1168           0 :                                          stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1169           0 :                                          stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
    1170             :             END IF
    1171             :          ELSE
    1172           0 :             IF (.NOT.l_xcExtended) THEN
    1173           0 :                READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
    1174           0 :                                          stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1175           0 :                                          stars%igfft2,stars%pgfft2
    1176           0 :                stars%ft2_gfx = 0.0
    1177           0 :                stars%ft2_gfy = 0.0
    1178             :             ELSE
    1179           0 :                READ (51,IOSTAT=ioStatus) stars%nstr,stars%nstr2,stars%rgphs,stars%sk3,stars%sk2,stars%kv3,stars%kv2,&
    1180           0 :                                          stars%ig,stars%ig2,igz,stars%kimax,stars%igfft,stars%pgfft,stars%kimax2,&
    1181           0 :                                          stars%igfft2,stars%pgfft2,stars%ft2_gfx,stars%ft2_gfy
    1182             :             END IF
    1183             :          END IF
    1184             : 
    1185           0 :          IF (ioStatus.NE.0) THEN
    1186           0 :             l_error = .TRUE.
    1187           0 :             RETURN
    1188             :          END IF
    1189             : 
    1190           0 :          CLOSE (51)
    1191             :       END IF
    1192             : 
    1193          38 :    END SUBROUTINE readStars
    1194             : 
    1195          27 :    SUBROUTINE writeStepfunction(stars)
    1196             : 
    1197             :       TYPE(t_stars),INTENT(IN) :: stars
    1198             : 
    1199             :       INTEGER                  :: mode, ifftd, i
    1200             : 
    1201             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
    1202             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
    1203             : #ifdef CPP_HDF
    1204             :       INTEGER(HID_T) :: fileID
    1205             : #endif
    1206             : 
    1207          27 :       ifftd=size(stars%ufft)
    1208             : 
    1209          27 :       CALL getIOMode(mode)
    1210             : 
    1211             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1212             : #ifdef CPP_HDF
    1213             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1214          27 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1215             : 
    1216          27 :          currentStepfunctionIndex = currentStepfunctionIndex + 1
    1217          27 :          CALL writeStepfunctionHDF(fileID, currentStepfunctionIndex, currentStarsIndex, currentStructureIndex, stars,.TRUE.)
    1218             :          CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1219          27 :                                  currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1220             : 
    1221          27 :          CALL closeCDNPOT_HDF(fileID)
    1222             : #endif
    1223             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
    1224             :          ! Write stars to stars file
    1225             :          STOP 'CDN_STREAM_MODE not yet implemented!'
    1226             :       ELSE
    1227           0 :          OPEN (14,file='wkf2',form='unformatted',status='unknown')
    1228             : 
    1229           0 :          WRITE (14) stars%ng3,ifftd
    1230           0 :          WRITE (14) (stars%ustep(i),i=1,stars%ng3)
    1231           0 :          WRITE (14) (stars%ufft(i),i=0,ifftd-1)
    1232             : 
    1233           0 :          CLOSE (14)
    1234             :       END IF
    1235             : 
    1236          27 :    END SUBROUTINE writeStepfunction
    1237             : 
    1238          38 :    SUBROUTINE readStepfunction(stars, atoms, cell, vacuum, l_error)
    1239             : 
    1240             :       TYPE(t_stars),INTENT(INOUT)   :: stars
    1241             :       TYPE(t_atoms), INTENT(IN)     :: atoms
    1242             :       TYPE(t_cell), INTENT(IN)      :: cell
    1243             :       TYPE(t_vacuum), INTENT(IN)    :: vacuum
    1244             :       LOGICAL, INTENT(OUT)          :: l_error
    1245             : 
    1246             :       TYPE(t_stars)                 :: starsTemp
    1247             :       TYPE(t_input)                 :: inputTemp
    1248             :       TYPE(t_atoms)                 :: atomsTemp
    1249             :       TYPE(t_cell)                  :: cellTemp
    1250             :       TYPE(t_vacuum)                :: vacuumTemp
    1251             :       TYPE(t_oneD)                  :: oneDTemp
    1252             : 
    1253             :       INTEGER        :: mode
    1254             :       INTEGER        :: ifftd, ng3Temp, ifftdTemp, ioStatus, i, starsIndexTemp
    1255             :       INTEGER        :: structureIndexTemp
    1256             :       LOGICAL        :: l_exist, l_same, l_sameTemp
    1257             : 
    1258             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
    1259             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
    1260             : #ifdef CPP_HDF
    1261             :       INTEGER(HID_T) :: fileID
    1262             : #endif
    1263             : 
    1264          38 :       l_error = .FALSE.
    1265          38 :       ioStatus = 0
    1266          38 :       ifftd = 27*stars%mx1*stars%mx2*stars%mx3
    1267             : 
    1268          38 :       CALL getIOMode(mode)
    1269             : 
    1270             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1271          38 :          INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
    1272          38 :          IF (l_exist) THEN
    1273             : #ifdef CPP_HDF
    1274             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1275          38 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1276          38 :             IF(currentStepfunctionIndex.GT.0) THEN
    1277          11 :                CALL peekStepfunctionHDF(fileID, currentStepfunctionIndex, starsIndexTemp, structureIndexTemp)
    1278          11 :                IF((starsIndexTemp.EQ.currentStarsIndex).AND.(structureIndexTemp.EQ.currentStructureIndex)) THEN
    1279          11 :                   CALL readStepfunctionHDF(fileID, currentStepfunctionIndex, stars)
    1280             :                ELSE
    1281             :                   mode = CDN_STREAM_MODE ! No adequate stepfunction entry found. Fall back to other IO modes.
    1282             :                END IF
    1283             :             ELSE
    1284             :                mode = CDN_STREAM_MODE ! No adequate stepfunction entry found. Fall back to other IO modes.
    1285             :             END IF
    1286          38 :             CALL closeCDNPOT_HDF(fileID)
    1287             : #endif
    1288             :          END IF
    1289          38 :          IF(.NOT.l_exist) THEN
    1290             :             mode = CDN_STREAM_MODE
    1291             :          END IF
    1292             :       END IF
    1293             : 
    1294          38 :       IF(mode.EQ.CDN_STREAM_MODE) THEN
    1295          27 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
    1296          27 :          IF (l_exist) THEN
    1297           0 :             STOP 'cdn.str code path not yet implemented!'
    1298             :          END IF
    1299             :          IF (.NOT.l_exist) THEN
    1300             :             mode = CDN_DIRECT_MODE
    1301             :          END IF
    1302             :       END IF
    1303             : 
    1304          11 :       IF (mode.EQ.CDN_DIRECT_MODE) THEN
    1305          27 :          INQUIRE(FILE='wkf2',EXIST=l_exist)
    1306          27 :          IF(.NOT.l_exist) THEN
    1307          27 :             l_error = .TRUE.
    1308          54 :             RETURN
    1309             :          END IF
    1310           0 :          OPEN (14,file='wkf2',form='unformatted',status='unknown')
    1311           0 :          ng3temp=0;ifftdTemp=0
    1312           0 :          READ (14,IOSTAT=ioStatus) ng3Temp, ifftdTemp
    1313           0 :          IF (ng3Temp.NE.stars%ng3) ioStatus = 1
    1314           0 :          IF (ifftdTemp.NE.ifftd) ioStatus = 1
    1315           0 :          IF (ioStatus.NE.0) THEN
    1316           0 :             l_error = .TRUE.
    1317           0 :             CLOSE (14)
    1318           0 :             RETURN
    1319             :          END IF
    1320           0 :          READ (14) (stars%ustep(i),i=1,stars%ng3)
    1321           0 :          READ (14) (stars%ufft(i),i=0,ifftd-1)
    1322             : 
    1323           0 :          CLOSE (14)
    1324             :       END IF
    1325             : 
    1326             :    END SUBROUTINE readStepfunction
    1327             : 
    1328          38 :    SUBROUTINE setStartingDensity(l_noco)
    1329             : 
    1330             :       LOGICAL,INTENT(IN) :: l_noco
    1331             : 
    1332             : #ifdef CPP_HDF
    1333             :       INTEGER(HID_T)    :: fileID
    1334             : #endif
    1335             :       INTEGER           :: currentStarsIndex,currentLatharmsIndex
    1336             :       INTEGER           :: currentStructureIndex,currentStepfunctionIndex
    1337             :       INTEGER           :: readDensityIndex, lastDensityIndex
    1338             :       INTEGER           :: sdIndex, ioStatus, mode
    1339             :       INTEGER           :: densityType
    1340             :       CHARACTER(LEN=1000) :: numberString
    1341             :       CHARACTER(LEN=30) :: archiveName
    1342             :       LOGICAL           :: l_exist
    1343             : 
    1344          38 :       IF (.NOT.juDFT_was_argument("-sd")) THEN
    1345          38 :          RETURN
    1346             :       END IF
    1347           0 :       numberString = juDFT_string_for_argument("-sd")
    1348           0 :       IF (TRIM(ADJUSTL(numberString)).EQ.'') THEN
    1349           0 :          CALL juDFT_error("No number for starting density set in command line.",calledby ="setStartingDensity")
    1350             :       END IF
    1351           0 :       ioStatus = 0
    1352           0 :       READ(numberString,'(i8)',iostat=ioStatus) sdIndex
    1353           0 :       IF(ioStatus.NE.0) THEN
    1354           0 :          CALL juDFT_error("Could not convert starting density index string to number.",calledby ="setStartingDensity")
    1355             :       END IF
    1356             : 
    1357           0 :       WRITE(*,'(a,i0,a)') 'Using density ', sdIndex, ' as starting density.'
    1358             : 
    1359           0 :       CALL getIOMode(mode)
    1360             : 
    1361             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1362             : #ifdef CPP_HDF
    1363             :          CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1364           0 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1365           0 :          densityType = DENSITY_TYPE_IN_const
    1366           0 :          IF(l_noco) THEN
    1367           0 :             densityType = DENSITY_TYPE_NOCO_IN_const
    1368             :          END IF
    1369           0 :          archiveName = ''
    1370           0 :          WRITE(archiveName,'(a,i0)') '/cdn-', sdIndex
    1371           0 :          l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType)
    1372           0 :          IF(.NOT.l_exist) THEN
    1373           0 :             WRITE(*,*) 'archiveName: ', TRIM(ADJUSTL(archiveName))
    1374           0 :             CALL juDFT_error("For selected starting density index no in-density is present.",calledby ="setStartingDensity")
    1375             :          END IF
    1376             :          CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1377           0 :                                  currentStepfunctionIndex,sdIndex,lastDensityIndex)
    1378           0 :          CALL closeCDNPOT_HDF(fileID)
    1379             : #endif
    1380             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
    1381             :          STOP 'CDN_STREAM_MODE not yet implemented!'
    1382             :       ELSE      
    1383           0 :          WRITE(*,*) 'Explicit setting of starting density in direct access mode'
    1384           0 :          WRITE(*,*) 'not implemented.'
    1385           0 :          WRITE(*,*) ''
    1386           0 :          WRITE(*,*) 'Ignoring -sd command line argument.'
    1387             :       END IF
    1388             : 
    1389             :    END SUBROUTINE setStartingDensity
    1390             : 
    1391          38 :    SUBROUTINE deleteDensities()
    1392             : 
    1393             : #ifdef CPP_HDF
    1394             :       INTEGER(HID_T)    :: fileID
    1395             : #endif
    1396             :       INTEGER           :: currentStarsIndex,currentLatharmsIndex
    1397             :       INTEGER           :: currentStructureIndex,currentStepfunctionIndex
    1398             :       INTEGER           :: readDensityIndex, lastDensityIndex
    1399             :       INTEGER           :: ioStatus, mode, i
    1400             :       INTEGER           :: startNumber, endNumber, separatorIndex
    1401             :       CHARACTER(LEN=1000) :: ddString
    1402             :       CHARACTER(LEN=30) :: archiveName
    1403             :       LOGICAL           :: l_exist, l_deleted
    1404             : 
    1405          38 :       IF (.NOT.juDFT_was_argument("-delden")) THEN
    1406          38 :          RETURN
    1407             :       END IF
    1408           0 :       ddString = juDFT_string_for_argument("-delden")
    1409           0 :       IF (TRIM(ADJUSTL(ddString)).EQ.'') THEN
    1410           0 :          CALL juDFT_error("Densities to be deleted not specified.",calledby ="deleteDensities")
    1411             :       END IF
    1412             : 
    1413           0 :       separatorIndex = -1
    1414           0 :       startNumber = -1
    1415           0 :       endNumber = -1
    1416           0 :       IF (TRIM(ADJUSTL(ddString)).EQ.'allbutlast') THEN
    1417           0 :          CALL getIOMode(mode)
    1418             : 
    1419             :          IF(mode.EQ.CDN_HDF5_MODE) THEN
    1420           0 :             INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
    1421           0 :             IF (l_exist) THEN
    1422             : #ifdef CPP_HDF
    1423             :                CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1424           0 :                                 currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1425           0 :                CALL closeCDNPOT_HDF(fileID)
    1426           0 :                startNumber = 1
    1427           0 :                endNumber = lastDensityIndex - 1
    1428             : #endif
    1429             :             END IF
    1430             :          END IF
    1431             :       ELSE
    1432           0 :          DO i = 1, LEN(TRIM(ADJUSTL(ddString)))
    1433           0 :             IF(VERIFY(ddString(i:i),'1234567890').NE.0) THEN
    1434           0 :                IF ((ddString(i:i).EQ.'-').AND.(separatorIndex.EQ.-1)) THEN
    1435             :                   separatorIndex = i
    1436             :                ELSE
    1437           0 :                   CALL juDFT_error("density deletion string format error",calledby ="deleteDensities")
    1438             :                END IF
    1439             :             END IF
    1440             :          END DO
    1441             : 
    1442           0 :          IF(separatorIndex.NE.-1) THEN
    1443           0 :             READ(ddString(1:separatorIndex-1),'(i7)') startNumber
    1444           0 :             READ(ddString(separatorIndex+1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') endNumber
    1445             :          ELSE
    1446           0 :             READ(ddString(1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') startNumber
    1447           0 :             READ(ddString(1:LEN(TRIM(ADJUSTL(ddString)))),'(i7)') endNumber
    1448             :          END IF
    1449             :       END IF
    1450             : 
    1451           0 :       CALL getIOMode(mode)
    1452             : 
    1453             :       IF(mode.EQ.CDN_HDF5_MODE) THEN
    1454           0 :          INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
    1455           0 :          IF (l_exist) THEN
    1456             : #ifdef CPP_HDF
    1457             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1458           0 :                           currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1459             : 
    1460           0 :             DO i = startNumber, endNumber
    1461           0 :                archiveName = ''
    1462           0 :                WRITE(archiveName,'(a,i0)') '/cdn-', i
    1463             : 
    1464           0 :                l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_UNDEFINED_const)
    1465           0 :                IF(.NOT.l_exist) THEN
    1466             :                   CYCLE
    1467             :                END IF
    1468             :                
    1469           0 :                l_deleted = deleteDensityEntryHDF(fileID,archiveName)
    1470           0 :                IF (l_deleted) THEN
    1471           0 :                   WRITE(*,*) 'deleted density entry ', TRIM(ADJUSTL(archiveName))
    1472             :                END IF
    1473             :             END DO
    1474             : 
    1475           0 :             CALL closeCDNPOT_HDF(fileID)
    1476             : #endif
    1477           0 :             WRITE(*,*) 'Please note:'
    1478           0 :             WRITE(*,*) 'The deletion of the densities does not free the associated disk space.'
    1479           0 :             WRITE(*,*) 'To do this you have to repack the cdn.hdf file.'
    1480           0 :             WRITE(*,*) 'It can be done by using the tool h5repack, e.g., by invoking'
    1481           0 :             WRITE(*,*) 'h5repack -i cdn.hdf -o cdn-packed.hdf'
    1482           0 :             WRITE(*,*) 'mv cdn-packed.hdf cdn.hdf'
    1483             :          ELSE
    1484           0 :             WRITE(*,*) "No cdn.hdf file found. No density entry deleted."
    1485             :          END IF
    1486             :       ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
    1487             :          STOP 'CDN_STREAM_MODE not yet implemented!'
    1488             :       ELSE      
    1489           0 :          WRITE(*,*) 'Explicit deletion of densities in direct access mode'
    1490           0 :          WRITE(*,*) 'not implemented.'
    1491           0 :          WRITE(*,*) ''
    1492           0 :          WRITE(*,*) 'Ignoring -delden command line argument.'
    1493             :       END IF
    1494             : 
    1495           0 :       CALL juDFT_end("Selected densities deleted.")
    1496             :       
    1497             :    END SUBROUTINE deleteDensities
    1498             : 
    1499           0 :    SUBROUTINE getIOMode(mode)
    1500             :       INTEGER, INTENT(OUT) :: mode
    1501             : 
    1502         671 :       mode = CDN_DIRECT_MODE
    1503             :       !IF (juDFT_was_argument("-stream_cdn")) THEN
    1504             :       !   mode=CDN_STREAM_MODE
    1505             :       !END IF
    1506         671 :       IF (.NOT.juDFT_was_argument("-no_cdn_hdf")) THEN !juDFT_was_argument("-hdf_cdn")) THEN
    1507             : #ifdef CPP_HDF
    1508         671 :          mode=CDN_HDF5_MODE
    1509             : #else
    1510             : !         WRITE(*,*) 'Code not compiled with HDF5 support.'
    1511             : !         WRITE(*,*) 'Falling back to direct access.'
    1512             : #endif
    1513             :       END IF
    1514           0 :    END SUBROUTINE getIOMode
    1515             : 
    1516          33 :    LOGICAL FUNCTION isDensityFilePresent(archiveType)
    1517             : 
    1518             :       INTEGER, INTENT(IN) :: archiveType
    1519             : 
    1520             :       LOGICAL             :: l_exist
    1521             :       INTEGER             :: mode
    1522             : 
    1523             :       INTEGER        :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
    1524             :       INTEGER        :: currentStepfunctionIndex,readDensityIndex,lastDensityIndex
    1525             : #ifdef CPP_HDF
    1526             :       INTEGER(HID_T) :: fileID
    1527             : #endif
    1528             : 
    1529          33 :       CALL getIOMode(mode)
    1530             : 
    1531             :       IF (mode.EQ.CDN_HDF5_MODE) THEN
    1532          33 :          INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
    1533          33 :          IF(l_exist) THEN
    1534             : #ifdef CPP_HDF
    1535             :             CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
    1536          33 :                              currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
    1537          33 :             CALL closeCDNPOT_HDF(fileID)
    1538          33 :             IF(readDensityIndex.GT.0) THEN
    1539             :                isDensityFilePresent = .TRUE.
    1540             :                RETURN
    1541             :             END IF
    1542             : #endif
    1543             :          END IF
    1544             :       END IF
    1545             : 
    1546          25 :       IF ((mode.EQ.CDN_STREAM_MODE).OR.(mode.EQ.CDN_HDF5_MODE)) THEN
    1547          25 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
    1548          25 :          IF(l_exist) THEN
    1549             :             isDensityFilePresent = l_exist
    1550             :             RETURN
    1551             :          END IF
    1552             :       END IF
    1553             : 
    1554             :       !cdn1 or rhomat_inp should be enough for any mode...
    1555          25 :       INQUIRE(FILE='cdn1',EXIST=l_exist)
    1556          25 :       IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const) THEN
    1557          19 :          isDensityFilePresent = l_exist
    1558          19 :          RETURN
    1559             :       END IF
    1560           6 :       IF (archiveType.NE.CDN_ARCHIVE_TYPE_NOCO_const) THEN
    1561           0 :          CALL juDFT_error("Illegal archive type selected.",calledby ="isDensityFilePresent")
    1562             :       END IF
    1563           6 :       IF (l_exist) THEN
    1564             :          isDensityFilePresent = l_exist
    1565             :          RETURN
    1566             :       END IF
    1567           6 :       INQUIRE(FILE='rhomat_inp',EXIST=l_exist)
    1568           6 :       isDensityFilePresent = l_exist
    1569           6 :    END FUNCTION isDensityFilePresent
    1570             : 
    1571           0 :    LOGICAL FUNCTION isCoreDensityPresent()
    1572             : 
    1573             :       LOGICAL             :: l_exist
    1574             :       INTEGER             :: mode
    1575             : 
    1576           0 :       CALL getIOMode(mode)
    1577             : 
    1578             :       IF (mode.EQ.CDN_HDF5_MODE) THEN
    1579             : #ifdef CPP_HDF
    1580           0 :          l_exist = isCoreDensityPresentHDF()
    1581           0 :          IF(l_exist) THEN
    1582             :             isCoreDensityPresent = l_exist
    1583             :             RETURN
    1584             :          END IF
    1585             : #endif
    1586             :       END IF
    1587             : 
    1588           0 :       IF ((mode.EQ.CDN_STREAM_MODE).OR.(mode.EQ.CDN_HDF5_MODE)) THEN
    1589           0 :          INQUIRE(FILE='cdn.str',EXIST=l_exist)
    1590           0 :          IF(l_exist) THEN
    1591           0 :             STOP 'Not yet implemented!'
    1592             :             RETURN
    1593             :          END IF
    1594             :       END IF
    1595             : 
    1596             :       !cdnc should be enough for any mode...
    1597           0 :       INQUIRE(FILE='cdnc',EXIST=l_exist)
    1598           0 :       IF (l_exist) THEN
    1599           0 :          isCoreDensityPresent = l_exist
    1600           0 :          RETURN
    1601             :       END IF
    1602             :       isCoreDensityPresent = .FALSE.
    1603             :    END FUNCTION isCoreDensityPresent
    1604             : 
    1605             : END MODULE m_cdn_io

Generated by: LCOV version 1.13