Commit 290bd4d8 authored by Gregor Michalicek's avatar Gregor Michalicek

Implemented autocorrection of Fermi energy in bandstructures

Also in this commit: improved Htr to eV factor in evaldos and
                     adapted the related tests.
parent fbcfc476
......@@ -880,7 +880,7 @@ CONTAINS
sliceplot,noco,sym,&
cell,&
l_mcd,ncored,ncore,e_mcd,&
results%ef,nsld,oneD)
results%ef,results%bandgap,nsld,oneD)
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
CALL Ek_write_sl(&
eig_id,dimension,kpts,atoms,vacuum,&
......
......@@ -17,7 +17,7 @@ CONTAINS
& sliceplot,noco,sym,&
& cell,&
& l_mcd,ncored,ncore,e_mcd,&
& efermi,nsld,oneD)
& efermi,bandgap,nsld,oneD)
USE m_eig66_io,ONLY:read_dos,read_eig
USE m_evaldos
USE m_cdninf
......@@ -40,7 +40,7 @@ CONTAINS
INTEGER,PARAMETER :: n2max=13
INTEGER, INTENT (IN) :: nsld,eig_id
INTEGER, INTENT (IN) :: ncored
REAL, INTENT (IN) :: efermi
REAL, INTENT (IN) :: efermi, bandgap
LOGICAL, INTENT (IN) :: l_mcd
! ..
! .. Array Arguments ..
......@@ -137,7 +137,7 @@ CONTAINS
IF (banddos%dos.AND.(banddos%ndir.LT.0)) THEN
CALL evaldos(&
& eig_id,input,banddos,vacuum,kpts,atoms,sym,noco,oneD,cell,&
& DIMENSION,efermi,&
& DIMENSION,efermi,bandgap,&
& l_mcd,ncored,ncore,e_mcd,nsld)
ENDIF
!
......
MODULE m_evaldos
CONTAINS
SUBROUTINE evaldos(eig_id,input,banddos,vacuum,kpts,atoms,sym,noco,oneD,cell,&
dimension,efermiarg, l_mcd,ncored,ncore,e_mcd,nsld)
dimension,efermiarg,bandgap,l_mcd,ncored,ncore,e_mcd,nsld)
!----------------------------------------------------------------------
!
! vk: k-vectors
......@@ -27,6 +27,8 @@
USE m_ptdos
USE m_smooth
USE m_types
USE m_constants
USE m_cdn_io
IMPLICIT NONE
INTEGER,INTENT(IN) :: eig_id
TYPE(t_dimension),INTENT(IN) :: dimension
......@@ -42,7 +44,7 @@
INTEGER, INTENT(IN) :: ncored
INTEGER, INTENT(IN) :: nsld
REAL, INTENT(IN) :: efermiarg
REAL, INTENT(IN) :: efermiarg, bandgap
LOGICAL, INTENT(IN) :: l_mcd
INTEGER, INTENT(IN) :: ncore(:)!(ntype)
......@@ -51,12 +53,11 @@
!+odim
! locals
INTEGER, PARAMETER :: lmax= 4, ned = 1301
REAL, PARAMETER :: factor = 27.2
INTEGER i,s,v,index,jspin,k,l,l1,l2,ln,n,nl,ntb,ntria,ntetra
INTEGER icore,qdim,n_orb
REAL as,de,efermi,emax,emin,qmt,sigma,totdos
REAL e_up,e_lo,e_test1,e_test2,fac,sumwei,dk
LOGICAL l_tria,l_orbcomp
REAL as,de,efermi,emax,emin,qmt,sigma,totdos,efermiPrev
REAL e_up,e_lo,e_test1,e_test2,fac,sumwei,dk,eFermiCorrection
LOGICAL l_tria,l_orbcomp,l_error
INTEGER itria(3,2*kpts%nkpt),nevk(kpts%nkpt),itetra(4,6*kpts%nkpt)
INTEGER, ALLOCATABLE :: ksym(:),jsym(:)
......@@ -97,10 +98,10 @@
ENDIF
!
! scale energies
sigma = banddos%sig_dos*factor
emin =min(banddos%e1_dos*factor,banddos%e2_dos*factor)
emax =max(banddos%e1_dos*factor,banddos%e2_dos*factor)
efermi = efermiarg*factor
sigma = banddos%sig_dos*hartree_to_ev_const
emin =min(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
emax =max(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
efermi = efermiarg*hartree_to_ev_const
WRITE (6,'(a)') 'DOS-Output is generated!'
......@@ -135,8 +136,8 @@
ENDDO
ENDDO
ENDDO
e_lo = e_lo*factor - efermi - emax
e_up = e_up*factor - efermi
e_lo = e_lo*hartree_to_ev_const - efermi - emax
e_up = e_up*hartree_to_ev_const - efermi
de = (e_up-e_lo)/(ned-1)
DO i=1,ned
e_grid(i) = e_lo + (i-1)*de
......@@ -237,7 +238,7 @@
!---- > convert eigenvalues to ev and shift them by efermi
!
DO i = 1 , nevk(k)
ev(i,k) = ev(i,k)*factor - efermi
ev(i,k) = ev(i,k)*hartree_to_ev_const - efermi
ENDDO
DO i = nevk(k) + 1, dimension%neigd
ev(i,k) = 9.9e+99
......@@ -418,8 +419,8 @@
DO icore = 1 , ncore(n)
DO i = 1 , ned-1
IF (e(i).GT.0) THEN ! take unoccupied part only
e_test1 = -e(i) - efermi +e_mcd(n,jspin,icore)*factor
e_test2 = -e(i+1)-efermi +e_mcd(n,jspin,icore)*factor
e_test1 = -e(i) - efermi +e_mcd(n,jspin,icore)*hartree_to_ev_const
e_test2 = -e(i+1)-efermi +e_mcd(n,jspin,icore)*hartree_to_ev_const
IF ((e_test2.LE.e_grid(l)).AND. (e_test1.GT.e_grid(l))) THEN
fac = (e_grid(l)-e_test1)/(e_test2-e_test1)
DO k = 3*(n-1)+1,3*(n-1)+3
......@@ -472,6 +473,20 @@
!------------------------------------------------------------------------------
IF (banddos%ndir == -4) THEN
IF(bandgap.LT.(8.0*input%tkb*hartree_to_ev_const)) THEN
CALL readPrevEFermi(eFermiPrev,l_error)
eFermiCorrection = 0.0
IF(.NOT.l_error) THEN
WRITE(*,*) 'Fermi energy is automatically corrected in bands.* files.'
WRITE(*,*) 'It is consistent with last calculated density!'
WRITE(*,*) 'No manual correction (e.g. in band.gnu file) required.'
eFermiCorrection = (eFermiPrev-efermiarg)*hartree_to_ev_const
ELSE
WRITE(*,*) 'Fermi energy in bands.* files may not be consistent with last density.'
WRITE(*,*) 'Please correct it manually (e.g. in band.gnu file).'
END IF
END IF
OPEN (18,FILE='bands'//spin12(jspin))
ntb = minval(nevk(:))
kx(1) = 0.0
......@@ -485,7 +500,7 @@
ENDDO
DO i = 1, ntb
DO k = 1, kpts%nkpt
write(18,'(2f15.9)') kx(k),ev(i,k)
write(18,'(2f15.9)') kx(k),ev(i,k)-eFermiCorrection
ENDDO
ENDDO
CLOSE (18)
......
......@@ -128,6 +128,7 @@ CONTAINS
efermi = results%ef
IF (nstef.LT.n) THEN
gap = e(INDEX(nstef+1)) - results%ef
results%bandgap = gap*hartree_to_ev_const
IF ( mpi%irank == 0 ) THEN
attributes = ''
WRITE(attributes(1),'(f20.10)') gap*hartree_to_ev_const
......
......@@ -212,6 +212,7 @@ CONTAINS
results%ts = 0.0
!-po
results%w_iks = 0.0
results%bandgap = 0.0
IF (input%gauss) THEN
CALL fergwt(kpts,input,mpi,ne, eig,results)
ELSE IF (input%tria) THEN
......
......@@ -706,6 +706,7 @@
REAL :: e_ldau !<total energy contribution of LDA+U
REAL :: tote
REAL :: last_distance
REAL :: bandgap
TYPE(t_energy_hf) :: te_hfex
REAL :: te_hfex_loc(2)
REAL, ALLOCATABLE :: w_iks(:,:,:)
......
......@@ -28,7 +28,7 @@ MODULE m_cdn_io
PUBLIC readDensity, writeDensity
PUBLIC isDensityFilePresent, isCoreDensityPresent
PUBLIC readCoreDensity, writeCoreDensity
PUBLIC setStartingDensity
PUBLIC setStartingDensity, readPrevEFermi
PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
PUBLIC CDN_ARCHIVE_TYPE_CDN_const
......@@ -87,7 +87,6 @@ MODULE m_cdn_io
fermiEnergy = 0.0
l_qfix = .FALSE.
WRITE(*,*) 'fermiEnergy and l_qfix set to default values in readDensity!'
CALL getMode(mode)
......@@ -495,6 +494,58 @@ MODULE m_cdn_io
END SUBROUTINE writeDensity
SUBROUTINE readPrevEFermi(eFermiPrev,l_error)
REAL, INTENT(OUT) :: eFermiPrev
LOGICAL, INTENT(OUT) :: l_error
INTEGER :: mode
#ifdef CPP_HDF
INTEGER(HID_T) :: fileID
#endif
INTEGER :: currentStarsIndex,currentLatharmsIndex
INTEGER :: currentStructureIndex
INTEGER :: readDensityIndex, lastDensityIndex
INTEGER :: starsIndex, latharmsIndex, structureIndex
INTEGER :: iter, jspins, previousDensityIndex
REAL :: fermiEnergy
LOGICAL :: l_qfix, l_exist
CHARACTER(LEN=30) :: archiveName
CALL getMode(mode)
eFermiPrev = 0.0
l_error = .FALSE.
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
CALL openCDN_HDF(fileID,currentStarsIndex,currentLatharmsIndex,currentStructureIndex,&
readDensityIndex,lastDensityIndex)
WRITE(archiveName,'(a,i0)') '/cdn-', readDensityIndex
CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_UNDEFINED_const,&
iter, starsIndex, latharmsIndex, structureIndex,&
previousDensityIndex, jspins, fermiEnergy, l_qfix)
archiveName = ''
WRITE(archiveName,'(a,i0)') '/cdn-', previousDensityIndex
l_exist = isDensityEntryPresentHDF(fileID,archiveName,DENSITY_TYPE_NOCO_OUT_const)
IF(l_exist) THEN
CALL peekDensityEntryHDF(fileID, archiveName, DENSITY_TYPE_NOCO_OUT_const,&
iter, starsIndex, latharmsIndex, structureIndex,&
previousDensityIndex, jspins, fermiEnergy, l_qfix)
eFermiPrev = fermiEnergy
ELSE
l_error = .TRUE.
END IF
CALL closeCDN_HDF(fileID)
#endif
ELSE IF(mode.EQ.CDN_STREAM_MODE) THEN
STOP 'cdn.str not yet implemented!'
ELSE
l_error = .TRUE.
END IF
END SUBROUTINE
SUBROUTINE readCoreDensity(input,atoms,dimension,rhcs,tecs,qints)
TYPE(t_atoms),INTENT(IN) :: atoms
......
......@@ -10,6 +10,6 @@ jt::testrun($executable,$workdir);
#now test output
$result=jt::test_fileexists("$workdir/band.1");
$result=jt::test_fileexists("$workdir/DOS.1");
$result+=jt::test_grepnumber("$workdir/DOS.1","9.40276","9.40276 (.....)",0.290,0.0001);
$result+=jt::test_grepnumber("$workdir/DOS.1","9.40670","9.40670 (.....)",0.290,0.0001);
jt::stageresult($workdir,$result,"2");
......@@ -8,6 +8,6 @@ jt::testrun("$executable -xmlInput",$workdir);
#now test output
$result=jt::test_fileexists("$workdir/band.1");
$result=jt::test_fileexists("$workdir/DOS.1");
$result+=jt::test_grepnumber("$workdir/DOS.1","9.40276","9.40276 (.....)",0.290,0.0001);
$result+=jt::test_grepnumber("$workdir/DOS.1","9.40670","9.40670 (.....)",0.290,0.0001);
jt::stageresult($workdir,$result,"1");
......@@ -9,6 +9,6 @@ jt::testrun($executable,$workdir);
#now test output
$result=jt::test_fileexists("$workdir/DOS.1");
$result+=jt::test_grepnumber("$workdir/DOS.1","10.88235","10.88235 (.....)",0.105,0.0001);
$result+=jt::test_grepnumber("$workdir/DOS.1","10.88690","10.88690 (.....)",0.105,0.0001);
jt::stageresult($workdir,$result,"2");
......@@ -7,6 +7,6 @@ jt::testrun("$executable -xmlInput",$workdir);
#now test output
$result=jt::test_fileexists("$workdir/DOS.1");
$result+=jt::test_grepnumber("$workdir/DOS.1","10.88486","10.88486 (.....)",0.105,0.0001);
$result+=jt::test_grepnumber("$workdir/DOS.1","10.88842","10.88842 (.....)",0.107,0.0001);
jt::stageresult($workdir,$result,"1");
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