find_enpara.f90 10.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
!--------------------------------------------------------------------------------
! 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) )
75 76
    f = 0.0
    vrd = 0.0
77 78 79 80 81 82 83 84 85 86 87
    ! 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)
88
    IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states")
89 90 91 92 93 94 95 96
    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
97
          IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
98 99 100
          nodeu = -1
       ELSE
          e = e + 0.01
101
          IF (e>1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
102 103 104 105 106 107 108 109 110 111 112 113
          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
114
          IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
       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
151
       WRITE(6,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') '  Atom',n,nqn,ch(l),' branch from',&
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
            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) )
194 195
    f = 0.0
    vrd = 0.0
196 197 198 199 200 201 202 203 204 205 206
    ! 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)
207
    IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states")
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
    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