Commit 075d7301 authored by Gregor Michalicek's avatar Gregor Michalicek

Use density matrix stored in cdn.hdf file

Note: Also in this commit:

1. Changed test values in NiO_ldauXML test since now the broyden files get deleted
   when a new density matrix is generated
2. The NiO_ldau test should not work any more since linear mixing for the density
   matrix is removed for the old inp file. (It is not specified in the n_mmp_mat
   file any more)
parent db9cf2ab
...@@ -30,9 +30,7 @@ MODULE m_cdn_io ...@@ -30,9 +30,7 @@ MODULE m_cdn_io
PUBLIC printDensityFileInfo PUBLIC printDensityFileInfo
PUBLIC readDensity, writeDensity PUBLIC readDensity, writeDensity
PUBLIC isDensityFilePresent, isCoreDensityPresent PUBLIC isDensityFilePresent, isCoreDensityPresent
PUBLIC isDensityMatrixPresent
PUBLIC readCoreDensity, writeCoreDensity PUBLIC readCoreDensity, writeCoreDensity
PUBLIC readDensityMatrix
PUBLIC readStars, writeStars PUBLIC readStars, writeStars
PUBLIC readStepfunction, writeStepfunction PUBLIC readStepfunction, writeStepfunction
PUBLIC setStartingDensity, readPrevEFermi, deleteDensities PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
...@@ -153,7 +151,7 @@ MODULE m_cdn_io ...@@ -153,7 +151,7 @@ MODULE m_cdn_io
LOGICAL, INTENT (OUT) :: l_qfix LOGICAL, INTENT (OUT) :: l_qfix
! local variables ! local variables
INTEGER :: mode, datend, k, i, iVac, j, iUnit, l INTEGER :: mode, datend, k, i, iVac, j, iUnit, l, numLines, ioStatus, iofl
LOGICAL :: l_exist, l_rhomatFile, l_DimChange LOGICAL :: l_exist, l_rhomatFile, l_DimChange
CHARACTER(LEN=30) :: filename CHARACTER(LEN=30) :: filename
...@@ -308,6 +306,23 @@ MODULE m_cdn_io ...@@ -308,6 +306,23 @@ MODULE m_cdn_io
CLOSE(iUnit) CLOSE(iUnit)
END IF END IF
INQUIRE(FILE='n_mmp_mat',EXIST=l_exist)
IF(l_exist.AND.atoms%n_u.GT.0) THEN
OPEN (69,file='n_mmp_mat',status='unknown',form='formatted')
READ (69,'(7f20.13)',IOSTAT=ioStatus) den%mmpMat
REWIND(69)
numLines = 0
DO
READ (69,*,iostat=iofl)
IF (iofl < 0) EXIT
numLines = numLines + 1
END DO
IF (MOD(numLines,14*input%jspins).NE.0) THEN
CALL juDFT_error("strange n_mmp_mat-file...", calledby = "readDensity")
END IF
CLOSE(69)
END IF
END SUBROUTINE readDensity END SUBROUTINE readDensity
SUBROUTINE writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,& SUBROUTINE writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,inOrOutCDN,&
...@@ -588,13 +603,28 @@ MODULE m_cdn_io ...@@ -588,13 +603,28 @@ MODULE m_cdn_io
END IF END IF
CLOSE(iUnit) CLOSE(iUnit)
END IF
!write density matrix to n_mmp_mat file !write density matrix to n_mmp_mat_out file
IF((archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).AND.(relCdnIndex.EQ.1)) THEN IF((densityType.EQ.DENSITY_TYPE_IN_const).AND.(relCdnIndex.EQ.1).AND.&
IF(atoms%n_u.GT.0) THEN ((archiveType.EQ.CDN_ARCHIVE_TYPE_CDN1_const).OR.(archiveType.EQ.CDN_ARCHIVE_TYPE_NOCO_const))) THEN
IF(atoms%n_u.GT.0) THEN
OPEN (69,file='n_mmp_mat_out',status='replace',form='formatted')
WRITE (69,'(7f20.13)') den%mmpMat(:,:,:,:)
CLOSE (69)
IF ((mode.EQ.CDN_HDF5_MODE).AND..NOT.(input%ldauLinMix.AND.(input%ldauMixParam.EQ.0.0))) THEN
INQUIRE(file='n_mmp_mat',exist=l_exist)
IF(l_exist) THEN
CALL system('mv n_mmp_mat n_mmp_mat_old')
PRINT *,"n_mmp_mat moved to n_mmp_mat_old"
END IF
ELSE
OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
WRITE (69,'(7f20.13)') den%mmpMat(:,:,:,:)
CLOSE (69)
END IF END IF
END IF END IF
END IF END IF
END SUBROUTINE writeDensity END SUBROUTINE writeDensity
...@@ -808,7 +838,7 @@ MODULE m_cdn_io ...@@ -808,7 +838,7 @@ MODULE m_cdn_io
l_writeStructure = .TRUE. l_writeStructure = .TRUE.
ELSE ELSE
CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, currentStructureIndex) CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, currentStructureIndex)
CALL compareStructure(atoms, vacuum, cell, sym, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same) CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same)
IF(.NOT.l_same) THEN IF(.NOT.l_same) THEN
currentStructureIndex = currentStructureIndex + 1 currentStructureIndex = currentStructureIndex + 1
l_writeStructure = .TRUE. l_writeStructure = .TRUE.
...@@ -1148,84 +1178,6 @@ MODULE m_cdn_io ...@@ -1148,84 +1178,6 @@ MODULE m_cdn_io
END SUBROUTINE readStepfunction END SUBROUTINE readStepfunction
SUBROUTINE readDensityMatrix(input, atoms, n_mmp, l_error)
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_atoms), INTENT(IN) :: atoms
COMPLEX, INTENT(OUT) :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
atoms%n_u,input%jspins)
LOGICAL, INTENT(OUT) :: l_error
INTEGER :: mode, ioStatus, numLines, iofl, i
LOGICAL :: l_exist
CHARACTER(LEN=30) :: filename
REAL :: alpha, spinf
l_error = .FALSE.
ioStatus = 0
CALL getMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
l_exist = .FALSE.
IF(.NOT.l_exist) THEN
mode = CDN_STREAM_MODE
END IF
END IF
IF(mode.EQ.CDN_STREAM_MODE) THEN
INQUIRE(FILE='cdn.str',EXIST=l_exist)
IF (l_exist) THEN
STOP 'cdn.str code path not yet implemented!'
END IF
IF (.NOT.l_exist) THEN
mode = CDN_DIRECT_MODE
END IF
END IF
IF (mode.EQ.CDN_DIRECT_MODE) THEN
filename = 'n_mmp_mat'
INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
IF(.NOT.l_exist) THEN
l_error = .TRUE.
RETURN
END IF
OPEN (69,file=TRIM(ADJUSTL(filename)),status='unknown',form='formatted')
READ (69,'(7f20.13)',IOSTAT=ioStatus) n_mmp
REWIND(69)
numLines = 0
DO
READ (69,*,iostat=iofl)
IF (iofl < 0) EXIT
numLines = numLines + 1
END DO
IF (MOD(numLines,14*input%jspins) == 1) THEN !read in alpha, spinf in this case
REWIND(69)
DO i = 1, numLines-1
READ (69,*,iostat=iofl)
END DO
READ (69,'(2(6x,f5.3))',IOSTAT=ioStatus) alpha, spinf
IF(input%ldauLinMix.AND.(input%ldauMixParam.LT.0.0)) THEN
input%ldauMixParam = alpha
input%ldauSpinf = spinf
END IF
ELSE IF (MOD(numLines,14*input%jspins).NE.0) THEN
CALL juDFT_error("strange n_mmp_mat-file...",calledby ="readDensityMatrix")
ELSE IF(input%ldauMixParam.LT.0.0) THEN
input%ldauLinMix = .FALSE.
END IF
CLOSE(69)
! IF (ioStatus.NE.0) THEN
! l_error = .TRUE.
! RETURN
! END IF
END IF
END SUBROUTINE readDensityMatrix
SUBROUTINE setStartingDensity(l_noco) SUBROUTINE setStartingDensity(l_noco)
LOGICAL,INTENT(IN) :: l_noco LOGICAL,INTENT(IN) :: l_noco
...@@ -1503,38 +1455,4 @@ MODULE m_cdn_io ...@@ -1503,38 +1455,4 @@ MODULE m_cdn_io
isCoreDensityPresent = .FALSE. isCoreDensityPresent = .FALSE.
END FUNCTION isCoreDensityPresent END FUNCTION isCoreDensityPresent
LOGICAL FUNCTION isDensityMatrixPresent()
LOGICAL :: l_exist
INTEGER :: mode
CHARACTER(LEN=12) :: filename
CALL getMode(mode)
filename="n_mmp_mat"
INQUIRE (file=filename,exist=l_exist)
IF(l_exist) THEN
isDensityMatrixPresent = l_exist
RETURN
END IF
! IF (mode.EQ.CDN_STREAM_MODE) THEN
! STOP 'Not yet implemented! (CDN_STREAM_MODE in isDensityMatrixPresent)'
! END IF
IF (mode.EQ.CDN_HDF5_MODE) THEN
!#ifdef CPP_HDF
! STOP 'Not yet implemented (uncomment line)!'
! l_exist = isDensityMatrixPresentHDF()
! IF(l_exist) THEN
! isDensityMatrixPresent = l_exist
! RETURN
! END IF
!#endif
END IF
isDensityMatrixPresent = .FALSE.
END FUNCTION isDensityMatrixPresent
END MODULE m_cdn_io END MODULE m_cdn_io
...@@ -22,6 +22,8 @@ MODULE m_cdnpot_io_common ...@@ -22,6 +22,8 @@ MODULE m_cdnpot_io_common
USE hdf5 USE hdf5
#endif #endif
IMPLICIT NONE
CONTAINS CONTAINS
SUBROUTINE compareStars(stars, refStars, l_same) SUBROUTINE compareStars(stars, refStars, l_same)
...@@ -58,9 +60,10 @@ MODULE m_cdnpot_io_common ...@@ -58,9 +60,10 @@ MODULE m_cdnpot_io_common
END SUBROUTINE compareStepfunctions END SUBROUTINE compareStepfunctions
SUBROUTINE compareStructure(atoms, vacuum, cell, sym, refAtoms, refVacuum,& SUBROUTINE compareStructure(input, atoms, vacuum, cell, sym, refInput, refAtoms, refVacuum,&
refCell, refSym, l_same) refCell, refSym, l_same)
TYPE(t_input),INTENT(IN) :: input, refInput
TYPE(t_atoms),INTENT(IN) :: atoms, refAtoms TYPE(t_atoms),INTENT(IN) :: atoms, refAtoms
TYPE(t_vacuum),INTENT(IN) :: vacuum, refVacuum TYPE(t_vacuum),INTENT(IN) :: vacuum, refVacuum
TYPE(t_cell),INTENT(IN) :: cell, refCell TYPE(t_cell),INTENT(IN) :: cell, refCell
...@@ -68,16 +71,26 @@ MODULE m_cdnpot_io_common ...@@ -68,16 +71,26 @@ MODULE m_cdnpot_io_common
LOGICAL, INTENT(OUT) :: l_same LOGICAL, INTENT(OUT) :: l_same
INTEGER :: i
l_same = .TRUE. l_same = .TRUE.
IF(atoms%ntype.NE.refAtoms%ntype) l_same = .FALSE. IF(atoms%ntype.NE.refAtoms%ntype) l_same = .FALSE.
IF(atoms%nat.NE.refAtoms%nat) l_same = .FALSE. IF(atoms%nat.NE.refAtoms%nat) l_same = .FALSE.
IF(atoms%lmaxd.NE.refAtoms%lmaxd) l_same = .FALSE. IF(atoms%lmaxd.NE.refAtoms%lmaxd) l_same = .FALSE.
IF(atoms%jmtd.NE.refAtoms%jmtd) l_same = .FALSE. IF(atoms%jmtd.NE.refAtoms%jmtd) l_same = .FALSE.
IF(atoms%n_u.NE.refAtoms%n_u) l_same = .FALSE.
IF(vacuum%dvac.NE.refVacuum%dvac) l_same = .FALSE. IF(vacuum%dvac.NE.refVacuum%dvac) l_same = .FALSE.
IF(sym%nop.NE.refSym%nop) l_same = .FALSE. IF(sym%nop.NE.refSym%nop) l_same = .FALSE.
IF(sym%nop2.NE.refSym%nop2) l_same = .FALSE. IF(sym%nop2.NE.refSym%nop2) l_same = .FALSE.
IF(atoms%n_u.EQ.refAtoms%n_u) THEN
DO i = 1, atoms%n_u
IF (atoms%lda_u(i)%atomType.NE.refAtoms%lda_u(i)%atomType) l_same = .FALSE.
IF (atoms%lda_u(i)%l.NE.refAtoms%lda_u(i)%l) l_same = .FALSE.
END DO
END IF
IF(ANY(ABS(cell%amat(:,:)-refCell%amat(:,:)).GT.1e-10)) l_same = .FALSE. IF(ANY(ABS(cell%amat(:,:)-refCell%amat(:,:)).GT.1e-10)) l_same = .FALSE.
IF(l_same) THEN IF(l_same) THEN
IF(ANY(atoms%nz(:).NE.refAtoms%nz(:))) l_same = .FALSE. IF(ANY(atoms%nz(:).NE.refAtoms%nz(:))) l_same = .FALSE.
...@@ -138,7 +151,7 @@ MODULE m_cdnpot_io_common ...@@ -138,7 +151,7 @@ MODULE m_cdnpot_io_common
TYPE(t_sym) :: symTemp TYPE(t_sym) :: symTemp
INTEGER :: starsIndexTemp, structureIndexTemp INTEGER :: starsIndexTemp, structureIndexTemp
LOGICAL :: l_same, l_writeAll LOGICAL :: l_same, l_writeAll, l_exist
l_storeIndices = .FALSE. l_storeIndices = .FALSE.
l_writeAll = .FALSE. l_writeAll = .FALSE.
...@@ -149,7 +162,7 @@ MODULE m_cdnpot_io_common ...@@ -149,7 +162,7 @@ MODULE m_cdnpot_io_common
CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, sym, currentStructureIndex) CALL writeStructureHDF(fileID, input, atoms, cell, vacuum, oneD, sym, currentStructureIndex)
ELSE ELSE
CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, currentStructureIndex) CALL readStructureHDF(fileID, inputTemp, atomsTemp, cellTemp, vacuumTemp, oneDTemp, symTemp, currentStructureIndex)
CALL compareStructure(atoms, vacuum, cell, sym, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same) CALL compareStructure(input, atoms, vacuum, cell, sym, inputTemp, atomsTemp, vacuumTemp, cellTemp, symTemp, l_same)
IF(.NOT.l_same) THEN IF(.NOT.l_same) THEN
currentStructureIndex = currentStructureIndex + 1 currentStructureIndex = currentStructureIndex + 1
......
...@@ -1011,7 +1011,6 @@ MODULE m_cdnpot_io_hdf ...@@ -1011,7 +1011,6 @@ MODULE m_cdnpot_io_hdf
CALL h5gcreate_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError) CALL h5gcreate_f(fileID, TRIM(ADJUSTL(groupName)), groupID, hdfError)
CALL io_write_attlog0(groupID,'l_film',input%film) CALL io_write_attlog0(groupID,'l_film',input%film)
CALL io_write_attlog0(groupID,'ldauLinMix',input%ldauLinMix)
CALL io_write_attreal0(groupID,'omtil',cell%omtil) CALL io_write_attreal0(groupID,'omtil',cell%omtil)
CALL io_write_attreal0(groupID,'area',cell%area) CALL io_write_attreal0(groupID,'area',cell%area)
...@@ -1351,7 +1350,6 @@ MODULE m_cdnpot_io_hdf ...@@ -1351,7 +1350,6 @@ MODULE m_cdnpot_io_hdf
CALL io_read_attchar0(groupID,'namgrp',sym%namgrp) CALL io_read_attchar0(groupID,'namgrp',sym%namgrp)
IF(fileFormatVersion.GE.29) THEN IF(fileFormatVersion.GE.29) THEN
CALL io_read_attlog0(groupID,'ldauLinMix',input%ldauLinMix)
CALL io_read_attint0(groupID,'n_u',atoms%n_u) CALL io_read_attint0(groupID,'n_u',atoms%n_u)
IF(ALLOCATED(atoms%lda_u)) DEALLOCATE(atoms%lda_u) IF(ALLOCATED(atoms%lda_u)) DEALLOCATE(atoms%lda_u)
ALLOCATE(atoms%lda_u(atoms%n_u)) ALLOCATE(atoms%lda_u(atoms%n_u))
...@@ -2176,7 +2174,7 @@ MODULE m_cdnpot_io_hdf ...@@ -2176,7 +2174,7 @@ MODULE m_cdnpot_io_hdf
INTEGER :: ntype,jmtd,nmzd,nmzxyd,nlhd,ng3,ng2 INTEGER :: ntype,jmtd,nmzd,nmzxyd,nlhd,ng3,ng2
INTEGER :: nmz, nvac, od_nq2, nmzxy, n_u, i, j INTEGER :: nmz, nvac, od_nq2, nmzxy, n_u, i, j
INTEGER :: localDensityType INTEGER :: localDensityType
LOGICAL :: l_film, l_exist, ldauLinMix, l_mmpMatDimEquals LOGICAL :: l_film, l_exist, l_mmpMatDimEquals
INTEGER(HID_T) :: archiveID, groupID, groupBID, generalGroupID INTEGER(HID_T) :: archiveID, groupID, groupBID, generalGroupID
INTEGER :: hdfError, fileFormatVersion INTEGER :: hdfError, fileFormatVersion
CHARACTER(LEN=30) :: groupName, groupBName, densityTypeName CHARACTER(LEN=30) :: groupName, groupBName, densityTypeName
...@@ -2282,7 +2280,6 @@ MODULE m_cdnpot_io_hdf ...@@ -2282,7 +2280,6 @@ MODULE m_cdnpot_io_hdf
CALL io_read_attint0(groupBID,'od_nq2',od_nq2) CALL io_read_attint0(groupBID,'od_nq2',od_nq2)
IF(fileFormatVersion.GE.29) THEN IF(fileFormatVersion.GE.29) THEN
CALL io_read_attint0(groupBID,'n_u',n_u) CALL io_read_attint0(groupBID,'n_u',n_u)
CALL io_read_attlog0(groupBID,'ldauLinMix',ldauLinMix)
IF(n_u.GT.0) THEN IF(n_u.GT.0) THEN
ALLOCATE(ldau_AtomType(n_u), ldau_l(n_u), ldau_l_amf(n_u)) ALLOCATE(ldau_AtomType(n_u), ldau_l(n_u), ldau_l_amf(n_u))
ALLOCATE(ldau_U(n_u), ldau_J(n_u)) ALLOCATE(ldau_U(n_u), ldau_J(n_u))
...@@ -2366,7 +2363,6 @@ MODULE m_cdnpot_io_hdf ...@@ -2366,7 +2363,6 @@ MODULE m_cdnpot_io_hdf
l_mmpMatDimEquals = .TRUE. l_mmpMatDimEquals = .TRUE.
IF(fileFormatVersion.GE.29) THEN IF(fileFormatVersion.GE.29) THEN
IF(input%ldauLinMix.NE.ldauLinMix) l_DimChange = .TRUE. ! Am I sure about this? parameter change => dim change?
IF(atoms%n_u.NE.n_u) THEN IF(atoms%n_u.NE.n_u) THEN
l_DimChange = .TRUE. l_DimChange = .TRUE.
l_mmpMatDimEquals = .FALSE. l_mmpMatDimEquals = .FALSE.
......
...@@ -42,7 +42,6 @@ CONTAINS ...@@ -42,7 +42,6 @@ CONTAINS
REAL wronk REAL wronk
LOGICAL n_exist LOGICAL n_exist
CHARACTER*8 l_type*2,l_form*9 CHARACTER*8 l_type*2,l_form*9
CHARACTER*12 ::filename
REAL f(atoms%jmtd,2),g(atoms%jmtd,2),theta(atoms%n_u),phi(atoms%n_u),zero(atoms%n_u) REAL f(atoms%jmtd,2),g(atoms%jmtd,2),theta(atoms%n_u),phi(atoms%n_u),zero(atoms%n_u)
REAL f0(atoms%n_u,input%jspins),f2(atoms%n_u,input%jspins),f4(atoms%n_u,input%jspins),f6(atoms%n_u,input%jspins) REAL f0(atoms%n_u,input%jspins),f2(atoms%n_u,input%jspins),f4(atoms%n_u,input%jspins),f6(atoms%n_u,input%jspins)
REAL, ALLOCATABLE :: u(:,:,:,:,:,:) REAL, ALLOCATABLE :: u(:,:,:,:,:,:)
...@@ -50,7 +49,7 @@ CONTAINS ...@@ -50,7 +49,7 @@ CONTAINS
! !
! look, whether density matrix exists already: ! look, whether density matrix exists already:
! !
IF (isDensityMatrixPresent().AND.atoms%n_u>0) THEN IF (ANY(inDen%mmpMat(:,:,:,:).NE.0.0).AND.atoms%n_u>0) THEN
! calculate slater integrals from u and j ! calculate slater integrals from u and j
CALL uj2f(input%jspins,atoms,f0,f2,f4,f6) CALL uj2f(input%jspins,atoms,f0,f2,f4,f6)
...@@ -118,7 +117,6 @@ CONTAINS ...@@ -118,7 +117,6 @@ CONTAINS
ELSE ELSE
IF (mpi%irank.EQ.0) THEN IF (mpi%irank.EQ.0) THEN
WRITE (*,*) 'no density matrix found ... skipping LDA+U' WRITE (*,*) 'no density matrix found ... skipping LDA+U'
WRITE(*,*) "File:",filename
ENDIF ENDIF
pot%mmpMat(:,:,:,:) = CMPLX(0.0,0.0) pot%mmpMat(:,:,:,:) = CMPLX(0.0,0.0)
results%e_ldau = 0.0 results%e_ldau = 0.0
......
...@@ -182,6 +182,7 @@ CONTAINS ...@@ -182,6 +182,7 @@ CONTAINS
archiveType = CDN_ARCHIVE_TYPE_CDN1_const archiveType = CDN_ARCHIVE_TYPE_CDN1_const
END IF END IF
ALLOCATE (inDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,atoms%n_u),input%jspins)) ALLOCATE (inDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,atoms%n_u),input%jspins))
inDen%mmpMat = CMPLX(0.0,0.0)
IF(mpi%irank.EQ.0) THEN IF(mpi%irank.EQ.0) THEN
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,& CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,fermiEnergyTemp,l_qfix,inDen) 0,fermiEnergyTemp,l_qfix,inDen)
...@@ -191,21 +192,7 @@ CONTAINS ...@@ -191,21 +192,7 @@ CONTAINS
CALL timestop("Qfix") CALL timestop("Qfix")
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,& CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,inDen) 0,-1.0,0.0,.FALSE.,inDen)
IF (isDensityMatrixPresent().AND.atoms%n_u>0) THEN
CALL readDensityMatrix(input,atoms,inDen%mmpMat,l_error)
IF(l_error) CALL juDFT_error('Error in reading density matrix!',calledby='fleur')
ELSE IF (atoms%n_u>0) THEN
IF(input%ldauLinMix) THEN
input%ldauMixParam = 0.05
input%ldauSpinf = 1.0
END IF
inDen%mmpMat = CMPLX(0.0,0.0)
END IF
END IF END IF
#ifdef CPP_MPI
CALL MPI_BCAST(input%ldauMixParam,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
CALL MPI_BCAST(input%ldauSpinf,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
#endif
! Initialize and load inDen density (end) ! Initialize and load inDen density (end)
! Initialize potentials (start) ! Initialize potentials (start)
......
...@@ -143,8 +143,8 @@ ...@@ -143,8 +143,8 @@
CALL initWannierDefaults(wann) CALL initWannierDefaults(wann)
input%minDistance = 0.0 input%minDistance = 0.0
input%ldauLinMix = .TRUE. input%ldauLinMix = .FALSE.
input%ldauMixParam = -1.0 input%ldauMixParam = 0.05
input%ldauSpinf = 1.0 input%ldauSpinf = 1.0
kpts%ntet = 1 kpts%ntet = 1
kpts%numSpecialPoints = 1 kpts%numSpecialPoints = 1
......
...@@ -54,7 +54,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,& ...@@ -54,7 +54,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
INTEGER i,imap,js,mit,irecl INTEGER i,imap,js,mit,irecl
INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2 INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
INTEGER iofl,n_u_keep INTEGER iofl,n_u_keep
LOGICAL l_exist,l_ldaU LOGICAL l_exist,l_ldaU, l_densityMatrixPresent
!Local Arrays !Local Arrays
REAL dist(6) REAL dist(6)
...@@ -73,6 +73,8 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,& ...@@ -73,6 +73,8 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
! vol = omtil ! vol = omtil
! ENDIF ! ENDIF
l_densityMatrixPresent = ANY(inDen%mmpMat(:,:,:,:).NE.0.0)
!In systems without inversions symmetry the interstitial star- !In systems without inversions symmetry the interstitial star-
!coefficients are complex. Thus twice as many numbers have to be !coefficients are complex. Thus twice as many numbers have to be
!stored. !stored.
...@@ -104,7 +106,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,& ...@@ -104,7 +106,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
! LDA+U (start) ! LDA+U (start)
n_mmpTemp = inDen%mmpMat n_mmpTemp = inDen%mmpMat
n_u_keep=atoms%n_u n_u_keep=atoms%n_u
IF (isDensityMatrixPresent()) THEN IF (l_densityMatrixPresent) THEN
!In an LDA+U caclulation, also the density matrix is included in the !In an LDA+U caclulation, also the density matrix is included in the
!supervectors (sm,fsm) if no linear mixing is performed on it. !supervectors (sm,fsm) if no linear mixing is performed on it.
IF(input%ldauLinMix) THEN IF(input%ldauLinMix) THEN
...@@ -279,12 +281,9 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,& ...@@ -279,12 +281,9 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
END IF END IF
END IF END IF
!write out mixed density
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen)
IF (atoms%n_u > 0) THEN IF (atoms%n_u > 0) THEN
IF (.NOT.isDensityMatrixPresent()) THEN IF (.NOT.l_densityMatrixPresent) THEN
inDen%mmpMat(:,:,:,:) = outDen%mmpMat(:,:,:,:)
INQUIRE(file='broyd',exist=l_exist) INQUIRE(file='broyd',exist=l_exist)
IF (l_exist) THEN IF (l_exist) THEN
CALL system('rm broyd') CALL system('rm broyd')
...@@ -296,14 +295,12 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,& ...@@ -296,14 +295,12 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
PRINT *,"broyd.7 file deleted because new density matrix is created." PRINT *,"broyd.7 file deleted because new density matrix is created."
ENDIF ENDIF
END IF END IF
OPEN (69,file='n_mmp_mat',status='replace',form='formatted')
WRITE (69,'(7f20.13)') inDen%mmpMat(:,:,:,:)
IF (input%ldauLinMix) THEN
WRITE (69,'(2(a6,f5.3))') 'alpha=',input%ldauMixParam,'spinf=',input%ldauSpinf
END IF
CLOSE (69)
ENDIF ENDIF
!write out mixed density
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
1,results%last_distance,results%ef,.TRUE.,inDen)
inDen%iter = inDen%iter + 1 inDen%iter = inDen%iter + 1
IF (input%imix.GT.0) THEN IF (input%imix.GT.0) THEN
......
...@@ -79,9 +79,8 @@ CONTAINS ...@@ -79,9 +79,8 @@ CONTAINS
END DO END DO
CALL closeXMLElement('ldaUDensityMatrix') CALL closeXMLElement('ldaUDensityMatrix')
! exit subroutine if density matrix does nbot exist ! exit subroutine if density matrix does not exist
IF(.NOT.isDensityMatrixPresent()) THEN IF(.NOT.ANY(n_mmp_in(:,:,:,:).NE.0.0)) THEN
n_mmp_in = n_mmp_out
RETURN RETURN
END IF END IF
......
...@@ -92,8 +92,7 @@ SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,& ...@@ -92,8 +92,7 @@ SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,&
END DO END DO
! for LDA+U: flip density matrix ! for LDA+U: flip density matrix
IF (isDensityMatrixPresent().AND.atoms%n_u>0) THEN IF (ANY(den%mmpMat(:,:,:,:).NE.0.0).AND.atoms%n_u>0) THEN
CALL readDensityMatrix(input,atoms,den%mmpMat,l_error)
DO i_u = 1, atoms%n_u DO i_u = 1, atoms%n_u