LCOV - code coverage report
Current view: top level - wannier - wann_rad_twf.f (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 26.1 % 119 31
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 1 1

            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 2.0-1