Skip to content

GitLab

fleur
fleur
Commit 1583ab14 authored by Gregor Michalicek's avatar Gregor Michalicek
Browse files
Options

Some development for RDMFT

parent 8300a407
Hide whitespace changes
Inline Side-by-side
Showing with 164 additions and 185 deletions
eigen/eigen.F90
View file @ 1583ab14
......@@ -141,7 +141,7 @@ CONTAINS
stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
CALL timestop("Setup of H&S matrices")
IF(hybrid%l_hybrid) THEN
IF(hybrid%l_hybrid.OR.input%l_rdmft) THEN
DO i = 1, hmat%matsize1
DO j = 1, i
......@@ -160,6 +160,9 @@ CONTAINS
! Write overlap matrix smat to direct access file olap
print *,"Wrong overlap matrix used, fix this later"
CALL write_olap(smat,(jsp-1)*kpts%nkpt+nk) ! Note: At this moment this only works without MPI parallelization
END IF ! hybrid%l_hybrid.OR.input%l_rdmft
IF(hybrid%l_hybrid) THEN
PRINT *,"TODO"
! STOP "TODO"
PRINT *,"BASIS:", lapw%nv(jsp), atoms%nlotot
......
main/fleur.F90
View file @ 1583ab14
......@@ -57,6 +57,7 @@ CONTAINS
USE m_juDFT_time
USE m_calc_hybrid
USE m_rdmft
USE m_io_hybrid
USE m_wann_optional
USE m_wannier
USE m_bs_comfort
......@@ -214,6 +215,10 @@ CONTAINS
iter = 0
END IF
ENDIF
!RDMFT
IF(input%l_rdmft) THEN
CALL open_hybrid_io1(DIMENSION,sym%invs)
END IF
CALL reset_eig(eig_id,noco%l_soc) ! This has to be placed after the calc_hybrid call but before eigen
......@@ -349,7 +354,7 @@ CONTAINS
SELECT TYPE(xcpot)
TYPE IS(t_xcpot_inbuild)
CALL rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,dimension,&
sphhar,sym,field,vTot,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
END SELECT
END IF
......
rdmft/bfgs_b2.f90
View file @ 1583ab14
......@@ -23,8 +23,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
gradientCorrections,iStep,mixParam,l_converged,convCrit)
USE m_juDFT
IMPLICIT NONE
INTEGER, INTENT(IN) :: vecLength
......@@ -72,19 +70,14 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
REAL, ALLOCATABLE :: subspaceBasis(:,:)
REAL, ALLOCATABLE :: subspaceInvHessMat(:,:)
WRITE(*,*) 'bfgs_b2 - Point A'
WRITE(*,*) 'bfgs_b2 (input): parameters(:vecLength): ', parameters(:vecLength)
WRITE(*,*) 'bfgs_b2 (input): gradient(:vecLength): ', gradient(:vecLength)
l_converged = .FALSE.
! 0. Add last data pair to history iff it is not obviously destructive for the positive definiteness of the Hessian matrix
WRITE(2000,*) '=============================================='
WRITE(2000,'(a,8f15.10)') 'parameters(:): ', parameters(:)
WRITE(2000,'(a,8f15.10)') 'lastParameters(:): ', lastParameters(:)
WRITE(2000,'(a,8f15.10)') 'gradient(:): ', gradient(:)
WRITE(2000,'(a,8f15.10)') 'lastGradient(:): ', lastGradient(:)
IF(ANY(lastGradient(:).NE.0.0).OR.ANY(lastParameters(:).NE.0.0)) THEN
WRITE(*,*) 'bfgs_b2 - Point A - 1'
lastGradCorrection(:) = gradient(:) - lastGradient(:)
lastParamCorrection(:) = parameters(:) - lastParameters(:)
temp = ddot(vecLength,lastParamCorrection(:),1,lastGradCorrection(:),1)
......@@ -127,8 +120,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END DO
END IF
WRITE(*,*) 'bfgs_b2 - Point B'
DO i = 1, historyLength
stepIndex = MOD(iStep - historyLength + i - 1,maxHistoryLength) + 1
......@@ -140,95 +131,34 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
norm = ddot(vecLength,paramCorrections(:,stepIndex),1,tempVecA(:),1)
CALL dgemm('N','N',vecLength,vecLength,1,-1.0/norm,tempVecA(:),vecLength,tempVecA(:),1,0.0,tempMatB,vecLength)
WRITE(2002,*) '======================================='
WRITE(2002,'(a,2i7)') 'stepIndex, iStep: ', stepIndex, iStep
WRITE(2002,'(a,8f15.8)') 'paramCorrections(:,stepIndex): ', paramCorrections(:,stepIndex)
WRITE(2002,*) 'HessMat:'
DO j = 1, vecLength
WRITE(2002,'(8f15.10)') hessMat(:,j)
END DO
WRITE(2002,'(a,f15.8)') 'norm: ', norm
WRITE(2002,'(a,8f15.8)') 'tempVecA(:): ', tempVecA(:)
WRITE(2002,*) 'tempMatB:'
DO j = 1, vecLength
WRITE(2002,'(8f15.10)') tempMatB(:,j)
END DO
! first correction term
norm = ddot(vecLength,gradientCorrections(:,stepIndex),1,paramCorrections(:,stepIndex),1)
CALL dgemm('N','N',vecLength,vecLength,1,1.0/norm,gradientCorrections(:,stepIndex),vecLength,gradientCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
hessMat(:,:) = hessMat(:,:) + tempMatB(:,:)
WRITE(2002,*) 'HessMat:'
DO j = 1, vecLength
WRITE(2002,'(8f15.10)') hessMat(:,j)
END DO
! Now the inverse of the Hessian matrix
tempMatB(:,:) = 0.0
tempVecA(:) = 0.0
CALL dgemv('N',vecLength,vecLength,1.0,invHessMat(:,:),vecLength,gradientCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
norm = ddot(vecLength,paramCorrections(:,stepIndex),1,gradientCorrections(:,stepIndex),1)
WRITE(2003,*) '======================================='
WRITE(2003,'(a,2i7)') 'stepIndex, iStep: ', stepIndex, iStep
WRITE(2003,'(a,8f15.8)') 'paramCorrections(:,stepIndex): ', paramCorrections(:,stepIndex)
WRITE(2003,'(a,8f15.8)') 'gradientCorrections(:,stepIndex): ', gradientCorrections(:,stepIndex)
WRITE(2003,*) 'invHessMat:'
DO j = 1, vecLength
WRITE(2003,'(8f15.10)') invHessMat(:,j)
END DO
WRITE(2003,'(a,f15.8)') 'norm: ', norm
WRITE(2003,'(a,8f15.8)') 'tempVecA(:): ', tempVecA(:)
! second correction term
CALL dgemm('N','N',vecLength,vecLength,1,1.0,tempVecA(:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempMatB,vecLength)
WRITE(2003,*) 'tempMatB:'
DO j = 1, vecLength
WRITE(2003,'(8f15.10)') tempMatB(:,j)
END DO
tempMatB(:,:) = tempMatB(:,:) + TRANSPOSE(tempMatB(:,:))
tempMatB(:,:) = -tempMatB(:,:) / norm
WRITE(2003,*) 'tempMatB:'
DO j = 1, vecLength
WRITE(2003,'(8f15.10)') tempMatB(:,j)
END DO
! first correction term
temp = norm + ddot(vecLength,gradientCorrections(:,stepIndex),1,tempVecA(:),1)
CALL dgemm('N','N',vecLength,vecLength,1,temp/(norm*norm),paramCorrections(:,stepIndex),vecLength,paramCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
invHessMat(:,:) = invHessMat(:,:) + tempMatB(:,:)
WRITE(2003,*) 'invHessMat:'
DO j = 1, vecLength
WRITE(2003,'(8f15.10)') invHessMat(:,j)
END DO
!Check whether invHessMat is the inverse of hessMat (for debugging only)
tempMatB(:,:) = 0.0
CALL dgemm('N','N',vecLength,vecLength,vecLength,1.0,hessMat(:,:),vecLength,invHessMat(:,:),vecLength,0.0,tempMatB,vecLength)
WRITE(2003,*) 'identityMatrix (hessMat*invHessMat):'
DO j = 1, vecLength
WRITE(2003,'(8f15.10)') tempMatB(:,j)
END DO
END DO
WRITE(2001,*) '========================================================='
WRITE(2001,*) 'HessMat:'
DO i = 1, vecLength
WRITE(2001,'(8f15.10)') hessMat(:,i)
END DO
WRITE(2001,*) 'invHessMat:'
DO i = 1, vecLength
WRITE(2001,'(8f15.10)') invHessMat(:,i)
END DO
WRITE(*,*) 'bfgs_b2 - Point C'
! 2. Determine next point based on the gradient, the inverse of the Hessian, and the box constraints
freeParams(:) = .TRUE.
......@@ -238,10 +168,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
converged = .FALSE.
DO WHILE (.NOT.converged)
WRITE(*,*) 'bfgs_b2 - Point D'
WRITE(1999,'(a,8f15.10)') 'prelimNextPoint(:vecLength): ', prelimNextPoint(:vecLength)
WRITE(1999,'(a,8f15.10)') 'prelimGradient(:vecLength): ', prelimGradient(:vecLength)
! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
numFixedParams = 0
DO i = 1, vecLength
......@@ -258,8 +184,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END DO
numFreeParams = vecLength - numFixedParams
WRITE(*,*) 'freeParams(:vecLength): ', freeParams(:vecLength)
! Perform basis transformation to obtain inverse Hessian matrix for the free parameter subspace
ALLOCATE (subspaceBasis(vecLength,numFreeParams))
......@@ -274,51 +198,22 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END IF
END DO
WRITE(*,*) 'subspaceBasis:'
DO i = 1, numFreeParams
WRITE(*,'(8f10.5)') subspaceBasis(:,i)
END DO
WRITE(*,*) 'bfgs_b2 - Point E'
ALLOCATE (tempMatC(vecLength,numFreeParams))
ALLOCATE (subspaceInvHessMat(numFreeParams,numFreeParams))
tempMatC(:,:) = 0.0
subspaceInvHessMat(:,:) = 0.0
CALL dgemm('N','N',vecLength,numFreeParams,vecLength,1.0,invHessMat(:,:),vecLength,subspaceBasis(:,:),vecLength,0.0,tempMatC(:,:),vecLength)
WRITE(*,*) 'tempMatC:'
DO i = 1, numFreeParams
WRITE(*,'(8f10.5)') tempMatC(:,i)
END DO
CALL dgemm('N','N',vecLength,numFreeParams,vecLength,1.0,invHessMat(:,:),vecLength,subspaceBasis(:,:),vecLength,0.0,tempMatC(:,:),vecLength)
CALL dgemm('N','N',numFreeParams,numFreeParams,vecLength,1.0,TRANSPOSE(subspaceBasis(:,:)),numFreeParams,tempMatC(:,:),vecLength,0.0,subspaceInvHessMat(:,:),numFreeParams)
DEALLOCATE(tempMatC)
DEALLOCATE(subspaceBasis)
WRITE(*,*) 'bfgs_b2 - Point F'
! Calculate next point (ignoring constraints)
WRITE(*,*) 'vecLength: ', vecLength
WRITE(*,*) 'numFreeParams: ', numFreeParams
WRITE(*,*) 'SHAPE(subspaceInvHessMat): ', SHAPE(subspaceInvHessMat)
WRITE(*,*) 'SHAPE(subspaceGrad): ', SHAPE(subspaceGrad)
WRITE(*,*) 'SHAPE(tempVecA): ', SHAPE(tempVecA)
WRITE(*,*) 'subspaceGrad(:numFreeParams): ', subspaceGrad(:numFreeParams)
WRITE(*,*) 'subspaceInvHessMat:'
DO i = 1, numFreeParams
WRITE(*,'(8f10.5)') subspaceInvHessMat(:,i)
END DO
tempVecA(:) = 0.0
CALL dgemv('N',numFreeParams,numFreeParams,1.0,subspaceInvHessMat(:,:),numFreeParams,-subspaceGrad(:numFreeParams),1,0.0,tempVecA(:numFreeParams),1)
WRITE(*,*) 'tempVecA(:numFreeParams): ', tempVecA(:numFreeParams)
WRITE(*,*) 'bfgs_b2 - Point F - 0'
j = 1
DO i = 1, vecLength
IF(freeParams(i)) THEN
......@@ -327,12 +222,8 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END IF
END DO
WRITE(*,*) 'bfgs_b2 - Point F - 1'
DEALLOCATE(subspaceInvHessMat)
WRITE(*,*) 'bfgs_b2 - Point F - 2'
! Constrain next point to box and check convergence
converged = .TRUE.
DO i = 1, vecLength
......@@ -352,7 +243,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END IF
END DO
WRITE(*,*) 'bfgs_b2 - Point F - 3'
! If not converged calculate prelimGradient for next point according to quadratic model
IF(.NOT.converged) THEN
......@@ -362,11 +252,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
prelimGradient(:) = prelimGradient(:) + diffGrad(:)
END IF
WRITE(*,*) 'prelimNextPoint(:vecLength): ', prelimNextPoint(:vecLength)
WRITE(*,*) 'prelimGradient(:vecLength): ', prelimGradient(:vecLength)
WRITE(*,*) 'bfgs_b2 - Point G'
END DO
! 3. Project prelimNextPoint onto the allowed subspace given by the equality constraints
......@@ -374,8 +259,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
WRITE(*,*) 'bfgs_b2 - Point H'
converged = .FALSE.
DO WHILE (.NOT.converged)
......@@ -429,8 +312,6 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
END DO
END DO
WRITE(*,*) 'bfgs_b2 - Point I'
! 4. Update lastParameters, lastGradient, parameters
lastParameters(:) = parameters(:)
......@@ -441,7 +322,7 @@ SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints
DEALLOCATE(invHessMat)
DEALLOCATE(hessMat)
WRITE(*,*) 'parameters(:vecLength): ', parameters(:vecLength)
WRITE(*,*) 'bfgs_b2 (output): parameters(:vecLength): ', parameters(:vecLength)
END SUBROUTINE bfgs_b2
......
rdmft/rdmft.F90
View file @ 1583ab14
......@@ -9,11 +9,12 @@ MODULE m_rdmft
CONTAINS
SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,dimension,&
sphhar,sym,field,vTot,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
USE m_types
USE m_juDFT
USE m_constants
USE m_eig66_io
#ifndef CPP_OLDINTEL
USE m_cdnval
USE m_cdngen
......@@ -32,6 +33,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
USE m_symm_hf
USE m_exchange_valence_hf
USE m_exchange_core
USE m_symmetrizeh
USE m_bfgs_b2
#ifdef CPP_MPI
......@@ -56,6 +58,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_field), INTENT(INOUT) :: field
TYPE(t_potden), INTENT(INOUT) :: vTot
TYPE(t_potden), INTENT(INOUT) :: vCoul
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_noco), INTENT(INOUT) :: noco
TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
......@@ -69,19 +72,19 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
#ifndef CPP_OLDINTEL
TYPE(t_cdnvalJob) :: cdnvalJob
TYPE(t_potden) :: singleStateDen, overallDen, overallVCoul
TYPE(t_potden) :: singleStateDen, overallDen, overallVCoul, vTotTemp
TYPE(t_regionCharges) :: regCharges
TYPE(t_dos) :: dos
TYPE(t_moments) :: moments
TYPE(t_mat) :: exMat
TYPE(t_mat) :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
TYPE(t_lapw) :: lapw
TYPE(t_hybdat) :: hybdat
INTEGER :: ikpt, iBand, jkpt, jBand
INTEGER :: jspin, jspmax, jsp, isp, ispin
INTEGER :: ikpt, iBand, jkpt, jBand, iAtom, i, na, itype, lh, j
INTEGER :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
INTEGER :: nsymop, nkpt_EIBZ, ikptf, iterHF, mnobd
INTEGER :: iState, iStep, numStates, maxHistoryLength, numRelevantStates
REAL :: fix, potDenInt, fermiEnergyTemp
REAL :: rdmftFunctionalValue, occStateI
REAL :: fix, potDenInt, fermiEnergyTemp, spinDegenFac
REAL :: rdmftFunctionalValue, occStateI, gradSum
REAL :: exchangeTerm, lagrangeMultiplier, equalityCriterion
REAL :: mixParam, convCrit, rdmftEnergy
REAL :: sumOcc, tempOcc, addCharge, subCharge, addChargeWeight, subChargeWeight
......@@ -125,10 +128,14 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
#endif
#ifndef CPP_OLDINTEL
WRITE(*,*) 'entered RDMFT subroutine'
! General initializations
mixParam = 0.0001
convCrit = 1.0e-8
convCrit = 1.0e-6
lagrangeMultiplier = 0.1 !results%ef
spinDegenFac = 2.0 / input%jspins ! This factor is used to compensate the missing second spin in non-spinpolarized calculations
neigTemp(:,:) = results%neig(:,:)
......@@ -229,12 +236,10 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
END IF
results%w_iks(iBand,ikpt,jsp) = occupationVec(iState) * kpts%wtkpt(ikpt)
tempOcc = tempOcc + occupationVec(iState) * kpts%wtkpt(ikpt)
WRITE(*,*) 'occ: ', iState, occupationVec(iState)
END DO
END DO
END DO
DEALLOCATE(occupationVec)
WRITE(*,*) 'sumOcc, tempOcc: ', sumOcc, tempOcc
! Some more initializations
......@@ -250,11 +255,25 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL overallDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL overallVCoul%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),1))
ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),input%jspins))
DO jspin = 1, input%jspins
CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
END DO
CALL vTotTemp%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTTOT)
vTotTemp = vTot
DO jsp = 1,SIZE(vTot%mt,4)
DO iAtom = 1,atoms%ntype
vTotTemp%mt(:atoms%jri(iAtom),0,iAtom,jsp) = sfp_const * vTot%mt(:atoms%jri(iAtom),0,iAtom,jsp) / atoms%rmsh(:atoms%jri(iAtom),iAtom)
END DO
END DO
vCoul%pw_w = CMPLX(0.0,0.0)
DO jspin = 1, input%jspins
CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin),stars%ufft)
END DO
vTotSSDen = 0.0
! Calculate all single state densities
......@@ -279,7 +298,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
cdnvalJob%nStart = iBand
cdnvalJob%nEnd = iBand
cdnvalJob%weights = 0.0
cdnvalJob%weights(1,ikpt) = 1.0
cdnvalJob%weights(1,ikpt) = spinDegenFac ! Note the doubling of the weight for non-spinpolarized calculations.
ELSE
cdnvalJob%noccbd = 0
cdnvalJob%nStart = 1
......@@ -304,19 +323,18 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
#ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,singleStateDen)
#endif
! For each state calculate Integral over KS effective potential times single state density
potDenInt = 0.0
CALL int_nv(jsp,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vTot,singleStateDen,potDenInt)
CALL int_nv(jsp,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vTotTemp,singleStateDen,potDenInt)
vTotSSDen(iBand,ikpt,jsp) = potDenInt
END DO
END DO
END DO
WRITE(*,*) 'Point A reached'
! Initializations for exchange contributions
WRITE(*,*) 'RDMFT: HF initializations start'
IF(ALLOCATED(hybrid%ne_eig)) DEALLOCATE(hybrid%ne_eig)
IF(ALLOCATED(hybrid%nbands)) DEALLOCATE(hybrid%nbands)
IF(ALLOCATED(hybrid%nobd)) DEALLOCATE(hybrid%nobd)
......@@ -330,21 +348,17 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
iterHF = 0
hybrid%l_calhf = .TRUE.
CALL open_hybrid_io1(DIMENSION,sym%invs)
! CALL open_hybrid_io1(DIMENSION,sym%invs)
CALL mixedbasis(atoms,kpts,dimension,input,cell,sym,xcpot,hybrid,enpara,mpi,vTot,l_restart)
WRITE(*,*) 'Point B reached'
CALL open_hybrid_io2(hybrid,DIMENSION,atoms,sym%invs)
CALL coulombmatrix(mpi,atoms,kpts,cell,sym,hybrid,xcpot,l_restart)
WRITE(*,*) 'Point C reached'
CALL hf_init(hybrid,kpts,atoms,input,DIMENSION,hybdat,sym%invs)
WRITE(*,*) 'Point D reached'
WRITE(*,*) 'RDMFT: HF initializations end'
ALLOCATE(parent(kpts%nkptf))
ALLOCATE(exDiag(dimension%neigd,ikpt,input%jspins))
......@@ -352,7 +366,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
ALLOCATE(lastParameters(numStates+1))
lastGradient = 0.0
lastParameters = 0.0
maxHistoryLength = 5!7
maxHistoryLength = 17!7
ALLOCATE(gradientCorrections(numStates+1,maxHistoryLength))
ALLOCATE(paramCorrections(numStates+1,maxHistoryLength))
gradientCorrections = 0.0
......@@ -364,7 +378,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
converged = .FALSE.
DO WHILE (.NOT.converged)
WRITE(*,*) 'Point E reached'
WRITE(*,*) 'RDMFT: convergence loop start'
DO jspin = 1, input%jspins
DO ikpt = 1,kpts%nkpt
......@@ -383,8 +397,6 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
sphhar,sym,vTot,oneD,cdnvalJob,overallDen,regCharges,dos,results,moments)
END DO
WRITE(*,*) 'overallDen%pw(1,1)', overallDen%pw(1,1)
CALL cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,vTot,overallDen,moments,results)
IF (mpi%irank.EQ.0) THEN
......@@ -394,8 +406,6 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallDen)
#endif
WRITE(*,*) 'Point E-1 reached'
! Calculate Coulomb potential for overall density (+including external potential)
CALL overallDen%sum_both_spin()!workden)
CALL overallVCoul%resetPotDen()
......@@ -407,8 +417,6 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
#endif
WRITE(*,*) 'Point E-2 reached'
overallVCoulSSDen = 0.0
DO jspin = 1, input%jspins
jsp = MERGE(1,jspin,noco%l_noco)
......@@ -428,14 +436,12 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
! For each state calculate integral over Coulomb potential times single state density
potDenInt = 0.0
CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,overallVCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
END DO
END DO
END DO
WRITE(*,*) 'Point E-3 reached'
! Construct exchange matrix diagonal
exDiag = 0.0
......@@ -460,9 +466,13 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
IF(ALLOCATED(hybdat%prod)) DEALLOCATE (hybdat%prod)
IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)
results%neig(:,:) = neigTemp(:,:)
CALL HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jspin,enpara,eig_id,&
hybdat,iterHF,sym%invs,vTot%mt(:,0,:,:),eig_irr)
results%neig(:,:) = highestState(:,:) + 1
mnobd = MAXVAL(hybrid%nobd)
DO ikpt = 1,kpts%nkpt
......@@ -479,40 +489,117 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
eig_irr,results,parent,pointer_EIBZ,n_q,wl_iks,iterHF,xcpot,noco,nsest,indx_sest,&
mpi,exMat)
CALL exchange_vccv1(ikpt,atoms,hybrid,hybdat,dimension,jspin,lapw,nsymop,nsest,indx_sest,mpi,1.0,results,exMat)
!Start of workaround for increased functionality of symmetrizeh (call it))
nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
CALL olap%alloc(sym%invs,nbasfcn)
CALL read_olap(olap, kpts%nkpt*(jspin-1)+ikpt)
IF (olap%l_real) THEN
DO i = 1, nbasfcn
DO j = 1, i
olap%data_r(i,j) = olap%data_r(j,i)
END DO
END DO
ELSE
DO i = 1, nbasfcn
DO j = 1, i
olap%data_c(i,j) = CONJG(olap%data_c(j,i))
END DO
END DO
olap%data_c = conjg(olap%data_c)
END IF
CALL zMat%init(olap%l_real,nbasfcn,dimension%neigd)
CALL read_eig(eig_id,ikpt,jspin,n_start=1,n_end=hybrid%nbands(ikpt),neig=nbands,zmat=zMat)
! CALL read_z(zMat,kpts%nkpt*(jspin-1)+ikpt)
zMat%matsize2 = hybrid%nbands(ikpt) ! reduce "visible matsize" for the following computations
CALL olap%multiply(zMat,trafo)
CALL invtrafo%alloc(olap%l_real,hybrid%nbands(ikpt),nbasfcn)
CALL trafo%TRANSPOSE(invtrafo)
IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)
DO i = 1, hybrid%nbands(ikpt)
DO j = 1, i-1
IF (exMat%l_real) THEN
exMat%data_r(i,j)=exMat%data_r(j,i)
ELSE
exMat%data_c(i,j)=conjg(exMat%data_c(j,i))
END IF
END DO
END DO
CALL exMat%multiply(invtrafo,tmpMat)
CALL trafo%multiply(tmpMat,exMatLAPW)
CALL symmetrizeh(atoms,kpts%bkf(:,ikpt),dimension,jspin,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,exMatLAPW)
IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c)
zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)
CALL exMatLAPW%multiply(zMat,tmpMat)
DO iBand = 1, highestState(ikpt,jsp)
IF (exMat%l_real) THEN
exDiag(iBand,ikpt,jspin) = exMat%data_r(iBand,iBand)
IF (zMat%l_real) THEN
exDiag(iBand,ikpt,jspin) = dot_product(zMat%data_r(:zMat%matsize1,iband),tmpMat%data_r(:,iband))
ELSE
exDiag(iBand,ikpt,jspin) = REAL(exMat%data_c(iBand,iBand))
exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
END IF
END DO
!End of workaround for increased functionality of symmetrizeh (call it))
! DO iBand = 1, highestState(ikpt,jsp)