LCOV - code coverage report
Current view: top level - wannier - wann_rad_twf.f (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 31 119 26.1 %
Date: 2024-03-29 04:21:46 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_wann_rad_twf
       8             :       use m_juDFT
       9             :       contains
      10           8 :       subroutine wann_rad_twf(
      11           8 :      >               nwfs,jmtd,natd,ind,rwf,zona,regio,
      12           8 :      >               us,dus,uds,duds,ff,gg,lmaxd,ikpt,
      13           8 :      >               ntypd,ntp,jri,rmsh,dx,
      14           8 :      >               nlod,flo,llo,nlo,
      15           8 :      <               rads)
      16             : c*********************************************************************
      17             : c..this subroutine calculates on the radial grid the radial function
      18             : c..corresponding to the twf, its mt sphere parameters, diffusivity
      19             : c..and region
      20             : c
      21             : c..Y. Mokrousov
      22             : c
      23             : c..tidied-up && extension of rwf-parameter
      24             : c..Frank Freimuth
      25             : c*********************************************************************
      26             : 
      27             :       USE m_constants
      28             :       use m_intgr, only : intgr3
      29             : 
      30             :       implicit none
      31             : 
      32             :       integer, intent (in)  :: nwfs,natd,jmtd,ntypd,lmaxd,ikpt,nlod
      33             :       integer, intent (in)  :: rwf(nwfs),ind(nwfs),ntp(natd),jri(ntypd)
      34             :       real,    intent (in)  :: zona(nwfs),regio(nwfs),rmsh(jmtd,ntypd)
      35             :       real,    intent (in)  :: us(0:lmaxd,ntypd),dus(0:lmaxd,ntypd)
      36             :       real,    intent (in)  :: uds(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
      37             :       real,    intent (in)  :: ff(ntypd,jmtd,2,0:lmaxd) 
      38             :       real,    intent (in)  :: gg(ntypd,jmtd,2,0:lmaxd),dx(ntypd) 
      39             :       real,    intent (in)  :: flo(ntypd,jmtd,2,nlod)
      40             :       integer, intent (in)  :: llo(nlod,ntypd)
      41             :       integer, intent (in)  :: nlo(ntypd)
      42             :       real,    intent (out) :: rads(nwfs,0:3,jmtd,2)
      43             : 
      44             :       integer :: nwf,nat,ntyp,j,n,l,lo
      45             :       real    :: rr,alpha,const,fact,wronk,radi
      46           8 :       real    :: acft,bcft,aa,bb,a1,b1,rho,radf(jmtd)
      47             : 
      48           8 :       call timestart("wann_rad_twf")
      49             : 
      50          72 :       do nwf = 1,nwfs
      51             : 
      52          64 :        if (abs(regio(nwf)-1.00000).ge.0.00001) then
      53             :           CALL juDFT_error("no projections outside muffin tins",calledby
      54           0 :      +         ="wann_rad_twf")
      55             :        endif
      56             : 
      57          64 :        nat = ind(nwf)
      58          64 :        ntyp = ntp(nat)
      59          64 :        alpha = zona(nwf)
      60             : 
      61      301632 :        rads (nwf,:,:,:) = 0.
      62             :  
      63          64 :        if (rwf(nwf).eq.-5) then 
      64             : c*********************************************************************
      65             : c..For rwf.eq.-5 the following way of implementing is realized.
      66             : c..Suppose, the twf has its center on some atom. The 
      67             : c..self-consistent potential V_l(r) and the energy parameters
      68             : c..E_l are known. A linear combination of the u_l and \dot{u}_l
      69             : c..is constructed, with the coefficients A and B such that 
      70             : c..on the mt boundary A*u_l + B*\dot{u}_l is smoothly continuous
      71             : c..with the 1s function 2*((zona)**(1.5))*exp(-r*zona)
      72             : c*********************************************************************
      73           0 :          aa = 2.*(zona(nwf)**(1.5))
      74           0 :          bb = zona(nwf)
      75             :         
      76           0 :          a1 = aa*exp(-bb*rmsh(jri(ntyp),ntyp))
      77           0 :          b1 = -aa*bb*exp(-bb*rmsh(jri(ntyp),ntyp))
      78             : 
      79           0 :          do l = 0,3
      80             : 
      81           0 :             wronk = us(l,ntyp)*duds(l,ntyp)-uds(l,ntyp)*dus(l,ntyp)
      82             : 
      83           0 :             acft = (a1*duds(l,ntyp) - b1*uds(l,ntyp))/wronk
      84           0 :             bcft = (b1*us(l,ntyp) - a1*dus(l,ntyp))/wronk
      85             :  
      86           0 :             do j = 1,jri(ntyp)
      87             :                rads(nwf,l,j,:) = acft*ff(ntyp,j,:,l) +
      88           0 :      +                           bcft*gg(ntyp,j,:,l)
      89             :             enddo 
      90             : 
      91             :          enddo
      92             :        elseif(rwf(nwf).eq.0) then
      93             : c***********************************************************
      94             : c..For rwf.eq.0 the zona parameter mixes the radial function
      95             : c..and its energy derivative linearly.
      96             : c***********************************************************
      97         320 :          do l = 0,3
      98      120896 :             do j = 1,jri(ntyp)
      99             :                rads(nwf,l,j,:) = ff(ntyp,j,:,l)*(1.-abs(zona(nwf))) +
     100      361984 :      +                         zona(nwf)*gg(ntyp,j,:,l)
     101             :             enddo 
     102             :          enddo   
     103             : c**************************************
     104             : c..project onto local orbitals
     105             : c**************************************
     106             :        elseif(rwf(nwf).eq.-6)then
     107           0 :          IF(nlod<1) CALL juDFT_error("nlod<1",calledby="wann_rad_twf"
     108           0 :      +         )
     109           0 :           do l=0,3
     110           0 :            do j=1,jri(ntyp)
     111           0 :             rads(nwf,l,j,:)=flo(ntyp,j,:,1)
     112             :            enddo!j
     113             :           enddo!l
     114             :        elseif(rwf(nwf).eq.-7)then
     115           0 :          IF(nlod<2) CALL juDFT_error("nlod<2",calledby="wann_rad_twf"
     116           0 :      +         )
     117           0 :           do l=0,3
     118           0 :            do j=1,jri(ntyp)
     119           0 :             rads(nwf,l,j,:)=flo(ntyp,j,:,2)
     120             :            enddo!j
     121             :           enddo!l
     122             :        elseif(rwf(nwf).eq.-8)then
     123           0 :           IF(nlod<3) CALL juDFT_error("nlod<3",calledby ="wann_rad_twf"
     124           0 :      +         )
     125           0 :           do l=0,3
     126           0 :            do j=1,jri(ntyp)
     127           0 :             rads(nwf,l,j,:)=flo(ntyp,j,:,3)
     128             :            enddo!j
     129             :           enddo!l
     130             :        elseif(rwf(nwf).eq.-9)then
     131           0 :           rads(nwf,:,:,:)=0
     132           0 :           do lo=1,nlo(ntyp)
     133           0 :            l=llo(lo,ntyp)
     134           0 :            do j=1,jri(ntyp)
     135           0 :             rads(nwf,l,j,:)=flo(ntyp,j,:,lo)     
     136             :            enddo!j
     137             :           enddo !lo
     138             :         
     139           0 :        elseif((-5.lt.rwf(nwf)).and.(rwf(nwf).lt.0)) then
     140             : c************************************************************
     141             : c..use the radial with l=abs(rwf)-1 for all components of the
     142             : c..trial wave function
     143             : c************************************************************
     144           0 :            do l = 0,3
     145           0 :             do j = 1,jri(ntyp)
     146             :                rads(nwf,l,j,:) = 
     147             :      =             ff(ntyp,j,:,abs(rwf(nwf))-1)*(1.-abs(zona(nwf))) +
     148           0 :      +             zona(nwf)*gg(ntyp,j,:,abs(rwf(nwf))-1)
     149             :             enddo 
     150             :            enddo   
     151             : 
     152             : c**************************************************************
     153             : c..Hydrogen-like test orbitals.
     154             : c..here we have tabulated the radial functions from orbitron.
     155             : c..for 1 <= rwf <= 6 different components of the test orbital
     156             : c..are constructed from different radials.
     157             : c  For Hybrides (e.g. sp3) it might be better to use the same
     158             : c  radials for all components. See below (7 <= rwf <= 12)
     159             : c**************************************************************
     160             :        elseif (rwf(nwf).eq.1) then
     161             : c..n=1
     162           0 :           do j = 1,jri(ntyp)
     163           0 :              rho = 2.*alpha*rmsh(j,ntyp)
     164             :              rads(nwf,0,j,1) = 2.*((alpha)**(1.5))*
     165           0 :      *            exp(-rho/2.)            
     166             :           enddo  
     167             : 
     168             :        elseif (rwf(nwf).eq.2) then
     169             : c..n=2
     170           0 :           do j = 1,jri(ntyp)
     171           0 :              rho = alpha*(rmsh(j,ntyp))
     172             :              rads(nwf,0,j,1) = (1./(2.*sqrt(2.)))*
     173             :      *              (2.-rho)*
     174           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     175             :              rads(nwf,1,j,1) = (1./(2.*sqrt(6.)))*
     176             :      *              rho*
     177           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     178             :           enddo          
     179             : 
     180             :        elseif (rwf(nwf).eq.3) then
     181             : c..n=3
     182           0 :           do j = 1,jri(ntyp)
     183           0 :              rho = 2*alpha*rmsh(j,ntyp)/3.
     184             :              rads(nwf,0,j,1) = (1./(9.*sqrt(3.)))*
     185             :      *              (6. - 6.*rho + rho**2)*
     186           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     187             :              rads(nwf,1,j,1) = (1./(9.*sqrt(6.)))*
     188             :      *              (4. - rho)*rho*
     189           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     190             :              rads(nwf,2,j,1) = (1./(9.*sqrt(30.)))*
     191             :      *               rho*rho*
     192           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     193             :           enddo
     194             : 
     195             :        elseif (rwf(nwf).eq.4) then
     196             : c..n=4
     197           0 :           do j = 1,jri(ntyp)
     198           0 :              rho = 2*alpha*rmsh(j,ntyp)/4.
     199             :              rads(nwf,0,j,1) = (1./96.)*
     200             :      *              (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
     201           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     202             :              rads(nwf,1,j,1) = (1./(32.*sqrt(15.)))*
     203             :      *              rho*(20. - 10.*rho + rho**2)*
     204           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     205             :              rads(nwf,2,j,1) = (1./(96.*sqrt(5.)))*
     206             :      *              rho*rho*(6. - rho)*
     207           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     208             :              rads(nwf,3,j,1) = (1./(96.*sqrt(35.)))*
     209             :      *              rho*rho*rho*
     210           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     211             :           enddo
     212             : 
     213             :        elseif (rwf(nwf).eq.5) then
     214           0 :           do j = 1,jri(ntyp)
     215           0 :              rho = 2*alpha*rmsh(j,ntyp)/5.
     216             :              rads(nwf,0,j,1) = (1./(300.*sqrt(5.)))*
     217             :      *              (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
     218             :      +                            +(rho**4))*
     219           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     220             :              rads(nwf,1,j,1) = (1./(150.*sqrt(30.)))*
     221             :      *              rho*(120. - 90.*rho + 18.*(rho**2) - (rho**3))*
     222           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     223             :              rads(nwf,2,j,1) = (1./(150.*sqrt(70.)))*
     224             :      *              rho*rho*(42. - 14.*rho + (rho*rho))*
     225           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     226             :              rads(nwf,3,j,1) = (1./(300.*sqrt(70.)))*
     227             :      *              rho*rho*rho*(8. - rho)*
     228           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     229             :           enddo
     230             :         elseif (rwf(nwf).eq.6) then
     231             : 
     232           0 :           do j = 1,jri(ntyp)
     233           0 :              rho = 2*alpha*rmsh(j,ntyp)/6.
     234             :              rads(nwf,0,j,1) = (1./2160.*(sqrt(6.)))*
     235             :      *              (720. - 1800.*rho + 1200.*rho*rho-
     236             :      -               300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
     237           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     238             :              rads(nwf,1,j,1) = (1./432.*(sqrt(210.)))*
     239             :      *              rho*(840. - 840.*rho + 252.*rho*rho-
     240             :      -               28.*rho*rho*rho + (rho**4))*
     241           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     242             :           enddo
     243             : 
     244             :         elseif (rwf(nwf).eq.7) then
     245             : c..n=1
     246           0 :           do j = 1,jri(ntyp)
     247           0 :              rho = 2.*alpha*rmsh(j,ntyp)
     248           0 :              do l=0,3
     249             :              rads(nwf,l,j,1) = 2.*((alpha)**(1.5))*
     250           0 :      *            exp(-rho/2.)            
     251             :              enddo
     252             :           enddo  
     253             : 
     254             :         elseif (rwf(nwf).eq.8) then
     255             : c..n=2
     256           0 :           do j = 1,jri(ntyp)
     257           0 :              rho = alpha*(rmsh(j,ntyp))
     258           0 :              do l=0,3
     259             :              rads(nwf,l,j,1) = (1./(2.*sqrt(2.)))*
     260             :      *              (2.-rho)*
     261           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     262             :              enddo
     263             :           enddo          
     264             : 
     265             :         elseif (rwf(nwf).eq.9) then
     266             : c..n=3
     267           0 :           do j = 1,jri(ntyp)
     268           0 :              rho = 2*alpha*rmsh(j,ntyp)/3.
     269           0 :              do l=0,3
     270             :              rads(nwf,l,j,1) = (1./(9.*sqrt(3.)))*
     271             :      *              (6. - 6.*rho + rho**2)*
     272           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     273             :              enddo
     274             :           enddo
     275             : 
     276             :         elseif (rwf(nwf).eq.10) then
     277             : c..n=4
     278           0 :           do j = 1,jri(ntyp)
     279           0 :              rho = 2*alpha*rmsh(j,ntyp)/4.
     280           0 :              do l=0,3
     281             :              rads(nwf,l,j,1) = (1./96.)*
     282             :      *              (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
     283           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     284             :              enddo
     285             :           enddo
     286             :         elseif (rwf(nwf).eq.11) then
     287             : 
     288           0 :           do j = 1,jri(ntyp)
     289           0 :              rho = 2*alpha*rmsh(j,ntyp)/5.
     290           0 :              do l=0,3
     291             :              rads(nwf,l,j,1) = (1./(300.*sqrt(5.)))*
     292             :      *              (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
     293             :      +                            +(rho**4))*
     294           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     295             :              enddo
     296             :           enddo
     297             : 
     298             :         elseif (rwf(nwf).eq.12) then
     299             : 
     300           0 :           do j = 1,jri(ntyp)
     301           0 :              rho = 2*alpha*rmsh(j,ntyp)/6.
     302           0 :              do l=0,3
     303             :              rads(nwf,0,j,1) = 0*(1./2160.*(sqrt(6.)))*
     304             :      *              (720. - 1800.*rho + 1200.*rho*rho-
     305             :      -               300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
     306           0 :      *              ((alpha)**(1.5))*exp(-rho/2.)
     307             :              enddo
     308             :           enddo
     309             : 
     310             :         else 
     311             : 
     312             :            CALL juDFT_error("radial function is not tabulated" ,calledby
     313           0 :      +          ="wann_rad_twf")
     314             : 
     315             : 
     316             :        endif
     317             : 
     318          72 :        if (ikpt.eq.1) then
     319          40 :         do l = 0,3
     320       15104 :          do j = 1,jri(ntyp)
     321             :           radf(j) = rmsh(j,ntyp)*rmsh(j,ntyp)*rads(nwf,l,j,1)*
     322       15104 :      *                rads(nwf,l,j,1)
     323             : c         radf(j) = rads(nwf,l,j,1)*rads(nwf,l,j,1)
     324             :          enddo
     325          32 :          call intgr3(radf,rmsh(1,ntyp),dx(ntyp),jri(ntyp),radi)
     326          32 :          write (oUnit,*)
     327          32 :          write (oUnit,*) 'Wannier Function N:',nwf
     328          32 :          write (oUnit,*) 'angular momentum',l
     329          32 :          write (oUnit,*) 'radial function at the MT boundary:',
     330          64 :      &                rads(nwf,l,jri(ntyp),1)
     331          72 :          write (oUnit,*) 'norma =',radi
     332             :         enddo
     333             :        endif 
     334             :  
     335             :       enddo ! by twfs
     336             : 
     337           8 :       call timestop("wann_rad_twf")
     338           8 :       return
     339             : 
     340             :       end subroutine wann_rad_twf
     341             :       end module m_wann_rad_twf

Generated by: LCOV version 1.14