LCOV - code coverage report
Current view: top level - global - find_enpara.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 54 119 45.4 %
Date: 2019-09-08 04:53:50 Functions: 2 3 66.7 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_find_enpara
       8             :   USE m_judft
       9             :   IMPLICIT NONE
      10             :   PRIVATE
      11             :   PUBLIC:: find_enpara
      12             : 
      13             : CONTAINS
      14             :  
      15             :   !> Function to determine the energy parameter given the quantum number and the potential
      16             :   !! Different schemes are implemented. Nqn (main quantum number) is used as a switch.
      17             :   !! This code was previously in lodpot.f
      18        3996 :   REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,mpi,vr,e_lo,e_up)RESULT(e)
      19             :     USE m_types_setup
      20             :     USE m_types_mpi
      21             :     USE m_radsra
      22             :     USE m_differ
      23             :     USE m_constants
      24             :     IMPLICIT NONE
      25             :     LOGICAL,INTENT(IN):: lo
      26             :     INTEGER,INTENT(IN):: l,n,nqn,jsp
      27             :     REAL,INTENT(OUT) :: e_lo,e_up
      28             :     TYPE(t_atoms),INTENT(IN)::atoms
      29             :     TYPE(t_mpi),INTENT(IN)  ::mpi
      30             :     REAL,INTENT(IN):: vr(:)
      31             : 
      32        3996 :     IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr,e_lo,e_up)
      33        3996 :     IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr,e_lo,e_up)
      34        3996 :   END FUNCTION find_enpara
      35             : 
      36             : 
      37        3996 :   REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr,e_lo,e_up)RESULT(e)
      38             :     USE m_types_setup
      39             :     USE m_types_mpi
      40             :     USE m_radsra
      41             :     USE m_differ
      42             :     USE m_constants
      43             :     IMPLICIT NONE
      44             :     LOGICAL,INTENT(IN):: lo
      45             :     INTEGER,INTENT(IN):: l,n,nqn,jsp
      46             :     REAL,INTENT(OUT)  :: e_lo,e_up
      47             :     TYPE(t_atoms),INTENT(IN)::atoms
      48             :     TYPE(t_mpi),INTENT(IN)  ::mpi
      49             :     REAL,INTENT(IN):: vr(:)
      50             : 
      51             : 
      52             :     INTEGER j,ilo,i
      53             :     INTEGER nodeu,node,ierr,msh
      54             :     REAL   lnd,e1
      55             :     REAL   d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
      56             :     LOGICAL start
      57             :     !     ..
      58             :     !     .. Local Arrays .. 
      59        3996 :     REAL, ALLOCATABLE :: f(:,:),vrd(:)
      60             :     CHARACTER(LEN=20)    :: attributes(6)
      61        3996 :     c=c_light(1.0)
      62             : 
      63             :     !Core potential setup done for each n,l now 
      64        3996 :     d = EXP(atoms%dx(n))
      65             :     ! set up core-mesh
      66        3996 :     rn = atoms%rmt(n)
      67        3996 :     msh = atoms%jri(n)
      68      922108 :     DO WHILE (rn < atoms%rmt(n) + 20.0)
      69      459056 :        msh = msh + 1
      70      463052 :        rn = rn*d
      71             :     ENDDO
      72        3996 :     rn = atoms%rmsh(1,n)*( d**(msh-1) )
      73        3996 :     ALLOCATE ( f(msh,2),vrd(msh) )
      74       11988 :     f = 0.0
      75     2865688 :     vrd = 0.0
      76             :     ! extend core potential (linear with slope t1 / a.u.)
      77        3996 :     vrd(:atoms%jri(n))=vr(:atoms%jri(n))
      78        3996 :     t1=0.125
      79        3996 :     t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
      80        3996 :     rr = atoms%rmt(n)
      81      463052 :     DO j = atoms%jri(n) + 1, msh
      82      459056 :        rr = d*rr
      83      463052 :        vrd(j) = rr*( t2 + rr*t1 )
      84             :     ENDDO
      85             : 
      86        3996 :     node = nqn - (l+1)
      87        3996 :     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")
      88        3996 :     e = 0.0 
      89             :     ! determine upper edge
      90        3996 :     nodeu = -1 ; start = .TRUE.
      91      775202 :     DO WHILE ( nodeu <= node ) 
      92             :        CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
      93      771206 :             atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
      94      775202 :        IF  ( ( nodeu > node ) .AND. start ) THEN
      95         586 :           e = e - 1.0
      96         586 :           IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
      97         586 :           nodeu = -1
      98             :        ELSE
      99      770620 :           e = e + 0.01
     100      770620 :           IF (e>1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
     101             :           start = .FALSE.
     102             :        ENDIF
     103             :     ENDDO
     104             : 
     105        3996 :     e_up = e
     106        3996 :     IF (node /= 0) THEN
     107             :        ! determine lower edge
     108        2080 :        nodeu = node + 1
     109     1335920 :        DO WHILE ( nodeu >= node ) 
     110             :           CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     111     1333840 :                atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     112     1333840 :           e = e - 0.01
     113     1333840 :           IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
     114             :        ENDDO
     115        2080 :        e_lo = e
     116             :     ELSE
     117        1916 :        e_lo = -9.99 
     118             :     ENDIF
     119             :     ! calculate core
     120        3996 :     e  = (e_up+e_lo)/2
     121        3996 :     fn = REAL(nqn) ; fl = REAL(l) ; fj = fl + 0.5
     122             :     CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
     123        3996 :          rn,d,msh,vrd, e, f(:,1),f(:,2),ierr)
     124        3996 :     IF (lo.AND.l>0) THEN
     125         116 :        e1  = (e_up+e_lo)/2
     126         116 :        fn = REAL(nqn) ; fl = REAL(l) ; fj = fl-0.5
     127             :        CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
     128         116 :             rn,d,msh,vrd, e1, f(:,1),f(:,2),ierr)
     129         116 :        e = (2.0*e + e1 ) / 3.0      
     130             :     ENDIF
     131             :     
     132             : 
     133             :    
     134        3996 :   END FUNCTION priv_method1
     135             : 
     136           0 :   REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr,e_lo,e_up)RESULT(e)
     137             :     USE m_types_setup
     138             :     USE m_types_mpi
     139             :     USE m_radsra
     140             :     USE m_differ
     141             :     USE m_constants
     142             :     IMPLICIT NONE
     143             :     LOGICAL,INTENT(IN):: lo
     144             :     INTEGER,INTENT(IN):: l,n,nqn,jsp
     145             :     REAL,INTENT(OUT)  :: e_lo,e_up
     146             :     TYPE(t_atoms),INTENT(IN)::atoms
     147             :     TYPE(t_mpi),INTENT(IN)  ::mpi
     148             :     REAL,INTENT(IN):: vr(:)
     149             : 
     150             :     INTEGER j,ilo,i
     151             :     INTEGER nodeu,node,ierr,msh
     152             :     REAL   lnd,e_up_temp,e_lo_temp,large_e_step
     153             :     REAL   d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
     154             :     LOGICAL start
     155             :     !     ..
     156             :     !     .. Local Arrays .. 
     157             : 
     158           0 :     REAL, ALLOCATABLE :: f(:,:),vrd(:)
     159             :     CHARACTER(LEN=20)    :: attributes(6)
     160             :     
     161           0 :     c=c_light(1.0)
     162             : 
     163             :     !Core potential setup done for each n,l now 
     164           0 :     d = EXP(atoms%dx(n))
     165             :     ! set up core-mesh
     166           0 :     rn = atoms%rmt(n)
     167           0 :     msh = atoms%jri(n)
     168           0 :     DO WHILE (rn < atoms%rmt(n) + 20.0)
     169           0 :        msh = msh + 1
     170           0 :        rn = rn*d
     171             :     ENDDO
     172           0 :     rn = atoms%rmsh(1,n)*( d**(msh-1) )
     173           0 :     ALLOCATE ( f(msh,2),vrd(msh) )
     174           0 :     f = 0.0
     175           0 :     vrd = 0.0
     176             :     ! extend core potential (linear with slope t1 / a.u.)
     177           0 :     vrd(:atoms%jri(n))=vr(:atoms%jri(n))
     178           0 :     t1=0.125
     179           0 :     t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
     180           0 :     rr = atoms%rmt(n)
     181           0 :     DO j = atoms%jri(n) + 1, msh
     182           0 :        rr = d*rr
     183           0 :        vrd(j) = rr*( t2 + rr*t1 )
     184             :     ENDDO
     185             :     ! search for branches
     186           0 :     node = ABS(nqn) - (l+1)
     187           0 :     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")
     188           0 :     e = 0.0 ! The initial value of e is arbitrary.
     189           0 :     large_e_step = 5.0 ! 5.0 Htr steps for coarse energy searches
     190             : 
     191             :     ! determine upper band edge
     192             :     ! Step 1: Coarse search for the band edge
     193             :     CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     194           0 :          atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     195           0 :     DO WHILE ( nodeu > node )
     196           0 :        e = e - large_e_step
     197             :        CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     198           0 :             atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     199             :     END DO
     200           0 :     DO WHILE ( nodeu <= node )
     201           0 :        e = e + large_e_step
     202             :        CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     203           0 :             atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     204             :     END DO
     205           0 :     e_up_temp = e
     206           0 :     e_lo_temp = e - large_e_step
     207             :     ! Step 2: Fine band edge determination by bisection search
     208           0 :     DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
     209           0 :        e = (e_up_temp + e_lo_temp) / 2.0
     210             :        CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     211           0 :             atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     212           0 :        IF (nodeu > node) THEN
     213             :           e_up_temp = e
     214             :        ELSE
     215           0 :           e_lo_temp = e
     216             :        END IF
     217             :     END DO
     218           0 :     e_up = (e_up_temp + e_lo_temp) / 2.0
     219           0 :     e    = e_up
     220             : 
     221             :     ! determine lower band edge
     222           0 :     IF (node == 0) THEN
     223           0 :        e_lo = -49.99
     224             :     ELSE
     225             :        ! Step 1: Coarse search for the band edge
     226           0 :        nodeu = node
     227           0 :        DO WHILE ( nodeu >= node )
     228           0 :           e = e - large_e_step
     229             :           CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     230           0 :                atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     231             :        ENDDO
     232           0 :        e_up_temp = e + large_e_step
     233           0 :        e_lo_temp = e
     234             :        ! Step 2: Fine band edge determination by bisection search
     235           0 :        DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
     236           0 :           e = (e_up_temp + e_lo_temp) / 2.0
     237             :           CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     238           0 :                atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     239           0 :           IF (nodeu < node) THEN
     240             :              e_lo_temp = e
     241             :           ELSE
     242           0 :              e_up_temp = e
     243             :           END IF
     244             :        END DO
     245           0 :        e_lo = (e_up_temp + e_lo_temp) / 2.0
     246             :     END IF
     247             : 
     248             : 
     249             :     ! determince notches by intersection
     250           0 :     ldmt= -99.0 !ldmt = logarithmic derivative @ MT boundary
     251           0 :     lnd = -l-1
     252           0 :     DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07) 
     253           0 :        e = (e_up+e_lo)/2
     254             :        CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
     255           0 :             atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
     256           0 :        ldmt = dus/us
     257           0 :        IF( ldmt .GT. lnd) THEN
     258           0 :           e_lo = e
     259           0 :        ELSE IF( ldmt .LT. lnd ) THEN
     260           0 :           e_up = e
     261             :           e_lo = e_lo
     262             :        END IF
     263             :     END DO
     264             : 
     265           0 :      END FUNCTION priv_method2
     266             :   
     267             : END MODULE m_find_enpara

Generated by: LCOV version 1.13