Commit 0720c364 authored by Daniel Wortmann's avatar Daniel Wortmann

Added functionality to do a fixed total moment calculation.

parent 62f0cb52
......@@ -6,7 +6,7 @@ MODULE m_fergwt
! c.l.fu
!*****************************************************************
CONTAINS
SUBROUTINE fergwt(kpts,input,mpi, ne,eig, results)
SUBROUTINE fergwt(kpts,input,mpi, ne,eig, ef,w_iks,seigv)
USE m_constants
USE m_types
......@@ -15,7 +15,8 @@ CONTAINS
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_input),INTENT(IN) :: input
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_results),INTENT(INOUT):: results
REAL, INTENT(INOUT) :: ef,seigv
REAL,INTENT(INOUT) :: w_iks(:,:,:)
! ..
! ..
! .. Array Arguments ..
......@@ -45,7 +46,7 @@ CONTAINS
nbnd = ne(k,jspin)
DO i = 1,nbnd
en = eig(i,k,jspin)
de = (en-results%ef)/input%delgau
de = (en-ef)/input%delgau
wt = 2.0
IF (de.GT.eup) wt = 0.0
IF (de.GE.elow .AND. de.LE.eup) THEN
......@@ -56,7 +57,7 @@ CONTAINS
END IF
END IF
s = s + wt*wtk
results%w_iks(i,k,jspin) = wt/2.
w_iks(i,k,jspin) = wt/2.
ENDDO
ENDDO
ENDDO
......@@ -65,22 +66,22 @@ CONTAINS
IF (ABS(zcdiff).LT.eps) EXIT conv_loop
IF (ifl.EQ.0) THEN
ifl = 1
ef0 = results%ef
results%ef = results%ef + 0.003
ef0 = ef
ef = ef + 0.003
s0 = s
ELSE
fac = (s0-s)/ (input%zelec-s)
IF (ABS(fac).LT.1.0e-1) THEN
ef0 = results%ef
ef0 = ef
s0 = s
IF (zcdiff.GE.zero) THEN
results%ef = results%ef + 0.003
ef = ef + 0.003
ELSE
results%ef = results%ef - 0.003
ef = ef - 0.003
END IF
ELSE
ef1 = results%ef
results%ef = results%ef + (ef0-results%ef)/fac
ef1 = ef
ef = ef + (ef0-ef)/fac
ef0 = ef1
s0 = s
END IF
......@@ -90,10 +91,10 @@ CONTAINS
IF ( mpi%irank == 0 ) WRITE (6,FMT=8000) eps
8000 FORMAT (10x,'warning: eps has been increased to',e12.5)
ENDDO conv_loop
workf = -hartree_to_ev_const*results%ef
workf = -hartree_to_ev_const*ef
IF ( mpi%irank == 0 ) THEN
WRITE (16,FMT=8010) results%ef,workf,s
WRITE (6,FMT=8010) results%ef,workf,s
WRITE (16,FMT=8010) ef,workf,s
WRITE (6,FMT=8010) ef,workf,s
END IF
8010 FORMAT (/,10x,'fermi energy=',f10.5,' har',3x,'work function=',&
f10.5,' ev',/,10x,'number of valence electrons=',f10.5)
......@@ -108,22 +109,21 @@ CONTAINS
nbnd = ne(k,jspin)
IF ( mpi%irank == 0 ) WRITE (6,FMT=8030) k
8030 FORMAT (/,5x,'k-point=',i5,/)
results%w_iks(:,k,jspin) = kpts%wtkpt(k)*results%w_iks(:,k,jspin)
IF ( mpi%irank == 0) WRITE (6,FMT=8040) (results%w_iks(i,k,jspin),i=1,nbnd)
w_iks(:,k,jspin) = kpts%wtkpt(k)*w_iks(:,k,jspin)
IF ( mpi%irank == 0) WRITE (6,FMT=8040) (w_iks(i,k,jspin),i=1,nbnd)
8040 FORMAT (5x,16f6.3)
ENDDO
ENDDO
s1 = 0.
s2 = 0.
results%seigv = 0.
DO jspin = 1,input%jspins
s = 0.
DO k = 1,kpts%nkpt
DO i = 1,ne(k,jspin)
s = s + results%w_iks(i,k,jspin)
results%seigv = results%seigv + results%w_iks(i,k,jspin)*eig(i,k,jspin)
s = s + w_iks(i,k,jspin)
seigv = seigv + w_iks(i,k,jspin)*eig(i,k,jspin)
en = eig(i,k,jspin)
de = (en-results%ef)/input%delgau
de = (en-ef)/input%delgau
! ---> correction term
IF (ABS(de).LT.3.) THEN
de = de*de
......@@ -133,12 +133,12 @@ CONTAINS
ENDDO
s1 = s1 + s
ENDDO
results%seigv = (2/input%jspins)*results%seigv
seigv = (2/input%jspins)*seigv
seigv1 = (1/input%jspins)*fact1*s2
chmom = s1 - input%jspins*s
IF ( mpi%irank == 0 ) THEN
WRITE (6,FMT=8050) results%seigv - seigv1,s1,chmom
WRITE (16,FMT=8050) results%seigv - seigv1,s1,chmom
WRITE (6,FMT=8050) seigv - seigv1,s1,chmom
WRITE (16,FMT=8050) seigv - seigv1,s1,chmom
END IF
8050 FORMAT (/,10x,'sum of eigenvalues-correction=',f12.5,/,10x,&
'sum of weight =',f12.5,/,10x,&
......
......@@ -6,8 +6,8 @@
MODULE m_ferhis
CONTAINS
SUBROUTINE ferhis(input,kpts,mpi,results, index,idxeig,idxkpt,idxjsp,n,&
nstef,ws,spindg,weight, e,ne,we, noco,cell)
SUBROUTINE ferhis(input,kpts,mpi, index,idxeig,idxkpt,idxjsp,n,&
nstef,ws,spindg,weight, e,ne,we, noco,cell,ef,seigv,w_iks,results)
!***********************************************************************
!
! This subroutine determines the fermi energy and the sum of the
......@@ -65,6 +65,8 @@ CONTAINS
! .. Scalar Arguments ..
INTEGER,INTENT(IN) :: n ,nstef
REAL,INTENT(IN) :: spindg,ws,weight
REAL,INTENT(INOUT) :: ef,seigv
REAL,INTENT(OUT) :: w_iks(:,:,:)
! ..
! .. Array Arguments ..
INTEGER, INTENT (IN) :: idxeig(:)!(dimension%neigd*kpts%nkpt*dimension%jspd)
......@@ -124,9 +126,9 @@ CONTAINS
WRITE (6,FMT='(''FERHIS: Fermi-Energy by histogram:'')')
END IF
efermi = results%ef
efermi = ef
IF (nstef.LT.n) THEN
gap = e(INDEX(nstef+1)) - results%ef
gap = e(INDEX(nstef+1)) - ef
results%bandgap = gap*hartree_to_ev_const
IF ( mpi%irank == 0 ) THEN
attributes = ''
......@@ -156,9 +158,9 @@ CONTAINS
!
!---> STATES ABOVE EF AVAILABLE
!
results%ef = 0.5* (e(INDEX(nstef+1))+results%ef)
emax = results%ef + 8.0*tkb
emin = results%ef - 8.0*tkb
ef = 0.5* (e(INDEX(nstef+1))+ef)
emax = ef + 8.0*tkb
emin = ef - 8.0*tkb
w_near_ef = 0.0
w_below_emin = 0.0
inkem = 0
......@@ -187,11 +189,11 @@ CONTAINS
!---> ADJUST FERMI-ENERGY BY NEWTON-METHOD
!
nocst = ink - 1
CALL ef_newton(n,mpi%irank, inkem,nocst,index,tkb,e, w_near_ef,results%ef,we)
CALL ef_newton(n,mpi%irank, inkem,nocst,index,tkb,e, w_near_ef,ef,we)
!
IF ( mpi%irank == 0 ) THEN
WRITE (16,FMT=8030) results%ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef)
WRITE (6,FMT=8030) results%ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef)
WRITE (16,FMT=8030) ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef)
WRITE (6,FMT=8030) ef,spindg*weight, spindg*w_below_emin,spindg* (w_below_emin+w_near_ef)
END IF
ELSE
......@@ -201,7 +203,7 @@ CONTAINS
IF ( mpi%irank == 0 ) WRITE (6,FMT=8020)
nocst = nstef
we(INDEX(nocst)) = we(INDEX(nocst)) - wfermi
results%ef = efermi
ef = efermi
tkb = 0.0
END IF
ELSE
......@@ -231,11 +233,11 @@ CONTAINS
!=======> DETERMINE OCCUPATION NUMBER AND WEIGHT OF EIGENVALUES
! FOR EACH K_POINT
!
results%w_iks(:,:,:) = 0.0
w_iks(:,:,:) = 0.0
IF ( mpi%irank == 0 ) WRITE (6,FMT=8080) nocst
DO i=1,nocst
results%w_iks(idxeig(INDEX(i)),idxkpt(INDEX(i)),idxjsp(INDEX(i))) = we(INDEX(i))
w_iks(idxeig(INDEX(i)),idxkpt(INDEX(i)),idxjsp(INDEX(i))) = we(INDEX(i))
ENDDO
!
!======> CHECK SUM OF VALENCE WEIGHTS
......@@ -244,7 +246,7 @@ CONTAINS
wvals = 0.0
DO js = 1,nspins
DO k = 1,kpts%nkpt
wvals = wvals + SUM(results%w_iks(:ne(k,js),k,js))
wvals = wvals + SUM(w_iks(:ne(k,js),k,js))
ENDDO
ENDDO
......@@ -268,7 +270,7 @@ CONTAINS
DO js = 1,nspins
DO kpt = 1 , kpts%nkpt
DO nocc=1,ne(kpt,js)
fermikn = results%w_iks(nocc,kpt,js)/kpts%wtkpt(kpt)
fermikn = w_iks(nocc,kpt,js)/kpts%wtkpt(kpt)
IF ( fermikn .GT. 0.0 .AND. fermikn .LT. 1.0 ) &
entropy = entropy + kpts%wtkpt(kpt) * ( fermikn * LOG( fermikn) + ( 1.0 - fermikn) * LOG( 1.0 - fermikn) )
END DO
......@@ -285,13 +287,13 @@ CONTAINS
!
!
results%seigv = spindg*DOT_PRODUCT(e(INDEX(:nocst)),we(INDEX(:nocst)))
seigv = seigv+spindg*DOT_PRODUCT(e(INDEX(:nocst)),we(INDEX(:nocst)))
IF (mpi%irank == 0) THEN
attributes = ''
WRITE(attributes(1),'(f20.10)') results%seigv
WRITE(attributes(1),'(f20.10)') seigv
WRITE(attributes(2),'(a)') 'Htr'
CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes)
WRITE (6,FMT=8040) results%seigv
WRITE (6,FMT=8040) seigv
END IF
......
......@@ -55,9 +55,9 @@ CONTAINS
!REAL, INTENT (OUT):: w(:,:,:) !(dimension%neigd,kpts%nkpt,dimension%jspd)
! ..
! .. Local Scalars ..
REAL del ,spindg,ssc ,ws,zc,tkb_1,weight
REAL del ,spindg,ssc ,ws,zc,weight,efermi
INTEGER i,idummy,j,jsp,k,l,n,nbands,nstef,nv,nmat,nspins
INTEGER n_help
INTEGER n_help,m_spins,mspin,sslice(2)
! ..
! .. Local Arrays ..
!
......@@ -138,7 +138,30 @@ CONTAINS
WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,eig(1:ne(k,jsp),k,jsp))
END IF
nv= -1
END DO
ENDDO
!finished reading of eigenvalues
IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues')
if (abs(input%fixed_moment)<1E-6) THEN
!this is a standard calculation
m_spins=1
else
!total moment is fixed
m_spins=2
end if
do mspin=1,m_spins
IF (m_spins == 1) THEN
sslice = (/1,nspins/)
ELSE
sslice = (/mspin,mspin/)
nspins = 1
ENDIF
n = 0
DO jsp = sslice(1),sslice(2)
!Generate a list of energies
DO k = 1,kpts%nkpt
!
!---> STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE
!---> CONECTION TO THE ORIGINAL ARRAYS
......@@ -149,12 +172,11 @@ CONTAINS
idxeig(n+j) = j+n_help
idxkpt(n+j) = k
idxjsp(n+j) = jsp
END DO
END DO
!---> COUNT THE NUMBER OF EIGENVALUES
n = n + ne(k,jsp)
END DO
END DO
IF (mpi%irank == 0) CALL closeXMLElement('eigenvalues')
CALL sort(n,e,index)
......@@ -173,18 +195,17 @@ CONTAINS
!
weight = input%zelec/spindg
results%seigv = 0.0e0
IF(m_spins /= 1) weight = weight/2.0 -(mspin-1.5)*input%fixed_moment
ws = 0.0e0
l = 0
DO WHILE ((ws+del).LT.weight)
l = l + 1
IF (l.GT.n) THEN
IF ( mpi%irank == 0 ) THEN
WRITE (16,FMT=8010) n,ws,weight
WRITE (6,FMT=8010) n,ws,weight
END IF
CALL juDFT_error("fermi",calledby="fermie")
8010 FORMAT (/,10x,'error: not enough wavefunctions.',i10,&
& 2d20.10)
CALL juDFT_error("Not enough eavefunctions",calledby="fermie")
8010 FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10)
END IF
ws = ws + we(INDEX(l))
results%seigv = results%seigv + e(INDEX(l))*we(INDEX(l))*spindg
......@@ -193,28 +214,44 @@ CONTAINS
results%ef = e(INDEX(l))
nstef = l
zc = input%zelec
IF(m_spins /= 1) THEN
zc = zc/2.0-(mspin-1.5)*input%fixed_moment
idxjsp = 1 !assume single spin in following calculations
IF (mspin == 1) THEN
WRITE(6,*) "Fixed total moment calculation"
WRITE(6,*) "Moment:",input%fixed_moment
write(6,*) "First Spin:"
ELSE
WRITE(6,*) "Second Spin:"
ENDIF
ENDIF
IF ( mpi%irank == 0 ) WRITE (6,FMT=8020) results%ef,nstef,results%seigv,ws,results%seigsc,ssc
!+po
results%ts = 0.0
!-po
results%w_iks = 0.0
results%w_iks(:,:,sslice(1):sslice(2)) = 0.0
results%bandgap = 0.0
IF (input%gauss) THEN
CALL fergwt(kpts,input,mpi,ne, eig,results)
CALL fergwt(kpts,input,mpi,ne(:,sslice(1):sslice(2)), eig(:,:,sslice(1):sslice(2)),results%ef,results%w_iks(:,:,sslice(1):sslice(2)),results%seigv)
ELSE IF (input%tria) THEN
CALL fertri(input,kpts,mpi%irank, ne,kpts%nkpt,nspins,zc,eig,kpts%bk,spindg,&
results%ef,results%seigv,results%w_iks)
CALL fertri(input,kpts,mpi%irank, ne(:,sslice(1):sslice(2)),kpts%nkpt,nspins,zc,eig(:,:,sslice(1):sslice(2)),kpts%bk,spindg,&
results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)))
ELSE
nspins = input%jspins
IF (noco%l_noco) nspins = 1
tkb_1 = input%tkb
CALL ferhis(input,kpts,mpi,results,index,idxeig,idxkpt,idxjsp, n,&
nstef,ws,spindg,weight,e,ne,we, noco,cell)
CALL ferhis(input,kpts,mpi,index,idxeig,idxkpt,idxjsp, n,&
nstef,ws,spindg,weight,e,ne(:,sslice(1):sslice(2)),we, noco,cell,results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)),results)
END IF
! 7.12.95 r.pentcheva seigscv must be calculated outside if (gauss)
results%seigscv = results%seigsc + results%seigv
!
IF (mspin == 2) THEN
WRITE(6,*) "Different Fermi-energies for both spins:"
WRITE(6,"(a,f0.3,a,f0.4,a,f0.4,a,f0.4)") "Fixed Moment:" &
,input%fixed_moment," Difference(EF):",efermi," - ",results%ef,"="&
,efermi-results%ef
ENDIF
efermi = results%ef
enddo
DEALLOCATE ( idxeig,idxjsp,idxkpt,index,e,eig,we )
attributes = ''
......
......@@ -350,7 +350,9 @@ SUBROUTINE r_inpXML(&
noco%l_noco = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@l_noco'))
input%swsp = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@swsp'))
input%lflip = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@lflip'))
input%fixed_moment=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@fixed_moment'))
IF (ABS(input%fixed_moment)>1E-8.AND.(input%jspins==1.OR.noco%l_noco)) CALL judft_error("Fixed moment only in collinear calculations with two spins")
dimension%jspd = input%jspins
! Read in Brillouin zone integration parameters
......
......@@ -142,6 +142,7 @@
<xsd:attribute default="F" name="l_J" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="swsp" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="lflip" type="FleurBool" use="optional"/>
<xsd:attribute default="0.0" name="fixed_moment" type="xsd:string" use="optional"/>
</xsd:complexType>
<xsd:simpleType name="SpinNumberType">
......
......@@ -410,6 +410,7 @@ MODULE m_types_setup
REAL :: elup
REAL :: rkmax
REAL :: zelec
REAL :: fixed_moment=0.0
CHARACTER(LEN=8) :: comment(10)
REAL,POINTER :: sigma !this is the difference in charge due to the electric field it points to the value stored in t_efield
LOGICAL :: l_core_confpot
......
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