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

Generated by: LCOV version 1.13