diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index 1780f478f630a1719e5b187f0849c0866764bcf5..9624f1490022f062265e8bad73e03679e86fa196 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -16,7 +16,7 @@ io/inpnoco.F90 io/loddop.f90 io/cdnpot_io_hdf.F90 io/cdn_io.F90 -io/pot_io.f90 +io/pot_io.F90 io/rw_inp.f90 io/rw_noco.f90 io/r_inpXML.F90 diff --git a/io/cdn_io.F90 b/io/cdn_io.F90 index 401ee747a10de6f4e31701da64b8a5a346d9d450..75966ec474f283eb0ac65cb4754e66535a284fdb 100644 --- a/io/cdn_io.F90 +++ b/io/cdn_io.F90 @@ -125,7 +125,7 @@ MODULE m_cdn_io CALL juDFT_error("Invalid inOrOutCDN selected.",calledby ="readDensity") END SELECT l_exist = isDensityEntryPresentHDF(fileID,archiveName,densityType) - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) END IF IF (l_exist) THEN @@ -135,7 +135,7 @@ MODULE m_cdn_io CALL readDensityHDF(fileID, archiveName, densityType,& fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy) - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) RETURN ELSE WRITE(*,*) 'cdn.hdf file or relevant density entry not found.' @@ -350,11 +350,11 @@ MODULE m_cdn_io fermiEnergy,l_qfix,iter+relCdnIndex,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy) IF(l_storeIndices) THEN - CALL writeHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,& - currentStructureIndex,readDensityIndex,lastDensityIndex) + CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,& + currentStructureIndex,readDensityIndex,lastDensityIndex) END IF - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) #endif ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN ! Write density to cdn.str file @@ -536,7 +536,7 @@ MODULE m_cdn_io ELSE l_error = .TRUE. END IF - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) #endif ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN STOP 'cdn.str not yet implemented!' @@ -576,7 +576,7 @@ MODULE m_cdn_io CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,& readDensityIndex,lastDensityIndex) CALL readCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints) - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) RETURN ELSE WRITE(*,*) 'No core density is available in HDF5 format.' @@ -645,7 +645,7 @@ MODULE m_cdn_io CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,& readDensityIndex,lastDensityIndex) CALL writeCoreDensityHDF(fileID,input,atoms,dimension,rhcs,tecs,qints) - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) #endif ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN ! Write core density to cdn.str file @@ -713,9 +713,9 @@ MODULE m_cdn_io WRITE(*,*) 'archiveName: ', TRIM(ADJUSTL(archiveName)) CALL juDFT_error("For selected starting density index no in-density is present.",calledby ="setStartingDensity") END IF - CALL writeHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,& - currentStructureIndex,sdIndex,lastDensityIndex) - CALL closeCDN_HDF(fileID) + CALL writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,& + currentStructureIndex,sdIndex,lastDensityIndex) + CALL closeCDNPOT_HDF(fileID) #endif ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN STOP 'CDN_STREAM_MODE not yet implemented!' diff --git a/io/cdnpot_io_hdf.F90 b/io/cdnpot_io_hdf.F90 index ef1e98e88f19e3ff93f2ba1ef2ee3fabbd39e8f0..28358b7e4e1bff59d4c36e2fb50cd09e4c08bec1 100644 --- a/io/cdnpot_io_hdf.F90 +++ b/io/cdnpot_io_hdf.F90 @@ -16,14 +16,16 @@ MODULE m_cdnpot_io_hdf PRIVATE #ifdef CPP_HDF - PUBLIC openCDN_HDF, closeCDN_HDF + PUBLIC openCDN_HDF, openPOT_HDF, closeCDNPOT_HDF PUBLIC writeStarsHDF, readStarsHDF PUBLIC writeLatharmsHDF, readLatharmsHDF PUBLIC writeStructureHDF, readStructureHDF PUBLIC writeDensityHDF, readDensityHDF + PUBLIC writePotentialHDF, readPotentialHDF PUBLIC writeCoreDensityHDF, readCoreDensityHDF - PUBLIC writeHeaderData - PUBLIC isCoreDensityPresentHDF, isDensityEntryPresentHDF + PUBLIC writeCDNHeaderData, writePOTHeaderData + PUBLIC isCoreDensityPresentHDF + PUBLIC isDensityEntryPresentHDF, isPotentialEntryPresentHDF PUBLIC peekDensityEntryHDF #endif @@ -31,6 +33,7 @@ MODULE m_cdnpot_io_hdf PUBLIC DENSITY_TYPE_IN_const, DENSITY_TYPE_OUT_const PUBLIC DENSITY_TYPE_NOCO_IN_const, DENSITY_TYPE_NOCO_OUT_const PUBLIC DENSITY_TYPE_PRECOND_const + PUBLIC POTENTIAL_TYPE_IN_const, POTENTIAL_TYPE_OUT_const INTEGER, PARAMETER :: DENSITY_TYPE_UNDEFINED_const = 0 INTEGER, PARAMETER :: DENSITY_TYPE_IN_const = 1 @@ -39,6 +42,9 @@ MODULE m_cdnpot_io_hdf INTEGER, PARAMETER :: DENSITY_TYPE_NOCO_OUT_const = 4 INTEGER, PARAMETER :: DENSITY_TYPE_PRECOND_const = 5 + INTEGER, PARAMETER :: POTENTIAL_TYPE_IN_const = 1 + INTEGER, PARAMETER :: POTENTIAL_TYPE_OUT_const = 2 + CONTAINS #ifdef CPP_HDF @@ -87,10 +93,48 @@ MODULE m_cdnpot_io_hdf CALL h5gclose_f(generalGroupID, hdfError) END IF - END SUBROUTINE openCDN_HDF - SUBROUTINE closeCDN_HDF(fileID) + + SUBROUTINE openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex) + + INTEGER(HID_T), INTENT(OUT) :: fileID + INTEGER, INTENT(OUT) :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex + + INTEGER(HID_T) :: generalGroupID + INTEGER :: hdfError + LOGICAL :: l_exist + + currentStarsIndex = 0 + currentLatharmsIndex = 0 + currentStructureIndex = 0 + + INQUIRE(FILE='pot.hdf',EXIST=l_exist) + IF(l_exist) THEN ! only open file + CALL h5fopen_f('pot.hdf', H5F_ACC_RDWR_F, fileID, hdfError, H5P_DEFAULT_F) + + CALL h5gopen_f(fileID, '/general', generalGroupID, hdfError) + ! read in primary attributes from the header '/general' + CALL io_read_attint0(generalGroupID,'currentStarsIndex',currentStarsIndex) + CALL io_read_attint0(generalGroupID,'currentLatharmsIndex',currentLatharmsIndex) + CALL io_read_attint0(generalGroupID,'currentStructureIndex',currentStructureIndex) + + CALL h5gclose_f(generalGroupID, hdfError) + ELSE ! create file + CALL h5fcreate_f('pot.hdf', H5F_ACC_TRUNC_F, fileID, hdfError, H5P_DEFAULT_F, H5P_DEFAULT_F) + + CALL h5gcreate_f(fileID, '/general', generalGroupID, hdfError) + ! write initial values to primary attributes in the header '/general' + CALL io_write_attint0(generalGroupID,'currentStarsIndex',currentStarsIndex) + CALL io_write_attint0(generalGroupID,'currentLatharmsIndex',currentLatharmsIndex) + CALL io_write_attint0(generalGroupID,'currentStructureIndex',currentStructureIndex) + + CALL h5gclose_f(generalGroupID, hdfError) + END IF + + END SUBROUTINE openPOT_HDF + + SUBROUTINE closeCDNPOT_HDF(fileID) INTEGER(HID_T), INTENT(IN) :: fileID @@ -98,9 +142,9 @@ MODULE m_cdnpot_io_hdf CALL h5fclose_f(fileID, hdfError) - END SUBROUTINE closeCDN_HDF + END SUBROUTINE closeCDNPOT_HDF - SUBROUTINE writeHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,& + SUBROUTINE writeCDNHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,& readDensityIndex,lastDensityIndex) INTEGER(HID_T), INTENT(IN) :: fileID @@ -123,7 +167,27 @@ MODULE m_cdnpot_io_hdf CALL h5gclose_f(generalGroupID, hdfError) - END SUBROUTINE writeHeaderData + END SUBROUTINE writeCDNHeaderData + + SUBROUTINE writePOTHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex) + + INTEGER(HID_T), INTENT(IN) :: fileID + INTEGER, INTENT(IN) :: currentStarsIndex + INTEGER, INTENT(IN) :: currentLatharmsIndex + INTEGER, INTENT(IN) :: currentStructureIndex + + INTEGER(HID_T) :: generalGroupID + INTEGER :: hdfError + + CALL h5gopen_f(fileID, '/general', generalGroupID, hdfError) + + CALL io_write_attint0(generalGroupID,'currentStarsIndex',currentStarsIndex) + CALL io_write_attint0(generalGroupID,'currentLatharmsIndex',currentLatharmsIndex) + CALL io_write_attint0(generalGroupID,'currentStructureIndex',currentStructureIndex) + + CALL h5gclose_f(generalGroupID, hdfError) + + END SUBROUTINE writePOTHeaderData SUBROUTINE writeStarsHDF(fileID, starsIndex, stars) @@ -1345,6 +1409,216 @@ MODULE m_cdnpot_io_hdf END SUBROUTINE writeDensityHDF + SUBROUTINE writePotentialHDF(input, fileID, archiveName, potentialType,& + starsIndex, latharmsIndex, structureIndex,& + iter,fr,fpw,fz,fzxy) + + TYPE(t_input), INTENT(IN) :: input + INTEGER(HID_T), INTENT(IN) :: fileID + INTEGER, INTENT(IN) :: potentialType + INTEGER, INTENT(IN) :: starsIndex, latharmsIndex, structureIndex + CHARACTER(LEN=*), INTENT(IN) :: archiveName + + INTEGER, INTENT (IN) :: iter + + REAL, INTENT (IN) :: fr(:,:,:,:) + REAL, INTENT (IN) :: fz(:,:,:) + COMPLEX, INTENT (IN) :: fpw(:,:) + COMPLEX, INTENT (IN) :: fzxy(:,:,:,:) + + INTEGER :: ntype,jmtd,nmzd,nmzxyd,nlhd,ng3,ng2 + INTEGER :: nmz, nvac, od_nq2, nmzxy + INTEGER :: hdfError + LOGICAL :: l_film, l_exist + INTEGER(HID_T) :: archiveID, groupID + CHARACTER(LEN=30) :: groupName, potentialTypeName + INTEGER(HSIZE_T) :: dims(7) + INTEGER :: dimsInt(7) + + INTEGER(HID_T) :: frSpaceID, frSetID + INTEGER(HID_T) :: fpwSpaceID, fpwSetID + INTEGER(HID_T) :: fzSpaceID, fzSetID + INTEGER(HID_T) :: fzxySpaceID, fzxySetID + + WRITE(groupname,'(a,i0)') '/structure-', structureIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Structure entry "//TRIM(ADJUSTL(groupName))//" does not exist",calledby ="writePotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + CALL io_read_attlog0(groupID,'l_film',l_film) + CALL io_read_attint0(groupID,'ntype',ntype) + CALL io_read_attint0(groupID,'jmtd',jmtd) + CALL io_read_attint0(groupID,'nmzd',nmzd) + CALL io_read_attint0(groupID,'nmzxyd',nmzxyd) + CALL io_read_attint0(groupID,'nmzxy',nmzxy) + CALL io_read_attint0(groupID,'nmz',nmz) + CALL io_read_attint0(groupID,'nvac',nvac) + CALL io_read_attint0(groupID,'od_nq2',od_nq2) + CALL h5gclose_f(groupID, hdfError) + + WRITE(groupname,'(a,i0)') '/latharms-', latharmsIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Latharms entry "//TRIM(ADJUSTL(groupName))//" does not exist",calledby ="writePotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + CALL io_read_attint0(groupID,'nlhd',nlhd) + CALL h5gclose_f(groupID, hdfError) + + WRITE(groupname,'(a,i0)') '/stars-', starsIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Stars entry "//TRIM(ADJUSTL(groupName))//" does not exist",calledby ="writePotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + CALL io_read_attint0(groupID,'ng3',ng3) + CALL io_read_attint0(groupID,'ng2',ng2) + CALL h5gclose_f(groupID, hdfError) + + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(archiveName))) + + SELECT CASE (potentialType) + CASE(POTENTIAL_TYPE_IN_const) + potentialTypeName = '/in' + CASE(POTENTIAL_TYPE_OUT_const) + potentialTypeName = '/out' + CASE DEFAULT + CALL juDFT_error("Unknown potential type selected",calledby ="writePotentialHDF") + END SELECT + + groupName = TRIM(ADJUSTL(archiveName))//TRIM(ADJUSTL(potentialTypeName)) + + IF(l_exist) THEN + CALL h5gopen_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError) + + CALL io_write_attint0(archiveID,'starsIndex',starsIndex) + CALL io_write_attint0(archiveID,'latharmsIndex',latharmsIndex) + CALL io_write_attint0(archiveID,'structureIndex',structureIndex) + CALL io_write_attint0(archiveID,'spins',input%jspins) + CALL io_write_attint0(archiveID,'iter',iter) + + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + + IF(l_exist) THEN + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + + dimsInt(:4)=(/jmtd,nlhd+1,ntype,input%jspins/) + CALL h5dopen_f(groupID, 'fr', frSetID, hdfError) + CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),fr) + CALL h5dclose_f(frSetID, hdfError) + + dimsInt(:3)=(/2,ng3,input%jspins/) + CALL h5dopen_f(groupID, 'fpw', fpwSetID, hdfError) + CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),fpw) + CALL h5dclose_f(fpwSetID, hdfError) + + IF (l_film) THEN + dimsInt(:3)=(/nmzd,2,input%jspins/) + CALL h5dopen_f(groupID, 'fz', fzSetID, hdfError) + CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),fz) + CALL h5dclose_f(fzSetID, hdfError) + + dimsInt(:5)=(/2,nmzxyd,ng2-1,2,input%jspins/) + CALL h5dopen_f(groupID, 'fzxy', fzxySetID, hdfError) + CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),fzxy) + CALL h5dclose_f(fzxySetID, hdfError) + END IF + + CALL h5gclose_f(groupID, hdfError) + ELSE + CALL h5gcreate_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + + dims(:4)=(/jmtd,nlhd+1,ntype,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(4,dims(:4),frSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fr", H5T_NATIVE_DOUBLE, frSpaceID, frSetID, hdfError) + CALL h5sclose_f(frSpaceID,hdfError) + CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),fr) + CALL h5dclose_f(frSetID, hdfError) + + dims(:3)=(/2,ng3,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(3,dims(:3),fpwSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fpw", H5T_NATIVE_DOUBLE, fpwSpaceID, fpwSetID, hdfError) + CALL h5sclose_f(fpwSpaceID,hdfError) + CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),fpw) + CALL h5dclose_f(fpwSetID, hdfError) + + IF (l_film) THEN + dims(:3)=(/nmzd,2,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(3,dims(:3),fzSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fz", H5T_NATIVE_DOUBLE, fzSpaceID, fzSetID, hdfError) + CALL h5sclose_f(fzSpaceID,hdfError) + CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),fz) + CALL h5dclose_f(fzSetID, hdfError) + + dims(:5)=(/2,nmzxyd,ng2-1,2,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(5,dims(:5),fzxySpaceID,hdfError) + CALL h5dcreate_f(groupID, "fzxy", H5T_NATIVE_DOUBLE, fzxySpaceID, fzxySetID, hdfError) + CALL h5sclose_f(fzxySpaceID,hdfError) + CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),fzxy) + CALL h5dclose_f(fzxySetID, hdfError) + END IF + + CALL h5gclose_f(groupID, hdfError) + END IF + + CALL h5gclose_f(archiveID, hdfError) + ELSE + CALL h5gcreate_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError) + + CALL io_write_attint0(archiveID,'starsIndex',starsIndex) + CALL io_write_attint0(archiveID,'latharmsIndex',latharmsIndex) + CALL io_write_attint0(archiveID,'structureIndex',structureIndex) + CALL io_write_attint0(archiveID,'spins',input%jspins) + CALL io_write_attint0(archiveID,'iter',iter) + + CALL h5gcreate_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + + dims(:4)=(/jmtd,nlhd+1,ntype,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(4,dims(:4),frSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fr", H5T_NATIVE_DOUBLE, frSpaceID, frSetID, hdfError) + CALL h5sclose_f(frSpaceID,hdfError) + CALL io_write_real4(frSetID,(/1,1,1,1/),dimsInt(:4),fr) + CALL h5dclose_f(frSetID, hdfError) + + dims(:3)=(/2,ng3,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(3,dims(:3),fpwSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fpw", H5T_NATIVE_DOUBLE, fpwSpaceID, fpwSetID, hdfError) + CALL h5sclose_f(fpwSpaceID,hdfError) + CALL io_write_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),fpw) + CALL h5dclose_f(fpwSetID, hdfError) + + IF (l_film) THEN + dims(:3)=(/nmzd,2,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(3,dims(:3),fzSpaceID,hdfError) + CALL h5dcreate_f(groupID, "fz", H5T_NATIVE_DOUBLE, fzSpaceID, fzSetID, hdfError) + CALL h5sclose_f(fzSpaceID,hdfError) + CALL io_write_real3(fzSetID,(/1,1,1/),dimsInt(:3),fz) + CALL h5dclose_f(fzSetID, hdfError) + + dims(:5)=(/2,nmzxyd,ng2-1,2,input%jspins/) + dimsInt = dims + CALL h5screate_simple_f(5,dims(:5),fzxySpaceID,hdfError) + CALL h5dcreate_f(groupID, "fzxy", H5T_NATIVE_DOUBLE, fzxySpaceID, fzxySetID, hdfError) + CALL h5sclose_f(fzxySpaceID,hdfError) + CALL io_write_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),fzxy) + CALL h5dclose_f(fzxySetID, hdfError) + END IF + + CALL h5gclose_f(groupID, hdfError) + + CALL h5gclose_f(archiveID, hdfError) + END IF + + END SUBROUTINE writePotentialHDF + SUBROUTINE readDensityHDF(fileID, archiveName, densityType,& fermiEnergy,l_qfix,iter,fr,fpw,fz,fzxy,cdom,cdomvz,cdomvxy) @@ -1411,13 +1685,13 @@ MODULE m_cdnpot_io_hdf CASE(DENSITY_TYPE_PRECOND_const) densityTypeName = '/precond' CASE DEFAULT - CALL juDFT_error("Unknown density type selected",calledby ="writeDensityHDF") + CALL juDFT_error("Unknown density type selected",calledby ="readDensityHDF") END SELECT groupName = TRIM(ADJUSTL(archiveName))//TRIM(ADJUSTL(densityTypeName)) l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) IF(.NOT.l_exist) THEN - CALL juDFT_error('density entry '//TRIM(ADJUSTL(groupName))//' does not exist.' ,calledby ="readDensityHDF") + CALL juDFT_error('Density entry '//TRIM(ADJUSTL(groupName))//' does not exist.' ,calledby ="readDensityHDF") END IF CALL h5gopen_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError) @@ -1517,6 +1791,129 @@ MODULE m_cdnpot_io_hdf END SUBROUTINE readDensityHDF + SUBROUTINE readPotentialHDF(fileID, archiveName, potentialType,& + iter,fr,fpw,fz,fzxy) + + INTEGER(HID_T), INTENT(IN) :: fileID + INTEGER, INTENT(IN) :: potentialType + CHARACTER(LEN=*), INTENT(IN) :: archiveName + + INTEGER, INTENT (OUT) :: iter + + REAL, INTENT (OUT) :: fr(:,:,:,:) + REAL, INTENT (OUT) :: fz(:,:,:) + COMPLEX, INTENT (OUT) :: fpw(:,:) + COMPLEX, INTENT (OUT) :: fzxy(:,:,:,:) + + INTEGER :: starsIndex, latharmsIndex, structureIndex + INTEGER :: jspins + INTEGER :: ntype,jmtd,nmzd,nmzxyd,nlhd,ng3,ng2 + INTEGER :: nmz, nvac, od_nq2, nmzxy + LOGICAL :: l_film, l_exist + INTEGER(HID_T) :: archiveID, groupID, groupBID + INTEGER :: hdfError + CHARACTER(LEN=30) :: groupName, groupBName, potentialTypeName + INTEGER :: dimsInt(7) + + INTEGER(HID_T) :: frSetID + INTEGER(HID_T) :: fpwSetID + INTEGER(HID_T) :: fzSetID + INTEGER(HID_T) :: fzxySetID + + + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(archiveName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error('density archive '//TRIM(ADJUSTL(archiveName))//' does not exist.' ,calledby ="readPotentialHDF") + END IF + + SELECT CASE (potentialType) + CASE(POTENTIAL_TYPE_IN_const) + potentialTypeName = '/in' + CASE(POTENTIAL_TYPE_OUT_const) + potentialTypeName = '/out' + CASE DEFAULT + CALL juDFT_error("Unknown potential type selected",calledby ="readPotentialHDF") + END SELECT + + groupName = TRIM(ADJUSTL(archiveName))//TRIM(ADJUSTL(potentialTypeName)) + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error('Potential entry '//TRIM(ADJUSTL(groupName))//' does not exist.' ,calledby ="readPotentialHDF") + END IF + + CALL h5gopen_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError) + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) + + CALL io_read_attint0(archiveID,'starsIndex',starsIndex) + CALL io_read_attint0(archiveID,'latharmsIndex',latharmsIndex) + CALL io_read_attint0(archiveID,'structureIndex',structureIndex) + CALL io_read_attint0(archiveID,'spins',jspins) + CALL io_read_attint0(archiveID,'iter',iter) + + WRITE(groupBName,'(a,i0)') '/structure-', structureIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupBName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Structure entry "//TRIM(ADJUSTL(groupBName))//" does not exist",calledby ="readPotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupBName)), groupBID, hdfError) + CALL io_read_attlog0(groupBID,'l_film',l_film) + CALL io_read_attint0(groupBID,'ntype',ntype) + CALL io_read_attint0(groupBID,'jmtd',jmtd) + CALL io_read_attint0(groupBID,'nmzd',nmzd) + CALL io_read_attint0(groupBID,'nmzxyd',nmzxyd) + CALL io_read_attint0(groupBID,'nmzxy',nmzxy) + CALL io_read_attint0(groupBID,'nmz',nmz) + CALL io_read_attint0(groupBID,'nvac',nvac) + CALL io_read_attint0(groupBID,'od_nq2',od_nq2) + CALL h5gclose_f(groupBID, hdfError) + + WRITE(groupBName,'(a,i0)') '/latharms-', latharmsIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupBName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Latharms entry "//TRIM(ADJUSTL(groupBName))//" does not exist",calledby ="readPotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupBName)), groupBID, hdfError) + CALL io_read_attint0(groupBID,'nlhd',nlhd) + CALL h5gclose_f(groupBID, hdfError) + + WRITE(groupBName,'(a,i0)') '/stars-', starsIndex + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupBName))) + IF(.NOT.l_exist) THEN + CALL juDFT_error("Stars entry "//TRIM(ADJUSTL(groupBName))//" does not exist",calledby ="readPotentialHDF") + END IF + CALL h5gopen_f(fileID, TRIM(ADJUSTL(groupBName)), groupBID, hdfError) + CALL io_read_attint0(groupBID,'ng3',ng3) + CALL io_read_attint0(groupBID,'ng2',ng2) + CALL h5gclose_f(groupBID, hdfError) + + dimsInt(:4)=(/jmtd,nlhd+1,ntype,jspins/) + CALL h5dopen_f(groupID, 'fr', frSetID, hdfError) + CALL io_read_real4(frSetID,(/1,1,1,1/),dimsInt(:4),fr) + CALL h5dclose_f(frSetID, hdfError) + + dimsInt(:3)=(/2,ng3,jspins/) + CALL h5dopen_f(groupID, 'fpw', fpwSetID, hdfError) + CALL io_read_complex2(fpwSetID,(/-1,1,1/),dimsInt(:3),fpw) + CALL h5dclose_f(fpwSetID, hdfError) + + IF (l_film) THEN + dimsInt(:3)=(/nmzd,2,jspins/) + CALL h5dopen_f(groupID, 'fz', fzSetID, hdfError) + CALL io_read_real3(fzSetID,(/1,1,1/),dimsInt(:3),fz) + CALL h5dclose_f(fzSetID, hdfError) + + dimsInt(:5)=(/2,nmzxyd,ng2-1,2,jspins/) + CALL h5dopen_f(groupID, 'fzxy', fzxySetID, hdfError) + CALL io_read_complex4(fzxySetID,(/-1,1,1,1,1/),dimsInt(:5),fzxy) + CALL h5dclose_f(fzxySetID, hdfError) + END IF + + CALL h5gclose_f(groupID, hdfError) + CALL h5gclose_f(archiveID, hdfError) + + END SUBROUTINE readPotentialHDF + + SUBROUTINE peekDensityEntryHDF(fileID, archiveName, densityType,& iter, starsIndex, latharmsIndex, structureIndex,& previousDensityIndex, jspins,& @@ -1540,7 +1937,7 @@ MODULE m_cdnpot_io_hdf l_exist = io_groupexists(fileID,TRIM(ADJUSTL(archiveName))) IF(.NOT.l_exist) THEN - CALL juDFT_error('density archive '//TRIM(ADJUSTL(archiveName))//' does not exist.' ,calledby ="readDensityHDF") + CALL juDFT_error('density archive '//TRIM(ADJUSTL(archiveName))//' does not exist.' ,calledby ="peekDensityHDF") END IF localDensityType = densityType @@ -1570,13 +1967,13 @@ MODULE m_cdnpot_io_hdf CASE(DENSITY_TYPE_PRECOND_const) densityTypeName = '/precond' CASE DEFAULT - CALL juDFT_error("Unknown density type selected",calledby ="writeDensityHDF") + CALL juDFT_error("Unknown density type selected",calledby ="peekDensityHDF") END SELECT groupName = TRIM(ADJUSTL(archiveName))//TRIM(ADJUSTL(densityTypeName)) l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) IF(.NOT.l_exist) THEN - CALL juDFT_error('density entry '//TRIM(ADJUSTL(groupName))//' does not exist.' ,calledby ="readDensityHDF") + CALL juDFT_error('density entry '//TRIM(ADJUSTL(groupName))//' does not exist.' ,calledby ="peekDensityHDF") END IF CALL h5gopen_f(fileID, TRIM(ADJUSTL(archiveName)), archiveID, hdfError) @@ -1745,7 +2142,7 @@ MODULE m_cdnpot_io_hdf CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,& readDensityIndex,lastDensityIndex) isCoreDensityPresentHDF = io_groupexists(fileID,'/cdnc') - CALL closeCDN_HDF(fileID) + CALL closeCDNPOT_HDF(fileID) END FUNCTION LOGICAL FUNCTION isDensityEntryPresentHDF(fileID,archiveName,densityType) @@ -1800,6 +2197,36 @@ MODULE m_cdnpot_io_hdf isDensityEntryPresentHDF = l_exist END FUNCTION + LOGICAL FUNCTION isPotentialEntryPresentHDF(fileID,archiveName,potentialType) + + INTEGER(HID_T), INTENT(IN) :: fileID + INTEGER , INTENT(IN) :: potentialType + CHARACTER(LEN=*), INTENT(IN) :: archiveName + + CHARACTER(LEN=30) :: groupName, potentialTypeName + LOGICAL :: l_exist + + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(archiveName))) + IF(.NOT.l_exist) THEN + isPotentialEntryPresentHDF = .FALSE. + RETURN + END IF + + SELECT CASE (potentialType) + CASE(POTENTIAL_TYPE_IN_const) + potentialTypeName = '/in' + CASE(POTENTIAL_TYPE_OUT_const) + potentialTypeName = '/out' + CASE DEFAULT + CALL juDFT_error("Unknown potential type selected",calledby ="isPotentialEntryPresentHDF") + END SELECT + + groupName = TRIM(ADJUSTL(archiveName))//TRIM(ADJUSTL(potentialTypeName)) + l_exist = io_groupexists(fileID,TRIM(ADJUSTL(groupName))) + + isPotentialEntryPresentHDF = l_exist + END FUNCTION + #endif END MODULE m_cdnpot_io_hdf diff --git a/io/pot_io.f90 b/io/pot_io.F90 similarity index 57% rename from io/pot_io.f90 rename to io/pot_io.F90 index bb76521aaf09d3e31fc69b9d5978a631b551fa24..ea5516ede4d7e5b20d9d82d6bd6a81d3f6246173 100644 --- a/io/pot_io.f90 +++ b/io/pot_io.F90 @@ -18,6 +18,10 @@ MODULE m_pot_io USE m_juDFT USE m_loddop USE m_wrtdop + USE m_cdnpot_io_hdf +#ifdef CPP_HDF + USE hdf5 +#endif IMPLICIT NONE PRIVATE @@ -58,19 +62,57 @@ MODULE m_pot_io LOGICAL :: l_exist CHARACTER(len=30) :: filename +#ifdef CPP_HDF + INTEGER(HID_T) :: fileID +#endif + INTEGER :: currentStarsIndex,currentLatharmsIndex + INTEGER :: currentStructureIndex + INTEGER :: potentialType + CHARACTER(LEN=30) :: archiveName + CALL getMode(mode) IF(mode.EQ.POT_HDF5_MODE) THEN +#ifdef CPP_HDF INQUIRE(FILE='pot.hdf',EXIST=l_exist) IF (l_exist) THEN - !load density from pot.hdf and exit subroutine + CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex) - RETURN + archiveName = 'illegalPotentialArchive' + IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN + archiveName = 'pottot' + END IF + IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN + archiveName = 'potcoul' + END IF + IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN + archiveName = 'potx' + END IF + + potentialType = POTENTIAL_TYPE_IN_const + + l_exist = isPotentialEntryPresentHDF(fileID,archiveName,potentialType) + + CALL closeCDNPOT_HDF(fileID) + END IF + + IF(l_exist) THEN + CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex) + + CALL readPotentialHDF(fileID, archiveName, potentialType,& + iter,fr,fpw,fz,fzxy) + + CALL closeCDNPOT_HDF(fileID) ELSE - WRITE(*,*) 'pot.hdf file not found.' + WRITE(*,*) 'Potential entry or pot.hdf file not found.' WRITE(*,*) 'Falling back to stream access file pot.str.' mode = POT_STREAM_MODE END IF +#else + WRITE(*,*) 'Not compiled for pot.hdf file usage.' + WRITE(*,*) 'Falling back to stream access file pot.str.' + mode = POT_STREAM_MODE +#endif END IF IF(mode.EQ.POT_STREAM_MODE) THEN @@ -114,15 +156,17 @@ MODULE m_pot_io END SUBROUTINE readPotential - SUBROUTINE writePotential(stars,vacuum,atoms,sphhar,input,sym,archiveType,& + SUBROUTINE writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,& iter,fr,fpw,fz,fzxy) TYPE(t_stars),INTENT(IN) :: stars TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_atoms),INTENT(IN) :: atoms + TYPE(t_cell), INTENT(IN) :: cell TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_input),INTENT(IN) :: input TYPE(t_sym),INTENT(IN) :: sym + TYPE(t_oneD),INTENT(IN) :: oneD INTEGER, INTENT (IN) :: iter INTEGER, INTENT (IN) :: archiveType @@ -133,14 +177,64 @@ MODULE m_pot_io ! local variables INTEGER :: mode, iUnit - LOGICAL :: l_exist + LOGICAL :: l_exist, l_storeIndices CHARACTER(len=30) :: filename +#ifdef CPP_HDF + INTEGER(HID_T) :: fileID +#endif + INTEGER :: currentStarsIndex,currentLatharmsIndex + INTEGER :: currentStructureIndex + INTEGER :: potentialType + CHARACTER(LEN=30) :: archiveName + CALL getMode(mode) IF(mode.EQ.POT_HDF5_MODE) THEN - ! Write potential to pot.hdf file - STOP 'POT_HDF5_MODE not yet implemented!' +#ifdef CPP_HDF + CALL openPOT_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex) + + l_storeIndices = .FALSE. + IF (currentStarsIndex.EQ.0) THEN + currentStarsIndex = 1 + l_storeIndices = .TRUE. + CALL writeStarsHDF(fileID, currentStarsIndex, stars) + END IF + IF (currentLatharmsIndex.EQ.0) THEN + currentLatharmsIndex = 1 + l_storeIndices = .TRUE. + CALL writeLatharmsHDF(fileID, currentLatharmsIndex, sphhar) + END IF + IF(currentStructureIndex.EQ.0) THEN + currentStructureIndex = 1 + l_storeIndices = .TRUE. + CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, currentStructureIndex) + END IF + + archiveName = 'illegalPotentialArchive' + IF (archiveType.EQ.POT_ARCHIVE_TYPE_TOT_const) THEN + archiveName = 'pottot' + END IF + IF (archiveType.EQ.POT_ARCHIVE_TYPE_COUL_const) THEN + archiveName = 'potcoul' + END IF + IF (archiveType.EQ.POT_ARCHIVE_TYPE_X_const) THEN + archiveName = 'potx' + END IF + + potentialType = POTENTIAL_TYPE_IN_const + + CALL writePotentialHDF(input, fileID, archiveName, potentialType,& + currentStarsIndex, currentLatharmsIndex, currentStructureIndex,& + iter,fr,fpw,fz,fzxy) + + IF(l_storeIndices) THEN + CALL writePOTHeaderData(fileID,currentStarsIndex,currentLatharmsIndex,& + currentStructureIndex) + END IF + + CALL closeCDNPOT_HDF(fileID) +#endif ELSE IF(mode.EQ.POT_STREAM_MODE) THEN ! Write potential to pot.str file STOP 'POT_STREAM_MODE not yet implemented!' @@ -170,8 +264,15 @@ MODULE m_pot_io INTEGER, INTENT(OUT) :: mode mode = POT_DIRECT_MODE - IF (juDFT_was_argument("-stream_pot")) mode=POT_STREAM_MODE - IF (juDFT_was_argument("-hdf_pot")) mode=POT_HDF5_MODE + IF (juDFT_was_argument("-stream_cdn")) mode=POT_STREAM_MODE + IF (juDFT_was_argument("-hdf_cdn")) THEN +#ifdef CPP_HDF + mode=POT_HDF5_MODE +#else + WRITE(*,*) 'Code not compiled with HDF5 support.' + WRITE(*,*) 'Falling back to direct access.' +#endif + END IF END SUBROUTINE getMode END MODULE m_pot_io diff --git a/main/vgen.F90 b/main/vgen.F90 index eb534937279986b60f9c2953dbdff07c5581fb43..fab7227b6cdad376ef8cb13222520303f09ee210 100644 --- a/main/vgen.F90 +++ b/main/vgen.F90 @@ -417,7 +417,7 @@ CONTAINS vpw_w(1:stars%ng3,js)=vpw_w(1:stars%ng3,js)/stars%nstr(1:stars%ng3) ! the PW-coulomb part is not ! used otherwise anyway. ENDDO - CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_COUL_const,& + CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_COUL_const,& iter,vr,vpw_w,vz,vxy) DO js = 1,input%jspins DO i = 1,stars%ng3 @@ -836,7 +836,7 @@ CONTAINS vpw_w(i,js)=vpw_w(i,js)/stars%nstr(i) ENDDO ENDDO - CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_TOT_const,& + CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_TOT_const,& iter,vr,vpw_w,vz,vxy) DO js=1,input%jspins @@ -845,7 +845,7 @@ CONTAINS ENDDO ENDDO - CALL writePotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_X_const,& + CALL writePotential(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,POT_ARCHIVE_TYPE_X_const,& iter,vxr,vxpw_w,vz,vxy) END IF