Commit b9c27918 authored by Daniel Wortmann's avatar Daniel Wortmann

Bugfixes in enpara handling...

parent e223c574
...@@ -915,9 +915,9 @@ CONTAINS ...@@ -915,9 +915,9 @@ CONTAINS
IF (.NOT.sliceplot%slice) THEN IF (.NOT.sliceplot%slice) THEN
DO n=1,atoms%ntype DO n=1,atoms%ntype
enpara%el1(0:3,n,ispin)=ener(0:3,n,ispin)/sqal(0:3,n,ispin) enpara%el1(0:3,n,ispin)=ener(0:3,n,ispin)/sqal(0:3,n,ispin)
enpara%ello1(:,n,ispin)=enerlo(:,n,ispin)/sqlo(:,n,ispin) IF (atoms%nlo(n)>0) enpara%ello1(:atoms%nlo(n),n,ispin)=enerlo(:atoms%nlo(n),n,ispin)/sqlo(:atoms%nlo(n),n,ispin)
ENDDO ENDDO
enpara%evac1(:,ispin)=pvac(:,ispin)/svac(:,ispin) IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=pvac(:vacuum%nvac,ispin)/svac(:vacuum%nvac,ispin)
END IF END IF
!---> check continuity of charge density !---> check continuity of charge density
......
...@@ -52,7 +52,7 @@ init/tetcon.f init/kvecon.f init/boxdim.f math/ylm4.f global/radsra.f math/intgr ...@@ -52,7 +52,7 @@ init/tetcon.f init/kvecon.f init/boxdim.f math/ylm4.f global/radsra.f math/intgr
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90 set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90 eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90
global/enpara.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/enpara.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 init/julia.f90 global/utility.F90 juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 init/julia.f90 global/utility.F90
init/compile_descr.F90 init/kpoints.f90 io/xmlOutput.F90 init/brzone2.f90) init/compile_descr.F90 init/kpoints.f90 io/xmlOutput.F90 init/brzone2.f90)
......
...@@ -27,4 +27,5 @@ global/phasy1.f90 ...@@ -27,4 +27,5 @@ global/phasy1.f90
global/qfix.f90 global/qfix.f90
global/radflo.F90 global/radflo.F90
global/utility.F90 global/utility.F90
global/find_enpara.f90
) )
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_find_enpara
USE m_judft
IMPLICIT NONE
PRIVATE
CHARACTER(len=1),PARAMETER,DIMENSION(0:9):: ch=(/'s','p','d','f','g','h','i','j','k','l'/)
PUBLIC:: find_enpara
CONTAINS
!> Function to determine the energy parameter given the quantum number and the potential
!! Different schemes are implemented. Nqn (main quantum number) is used as a switch.
!! This code was previously in lodpot.f
REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
USE m_types_setup
USE m_types_mpi
USE m_radsra
USE m_differ
USE m_xmlOutput
USE m_constants
IMPLICIT NONE
LOGICAL,INTENT(IN):: lo
INTEGER,INTENT(IN):: l,n,nqn,jsp
TYPE(t_atoms),INTENT(IN)::atoms
TYPE(t_mpi),INTENT(IN) ::mpi
REAL,INTENT(IN):: vr(:)
IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr)
IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr)
END FUNCTION find_enpara
REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
USE m_types_setup
USE m_types_mpi
USE m_radsra
USE m_differ
USE m_xmlOutput
USE m_constants
IMPLICIT NONE
LOGICAL,INTENT(IN):: lo
INTEGER,INTENT(IN):: l,n,nqn,jsp
TYPE(t_atoms),INTENT(IN)::atoms
TYPE(t_mpi),INTENT(IN) ::mpi
REAL,INTENT(IN):: vr(:)
INTEGER j,ilo,i
INTEGER nodeu,node,ierr,msh
REAL e_up,e_lo,lnd,e1
REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
LOGICAL start
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: f(:,:),vrd(:)
CHARACTER(LEN=20) :: attributes(6)
c=c_light(1.0)
!Core potential setup done for each n,l now
d = EXP(atoms%dx(n))
! set up core-mesh
rn = atoms%rmt(n)
msh = atoms%jri(n)
DO WHILE (rn < atoms%rmt(n) + 20.0)
msh = msh + 1
rn = rn*d
ENDDO
rn = atoms%rmsh(1,n)*( d**(msh-1) )
ALLOCATE ( f(msh,2),vrd(msh) )
! extend core potential (linear with slope t1 / a.u.)
vrd(:atoms%jri(n))=vr(:atoms%jri(n))
t1=0.125
t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
rr = atoms%rmt(n)
DO j = atoms%jri(n) + 1, msh
rr = d*rr
vrd(j) = rr*( t2 + rr*t1 )
ENDDO
node = nqn - (l+1)
e = 0.0
! determine upper edge
nodeu = -1 ; start = .TRUE.
DO WHILE ( nodeu <= node )
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
IF ( ( nodeu > node ) .AND. start ) THEN
e = e - 1.0
nodeu = -1
ELSE
e = e + 0.01
start = .FALSE.
ENDIF
ENDDO
e_up = e
IF (node /= 0) THEN
! determine lower edge
nodeu = node + 1
DO WHILE ( nodeu >= node )
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
e = e - 0.01
ENDDO
e_lo = e
ELSE
e_lo = -9.99
ENDIF
! calculate core
e = (e_up+e_lo)/2
fn = REAL(nqn) ; fl = REAL(l) ; fj = fl + 0.5
CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
rn,d,msh,vrd, e, f(:,1),f(:,2),ierr)
IF (lo.AND.l>0) THEN
e1 = (e_up+e_lo)/2
fn = REAL(nqn) ; fl = REAL(l) ; fj = fl-0.5
CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
rn,d,msh,vrd, e1, f(:,1),f(:,2),ierr)
e = (2.0*e + e1 ) / 3.0
ENDIF
IF (mpi%irank == 0) THEN
attributes = ''
WRITE(attributes(1),'(i0)') n
WRITE(attributes(2),'(i0)') jsp
WRITE(attributes(3),'(i0,a1)') nqn, ch(l)
WRITE(attributes(4),'(f8.2)') e_lo
WRITE(attributes(5),'(f8.2)') e_up
WRITE(attributes(6),'(f16.10)') e
IF (lo) THEN
CALL writeXMLElementForm('loAtomicEP',(/'atomType ','spin ','branch ',&
'branchLowest ','branchHighest','value '/),&
attributes,RESHAPE((/10,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
ELSE
CALL writeXMLElementForm('atomicEP',(/'atomType ','spin ','branch ',&
'branchLowest ','branchHighest','value '/),&
attributes,RESHAPE((/12,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
ENDIF
WRITE(6,'(a6,i3,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') ' Atom',n,nqn,ch(l),' branch from',&
e_lo, ' to',e_up,' htr. ; e_l =',e
ENDIF
END FUNCTION priv_method1
REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
USE m_types_setup
USE m_types_mpi
USE m_radsra
USE m_differ
USE m_xmlOutput
USE m_constants
IMPLICIT NONE
LOGICAL,INTENT(IN):: lo
INTEGER,INTENT(IN):: l,n,nqn,jsp
TYPE(t_atoms),INTENT(IN)::atoms
TYPE(t_mpi),INTENT(IN) ::mpi
REAL,INTENT(IN):: vr(:)
INTEGER j,ilo,i
INTEGER nodeu,node,ierr,msh
REAL e_up,e_lo,lnd,e_up_temp,e_lo_temp,large_e_step
REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
LOGICAL start
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: f(:,:),vrd(:)
CHARACTER(LEN=20) :: attributes(6)
c=c_light(1.0)
!Core potential setup done for each n,l now
d = EXP(atoms%dx(n))
! set up core-mesh
rn = atoms%rmt(n)
msh = atoms%jri(n)
DO WHILE (rn < atoms%rmt(n) + 20.0)
msh = msh + 1
rn = rn*d
ENDDO
rn = atoms%rmsh(1,n)*( d**(msh-1) )
ALLOCATE ( f(msh,2),vrd(msh) )
! extend core potential (linear with slope t1 / a.u.)
vrd(:atoms%jri(n))=vr(:atoms%jri(n))
t1=0.125
t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
rr = atoms%rmt(n)
DO j = atoms%jri(n) + 1, msh
rr = d*rr
vrd(j) = rr*( t2 + rr*t1 )
ENDDO
! search for branches
node = ABS(nqn) - (l+1)
e = 0.0 ! The initial value of e is arbitrary.
large_e_step = 5.0 ! 5.0 Htr steps for coarse energy searches
! determine upper band edge
! Step 1: Coarse search for the band edge
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
DO WHILE ( nodeu > node )
e = e - large_e_step
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
END DO
DO WHILE ( nodeu <= node )
e = e + large_e_step
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
END DO
e_up_temp = e
e_lo_temp = e - large_e_step
! Step 2: Fine band edge determination by bisection search
DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
e = (e_up_temp + e_lo_temp) / 2.0
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
IF (nodeu > node) THEN
e_up_temp = e
ELSE
e_lo_temp = e
END IF
END DO
e_up = (e_up_temp + e_lo_temp) / 2.0
e = e_up
! determine lower band edge
IF (node == 0) THEN
e_lo = -49.99
ELSE
! Step 1: Coarse search for the band edge
nodeu = node
DO WHILE ( nodeu >= node )
e = e - large_e_step
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
ENDDO
e_up_temp = e + large_e_step
e_lo_temp = e
! Step 2: Fine band edge determination by bisection search
DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
e = (e_up_temp + e_lo_temp) / 2.0
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
IF (nodeu < node) THEN
e_lo_temp = e
ELSE
e_up_temp = e
END IF
END DO
e_lo = (e_up_temp + e_lo_temp) / 2.0
END IF
! determince notches by intersection
ldmt= -99.0 !ldmt = logarithmic derivative @ MT boundary
lnd = -l-1
DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07)
e = (e_up+e_lo)/2
CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
ldmt = dus/us
IF( ldmt .GT. lnd) THEN
e_lo = e
ELSE IF( ldmt .LT. lnd ) THEN
e_up = e
e_lo = e_lo
END IF
END DO
IF (mpi%irank == 0) THEN
attributes = ''
WRITE(attributes(1),'(i0)') n
WRITE(attributes(2),'(i0)') jsp
WRITE(attributes(3),'(i0,a1)') ABS(nqn), ch(l)
WRITE(attributes(4),'(f16.10)') ldmt
WRITE(attributes(5),'(f16.10)') e
IF (lo) THEN
CALL writeXMLElementForm('heloAtomicEP',(/'atomType ','spin ','branch ',&
'logDerivMT ','value '/),&
attributes(1:5),reshape((/8,4,6,12,5+17,6,1,3,16,16/),(/5,2/)))
ELSE
CALL writeXMLElementForm('heAtomicEP',(/'atomType ','spin ','branch ',&
'logDerivMT ','value '/),&
attributes(1:5),reshape((/10,4,6,12,5+17,6,1,3,16,16/),(/5,2/)))
ENDIF
WRITE (6,'(a7,i3,i2,a1,a12,f7.2,a4,f7.2,a5)') " Atom ",n,nqn,ch(l)," branch, D = ",&
ldmt, " at ",e," htr."
ENDIF
END FUNCTION priv_method2
END MODULE m_find_enpara
...@@ -1263,7 +1263,7 @@ SUBROUTINE r_inpXML(& ...@@ -1263,7 +1263,7 @@ SUBROUTINE r_inpXML(&
atoms%nlod = max(atoms%nlod,atoms%nlo(iType)) atoms%nlod = max(atoms%nlod,atoms%nlo(iType))
END DO END DO
atoms%nlod = max(atoms%nlod,2) ! for chkmt atoms%nlod = max(atoms%nlod,2) ! for chkmt
ALLOCATE(atoms%llo(atoms%nlod,atoms%ntype)) ALLOCATE(atoms%llo(atoms%nlod,atoms%ntype));atoms%llo=-1
ALLOCATE(atoms%ulo_der(atoms%nlod,atoms%ntype)) ALLOCATE(atoms%ulo_der(atoms%nlod,atoms%ntype))
ALLOCATE(atoms%l_dulo(atoms%nlod,atoms%ntype)) ! For what is this? ALLOCATE(atoms%l_dulo(atoms%nlod,atoms%ntype)) ! For what is this?
...@@ -1483,12 +1483,25 @@ SUBROUTINE r_inpXML(& ...@@ -1483,12 +1483,25 @@ SUBROUTINE r_inpXML(&
atoms%llod = max(abs(atoms%llo(iLLO,iType)),atoms%llod) atoms%llod = max(abs(atoms%llo(iLLO,iType)),atoms%llod)
DO jsp = 1, input%jspins DO jsp = 1, input%jspins
enpara%ello0(iLLO,iType,jsp) = speciesLOeParams(loOrderList(iLLO)) enpara%ello0(iLLO,iType,jsp) = speciesLOeParams(loOrderList(iLLO))
IF (enpara%ello0(iLLO,iType,jsp)==NINT(enpara%ello0(iLLO,iType,jsp))) THEN
enpara%qn_ello(iLLO,iType,jsp)=NINT(enpara%ello0(iLLO,iType,jsp))
enpara%ello0(iLLO,iType,jsp)=0
ELSE
enpara%qn_ello(iLLO,iType,jsp)=0
ENDIF
enpara%skiplo(iType,jsp)=enpara%skiplo(iType,jsp)+(2*atoms%llo(iLLO,itype)+1)
END DO END DO
END DO END DO
! Energy parameters ! Energy parameters
DO jsp = 1, input%jspins DO jsp = 1, input%jspins
DO l = 0, 3 DO l = 0, 3
enpara%el0(l,iType,jsp) = speciesEParams(l) enpara%el0(l,iType,jsp) = speciesEParams(l)
IF (enpara%el0(l,iType,jsp)==NINT(enpara%el0(l,iType,jsp))) THEN
enpara%qn_el(l,iType,jsp)=nint(enpara%el0(l,iType,jsp))
enpara%el0(l,iType,jsp)=0
ELSE
enpara%qn_el(l,iType,jsp)=0
ENDIF
END DO END DO
DO l = 4,atoms%lmax(iType) DO l = 4,atoms%lmax(iType)
enpara%el0(l,iType,jsp) = enpara%el0(3,iType,jsp) enpara%el0(l,iType,jsp) = enpara%el0(3,iType,jsp)
......
...@@ -445,11 +445,10 @@ SUBROUTINE w_inpXML(& ...@@ -445,11 +445,10 @@ SUBROUTINE w_inpXML(&
320 FORMAT(' <atomicCutoffs lmax="',i0,'" lnonsphr="',i0,'"/>') 320 FORMAT(' <atomicCutoffs lmax="',i0,'" lnonsphr="',i0,'"/>')
WRITE (fileNum,320) atoms%lmax(iAtomType),atoms%lnonsph(iAtomType) WRITE (fileNum,320) atoms%lmax(iAtomType),atoms%lnonsph(iAtomType)
IF (ALL((enpara%el0(0:3,iAtomType,1)-INT(enpara%el0(0:3,iAtomType,1))).LE.0.00000001)) THEN IF (ALL(enpara%qn_el(0:3,iAtomType,1).ne.0)) THEN
! <energyParameters s="3" p="3" d="3" f="4"/> ! <energyParameters s="3" p="3" d="3" f="4"/>
321 FORMAT(' <energyParameters s="',i0,'" p="',i0,'" d="',i0,'" f="',i0,'"/>') 321 FORMAT(' <energyParameters s="',i0,'" p="',i0,'" d="',i0,'" f="',i0,'"/>')
WRITE (fileNum,321) INT(enpara%el0(0,iAtomType,1)),INT(enpara%el0(1,iAtomType,1)),& WRITE (fileNum,321) enpara%qn_el(0:3,iAtomType,1)
INT(enpara%el0(2,iAtomType,1)),INT(enpara%el0(3,iAtomType,1))
END IF END IF
IF(ANY(xmlElectronStates(:,iAtomType).NE.noState_const)) THEN IF(ANY(xmlElectronStates(:,iAtomType).NE.noState_const)) THEN
......
This diff is collapsed.
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