Commit 7784d712 authored by Gregor Michalicek's avatar Gregor Michalicek

Use density IO framewaork for more files

...not fully tested yet.
parent ac09716d
......@@ -133,7 +133,7 @@ MODULE m_cdn_io
SUBROUTINE readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
relCdnIndex,fermiEnergy,l_qfix,den)
relCdnIndex,fermiEnergy,l_qfix,den,inFilename)
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
......@@ -151,6 +151,8 @@ MODULE m_cdn_io
REAL, INTENT (OUT) :: fermiEnergy
LOGICAL, INTENT (OUT) :: l_qfix
CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: inFilename
! local variables
INTEGER :: mode, datend, k, i, iVac, j, iUnit, l, numLines, ioStatus, iofl
LOGICAL :: l_exist, l_rhomatFile, l_DimChange
......@@ -179,10 +181,13 @@ MODULE m_cdn_io
densityType = 0
archiveName = ''
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
filename = 'cdn.hdf'
IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF (l_exist) THEN
CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
IF (archiveType.EQ.CDN_ARCHIVE_TYPE_CDN_const) THEN
archiveName = 'cdn'
......@@ -213,7 +218,7 @@ MODULE m_cdn_io
IF (l_exist) THEN
CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
CALL readDensityHDF(fileID, input, stars, sphhar, atoms, vacuum, oneD, archiveName, densityType,&
fermiEnergy,l_qfix,l_DimChange,den)
......@@ -225,8 +230,8 @@ MODULE m_cdn_io
1,-1.0,fermiEnergy,l_qfix,den)
END IF
ELSE
WRITE(*,*) 'cdn.hdf file or relevant density entry not found.'
WRITE(*,*) 'Falling back to stream access file cdn.str.'
WRITE(*,*) TRIM(ADJUSTL(filename))//' file or relevant density entry not found.'
WRITE(*,*) 'Falling back to stream access.'
mode = CDN_STREAM_MODE
END IF
#endif
......@@ -260,6 +265,8 @@ MODULE m_cdn_io
filename = 'cdn'
END IF
IF(PRESENT(inFilename)) filename = inFilename
INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF(.NOT.l_exist) THEN
CALL juDFT_error("charge density file "//TRIM(ADJUSTL(filename))//" missing",calledby ="readDensity")
......@@ -341,7 +348,7 @@ MODULE m_cdn_io
END SUBROUTINE readDensity
SUBROUTINE writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
relCdnIndex,distance,fermiEnergy,l_qfix,den)
relCdnIndex,distance,fermiEnergy,l_qfix,den,inFilename)
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
......@@ -359,6 +366,8 @@ MODULE m_cdn_io
REAL, INTENT (IN) :: fermiEnergy, distance
LOGICAL, INTENT (IN) :: l_qfix
CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: inFilename
TYPE(t_stars) :: starsTemp
TYPE(t_vacuum) :: vacuumTemp
TYPE(t_atoms) :: atomsTemp
......@@ -406,7 +415,7 @@ MODULE m_cdn_io
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
CALL checkAndWriteMetadataHDF(fileID, input, atoms, cell, vacuum, oneD, stars, sphhar, sym,&
currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
......@@ -491,6 +500,8 @@ MODULE m_cdn_io
filename = 'cdn'
END IF
IF(PRESENT(inFilename)) filename = inFilename
IF ((relCdnIndex.EQ.1).AND.(archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(den%iter.EQ.0)) THEN
INQUIRE(file=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF(l_exist) THEN
......
......@@ -54,15 +54,18 @@ MODULE m_cdnpot_io_hdf
#ifdef CPP_HDF
SUBROUTINE openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
currentStepfunctionIndex,readDensityIndex,lastDensityIndex)
currentStepfunctionIndex,readDensityIndex,lastDensityIndex,inFilename)
INTEGER(HID_T), INTENT(OUT) :: fileID
INTEGER, INTENT(OUT) :: currentStarsIndex,currentLatharmsIndex,currentStructureIndex
INTEGER, INTENT(OUT) :: currentStepfunctionIndex, readDensityIndex,lastDensityIndex
CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: inFilename
INTEGER(HID_T) :: generalGroupID
INTEGER :: hdfError, fileFormatVersion
LOGICAL :: l_exist
CHARACTER(LEN=30) :: filename
currentStarsIndex = 0
currentLatharmsIndex = 0
......@@ -72,9 +75,12 @@ MODULE m_cdnpot_io_hdf
lastDensityIndex = 0
fileFormatVersion = 0
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
filename = 'cdn.hdf'
IF(PRESENT(inFilename)) filename = TRIM(ADJUSTL(inFilename))//'.hdf'
INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF(l_exist) THEN ! only open file
CALL h5fopen_f('cdn.hdf', H5F_ACC_RDWR_F, fileID, hdfError, H5P_DEFAULT_F)
CALL h5fopen_f(TRIM(ADJUSTL(filename)), H5F_ACC_RDWR_F, fileID, hdfError, H5P_DEFAULT_F)
CALL h5gopen_f(fileID, '/general', generalGroupID, hdfError)
! read in primary attributes from the header '/general'
......@@ -88,11 +94,11 @@ MODULE m_cdnpot_io_hdf
CALL h5gclose_f(generalGroupID, hdfError)
IF(fileFormatVersion.GT.FILE_FORMAT_VERSION_const) THEN
WRITE(*,'(a,i4)') 'cdn.hdf has file format version ', fileFormatVersion
CALL juDFT_error('cdn.hdf file format not readable.' ,calledby ="openCDN_HDF")
WRITE(*,'(a,i4)') TRIM(ADJUSTL(filename))//' has file format version ', fileFormatVersion
CALL juDFT_error(TRIM(ADJUSTL(filename))//' file format not readable.' ,calledby ="openCDN_HDF")
END IF
ELSE ! create file
CALL h5fcreate_f('cdn.hdf', H5F_ACC_TRUNC_F, fileID, hdfError, H5P_DEFAULT_F, H5P_DEFAULT_F)
CALL h5fcreate_f(TRIM(ADJUSTL(filename)), 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'
......
......@@ -12,7 +12,7 @@ CONTAINS
SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
dimension,kpts,atoms,sphhar,stars,sym,obsolete,&
enpara,cell,noco,jij,vTot,results,oneD,coreSpecInput,&
inDen,outDen)
archiveType,outDen)
!*****************************************************
! Charge density generator
......@@ -61,10 +61,10 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_coreSpecInput),INTENT(IN) :: coreSpecInput
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_potden),INTENT(INOUT) :: inDen,outDen
TYPE(t_potden),INTENT(INOUT) :: outDen
!Scalar Arguments
INTEGER, INTENT (IN) :: eig_id
INTEGER, INTENT (IN) :: eig_id, archiveType
! Local type instances
TYPE(t_noco) :: noco_new
......@@ -85,8 +85,6 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
COMPLEX,ALLOCATABLE :: qa21(:), cdomvz(:,:)
!pk non-collinear (end)
CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
IF (mpi%irank.EQ.0) THEN
INQUIRE(file='enpara',exist=l_enpara)
IF (l_enpara) OPEN (40,file ='enpara',form = 'formatted',status ='unknown')
......@@ -101,8 +99,6 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
qvac(:,:,:,:) = 0.0
qvlay(:,:,:,:,:) = 0.0
outDen%iter = inDen%iter
!Set up pointer for backtransformation of from g-vector in
!positive domain fof carge density fftibox into stars
!In principle this can also be done in main program once.
......@@ -187,23 +183,8 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
8010 FORMAT (i13,2x,3e20.8,5x,2e20.8)
IF (sliceplot%slice) THEN
OPEN (20,file='cdn_slice',form='unformatted',status='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym, 20, outDen%iter,outDen%mt,outDen%pw,outDen%vacz,outDen%vacxy)
IF (noco%l_noco) THEN
WRITE (20) (outDen%pw(k,3),k=1,stars%ng3)
IF (input%film) THEN
ALLOCATE(cdomvz(vacuum%nmz,vacuum%nvac))
DO ivac = 1, vacuum%nvac
DO j = 1, vacuum%nmz
cdomvz(j,ivac) = CMPLX(outDen%vacz(j,ivac,3),outDen%vacz(j,ivac,4))
END DO
END DO
WRITE (20) ((cdomvz(j,ivac),j=1,vacuum%nmz),ivac=1,vacuum%nvac)
WRITE (20) (((outDen%vacxy(j,k-1,ivac,3),j=1,vacuum%nmzxy),k=2,oneD%odi%nq2) ,ivac=1,vacuum%nvac)
DEALLOCATE(cdomvz)
END IF
END IF
CLOSE(20)
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,-1.0,0.0,.FALSE.,outDen,'cdn_slice')
CALL juDFT_end("slice OK")
END IF
......
......@@ -176,11 +176,10 @@ CONTAINS
! Initialize and load inDen density (start)
CALL inDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
CALL inDenRot%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
IF (noco%l_noco) THEN
archiveType = CDN_ARCHIVE_TYPE_NOCO_const
ELSE
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
END IF
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
IF(mpi%irank.EQ.0) THEN
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,inDen)
......@@ -569,10 +568,12 @@ CONTAINS
!-Wannier
CALL timestart("generation of new charge density (total)")
CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
outDen%iter = inDen%iter
CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
DIMENSION,kpts,atoms,sphhar,stars,sym,obsolete,&
enpara_out,cell,noco,jij,vTot,results,oneD,coreSpecInput,&
inDen,outDen)
archiveType,outDen)
IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) dimension%neigd=dimension%neigd/2
!+t3e
......
......@@ -69,9 +69,9 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
! Local Scalars
INTEGER :: nrhomfile=26
INTEGER iden,ivac,ifft2,ifft3
INTEGER iden,ivac,ifft2,ifft3,archiveType
INTEGER imz,ityp,iri,ilh,imesh,lh,iq2,iq3,iter
REAL cdnup,cdndown,chden,mgden,theta,phi,zero,rho_11,rziw
REAL cdnup,cdndown,chden,mgden,theta,phi,zero,rho_11,rziw,fermiEnergyTemp
REAL rho_22,rho_21r,rho_21i,rhotot,mx,my,mz,fix,vz_r,vz_i
COMPLEX czero
CHARACTER*8 dop,iop,name(10)
......@@ -86,9 +86,9 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
!---> for testing: output of offdiag. output density matrix. to plot the
!---> offdiag. part of the output density matrix, that part has to be
!---> written the file rhomt21 in cdnmt.
LOGICAL :: l_fmpl2
REAL :: cdn11, cdn22
COMPLEX :: cdn21
LOGICAL :: l_fmpl2, l_qfix
REAL :: cdn11, cdn22
COMPLEX :: cdn21
COMPLEX, ALLOCATABLE :: rho21(:,:,:)
!---> end of test part
......@@ -120,34 +120,28 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
END IF
!---> reload the density matrix from file rhomat_inp
IF (.NOT. sliceplot%slice) THEN
OPEN (nrhomfile,FILE='rhomat_inp',FORM='unformatted',STATUS='unknown')
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
IF (.NOT.sliceplot%slice) THEN
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,den)
ELSE
OPEN (nrhomfile,FILE='cdn_slice',FORM='unformatted',STATUS='unknown')
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,den,'cdn_slice')
END IF
!---> first the diagonal elements of the density matrix
CALL loddop(stars,vacuum,atoms,sphhar,input,sym,nrhomfile,&
iter,rho,qpw,rht,rhtxy)
!---> and then the off-diagonal part
READ (nrhomfile,END=100,ERR=50) (cdom(iq3),iq3=1,stars%ng3)
IF (input%film) THEN
READ (nrhomfile,END=75,ERR=50) ((cdomvz(imz,ivac),imz=1,vacuum%nmz),ivac=1,vacuum%nvac)
READ (nrhomfile,END=75,ERR=50) (((cdomvxy(imz,iq2-1,ivac),imz=1,vacuum%nmzxy),iq2=2,stars%ng2),ivac=1,vacuum%nvac)
rho(:,0:,1:,:input%jspins) = den%mt(:,0:,1:,:input%jspins)
qpw(1:,:input%jspins) = den%pw(1:,:input%jspins)
rht(1:,1:,:input%jspins) = den%vacz(1:,1:,:input%jspins)
rhtxy(1:,1:,1:,:input%jspins) = den%vacxy(1:,1:,1:,:input%jspins)
IF(noco%l_noco) THEN
cdom = den%pw(:,3)
cdomvz(:,:) = CMPLX(den%vacz(:,:,3),den%vacz(:,:,4))
cdomvxy = den%vacxy(:,:,:,3)
END IF
GOTO 150
50 WRITE(6,*)'rhodirgen: ERROR: Problems while reading density'
WRITE(6,*)'matrix from file rhomat_inp.'
CALL juDFT_error("rhomatdir: ERROR while reading file rhomat_inp",calledby ="pldngen")
75 WRITE(6,*)'rhomatdir: ERROR: reached end of file rhomat_inp'
WRITE(6,*)'while reading the vacuum part of the off-diagonal'
WRITE(6,*)'element of the desity matrix.'
CALL juDFT_error("rhomatdir: ERROR while reading file rhomat_inp",calledby ="pldngen")
100 WRITE(6,*)'rhodirgen: WARNING: The file rhomat_inp does not'
WRITE(6,*)'contain off-diagonal part of the density matrix.'
WRITE(6,*)'Assuming collinear magnetization.'
150 CLOSE (nrhomfile)
IF (.NOT. sliceplot%slice) THEN
CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,inp%jspins,noco%l_noco,POTDEN_TYPE_DEN)
IF (.NOT. sliceplot%slice) THEN
CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
den%iter = iter
den%mt(:,0:,1:,:input%jspins) = rho(:,0:,1:,:input%jspins)
den%pw(1:,:input%jspins) = qpw(1:,:input%jspins)
......@@ -340,25 +334,28 @@ SUBROUTINE pldngen(sym,stars,atoms,sphhar,vacuum,&
0,-1.0,0.0,.FALSE.,den)
!---> save mx to file mdnx
OPEN (72,FILE='mdnx',FORM='unformatted',STATUS='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar,&
inp,sym,72,iter,rho(:,0:,1:,2:2),qpw(1:,2:2),rht(1:,1:,2:2),&
rhtxy(1:,1:,1:,2:2))
CLOSE (72)
den%mt(:,0:,1:,1) = rho(:,0:,1:,2)
den%pw(1:,1) = qpw(1:,2)
den%vacz(1:,1:,1) = rht(1:,1:,2)
den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,2)
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,den,'mdnx')
!---> save my to file mdny
OPEN (72,FILE='mdny',FORM='unformatted',STATUS='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar,&
inp,sym,72,iter,rho(:,0:,1:,3:3),qpw(1:,3:3),rht(1:,1:,3:3),&
rhtxy(1:,1:,1:,3:3))
CLOSE (72)
den%mt(:,0:,1:,1) = rho(:,0:,1:,3)
den%pw(1:,1) = qpw(1:,3)
den%vacz(1:,1:,1) = rht(1:,1:,3)
den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,3)
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,den,'mdny')
!---> save mz to file mdnz
OPEN (72,FILE='mdnz',FORM='unformatted',STATUS='unknown')
CALL wrtdop(stars,vacuum,atoms,sphhar,&
inp,sym,72,iter,rho(:,0:,1:,4:4),qpw(1:,4:4),rht(1:,1:,4:4),&
rhtxy(1:,1:,1:,4:4))
CLOSE (72)
den%mt(:,0:,1:,1) = rho(:,0:,1:,4)
den%pw(1:,1) = qpw(1:,4)
den%vacz(1:,1:,1) = rht(1:,1:,4)
den%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,4)
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,inp,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,den,'mdnz')
DEALLOCATE (qpw,rhtxy,cdom,cdomvz,cdomvxy,ris,fftwork,rvacxy,rho,rht)
......
......@@ -142,10 +142,8 @@ SUBROUTINE plotdop(oneD,dimension,stars,vacuum,sphhar,atoms,&
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den(i))
ELSE
OPEN(20,file = cdnFilenames(i),form='unformatted',status='old')
CALL loddop(stars,vacuum,atoms,sphhar,input,sym,20,&
den(i)%iter,den(i)%mt,den(i)%pw,den(i)%vacz,den(i)%vacxy)
CLOSE(20)
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den(i),TRIM(ADJUSTL(cdnFilenames(i))))
END IF
! Subtract core charge if input%score is set
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment