LCOV - code coverage report
Current view: top level - global - radflo.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 35 86 40.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          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_radflo
       8             :   USE m_juDFT
       9             : CONTAINS
      10         496 :   SUBROUTINE radflo(atoms, ntyp,jsp,ello,vr, f,g,mpi, usdus,&
      11         496 :        uuilon,duilon,ulouilopn,flo,lout_all)
      12             :     !
      13             :     !***********************************************************************
      14             :     ! generates the scalar relativistic wavefunctions (flo) needed for the
      15             :     ! local orbitals at atom type n for angular momentum l.
      16             :     ! the values of the function and the radial derivative on the sphere
      17             :     ! (ulos,dulos) boundaries, the overlap with the other radial functions
      18             :     ! (uulon,dulon) and between the different local orbitals (uloulopn) are
      19             :     ! also calculated.
      20             :     !
      21             :     ! ulos    : the value of the radial function of a local orbital
      22             :     !           at the muffin tin radius
      23             :     ! dulos   : the value of the radial derivative of the radial function
      24             :     !           function of a local orbital at the muffin tin radius
      25             :     ! uulon   : overlap integral between the radial functions of a local
      26             :     !           obital and the flapw radial function with the same l
      27             :     ! dulon   : overlap integral between the radial functions of a local
      28             :     !           obital and the energy derivative of the flapw radial
      29             :     !           function with the same l
      30             :     ! uloulopn: overlap integral between the radial functions of two
      31             :     !           different local orbitals
      32             :     !           (only needed if they have the same l)
      33             :     ! l_dulo  : whether we use a dot u as local orbital (see setlomap.F)
      34             :     !
      35             :     ! p.kurz jul. 1996 gb 2001
      36             :     !
      37             :     ! ulo_der : specifies the order of the energy derivative to be used
      38             :     !           (0, 1, 2, ... for primitive, first, second derivatives etc.)
      39             :     ! uuilon  : overlap integral between the radial functions of the
      40             :     !           integral (multiplied by ulo_der) of a local orbital and the
      41             :     !           flapw radial function with the same l
      42             :     ! duilon  : overlap integral between the radial functions of the
      43             :     !           integral of a local orbital and the energy derivative of the
      44             :     !           flapw radial function with the same l
      45             :     ! ulouilopn: overlap integral between the radial functions of the
      46             :     !           integral of a local orbital and another local orbital with
      47             :     !           the same l.
      48             :     !
      49             :     ! C. Friedrich Feb. 2005
      50             :     !***********************************************************************
      51             :     !
      52             :     USE m_intgr, ONLY : intgr0
      53             :     USE m_constants, ONLY : c_light
      54             :     USE m_radsra
      55             :     USE m_radsrdn
      56             :     USE m_differ
      57             : #include "cpp_double.h"
      58             :     USE m_types
      59             :     IMPLICIT NONE
      60             :     TYPE(t_usdus),INTENT(INOUT):: usdus !lo part is calculated here
      61             :     TYPE(t_mpi),INTENT(IN)     :: mpi
      62             :     TYPE(t_atoms),INTENT(IN)   :: atoms
      63             :     !     ..
      64             :     !     .. Scalar Arguments ..
      65             :     INTEGER, INTENT (IN) :: jsp,ntyp 
      66             :     LOGICAL, INTENT (IN), OPTIONAL :: lout_all
      67             :     !     ..
      68             :     !     .. Array Arguments ..
      69             :     REAL,    INTENT (IN) :: ello(atoms%nlod,atoms%ntype),vr(atoms%jmtd)
      70             :     REAL,    INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
      71             :     REAL,    INTENT (OUT):: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
      72             :     REAL,    INTENT (OUT):: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
      73             :     REAL,    INTENT (OUT):: flo(atoms%jmtd,2,atoms%nlod)
      74             :     !     ..
      75             :     !     .. Local Scalars ..
      76             :     INTEGER i,j,k,l,ilo,jlo,nodelo,noded ,ierr,msh
      77             :     REAL    rn,t1,t2,d ,rr,fn,fl,fj,e,uds,duds,ddn,c,rmt
      78             :     LOGICAL ofdiag, loutput
      79             :     !     ..
      80             :     !     .. Local Arrays ..
      81         992 :     REAL dulo(atoms%jmtd),ulo(atoms%jmtd),glo(atoms%jmtd,2),filo(atoms%jmtd,2,atoms%nlod)
      82         496 :     REAL, ALLOCATABLE :: f_rel(:,:),vrd(:)
      83         248 :     REAL help(atoms%nlod+2,atoms%nlod+2)
      84             :     !     ..
      85         248 :     c = c_light(1.0)
      86             :     !
      87             :     IF ( PRESENT(lout_all) ) THEN
      88         248 :        loutput = ( mpi%irank == 0 ) .OR. lout_all
      89             :     ELSE
      90             :        loutput = ( mpi%irank == 0 )
      91             :     END IF
      92         248 :     !$    loutput=.false.
      93             :     IF (loutput) THEN
      94             :        WRITE (6,FMT=8000)
      95             :     END IF
      96         248 :     ofdiag = .FALSE.
      97             :     !---> calculate the radial wavefunction with the appropriate
      98             :     !---> energy parameter (ello)
      99         640 :     DO ilo = 1,atoms%nlo(ntyp)
     100             :        CALL radsra(ello(ilo,ntyp),atoms%llo(ilo,ntyp),vr(:),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c,&
     101         392 :             usdus%ulos(ilo,ntyp,jsp),usdus%dulos(ilo,ntyp,jsp),nodelo,flo(:,1,ilo), flo(:,2,ilo))
     102             :        !
     103             :        !+apw+lo
     104         392 :        IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN
     105           0 :           IF (atoms%ulo_der(ilo,ntyp).LE.8) THEN
     106             :              !--->    calculate orthogonal energy derivative at e
     107           0 :              i = atoms%ulo_der(ilo,ntyp)
     108           0 :              IF(atoms%l_dulo(ilo,ntyp)) i=1
     109             :              CALL radsrdn(ello(ilo,ntyp),atoms%llo(ilo,ntyp),vr(:),atoms%rmsh(1,ntyp), atoms%dx(ntyp),&
     110           0 :                   atoms%jri(ntyp),c, uds,duds,ddn,noded,glo,filo(:,:,ilo), flo(:,:,ilo),usdus%dulos(ilo,ntyp,jsp),i)
     111           0 :              DO i=1,atoms%jri(ntyp)
     112           0 :                 flo(i,1,ilo) = glo(i,1)
     113           0 :                 flo(i,2,ilo) = glo(i,2)
     114             :              ENDDO
     115           0 :              nodelo = noded
     116           0 :              ddn    = SQRT(ddn)
     117           0 :              IF(atoms%l_dulo(ilo,ntyp)) ddn=1.0
     118           0 :              flo (:,:,ilo)   = flo (:,:,ilo)/ddn ! Normalize ulo (flo) if APW+lo is not used
     119           0 :              filo(:,:,ilo)   = filo(:,:,ilo)/ddn ! and scale its integral (filo) accordingly
     120           0 :              usdus%dulos(ilo,ntyp,jsp) = duds/ddn          !   (setabc1lo and slomat assume <flo|flo>=1)
     121           0 :              usdus%ulos (ilo,ntyp,jsp) =  uds/ddn          ! 
     122             :           ELSE
     123             :              !
     124             :              !          test:
     125             :              !
     126             :              ! set up core-mesh
     127           0 :              d = EXP(atoms%dx(ntyp))
     128           0 :              rmt = atoms%rmsh(1,ntyp)
     129           0 :              DO i = 1, atoms%jri(ntyp) - 1
     130           0 :                 rmt = rmt * d
     131             :              ENDDO
     132           0 :              rn = rmt
     133           0 :              msh = atoms%jri(ntyp)
     134           0 :              DO WHILE (rn < rmt + 20.0)
     135           0 :                 msh = msh + 1
     136           0 :                 rn = rn*d
     137             :              ENDDO
     138           0 :              rn = atoms%rmsh(1,ntyp)*( d**(msh-1) )
     139           0 :              ALLOCATE ( f_rel(msh,2),vrd(msh) )
     140             : 
     141             :              ! extend core potential (linear with slope t1 / a.u.)
     142             : 
     143           0 :              DO j = 1, atoms%jri(ntyp)
     144           0 :                 vrd(j) = vr(j)
     145             :              ENDDO
     146           0 :              t1=0.125
     147           0 :              t2 = vrd(atoms%jri(ntyp))/rmt - rmt*t1
     148           0 :              rr = rmt
     149           0 :              DO j = atoms%jri(ntyp) + 1, msh
     150           0 :                 rr = d*rr
     151           0 :                 vrd(j) = rr*( t2 + rr*t1 )
     152             :              ENDDO
     153           0 :              e = ello(ilo,ntyp)
     154           0 :              fn = 6.0 ; fl = 1.0 ; fj = 0.5
     155             :              CALL differ(fn,fl,fj,c,82.0,atoms%dx(ntyp),atoms%rmsh(1,ntyp),&
     156           0 :                   rn,d,msh,vrd, e, f_rel(1,1),f_rel(1,2),ierr)
     157             : 
     158           0 :              f_rel(:,1) = 2 * f_rel(:,1)
     159           0 :              f_rel(:,2) = 2 * f_rel(:,2)
     160           0 :              rn = atoms%rmsh(1,ntyp)
     161           0 :              DO i = 1, atoms%jri(ntyp) 
     162           0 :                 rn = rn * d
     163           0 :                 WRITE(123,'(5f20.15)') rn,f_rel(i,1),f_rel(i,2), f(i,1,1),f(i,2,1)
     164             :              ENDDO
     165             : 
     166           0 :              STOP
     167             :           ENDIF
     168             : 
     169             :        ENDIF
     170             :        !-apw+lo
     171             :        !
     172             :        !--->    calculate the overlap between these fcn. and the radial functions
     173             :        !--->    of the flapw basis with the same l
     174      279088 :        DO i = 1,atoms%jri(ntyp)
     175      278696 :           ulo(i) = f(i,1,atoms%llo(ilo,ntyp))*flo(i,1,ilo) + f(i,2,atoms%llo(ilo,ntyp))*flo(i,2,ilo)
     176      279088 :           dulo(i) = g(i,1,atoms%llo(ilo,ntyp))*flo(i,1,ilo) + g(i,2,atoms%llo(ilo,ntyp))*flo(i,2,ilo)
     177             :        END DO
     178         392 :        CALL intgr0(ulo, atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,jsp))
     179         392 :        CALL intgr0(dulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,jsp))
     180         392 :        IF (atoms%l_dulo(ilo,ntyp)) usdus%dulon(ilo,ntyp,jsp) = 0.0
     181             :        IF (loutput) THEN
     182             :           WRITE (6,FMT=8010) ilo,atoms%llo(ilo,ntyp),ello(ilo,ntyp),&
     183             :                usdus%ulos(ilo,ntyp,jsp),usdus%dulos(ilo,ntyp,jsp),nodelo,usdus%uulon(ilo,ntyp,jsp),&
     184             :                usdus%dulon(ilo,ntyp,jsp)
     185             :        END IF
     186             :        !
     187             :        !--->   case LO = energy derivative (ulo_der>=1):
     188             :        !--->   calculate the overlap between the LO-integral (filo) and the radial functions
     189         392 :        IF(atoms%ulo_der(ilo,ntyp).GE.1) THEN
     190           0 :           DO i=1,atoms%jri(ntyp)
     191           0 :              ulo(i) = f(i,1,atoms%llo(ilo,ntyp))*filo(i,1,ilo) + f(i,2,atoms%llo(ilo,ntyp))*filo(i,2,ilo)
     192           0 :              dulo(i) = g(i,1,atoms%llo(ilo,ntyp))*filo(i,1,ilo) + g(i,2,atoms%llo(ilo,ntyp))*filo(i,2,ilo)
     193             :           ENDDO
     194           0 :           CALL intgr0(ulo, atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),uuilon(ilo,ntyp))
     195           0 :           CALL intgr0(dulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),duilon(ilo,ntyp))
     196             :        ELSE
     197         392 :           uuilon(ilo,ntyp)=0
     198         392 :           duilon(ilo,ntyp)=0
     199             :        ENDIF
     200             :        !--->   calculate overlap between radial fcn. of different local
     201             :        !--->   orbitals (only if both have the same l)
     202             :        !         uloulopn(ilo,ilo,ntyp) = 1.0
     203             :        !         DO jlo = 1, (ilo-1)
     204        1712 :        DO jlo = 1, ilo
     205         928 :           IF (atoms%llo(ilo,ntyp).EQ.atoms%llo(jlo,ntyp)) THEN
     206      279088 :              DO i = 1,atoms%jri(ntyp)
     207      279088 :                 ulo(i) = flo(i,1,ilo)*flo(i,1,jlo) + flo(i,2,ilo)*flo(i,2,jlo)
     208             :              END DO
     209         392 :              CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uloulopn(ilo,jlo,ntyp,jsp))
     210         392 :              usdus%uloulopn(jlo,ilo,ntyp,jsp) = usdus%uloulopn(ilo,jlo,ntyp,jsp)
     211         392 :              ofdiag = .TRUE.
     212             :           ELSE
     213         144 :              usdus%uloulopn(ilo,jlo,ntyp,jsp) = 0.0
     214         144 :              usdus%uloulopn(jlo,ilo,ntyp,jsp) = 0.0
     215             :           END IF
     216             :        END DO
     217             :     END DO
     218             :     !
     219             :     !---> case: one of LOs = energy derivative (ulo_der>=1):
     220             :     !---> calculate overlap between LOs and integrals of LOs
     221         640 :     DO ilo = 1,atoms%nlo(ntyp)
     222        1320 :        DO jlo = 1,atoms%nlo(ntyp)
     223        1072 :           IF(atoms%ulo_der(jlo,ntyp).GE.1.AND.atoms%llo(ilo,ntyp).EQ.atoms%llo(jlo,ntyp)) THEN
     224           0 :              DO i = 1,atoms%jri(ntyp)
     225           0 :                 ulo(i) = flo(i,1,ilo)*filo(i,1,jlo) + flo(i,2,ilo)*filo(i,2,jlo)
     226             :              ENDDO
     227           0 :              CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ulouilopn(ilo,jlo,ntyp))
     228             :           ELSE
     229         680 :              ulouilopn(ilo,jlo,ntyp)=0.0
     230             :           ENDIF
     231             :        ENDDO
     232             :     ENDDO
     233             :     ! 
     234             :     IF ( (ofdiag).AND.(loutput) ) THEN
     235             :        WRITE (6,FMT=*) 'overlap matrix between different local orbitals'
     236             :        WRITE (6,FMT=8020) (i,i=1,atoms%nlo(ntyp))
     237             :        DO ilo = 1,atoms%nlo(ntyp)
     238             :           WRITE (6,FMT='(i3,40e13.6)') ilo, (usdus%uloulopn(ilo,jlo,ntyp,jsp),jlo=1,atoms%nlo(ntyp))
     239             : 
     240             :        END DO
     241             :     END IF
     242             :     !
     243             :     !
     244             :     !     Diagonalize overlap matrix of normalized MT functions
     245             :     !       to check linear dependencies
     246             :     !       help is overlap matrix of all MT functions (flapw & LO)
     247             :     IF (.FALSE.) THEN
     248             :        DO l=0,MAXVAL(atoms%llo(1:atoms%nlo(ntyp),ntyp))
     249             :           IF(ALL(atoms%llo(1:atoms%nlo(ntyp),ntyp).NE.l)) CYCLE
     250             :           help(1,1)=1.0
     251             :           help(1,2)=0.0
     252             :           help(2,1)=0.0
     253             :           help(2,2)=1.0
     254             :           DO i=1,atoms%jri(ntyp)
     255             :              ulo(i)=g(i,1,l)**2+g(i,2,l)**2
     256             :           ENDDO
     257             :           CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ddn)
     258             :           ddn=SQRT(ddn)
     259             :           j=2
     260             :           k=2
     261             :           DO ilo=1,atoms%nlo(ntyp)
     262             :              IF(atoms%llo(ilo,ntyp).EQ.l) THEN
     263             :                 j=j+1
     264             :                 help(1,j)=usdus%uulon(ilo,ntyp,jsp)
     265             :                 help(2,j)=usdus%dulon(ilo,ntyp,jsp)/ddn
     266             :                 k=2
     267             :                 DO jlo=1,ilo
     268             :                    IF(atoms%llo(jlo,ntyp).EQ.l) THEN
     269             :                       k=k+1
     270             :                       help(k,j)=usdus%uloulopn(ilo,jlo,ntyp,jsp)
     271             :                    ENDIF
     272             :                 ENDDO
     273             :              ENDIF
     274             :           ENDDO
     275             :           IF ( loutput ) THEN
     276             :              WRITE(6,'(A,I2)')&
     277             :                   &      'Overlap matrix of normalized MT functions '//&
     278             :                   &      'and its eigenvalues for l =',l
     279             :              DO i=1,j
     280             :                 WRITE(6,'(30F11.7)') (help(k,i),k=1,i)
     281             :              ENDDO
     282             :           END IF
     283             :           CALL CPP_LAPACK_ssyev('N','U',j,help,atoms%nlod+2,ulo,dulo, i)
     284             :           IF(i/=0)  CALL juDFT_error("ssyev failed.",calledby ="radflo")
     285             :           IF ( loutput ) THEN
     286             :              WRITE(6,'(30F11.7)') (ulo(i),i=1,j)
     287             :           END IF
     288             :        ENDDO
     289             :     ENDIF ! .false.
     290             :     !
     291             : 
     292             : 8000 FORMAT (/,t20,'radial function for local orbitals',/,t2,'lo',t6,&
     293             :          &       'l',t11,'energy',t29,'value',t42,'derivative',t56,'nodes',&
     294             :          &       t63,'ovlp with u',t78,'ovlp with udot')
     295             : 8010 FORMAT (i3,i3,f10.5,5x,1p,2e16.7,i5,2e16.7)
     296             : 8020 FORMAT (' lo',i9,19i15)
     297         248 :     RETURN
     298             :   END SUBROUTINE radflo
     299             : END MODULE m_radflo

Generated by: LCOV version 1.13