LCOV - code coverage report
Current view: top level - wannier - wann_amn.f (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 67.3 % 165 111
Test Date: 2025-06-14 04:34:23 Functions: 50.0 % 2 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_amn
       8              :       use m_juDFT
       9              :       contains
      10            8 :       subroutine wann_amn(
      11            8 :      >               chi,nslibd,nwfsd,ntypd,nlod,llod,llo,nlo,
      12            8 :      >               lmaxd,jmtd,lmd,neq,natd,ikpt,nbnd,
      13            8 :      >               rmsh,rmt,jri,dx,lmax,
      14            8 :      >               us,dus,uds,duds,flo,
      15            8 :      >               ff,gg,acof,bcof,ccof,l_nocosoc,jspin,
      16              :      &               l_amn2_in,
      17            8 :      <               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
      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            8 :       integer          :: ind(nwfsd),ntp(natd),banddummy
      89            8 :       integer          :: lwf(nwfsd),mrwf(nwfsd),rwf(nwfsd),spi(nwfsd)
      90            8 :       real             :: alpha(nwfsd),beta(nwfsd),gamma(nwfsd)
      91            8 :       real             :: jwf(nwfsd),jmwf(nwfsd),pi
      92            8 :       real,allocatable :: posshifts(:,:)
      93            8 :       real,allocatable :: weights(:)
      94            8 :       real             :: zona(nwfsd),regio(nwfsd),amx(3,3,nwfsd)
      95              :       real             :: rr,vl,vld,r0,vlo
      96            8 :       real             :: rads(nwfsd,0:3,jmtd,2),vlpr(jmtd),vlprd(jmtd)
      97            8 :       complex          :: tlmwf(0:3,-3:3,nwfsd),tlmwft(0:3,-3:3,nwfsd)
      98              :       logical          :: l_oldproj,l_amn2,l_file
      99            8 :       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            8 :       call timestart("wann_amn")
     108           72 :       spi(:) = 0
     109              : 
     110            8 :       l_amn2=l_amn2_in
     111            8 :       if(.not.l_amn2_in)then
     112            8 :          inquire(file='pro2.1',exist=l_amn2_in)
     113              :       endif
     114            8 :       pi=pimach()
     115              : c..generates an array giving the atom type for each atom
     116            8 :       call timestart("gen ntp")
     117           24 :       ntp(:) = 0
     118              :       iatom = 0
     119           16 :       do ntyp = 1,ntypd
     120           32 :          do nat = 1,neq(ntyp)
     121           16 :             iatom = iatom + 1
     122           24 :             ntp(iatom) = ntyp
     123              :          enddo
     124              :       enddo 
     125            8 :       call timestop("gen ntp")
     126              : 
     127            8 :       call timestart("inquire pro2")
     128            8 :       if(l_amn2)then
     129            0 :        do j=jspin,1,-1
     130            0 :          inquire(file=trim('pro2'//spin12(j)),exist=l_file)
     131            0 :          if(l_file)then
     132            0 :             filename='pro2'//spin12(j)
     133            0 :             exit
     134              :          endif
     135              :        enddo
     136              :       else   
     137              : c..reading the proj.1 / proj.2 / proj file
     138           16 :        do j=jspin,0,-1
     139           16 :          inquire(file=trim('proj'//spin12(j)),exist=l_file)
     140           16 :          if(l_file)then
     141            8 :             filename='proj'//spin12(j)
     142            8 :             exit
     143              :          endif
     144              :        enddo
     145              :       endif
     146            8 :       call timestop("inquire pro2")
     147              : 
     148            8 :       if(l_file)then
     149            8 :         open (203,file=trim(filename),status='old')
     150            8 :         rewind (203)
     151              :       else
     152            0 :          CALL juDFT_error("no proj/proj.1/proj.2",calledby="wann_amn")
     153              :       endif  
     154              : 
     155            8 :       if(l_amn2)then
     156            0 :        allocate(posshifts(3,nwfsd))
     157            0 :        allocate(weights(nwfsd))
     158              :       endif    
     159              : 
     160            8 :       call timestart("read proj 203")
     161            8 :       if(l_nocosoc)then
     162            0 :        read (203,*)nwfs,banddummy,l_oldproj
     163            0 :        if(.not.l_oldproj)then
     164            0 :         do nwf=1,nwfs
     165            0 :          read(203,*)ind(nwf),lwf(nwf),jwf(nwf),jmwf(nwf),rwf(nwf)
     166            0 :          read(203,*)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              :        else
     170            0 :         do nwf = 1,nwfs
     171              :          read (203,*) 
     172            0 :      &            ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf),spi(nwf)
     173              :          read (203,*) 
     174            0 :      &            alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
     175            0 :          if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
     176              :         enddo !nwf
     177              :        endif !oldproj
     178              :       else
     179            8 :        read (203,*) nwfs
     180           72 :        do nwf = 1,nwfs
     181              :          read (203,*) 
     182           64 :      &            ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf)
     183              :          read (203,*) 
     184           64 :      &            alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
     185           72 :          if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
     186              :        enddo !nwf
     187              :       endif !l_nocosoc
     188            8 :       rewind (203)
     189            8 :       close (203)
     190            8 :       call timestop("read proj 203")
     191              : 
     192            8 :       call timestart("output trail WFs")
     193            8 :       if (ikpt.eq.1) then
     194            1 :       write (oUnit,*) 'Number of trial WFs:',nwfs
     195            1 :       write (oUnit,*)
     196            9 :       do nwf = 1,nwfs
     197            8 :         write (oUnit,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
     198            8 :         write (oUnit,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
     199            8 :         write (oUnit,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
     200            8 :         write (oUnit,*) 'alpha=',alpha(nwf),
     201           16 :      &          ' beta=',beta(nwf),' gamma=',gamma(nwf)
     202            9 :         write (oUnit,*)
     203              :       enddo 
     204              :       endif
     205            8 :       call timestop("output trail WFs")
     206              : 
     207              : c..generating the radial twf function
     208              : 
     209       278856 :       rads(:,:,:,:) = 0.
     210              : 
     211              :       call wann_rad_twf(
     212              :      >         nwfs,jmtd,natd,ind,rwf,zona,regio,
     213              :      >         us,dus,uds,duds,ff,gg,lmaxd,ikpt,
     214              :      >         ntypd,ntp,jri,rmsh,dx,
     215              :      >         nlod,flo,llo,nlo, 
     216            8 :      <         rads)
     217              : 
     218            8 :       call timestart("write rads")
     219            8 :       open (100,file='rads')
     220         3776 :       do i = 1,jmtd
     221              : c       write (100,'(i3,2x,4f10.6)') i,rads(1,0:3,i,1)
     222         3776 :         write (100,'(f10.6,2x,4f10.6)') rmsh(i,1),rads(1,0:3,i,1)
     223              :       enddo
     224            8 :       close(100)
     225            8 :       call timestop("write rads")
     226              : 
     227              : c..now generate the coefficients in the expansion in lattice 
     228              : c..harmonics of the angular part of the twf
     229            8 :       call timestart("generate coeffs in latham")
     230         2312 :       tlmwft(:,:,:) = cmplx(0.,0.)
     231         2312 :       tlmwf(:,:,:)  = cmplx(0.,0.)
     232              :       
     233            8 :       if(l_nocosoc.and..not.l_oldproj)then
     234            0 :          call soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
     235              :       else
     236              :          call wann_tlmw(
     237              :      >          nwfs,lwf,mrwf,
     238              :      >          l_nocosoc,spi,jspin,
     239            8 :      <          tlmwf)  
     240              :       endif
     241              : 
     242            8 :       call eulerrot(nwfs,alpha,beta,gamma,amx)
     243              : 
     244            8 :       imx(:,:) = 0. ; imx(1,1) = 1. ; imx(2,2) = 1. ; imx(3,3) = 1.
     245              : 
     246              : c..performing the wigner rotation
     247              : c..These rotations are specified in the proj file in terms of the 
     248              : c..euler angles. The wave functions inside the muffin-tins are represen-
     249              : c..-ted in terms of the global-frame-lattice-harmonics, therefore, 
     250              : c..the wigner transformation has to be applied to tlmwf, resulting
     251              : c..in the tlmwft matrix, which is suitable for further use  
     252              : 
     253              : c..given the euler angles, the following procedure has to be 
     254              : c..performed in order to match the two local frames:
     255              : c..1. Rotation around the z-axis by alpha
     256              : c..2. Rotation around the x-axis by beta
     257              : c..3. Rotation around the z-axis by gamma again.
     258            8 :       call d_wigner(nwfs,amx,imx,3,wign)
     259              : 
     260              : c..now we transform the tlmwf coefficients
     261              : 
     262           72 :       do nwf = 1,nwfs
     263          512 :          tlmwft(0,:,nwf) = tlmwf(0,:,nwf)
     264          264 :          do l = 1,3
     265         1216 :             do m = -l,l
     266         6464 :                do mp = -l,l
     267              :                   tlmwft(l,m,nwf) = tlmwft(l,m,nwf) + 
     268         6272 :      +                     wign(mp,m,l,nwf)*tlmwf(l,mp,nwf)
     269              :                enddo 
     270              :             enddo
     271              :          enddo         
     272              :       enddo
     273            8 :       call timestop("generate coeffs in latham")
     274              : 
     275              : c..calculating the amn matrix
     276              : 
     277         7544 :       vlpr(:) = 0. ; vlprd(:) = 0.
     278              : 
     279              : c...sum by wfs, each of them is localized at a certain mt
     280            8 :       call timestart("sum by wfs")
     281           72 :       do nwf = 1,nwfs
     282           64 :          if(l_amn2.and.present(bkpt))then
     283            0 :            arg=-bkpt(1)*posshifts(1,nwf)
     284            0 :            arg=arg-bkpt(2)*posshifts(2,nwf)
     285            0 :            arg=arg-bkpt(3)*posshifts(3,nwf)
     286            0 :            arg=2*pi*arg
     287            0 :            factor=cmplx(cos(arg),sin(arg))*weights(nwf)
     288              :          else
     289              :            factor=1.0
     290              :          endif
     291           64 :          factor = factor*chi
     292              : 
     293           64 :          nat = ind(nwf)
     294           64 :          ntyp = ntp(nat)
     295              : c...sum by bands
     296          584 :          do ne = 1,nslibd
     297              : c...sum by l,m
     298         2560 :             do l = 0,min(lmax(ntyp),3)
     299       966656 :                do j = 1,jri(ntyp)
     300              :                   vlpr(j) = ff(ntyp,j,1,l)*rads(nwf,l,j,1)+
     301       964608 :      +                      ff(ntyp,j,2,l)*rads(nwf,l,j,2) 
     302              :                   vlprd(j) = gg(ntyp,j,1,l)*rads(nwf,l,j,1) +
     303       964608 :      +                       gg(ntyp,j,2,l)*rads(nwf,l,j,2)
     304       966656 :                   if (rwf(nwf).gt.0)then
     305            0 :                      vlpr(j) = vlpr(j)*rmsh(j,ntyp)
     306            0 :                      vlprd(j) = vlprd(j)*rmsh(j,ntyp)
     307              :                   endif
     308              : 
     309              :                enddo
     310              : 
     311              : c..these integrations are not necessary if rads is the lin.comb. 
     312              : c..of the u_l and \dot{u}_l, but are necessary for other ways of
     313              : c..constructing the radial part, therefore, we do it anyway
     314              :  
     315         2048 :                r0 = rmsh(1,ntyp)
     316         2048 :                call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vl)
     317         2048 :                call intgr3(vlprd,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vld)
     318        10752 :                do m = -l,l
     319         8192 :                   lm = l*(l+1) + m
     320              :                   amn(ne,nwf) = amn(ne,nwf) +
     321              :      +               tlmwft(l,m,nwf)*conjg(( acof(ne,lm,nat)*vl + 
     322        10240 :      +               bcof(ne,lm,nat)*vld )*(ImagUnit)**l)*factor
     323              :      
     324              :                enddo
     325              :             enddo
     326              : 
     327              : c..local orbitals
     328          576 :             if (nlo(ntyp).ge.1) then
     329            0 :                do lo = 1,nlo(ntyp)
     330            0 :                   l = llo(lo,ntyp)
     331            0 :                   do j = 1,jri(ntyp)
     332              :                      vlpr(j) = flo(ntyp,j,1,lo)*rads(nwf,l,j,1) +
     333            0 :      +                         flo(ntyp,j,2,lo)*rads(nwf,l,j,2)
     334            0 :                      if (rwf(nwf).gt.0)then
     335            0 :                      vlpr(j) = vlpr(j)*rmsh(j,ntyp)
     336              :                      endif
     337              :                   enddo
     338              : 
     339              : 
     340            0 :                   r0 = rmsh(1,ntyp)
     341            0 :                   call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vlo)
     342            0 :                   do m = -l,l
     343              :                      amn(ne,nwf) = amn(ne,nwf) +
     344              :      +          tlmwft(l,m,nwf)*conjg((ccof(m,ne,lo,nat)*vlo)*
     345            0 :      *          (ImagUnit)**l)*factor
     346              :                   enddo
     347              :                enddo
     348              :             endif 
     349              : 
     350              :          enddo
     351              :       enddo
     352            8 :       call timestop("sum by wfs")
     353            8 :       call timestop("wann_amn")
     354            8 :       return
     355              : 
     356            8 :       end subroutine wann_amn
     357              : 
     358              : c*********************************************************************
     359              : c     Calculate the expansion of the angular part of the trial orbital
     360              : c     in terms of spherical harmonics.
     361              : c     Version for Soc: The trial orbital is a spinor.
     362              : c     Frank Freimuth, June 2007
     363              : c*********************************************************************
     364            0 :       subroutine soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
     365              :       use m_clebsch
     366              :       implicit none
     367              : 
     368              :       integer, intent (in)  :: nwfs
     369              :       integer, intent (in)  :: lwf(nwfs)
     370              :       real,    intent (in)  :: jwf(nwfs),jmwf(nwfs)
     371              :       integer, intent (in)  :: jspin
     372              :       complex, intent (out) :: tlmwf(0:3,-3:3,nwfs)
     373              : 
     374              :       integer nwf,l,m
     375              :       complex tlm(0:3,-3:3,1:7)
     376              :       real j,jm
     377              : 
     378            0 :       call timestart("soc_tlmw")
     379            0 :       do nwf=1,nwfs
     380            0 :          l=lwf(nwf)
     381            0 :          IF(l<0) CALL juDFT_error("not yet implemented",calledby
     382            0 :      +        ="wann_amn")
     383            0 :          j=jwf(nwf)
     384            0 :          jm=jmwf(nwf)
     385            0 :          IF(j<0) CALL juDFT_error("jwf",calledby ="wann_amn")
     386            0 :          IF( ABS(jm) - j  >1e-10) CALL juDFT_error("jmwf",calledby
     387            0 :      +        ="wann_amn")
     388            0 :          if( abs( l + 0.5 -j ).gt. 1e-10 .and.
     389              :      &       abs( l - 0.5 -j ).gt. 1e-10 )
     390              :      &     CALL juDFT_error ('regula trianguli violata',
     391            0 :      &                     calledby="wann_amn")
     392            0 :          tlmwf(0:3,-3:3,nwf)=cmplx(0.0,0.0)
     393            0 :          do m=-l,l
     394            0 :             tlmwf(l,m,nwf)=clebsch(real(l),0.5,real(m),1.5-jspin,j,jm)
     395              :          enddo
     396              :       enddo
     397              : 
     398            0 :       call timestop("soc_tlmw")
     399            0 :       end subroutine soc_tlmw
     400              : 
     401              :       end module m_wann_amn
        

Generated by: LCOV version 2.0-1