Commit 6e605703 authored by Gregor Michalicek's avatar Gregor Michalicek

Move reading of density matrix out of SCF loop

Up to now this is only a partial solution. There is still reading
of the n_mmp_mat file in the SCF loop to read in the optional mixing
parameters.
parent 7172bf70
......@@ -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,&
inIter,outDen)
inIter,inDen,outDen)
!*****************************************************
! Charge density generator
......@@ -63,7 +63,7 @@ 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) :: outDen
TYPE(t_potden),INTENT(INOUT) :: inDen,outDen
!Scalar Arguments
INTEGER, INTENT (IN) :: eig_id
......@@ -162,7 +162,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
END DO
! lda+u
IF ((atoms%n_u.GT.0).and.(mpi%irank.EQ.0)) CALL u_mix(atoms,input%jspins,outDen%mmpMat)
IF ((atoms%n_u.GT.0).and.(mpi%irank.EQ.0)) CALL u_mix(atoms,input%jspins,inDen%mmpMat,outDen%mmpMat)
!+t3e
IF (mpi%irank.EQ.0) THEN
......
......@@ -171,7 +171,7 @@ CONTAINS
results%last_distance = -1.0
IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')
! Initialize and load inDen density without density matrix(start)
! Initialize and load inDen density (start)
CALL inDen%init(stars,atoms,sphhar,vacuum,oneD,DIMENSION%jspd,.FALSE.,POTDEN_TYPE_DEN)
IF (noco%l_noco) THEN
ALLOCATE (inDen%cdom(stars%ng3),inDen%cdomvz(vacuum%nmzd,2))
......@@ -181,6 +181,7 @@ CONTAINS
ALLOCATE (inDen%cdom(1),inDen%cdomvz(1,1),inDen%cdomvxy(1,1,1))
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
END IF
ALLOCATE (inDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,atoms%n_u),input%jspins))
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%iter,inDen%mt,inDen%pw,inDen%vacz,inDen%vacxy,&
......@@ -192,9 +193,14 @@ CONTAINS
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,-1.0,0.0,.FALSE.,inDen%iter,inDen%mt,inDen%pw,inDen%vacz,inDen%vacxy,inDen%cdom,&
inDen%cdomvz,inDen%cdomvxy)
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
inDen%mmpMat = CMPLX(0.0,0.0)
END IF
END IF
ALLOCATE (inDen%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MAX(1,atoms%n_u),input%jspins))
! Initialize and load inDen density without density matrix(end)
! Initialize and load inDen density (end)
! Initialize potentials (start)
CALL vTot%init(stars,atoms,sphhar,vacuum,oneD,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
......@@ -267,15 +273,6 @@ CONTAINS
END IF
END IF
! Initialize and load inDen density matrix and broadcast inDen(start)
IF(mpi%irank.EQ.0) THEN
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
inDen%mmpMat = CMPLX(0.0,0.0)
END IF
END IF
#ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen)
#endif
......@@ -595,7 +592,7 @@ CONTAINS
CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
DIMENSION,kpts,atoms,sphhar,stars,sym,obsolete,&
enpara,cell,noco,jij,vTot,results,oneD,coreSpecInput,&
inDen%iter,outDen)
inDen%iter,inDen,outDen)
IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) dimension%neigd=dimension%neigd/2
!+t3e
......
......@@ -60,6 +60,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
REAL dist(6)
REAL, ALLOCATABLE :: sm(:), fsm(:)
CHARACTER(LEN=20) :: attributes(2)
COMPLEX :: n_mmpTemp(-3:3,-3:3,atoms%n_u,input%jspins)
!External functions
REAL CPP_BLAS_sdot
......@@ -100,6 +101,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
2*vacuum%nmzd*vacuum%nvac
END IF
n_mmpTemp = inDen%mmpMat
n_u_keep=atoms%n_u
INQUIRE (file='n_mmp_mat',exist=l_ldaU)
IF (l_ldaU) THEN
......@@ -118,8 +120,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
IF ( MOD(i,14*input%jspins) == 1 ) THEN ! was already mixed in u_mix
atoms%n_u = 0
ELSE IF ( MOD(i,28*input%jspins)== 0 ) THEN ! mix here
atoms%n_u = i / (28 * input%jspins ) ! calculate number of U parameters
mmap = mmap + 7 * i / 2 ! add 7*7 complex numbers per atoms%n_u
mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
ELSE
CALL juDFT_error("strange n_mmp_mat-file...",calledby ="mix")
END IF
......@@ -291,6 +292,10 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
CLOSE (59)
END IF
IF(atoms%n_u.NE.n_u_keep) THEN
inDen%mmpMat = n_mmpTemp
END IF
atoms%n_u=n_u_keep
7900 FORMAT (/,'----> distance of charge densities for spin ',i2,' it=',i5,':',f13.6,' me/bohr**3')
......
......@@ -12,7 +12,7 @@ MODULE m_umix
! --------------------------------------------------------
! Extension to multiple U per atom type by G.M. 2017
CONTAINS
SUBROUTINE u_mix(atoms,jspins,n_mmp_new)
SUBROUTINE u_mix(atoms,jspins,n_mmp_in,n_mmp_out)
USE m_types
USE m_nmat_rot
......@@ -23,7 +23,8 @@ CONTAINS
IMPLICIT NONE
TYPE(t_atoms),INTENT(IN) :: atoms
INTEGER, INTENT (IN) :: jspins
COMPLEX, INTENT (INOUT) :: n_mmp_new(-3:3,-3:3,atoms%n_u,jspins)
COMPLEX, INTENT (INOUT) :: n_mmp_out(-3:3,-3:3,atoms%n_u,jspins)
COMPLEX, INTENT (INOUT) :: n_mmp_in (-3:3,-3:3,atoms%n_u,jspins)
!
! ... Locals ...
INTEGER j,k,iofl,l,itype,ios,i_u,jsp,lty(atoms%n_u)
......@@ -31,7 +32,7 @@ CONTAINS
REAL theta(atoms%n_u),phi(atoms%n_u),zero(atoms%n_u)
LOGICAL n_exist
CHARACTER(LEN=20) :: attributes(6)
COMPLEX,ALLOCATABLE :: n_mmp(:,:,:,:),n_mmp_old(:,:,:,:)
COMPLEX,ALLOCATABLE :: n_mmp(:,:,:,:)
!
! check for possible rotation of n_mmp
!
......@@ -51,10 +52,10 @@ CONTAINS
END DO
CLOSE (68)
zero = 0.0
CALL nmat_rot(zero,-theta,-phi,3,atoms%n_u,jspins,lty,n_mmp_new)
CALL nmat_rot(zero,-theta,-phi,3,atoms%n_u,jspins,lty,n_mmp_out)
END IF
! Write out n_mmp_new to out.xml file
! Write out n_mmp_out to out.xml file
CALL openXMLElementNoAttributes('ldaUDensityMatrix')
DO jsp = 1, jspins
......@@ -72,7 +73,7 @@ CONTAINS
WRITE(attributes(6),'(f15.8)') jParam
CALL writeXMLElementMatrixPoly('densityMatrixFor',&
(/'spin ','atomType','uIndex ','l ','U ','J '/),&
attributes,n_mmp_new(-l:l,-l:l,i_u,jsp))
attributes,n_mmp_out(-l:l,-l:l,i_u,jsp))
END DO
END DO
CALL closeXMLElement('ldaUDensityMatrix')
......@@ -85,9 +86,9 @@ CONTAINS
IF (n_exist) THEN
ALLOCATE (n_mmp_old(-3:3,-3:3,atoms%n_u,jspins))
ALLOCATE ( n_mmp(-3:3,-3:3,atoms%n_u,jspins))
READ (69,9000) n_mmp_old(:,:,:,:)
ALLOCATE ( n_mmp(-3:3,-3:3,atoms%n_u,jspins))
READ (69,9000) n_mmp(:,:,:,:)
n_mmp = CMPLX(0.0,0.0)
READ (69,'(2(6x,f5.3))',IOSTAT=iofl) alpha,spinf
IF ( iofl == 0 ) THEN
......@@ -100,8 +101,8 @@ CONTAINS
DO i_u = 1, atoms%n_u
DO j = -3,3
DO k = -3,3
sum1 = sum1 + ABS(n_mmp_new(k,j,i_u,1) - n_mmp_old(k,j,i_u,1))
n_mmp(k,j,i_u,1) = alpha * n_mmp_new(k,j,i_u,1) + (1.0-alpha) * n_mmp_old(k,j,i_u,1)
sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
n_mmp(k,j,i_u,1) = alpha * n_mmp_out(k,j,i_u,1) + (1.0-alpha) * n_mmp_in(k,j,i_u,1)
END DO
END DO
END DO
......@@ -113,18 +114,18 @@ CONTAINS
DO i_u = 1,atoms%n_u
DO j = -3,3
DO k = -3,3
sum1 = sum1 + ABS(n_mmp_new(k,j,i_u,1) - n_mmp_old(k,j,i_u,1))
sum2 = sum2 + ABS(n_mmp_new(k,j,i_u,2) - n_mmp_old(k,j,i_u,2))
n_mmp(k,j,i_u,1) = gam * n_mmp_new(k,j,i_u,1) + &
(1.0-gam) * n_mmp_old(k,j,i_u,1) + &
del * n_mmp_new(k,j,i_u,2) - &
del * n_mmp_old(k,j,i_u,2)
n_mmp(k,j,i_u,2) = gam * n_mmp_new(k,j,i_u,2) + &
(1.0-gam) * n_mmp_old(k,j,i_u,2) + &
del * n_mmp_new(k,j,i_u,1) - &
del * n_mmp_old(k,j,i_u,1)
sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
sum2 = sum2 + ABS(n_mmp_out(k,j,i_u,2) - n_mmp_in(k,j,i_u,2))
n_mmp(k,j,i_u,1) = gam * n_mmp_out(k,j,i_u,1) + &
(1.0-gam) * n_mmp_in (k,j,i_u,1) + &
del * n_mmp_out(k,j,i_u,2) - &
del * n_mmp_in (k,j,i_u,2)
n_mmp(k,j,i_u,2) = gam * n_mmp_out(k,j,i_u,2) + &
(1.0-gam) * n_mmp_in (k,j,i_u,2) + &
del * n_mmp_out(k,j,i_u,1) - &
del * n_mmp_in (k,j,i_u,1)
END DO
END DO
END DO
......@@ -133,6 +134,7 @@ CONTAINS
ENDIF
WRITE (69,9000) n_mmp
WRITE (69,'(2(a6,f5.3))') 'alpha=',alpha,'spinf=',spinf
n_mmp_in = n_mmp
ELSEIF (iofl > 0 ) THEN
!
......@@ -149,7 +151,7 @@ CONTAINS
DO i_u = 1, atoms%n_u
DO j = -3,3
DO k = -3,3
sum1 = sum1 + ABS(n_mmp_new(k,j,i_u,1) - n_mmp_old(k,j,i_u,1))
sum1 = sum1 + ABS(n_mmp_out(k,j,i_u,1) - n_mmp_in(k,j,i_u,1))
END DO
END DO
END DO
......@@ -161,30 +163,31 @@ CONTAINS
DO i_u = 1, atoms%n_u
DO j = -3,3
DO k = -3,3
sum2 = sum2 + ABS(n_mmp_new(k,j,i_u,2) - n_mmp_old(k,j,i_u,2))
sum2 = sum2 + ABS(n_mmp_out(k,j,i_u,2) - n_mmp_in(k,j,i_u,2))
END DO
END DO
END DO
DO j=-3,3
WRITE(6,'(14f12.6)') (n_mmp_old(k,j,1,2),k=-3,3)
WRITE(6,'(14f12.6)') (n_mmp_in(k,j,1,2),k=-3,3)
END DO
WRITE (6,'(a23,f12.6)') 'n_mmp distance spin 2 =',sum2
DO j=-3,3
WRITE(6,'(14f12.6)') (n_mmp_new(k,j,1,2),k=-3,3)
WRITE(6,'(14f12.6)') (n_mmp_out(k,j,1,2),k=-3,3)
END DO
END IF
REWIND(69)
WRITE (69,9000) n_mmp_old
WRITE (69,9000) n_mmp_new
WRITE (69,9000) n_mmp_in
WRITE (69,9000) n_mmp_out
END IF ! iofl == 0
DEALLOCATE ( n_mmp_old,n_mmp )
DEALLOCATE (n_mmp)
ELSE
!
! first time with lda+u; write new n_mmp
!
WRITE (69,9000) n_mmp_new
WRITE (69,9000) n_mmp_out
WRITE (69,'(2(a6,f5.3))') 'alpha=',0.05,'spinf=',1.0
n_mmp_in = n_mmp_out
END IF
9000 FORMAT(7f20.13)
......
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