LCOV - code coverage report
Current view: top level - wannier - wann_amn.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 146 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 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_amn
       8             :       use m_juDFT
       9             :       contains
      10           0 :       subroutine wann_amn(
      11           0 :      >               chi,nslibd,nwfsd,ntypd,nlod,llod,llo,nlo,
      12           0 :      >               lmaxd,jmtd,lmd,neq,natd,ikpt,nbnd,
      13           0 :      >               rmsh,rmt,jri,dx,lmax,
      14           0 :      >               us,dus,uds,duds,flo,
      15           0 :      >               ff,gg,acof,bcof,ccof,l_nocosoc,jspin,
      16             :      &               l_amn2_in,
      17           0 :      <               amn,
      18             :      >               bkpt) 
      19             : c*********************************************************************
      20             : c...here: twf - trial wf                      Y. Mokrousov June'06
      21             : c..                                                  and August'06
      22             : 
      23             : c...For the moment calculates only the overlaps inside the spheres.
      24             : 
      25             : c...As an input file the routine reads the file 'proj', where
      26             : c...it reads the following parameters:
      27             : c...nwfs: total number of wfs
      28             : c...For each atom (not atom type) the following values
      29             : c...are specified:
      30             : c..   lwf, mrwf - specify the angular part of the twf (see tutorial)
      31             : c..   rwf       - specifies the radial index of the twf  
      32             : c..   alpha,beta,gamma - the euler angles of the twf 
      33             : c..   zona      - the diffusivity of the twf
      34             : c..   regio     - 2*Rmt for instance, will calculate the radial 
      35             : c..               integrals of the twf and the wavefunction in the
      36             : c..               sphere around the atom center with the radius of
      37             : c..               2*Rmt. Outside the MT the integrals with the   
      38             : c..               radial part of the planewave (in terms of bessels)
      39             : c..               will be calculated without seeing other MTs, vacuum
      40             : c..               boundaries whatsoever. But for the moment works
      41             : c..               only for 1*Rmt (so regio=1.0)
      42             : c*********************************************************************
      43             : c     l_amn2_in=.true. => Variant of standard wann_amn:
      44             : c     Projections are defined in file pro2.1/pro2.2 
      45             : c     Includes the possibility to define displacements 
      46             : c     of the MT-localized
      47             : c     trial orbitals onto atoms in neighboring unit cells. Otherwise
      48             : c     equal to standard case. May
      49             : c     be useful for covalent systems: The proj 
      50             : c     file specifies the localized
      51             : c     trial orbital for one partner of the covalent bond 
      52             : c     while proj2 does 
      53             : c     the same for the second. In many cases the second partner will be
      54             : c     located in a neighboring unit cell, however.
      55             : c     Frank Freimuth, March 2007
      56             : c************************************************************************* 
      57             : c     Inclusion of Noco&&soc; tidied up version
      58             : c     Frank Freimuth
      59             : c*********************************************************************
      60             :       use m_intgr, only : intgr3
      61             :       use m_constants, only: pimach, ImagUnit
      62             :       use m_dwigner
      63             :       use m_wann_tlmw
      64             :       use m_wann_rad_twf
      65             :       use m_eulerrot
      66             :       implicit none
      67             : 
      68             :       integer, intent (in)    :: nwfsd,nslibd,ntypd,nlod,llod
      69             :       integer, intent (in)    :: jmtd,lmaxd,lmd,natd,ikpt,nbnd
      70             :       integer, intent (in)    :: jri(ntypd),nlo(ntypd),llo(nlod,ntypd)
      71             :       integer, intent (in)    :: neq(ntypd),lmax(ntypd),jspin
      72             :       real,    intent (in)    :: rmsh(jmtd,ntypd),rmt(ntypd),dx(ntypd)
      73             :       real,    intent (in)    :: us(0:lmaxd,ntypd),dus(0:lmaxd,ntypd)
      74             :       real,    intent (in)    :: uds(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
      75             :       real,    intent (in)    :: ff(ntypd,jmtd,2,0:lmaxd)
      76             :       real,    intent (in)    :: gg(ntypd,jmtd,2,0:lmaxd)
      77             :       real,    intent (in)    :: flo(ntypd,jmtd,2,nlod)
      78             :       complex, intent (in)    :: acof(nslibd,0:lmd,natd)
      79             :       complex, intent (in)    :: bcof(nslibd,0:lmd,natd)
      80             :       complex, intent (in)    :: ccof(-llod:llod,nslibd,nlod,natd)
      81             :       logical, intent (in)    :: l_nocosoc
      82             :       logical, intent (inout) :: l_amn2_in
      83             :       complex, intent (inout) :: amn(nbnd,nwfsd)
      84             :       real,intent(in),optional:: bkpt(3)
      85             :       complex, intent (in)    :: chi
      86             : c...local
      87             :       integer          :: nwf,nwfs,nat,j,ntyp,ne,l,m,lm,iatom,i,mp,lo
      88           0 :       integer          :: ind(nwfsd),ntp(natd),banddummy
      89           0 :       integer          :: lwf(nwfsd),mrwf(nwfsd),rwf(nwfsd),spi(nwfsd)
      90           0 :       real             :: alpha(nwfsd),beta(nwfsd),gamma(nwfsd)
      91           0 :       real             :: jwf(nwfsd),jmwf(nwfsd),pi
      92           0 :       real,allocatable :: posshifts(:,:)
      93           0 :       real,allocatable :: weights(:)
      94           0 :       real             :: zona(nwfsd),regio(nwfsd),amx(3,3,nwfsd)
      95             :       real             :: rr,vl,vld,r0,vlo
      96           0 :       real             :: rads(nwfsd,0:3,jmtd,2),vlpr(jmtd),vlprd(jmtd)
      97           0 :       complex          :: tlmwf(0:3,-3:3,nwfsd),tlmwft(0:3,-3:3,nwfsd)
      98             :       logical          :: l_oldproj,l_amn2,l_file
      99           0 :       complex          :: d_wgn(-3:3,-3:3,1:3),wign(-3:3,-3:3,3,nwfsd)
     100             :       complex          :: factor
     101             :       real             :: arg
     102             :       real             :: mrott(3,3),bmatt(3,3),imx(3,3)
     103             :       character(len=2) :: spin12(0:2)
     104             :       data spin12/'  ', '.1', '.2'/
     105             :       character(len=6) :: filename
     106             : 
     107           0 :       l_amn2=l_amn2_in
     108           0 :       if(.not.l_amn2_in)then
     109           0 :          inquire(file='pro2.1',exist=l_amn2_in)
     110             :       endif
     111           0 :       pi=pimach()
     112             : c..generates an array giving the atom type for each atom
     113           0 :       ntp(:) = 0
     114             :       iatom = 0
     115           0 :       do ntyp = 1,ntypd
     116           0 :          do nat = 1,neq(ntyp)
     117           0 :             iatom = iatom + 1
     118           0 :             ntp(iatom) = ntyp
     119             :          enddo
     120             :       enddo 
     121             : 
     122           0 :       if(l_amn2)then
     123           0 :        do j=jspin,1,-1
     124           0 :          inquire(file=trim('pro2'//spin12(j)),exist=l_file)
     125           0 :          if(l_file)then
     126           0 :             filename='pro2'//spin12(j)
     127           0 :             exit
     128             :          endif
     129             :        enddo
     130             :       else   
     131             : c..reading the proj.1 / proj.2 / proj file
     132           0 :        do j=jspin,0,-1
     133           0 :          inquire(file=trim('proj'//spin12(j)),exist=l_file)
     134           0 :          if(l_file)then
     135           0 :             filename='proj'//spin12(j)
     136           0 :             exit
     137             :          endif
     138             :        enddo
     139             :       endif
     140             : 
     141           0 :       if(l_file)then
     142           0 :         open (203,file=trim(filename),status='old')
     143           0 :         rewind (203)
     144             :       else
     145           0 :          CALL juDFT_error("no proj/proj.1/proj.2",calledby="wann_amn")
     146             :       endif  
     147             : 
     148           0 :       if(l_amn2)then
     149           0 :        allocate(posshifts(3,nwfsd))
     150           0 :        allocate(weights(nwfsd))
     151             :       endif    
     152             : 
     153           0 :       if(l_nocosoc)then
     154           0 :        read (203,*)nwfs,banddummy,l_oldproj
     155           0 :        if(.not.l_oldproj)then
     156           0 :         do nwf=1,nwfs
     157           0 :          read(203,*)ind(nwf),lwf(nwf),jwf(nwf),jmwf(nwf),rwf(nwf)
     158           0 :          read(203,*)alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
     159           0 :          if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
     160             :         enddo !nwf
     161             :        else
     162           0 :         do nwf = 1,nwfs
     163             :          read (203,*) 
     164           0 :      &            ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf),spi(nwf)
     165             :          read (203,*) 
     166           0 :      &            alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
     167           0 :          if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
     168             :         enddo !nwf
     169             :        endif !oldproj
     170             :       else
     171           0 :        read (203,*) nwfs
     172           0 :        do nwf = 1,nwfs
     173             :          read (203,*) 
     174           0 :      &            ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf)
     175             :          read (203,*) 
     176           0 :      &            alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
     177           0 :          if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
     178             :        enddo !nwf
     179             :       endif !l_nocosoc
     180           0 :       rewind (203)
     181           0 :       close (203)
     182             : 
     183           0 :       if (ikpt.eq.1) then
     184           0 :       write (6,*) 'Number of trial WFs:',nwfs
     185           0 :       write (6,*)
     186           0 :       do nwf = 1,nwfs
     187           0 :         write (6,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
     188           0 :         write (6,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
     189           0 :         write (6,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
     190           0 :         write (6,*) 'alpha=',alpha(nwf),
     191           0 :      &          ' beta=',beta(nwf),' gamma=',gamma(nwf)
     192           0 :         write (6,*)
     193             :       enddo 
     194             :       endif
     195             : 
     196             : c..generating the radial twf function
     197             : 
     198           0 :       rads(:,:,:,:) = 0.
     199             : 
     200             :       call wann_rad_twf(
     201             :      >         nwfs,jmtd,natd,ind,rwf,zona,regio,
     202             :      >         us,dus,uds,duds,ff,gg,lmaxd,ikpt,
     203             :      >         ntypd,ntp,jri,rmsh,dx,
     204             :      >         nlod,flo,llo,nlo, 
     205           0 :      <         rads)
     206             : 
     207           0 :       open (100,file='rads')
     208           0 :       do i = 1,jmtd
     209             : c       write (100,'(i3,2x,4f10.6)') i,rads(1,0:3,i,1)
     210           0 :         write (100,'(f10.6,2x,4f10.6)') rmsh(i,1),rads(1,0:3,i,1)
     211             :       enddo
     212           0 :       close(100)
     213             : 
     214             : c..now generate the coefficients in the expansion in lattice 
     215             : c..harmonics of the angular part of the twf
     216             :  
     217           0 :       tlmwft(:,:,:) = cmplx(0.,0.)
     218           0 :       tlmwf(:,:,:)  = cmplx(0.,0.)
     219             :       
     220           0 :       if(l_nocosoc.and..not.l_oldproj)then
     221           0 :          call soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
     222             :       else
     223             :          call wann_tlmw(
     224             :      >          nwfs,lwf,mrwf,
     225             :      >          l_nocosoc,spi,jspin,
     226           0 :      <          tlmwf)  
     227             :       endif
     228             : 
     229           0 :       call eulerrot(nwfs,alpha,beta,gamma,amx)
     230             : 
     231           0 :       imx(:,:) = 0. ; imx(1,1) = 1. ; imx(2,2) = 1. ; imx(3,3) = 1.
     232             : 
     233             : c..performing the wigner rotation
     234             : c..These rotations are specified in the proj file in terms of the 
     235             : c..euler angles. The wave functions inside the muffin-tins are represen-
     236             : c..-ted in terms of the global-frame-lattice-harmonics, therefore, 
     237             : c..the wigner transformation has to be applied to tlmwf, resulting
     238             : c..in the tlmwft matrix, which is suitable for further use  
     239             : 
     240             : c..given the euler angles, the following procedure has to be 
     241             : c..performed in order to match the two local frames:
     242             : c..1. Rotation around the z-axis by alpha
     243             : c..2. Rotation around the x-axis by beta
     244             : c..3. Rotation around the z-axis by gamma again.
     245           0 :       call d_wigner(nwfs,amx,imx,3,wign)
     246             : 
     247             : c..now we transform the tlmwf coefficients
     248             : 
     249           0 :       do nwf = 1,nwfs
     250           0 :          tlmwft(0,:,nwf) = tlmwf(0,:,nwf)
     251           0 :          do l = 1,3
     252           0 :             do m = -l,l
     253           0 :                do mp = -l,l
     254             :                   tlmwft(l,m,nwf) = tlmwft(l,m,nwf) + 
     255           0 :      +                     wign(mp,m,l,nwf)*tlmwf(l,mp,nwf)
     256             :                enddo 
     257             :             enddo
     258             :          enddo         
     259             :       enddo
     260             : 
     261             : c..calculating the amn matrix
     262             : 
     263           0 :       vlpr(:) = 0. ; vlprd(:) = 0.
     264             : 
     265             : c...sum by wfs, each of them is localized at a certain mt
     266           0 :       do nwf = 1,nwfs
     267           0 :          if(l_amn2)then
     268           0 :            arg=-bkpt(1)*posshifts(1,nwf)
     269           0 :            arg=arg-bkpt(2)*posshifts(2,nwf)
     270           0 :            arg=arg-bkpt(3)*posshifts(3,nwf)
     271           0 :            arg=2*pi*arg
     272           0 :            factor=cmplx(cos(arg),sin(arg))*weights(nwf)
     273             :          else
     274             :            factor=1.0
     275             :          endif
     276           0 :          factor = factor*chi
     277             : 
     278           0 :          nat = ind(nwf)
     279           0 :          ntyp = ntp(nat)
     280             : c...sum by bands
     281           0 :          do ne = 1,nslibd
     282             : c...sum by l,m
     283           0 :             do l = 0,min(lmax(ntyp),3)
     284           0 :                do j = 1,jri(ntyp)
     285             :                   vlpr(j) = ff(ntyp,j,1,l)*rads(nwf,l,j,1)+
     286           0 :      +                      ff(ntyp,j,2,l)*rads(nwf,l,j,2) 
     287             :                   vlprd(j) = gg(ntyp,j,1,l)*rads(nwf,l,j,1) +
     288           0 :      +                       gg(ntyp,j,2,l)*rads(nwf,l,j,2)
     289           0 :                   if (rwf(nwf).gt.0)then
     290           0 :                      vlpr(j) = vlpr(j)*rmsh(j,ntyp)
     291           0 :                      vlprd(j) = vlprd(j)*rmsh(j,ntyp)
     292             :                   endif
     293             : 
     294             :                enddo
     295             : 
     296             : c..these integrations are not necessary if rads is the lin.comb. 
     297             : c..of the u_l and \dot{u}_l, but are necessary for other ways of
     298             : c..constructing the radial part, therefore, we do it anyway
     299             :  
     300           0 :                r0 = rmsh(1,ntyp)
     301           0 :                call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vl)
     302           0 :                call intgr3(vlprd,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vld)
     303           0 :                do m = -l,l
     304           0 :                   lm = l*(l+1) + m
     305             :                   amn(ne,nwf) = amn(ne,nwf) +
     306             :      +               tlmwft(l,m,nwf)*conjg(( acof(ne,lm,nat)*vl + 
     307           0 :      +               bcof(ne,lm,nat)*vld )*(ImagUnit)**l)*factor
     308             :      
     309             :                enddo
     310             :             enddo
     311             : 
     312             : c..local orbitals
     313           0 :             if (nlo(ntyp).ge.1) then
     314           0 :                do lo = 1,nlo(ntyp)
     315           0 :                   l = llo(lo,ntyp)
     316           0 :                   do j = 1,jri(ntyp)
     317             :                      vlpr(j) = flo(ntyp,j,1,lo)*rads(nwf,l,j,1) +
     318           0 :      +                         flo(ntyp,j,2,lo)*rads(nwf,l,j,2)
     319           0 :                      if (rwf(nwf).gt.0)then
     320           0 :                      vlpr(j) = vlpr(j)*rmsh(j,ntyp)
     321             :                      endif
     322             :                   enddo
     323             : 
     324             : 
     325           0 :                   r0 = rmsh(1,ntyp)
     326           0 :                   call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vlo)
     327           0 :                   do m = -l,l
     328             :                      amn(ne,nwf) = amn(ne,nwf) +
     329             :      +          tlmwft(l,m,nwf)*conjg((ccof(m,ne,lo,nat)*vlo)*
     330           0 :      *          (ImagUnit)**l)*factor
     331             :                   enddo
     332             :                enddo
     333             :             endif 
     334             : 
     335             :          enddo
     336             :       enddo
     337             : 
     338           0 :       return
     339             : 
     340           0 :       end subroutine wann_amn
     341             : 
     342             : c*********************************************************************
     343             : c     Calculate the expansion of the angular part of the trial orbital
     344             : c     in terms of spherical harmonics.
     345             : c     Version for Soc: The trial orbital is a spinor.
     346             : c     Frank Freimuth, June 2007
     347             : c*********************************************************************
     348           0 :       subroutine soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
     349             :       use m_clebsch
     350             :       implicit none
     351             : 
     352             :       integer, intent (in)  :: nwfs
     353             :       integer, intent (in)  :: lwf(nwfs)
     354             :       real,    intent (in)  :: jwf(nwfs),jmwf(nwfs)
     355             :       integer, intent (in)  :: jspin
     356             :       complex, intent (out) :: tlmwf(0:3,-3:3,nwfs)
     357             : 
     358             :       integer nwf,l,m
     359             :       complex tlm(0:3,-3:3,1:7)
     360             :       real j,jm
     361             : 
     362             : 
     363           0 :       do nwf=1,nwfs
     364           0 :          l=lwf(nwf)
     365           0 :          IF(l<0) CALL juDFT_error("not yet implemented",calledby
     366           0 :      +        ="wann_amn")
     367           0 :          j=jwf(nwf)
     368           0 :          jm=jmwf(nwf)
     369           0 :          IF(j<0) CALL juDFT_error("jwf",calledby ="wann_amn")
     370           0 :          IF( ABS(jm) - j  >1e-10) CALL juDFT_error("jmwf",calledby
     371           0 :      +        ="wann_amn")
     372           0 :          if( abs( l + 0.5 -j ).gt. 1e-10 .and.
     373             :      &       abs( l - 0.5 -j ).gt. 1e-10 )
     374             :      &     CALL juDFT_error ('regula trianguli violata',
     375           0 :      &                     calledby="wann_amn")
     376           0 :          tlmwf(0:3,-3:3,nwf)=cmplx(0.0,0.0)
     377           0 :          do m=-l,l
     378           0 :             tlmwf(l,m,nwf)=clebsch(real(l),0.5,real(m),1.5-jspin,j,jm)
     379             :          enddo
     380             :       enddo
     381             : 
     382           0 :       end subroutine soc_tlmw
     383             : 
     384             :       end module m_wann_amn

Generated by: LCOV version 1.13