LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_tlmplm2.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 271 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : c*******************************************c
       2             : c   Muffin tin contribution to uHu-dmi      c
       3             : c  <u_{k+b1} | H_{k}^{mt} | u_{\theta+b2}>  c
       4             : c*******************************************c
       5             : c   Routine to set up T(lm,lmp) for all     c
       6             : c   pairs (b1,b2) and every atom n          c
       7             : c                                           c
       8             : c   Includes spherical and non-sph.         c
       9             : c   contributions to Hamiltonian:           c
      10             : c   H^{mt} ~ H^{sph} + H^{non-sph}          c
      11             : c                                           c
      12             : c   possible SOC contribution is dealt      c
      13             : c   with in a separate routine              c
      14             : c*******************************************c
      15             : c               J.-P. Hanke, Feb. 2016      c
      16             : c*******************************************c
      17             : c  TODO: local orbital part: epar(jspin)    c
      18             : c*******************************************c
      19             :       MODULE m_wann_uHu_tlmplm2
      20             :       CONTAINS
      21           0 :       SUBROUTINE wann_uHu_tlmplm2(
      22             :      >                  memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd,
      23           0 :      >                  ntype,dx,rmsh,jri,lmax,ntypsy,natd,
      24           0 :      >                  lnonsph,lmd,lmplmd,clnu,mlh,nmem,llh,nlh,neq,
      25             :      >                  irank,mlotot,mlolotot,
      26           0 :      >                  vr,nlod,llod,loplod,ello,llo,nlo,lo1l,l_dulo,
      27           0 :      >                  ulo_der,f,g,flo,f_b,g_b,flo_b,
      28           0 :      >                  kdiff,kdiff2,nntot,nntot2,bmat,bbmat,vr0,epar,
      29           0 :      >                  invsat,jspin,jspin_b,
      30             :      >                  l_skip_sph,l_skip_non,l_skip_loc,
      31           0 :      <                  tuu,tud,tdu,tdd,
      32           0 :      >                  tuulo,tulou,tdulo,tulod,tuloulo)
      33             : 
      34             :       USE m_intgr, ONLY : intgr3
      35             :       USE m_radflo
      36             :       USE m_radfun
      37             :       USE m_tlo
      38             :       USE m_sphbes
      39             :       USE m_ylm
      40             :       USE m_gaunt, ONLY: gaunt1
      41             :       USE m_dujdr
      42             :       USE m_wann_uHu_radintsra5
      43             :       USE m_wann_uHu_tlo
      44             :       USE m_constants
      45             : 
      46             :       IMPLICIT NONE
      47             : C     ..
      48             : C     .. Scalar Arguments ..
      49             :       INTEGER, INTENT (IN) :: memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd
      50             :       INTEGER, INTENT (IN) :: lmd,lmplmd,ntype,irank
      51             :       INTEGER, INTENT (IN) :: nlod,llod,loplod,natd,mlotot,mlolotot
      52             :       INTEGER, INTENT (IN) :: nntot,nntot2,jspin,jspin_b
      53             :       LOGICAL, INTENT (IN) :: l_skip_sph,l_skip_non,l_skip_loc
      54             : C     ..
      55             : C     .. Array Arguments ..
      56             :       COMPLEX, INTENT (OUT):: 
      57             :      >       tdd(0:lmd,0:lmd,ntypd,nntot*nntot2),
      58             :      >       tdu(0:lmd,0:lmd,ntypd,nntot*nntot2),
      59             :      >       tud(0:lmd,0:lmd,ntypd,nntot*nntot2),
      60             :      >       tuu(0:lmd,0:lmd,ntypd,nntot*nntot2)
      61             :       COMPLEX, INTENT (OUT):: 
      62             :      >       tdulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
      63             :      >       tuulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
      64             :      >       tulou(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
      65             :      >       tulod(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot2),
      66             :      >       tuloulo(nlod,-llod:llod,nlod,-llod:llod,ntypd,nntot*nntot2)
      67             :       COMPLEX, INTENT (IN) :: clnu(memd,0:nlhd,ntypsd)
      68             :       INTEGER, INTENT (IN) :: llo(nlod,ntypd),nlo(ntypd)
      69             :       INTEGER, INTENT (IN) :: lo1l(0:llod,ntypd),neq(ntypd)
      70             :       INTEGER, INTENT (IN) :: mlh(memd,0:nlhd,ntypsd),nlh(ntypsd)
      71             :       INTEGER, INTENT (IN) :: nmem(0:nlhd,ntypsd),llh(0:nlhd,ntypsd)
      72             :       INTEGER, INTENT (IN) :: jri(ntypd),lmax(ntypd),ntypsy(natd)
      73             :       INTEGER, INTENT (IN) :: lnonsph(ntypd),ulo_der(nlod,ntypd)
      74             :       INTEGER, INTENT (IN) :: invsat(natd)
      75             :       REAL,    INTENT (IN) :: dx(ntypd),rmsh(jmtd,ntypd)
      76             :       REAL,    INTENT (IN) :: vr(jmtd,0:nlhd,ntypd)
      77             :       REAL,    INTENT (IN) :: ello(nlod,ntypd)
      78             :       REAL,    INTENT (IN) :: epar(0:lmaxd,ntypd,2)
      79             :       REAL,    INTENT (IN) :: vr0(jmtd,ntypd)
      80             :       REAL,    INTENT (IN) :: f  (ntypd,jmtd,2,0:lmaxd),
      81             :      >                        g  (ntypd,jmtd,2,0:lmaxd),
      82             :      >                        f_b(ntypd,jmtd,2,0:lmaxd),
      83             :      >                        g_b(ntypd,jmtd,2,0:lmaxd),
      84             :      >                        flo  (ntypd,jmtd,2,nlod),
      85             :      >                        flo_b(ntypd,jmtd,2,nlod)
      86             :       REAL,    INTENT (IN) :: bmat(3,3),bbmat(3,3)
      87             :       REAL,    INTENT (IN) :: kdiff (3,nntot)
      88             :       REAL,    INTENT (IN) :: kdiff2(3,nntot2)
      89             :       LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd)
      90             : C     ..
      91             : C     .. Local Scalars ..
      92             :       COMPLEX cil,cci,cil1,fac1,fac2,fac3,cilp,cciljj1,ciljj2
      93             :       INTEGER i,ii,l,l2,lamda,lh,lm,lmin,lmp,lp
      94             :       INTEGER lp1,lpl,m,mem,mems,mp,mu,n,nh,noded,nodeu,nrec,nsym,na,ne
      95             :       INTEGER ljj1,ljj2,mjj1,mjj2,lo,lop
      96             :       INTEGER ltil1,ltil2,mtil1,mtil2
      97             :       INTEGER mlo, mlolo, iu
      98             :       INTEGER :: lmini,lmini2,lmini3
      99             :       INTEGER :: lmaxi,lmaxi2,lmaxi3
     100             :       INTEGER :: lmx,lmx2,lmx3
     101             :       INTEGER :: lzb,ll,llp,lljj1,lljj2,lmjj1,lmjj2
     102             :       INTEGER :: ikpt_b,ikpt_b2
     103             :       INTEGER :: lwn,indexx
     104             :       REAL gs,rk1,rk2,sign
     105             :       REAL gc1,gc2,gc3,invsfct,e,temp_val1,temp_val2
     106             :       REAL t1,t0,t_sphbes,t_ylm,t_radint,t_gggint,t_sphint,t_lo
     107             :       REAL tt1,tt2,tt3
     108             :       LOGICAL :: l_nns
     109             : C     ..
     110             : C     .. Local Arrays ..
     111           0 :       COMPLEX, ALLOCATABLE :: yl1(:),yl2(:)
     112           0 :       REAL, ALLOCATABLE :: dvd_non(:,:,:,:,:),dvu_non(:,:,:,:,:)
     113           0 :       REAL, ALLOCATABLE :: uvd_non(:,:,:,:,:),uvu_non(:,:,:,:,:)
     114           0 :       REAL, ALLOCATABLE :: dvd_sph(:,:,:,:,:),dvu_sph(:,:,:,:,:)
     115           0 :       REAL, ALLOCATABLE :: uvd_sph(:,:,:,:,:),uvu_sph(:,:,:,:,:)
     116           0 :       REAL, ALLOCATABLE :: p(:,:,:,:),p_b(:,:,:,:),dp_b(:,:,:,:)
     117           0 :       REAL, ALLOCATABLE :: q(:,:,:,:),q_b(:,:,:,:),dq_b(:,:,:,:)
     118           0 :       REAL, ALLOCATABLE :: dp(:,:,:,:),dq(:,:,:,:)
     119           0 :       REAL, ALLOCATABLE :: x(:),vv(:),vm(:)
     120           0 :       REAL, ALLOCATABLE :: jlpp(:),jj1(:,:),jj2(:,:)
     121             :       REAL :: bpt(3),bpt2(3),bkrot(3)
     122             : 
     123             : #if (defined(CPP_MPI) && !defined(CPP_T90))
     124             :       INCLUDE 'mpif.h'
     125             :       INTEGER ierr(3)
     126             : #endif
     127             : C     .. Intrinsic Functions ..
     128             :       INTRINSIC abs,cmplx,max,mod
     129             : C     ..
     130             : 
     131           0 :       t_sphbes = 0.
     132           0 :       t_ylm    = 0.
     133           0 :       t_radint = 0.
     134           0 :       t_sphint = 0.
     135           0 :       t_gggint = 0.
     136           0 :       t_lo = 0.
     137           0 :       tt1 = 0.
     138           0 :       tt2 = 0.
     139           0 :       tt3 = 0.
     140             : 
     141           0 :       cci = cmplx(0.,-1.)
     142           0 :       tulou = cmplx(0.,0.); tulod = cmplx(0.,0.)
     143           0 :       tdulo = cmplx(0.,0.); tuulo = cmplx(0.,0.); tuloulo = cmplx(0.,0.)
     144           0 :       tuu = cmplx(0.,0.); tdu = cmplx(0.,0.)
     145           0 :       tud = cmplx(0.,0.); tdd = cmplx(0.,0.)
     146             : 
     147           0 :       allocate( yl1((lmaxd+1)**2), yl2((lmaxd+1)**2) )
     148           0 :       allocate( dvd_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     149           0 :       allocate( dvu_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     150           0 :       allocate( uvd_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     151           0 :       allocate( uvu_non(1:nlhd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     152           0 :       allocate( dvd_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     153           0 :       allocate( dvu_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     154           0 :       allocate( uvd_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     155           0 :       allocate( uvu_sph(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd) )
     156           0 :       allocate( p(jmtd,2,0:lmaxd,0:lmaxd) )
     157           0 :       allocate( p_b(jmtd,2,0:lmaxd,0:lmaxd) )
     158           0 :       allocate( q(jmtd,2,0:lmaxd,0:lmaxd) )
     159           0 :       allocate( q_b(jmtd,2,0:lmaxd,0:lmaxd) )
     160           0 :       allocate( dp(jmtd,2,0:lmaxd,0:lmaxd) )
     161           0 :       allocate( dp_b(jmtd,2,0:lmaxd,0:lmaxd) )
     162           0 :       allocate( dq(jmtd,2,0:lmaxd,0:lmaxd) )
     163           0 :       allocate( dq_b(jmtd,2,0:lmaxd,0:lmaxd) )
     164           0 :       allocate( x(jmtd),vv(jmtd) )
     165           0 :       allocate( vm(jmtd) )
     166           0 :       allocate( jlpp(0:lmaxd), jj1(0:lmaxd,jmtd), jj2(0:lmaxd,jmtd) )
     167             : 
     168             : c***************************************c
     169             : c  begin 1st neighbor loop:   <k+b1|    c
     170             : c***************************************c
     171           0 :       do ikpt_b=1,nntot
     172           0 :        bpt = kdiff(:,ikpt_b)
     173           0 :        rk1 = sqrt(dot_product(bpt,matmul(bbmat,bpt)))
     174             : !       rk1= sqrt(dotirp(bpt,bpt,bbmat))
     175             :        !if(irank.eq.0) write(*,*)'rk1',rk1
     176             : 
     177             : c***************************************c
     178             : c  begin 2nd neighbor loop:   |k+b2>    c
     179             : c***************************************c
     180           0 :        do ikpt_b2=1,nntot2
     181           0 :         indexx = ikpt_b2 + (ikpt_b-1)*nntot2
     182           0 :         bpt2 = kdiff2(:,ikpt_b2)
     183           0 :        rk2 = sqrt(dot_product(bpt2,matmul(bbmat,bpt2)))
     184             : !        rk2= sqrt(dotirp(bpt2,bpt2,bbmat))
     185           0 :         if(rk2.ne.0.0) stop 'dmi but rk2.ne.0'
     186             : 
     187           0 :         na = 1
     188           0 :         mlo = 1 ; mlolo = 1
     189             : 
     190             :       ! loop over atoms
     191           0 :       DO 210 n = 1,ntype
     192           0 :        lwn = lmax(n)
     193           0 :        nsym = ntypsy(na)
     194           0 :        nh = nlh(nsym)
     195             : 
     196             :        !if(indexx.eq.1 .and. irank.eq.0) then
     197             :        ! write(*,*)'atom n=',n
     198             :        ! write(*,*)'spins=',jspin,jspin_b
     199             :        ! do l=0,lwn
     200             :        !  write(*,*)'l=',l,epar(l,n,jspin),epar(l,n,jspin_b)
     201             :        ! enddo
     202             :        !endif
     203             : 
     204             :        l_nns = l_skip_non .OR. (lnonsph(n).LT.0) .OR.
     205           0 :      >         ( (invsat(na).ne.0).and.(invsat(na).ne.1) )
     206           0 :        if(invsat(na).eq.0) invsfct = 1.0
     207           0 :        if(invsat(na).eq.1) invsfct = 2.0
     208           0 :        if(indexx.eq.1 .and. irank.eq.0) write(*,*)'invsfct=',invsfct
     209             : 
     210             :        ! set up spherical Bessel functions
     211             :        ! jj1(l,b1*r) and jj2(l,b2*r)
     212           0 :        call cpu_time(t0)
     213           0 :        do i=1,jri(n)
     214           0 :           gs = rk1*rmsh(i,n)
     215           0 :           call sphbes(lwn,gs,jlpp)
     216           0 :           jj1(:,i) = jlpp(:)
     217             : 
     218           0 :           gs = rk2*rmsh(i,n)
     219           0 :           call sphbes(lwn,gs,jlpp)
     220           0 :           jj2(:,i) = jlpp(:)
     221             : 
     222             : c          if(rk1.eq.0.0) then
     223             : c           if(jj1(0,i).ne.1.0) write(*,*)'jj1.ne.1',i,rmsh(i,n),jj1(0,i)
     224             : c           if(ANY(jj1(1:lwn,i).ne.0.0)) write(*,*)'jj1.ne.0',i,rmsh(i,n)
     225             : c          endif
     226           0 :           if(jj2(0,i).ne.1.0) stop 'dmi but jj2.ne.1'
     227           0 :           if(ANY(jj2(1:lwn,i).ne.0.0)) stop 'dmi but jj2.ne.0'
     228             :          
     229             :        enddo
     230           0 :        call cpu_time(t1)
     231           0 :        t_sphbes=t_sphbes + t1-t0
     232             : 
     233             :        ! set up spherical harmonics
     234             :        ! yl1(lm) for b1 and yl2(lm) for b2
     235           0 :        yl1 = cmplx(0.,0.)
     236           0 :        bkrot=MATMUL(bpt,bmat)
     237             : !       call cotra3(bpt,bkrot,bmat)
     238           0 :        call ylm4(lwn,bkrot,yl1)
     239             : 
     240           0 :        yl2 = cmplx(0.,0.)
     241           0 :        bkrot=MATMUL(bpt2,bmat)
     242             : !       call cotra3(bpt2,bkrot,bmat)
     243           0 :        if(ANY(bkrot.ne.0.0)) stop 'dmi but bkrot.ne.0'
     244           0 :        call ylm4(lwn,bkrot,yl2)
     245             :        !if(indexx.eq.1 .and. irank.eq.0) write(*,*)'yl2',yl2(1)
     246             : 
     247           0 :        call cpu_time(t0)
     248           0 :        t_ylm=t_ylm + t0-t1
     249             : 
     250             :        ! set up products of radial functions and sph. Bessel
     251           0 :        DO ljj1 = 0, lwn
     252           0 :         DO l = 0, lwn
     253           0 :          DO i=1,jri(n)
     254           0 :           p(i,:,l,ljj1) = f(n,i,:,l)*jj1(ljj1,i)
     255           0 :           q(i,:,l,ljj1) = g(n,i,:,l)*jj1(ljj1,i)
     256           0 :           p_b(i,:,l,ljj1) = f_b(n,i,:,l)*jj2(ljj1,i)
     257           0 :           q_b(i,:,l,ljj1) = g_b(n,i,:,l)*jj2(ljj1,i)
     258             :          ENDDO
     259             :          CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),f(n,:,:,l),
     260           0 :      >              jj1,rk1,ljj1,lmaxd,dp(:,:,l,ljj1))
     261             :          CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),g(n,:,:,l),
     262           0 :      >              jj1,rk1,ljj1,lmaxd,dq(:,:,l,ljj1))
     263             :          CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),f_b(n,:,:,l),
     264           0 :      >              jj2,rk2,ljj1,lmaxd,dp_b(:,:,l,ljj1))
     265             :          CALL dujdr(jmtd,jri(n),rmsh(1,n),dx(n),g_b(n,:,:,l),
     266           0 :      >              jj2,rk2,ljj1,lmaxd,dq_b(:,:,l,ljj1))
     267             :         ENDDO
     268             :        ENDDO
     269             : 
     270             :        if(.false. .and. irank.eq.0 .and. n.eq.1) then
     271             :         open(876,file='out_radff')
     272             :         open(877,file='out_radgg')
     273             :         open(878,file='out_raddf')
     274             :         open(879,file='out_raddg')
     275             :         open(880,file='out_radjj')
     276             :         do i=1,jri(n)
     277             :          do l=0,lwn
     278             :           write(876,*)rmsh(i,n),f(n,i,1,l),f(n,i,2,l)
     279             :           write(877,*)rmsh(i,n),g(n,i,1,l),g(n,i,2,l)
     280             :           write(878,*)rmsh(i,n),dp(i,1,l,1),dp(i,2,l,1)
     281             :           write(879,*)rmsh(i,n),dq(i,1,l,1),dq(i,2,l,1)
     282             :           write(880,*)rmsh(i,n),jj1(l,i),jj2(l,i)
     283             :          enddo
     284             :         enddo
     285             :         close(876)
     286             :         close(877)
     287             :         close(878)
     288             :         close(879)
     289             :         close(880)
     290             :         stop
     291             :        endif
     292             : 
     293             : c****************************************************c
     294             : c   compute radial integrals                         c
     295             : c   <u(l')jj1(ljj1) | v(lamda,nu) | u(l)jj2(ljj2)>   c
     296             : c   where v(lamda,nu) ~ V^{sph} + V^{non-sph}        c
     297             : c****************************************************c
     298             :        ! spherical part
     299           0 :        IF(.not.l_skip_sph) THEN
     300           0 :        uvu_sph = 0.
     301           0 :        uvd_sph = 0.
     302           0 :        dvu_sph = 0.
     303           0 :        dvd_sph = 0.
     304           0 :        DO ljj1 = 0, lwn
     305           0 :         DO lp = 0, lwn
     306           0 :          lmini = abs(lp-ljj1)
     307           0 :          lmaxi = lp+ljj1
     308           0 :          lmx = min(lmaxi,lwn)
     309           0 :          DO lh=lmini,lmx
     310           0 :           if(mod(lmaxi+lh,2).eq.1) cycle
     311           0 :           DO ljj2 = 0,0!lwn
     312           0 :              lmini2= abs(lh-ljj2)
     313           0 :              lmaxi2= lh+ljj2
     314           0 :              lmx2 = min(lmaxi2,lwn)
     315           0 :              DO l = lmini2, lmx2
     316           0 :                 if(mod(lmaxi2+l,2).eq.1) cycle
     317           0 :                 e=(epar(lp,n,jspin)+epar(l,n,jspin_b))/2.0!epar(lh,n,jspin)
     318           0 :                 lzb = lh
     319             : 
     320             :                 !if(irank.eq.0) write(*,*)lh,l,lp,ljj2,ljj1
     321             : 
     322             :                 !if(irank.eq.0) write(*,*)'uvu',epar(lp,n,jspin)
     323             :                 call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
     324             :      >            e,vr0(1,n),
     325             :      >            p(:,:,lp,ljj1),p_b(:,:,l,ljj2),
     326             :      >            dp(:,:,lp,ljj1),dp_b(:,:,l,ljj2),
     327           0 :      >            lmaxd,lzb,uvu_sph(lh,l,lp,ljj2,ljj1),irank)
     328             : 
     329             :                 !if(irank.eq.0) write(*,*)'uvd',1.0
     330             :                 call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
     331             :      >            e,vr0(1,n),
     332             :      >            p(:,:,lp,ljj1),q_b(:,:,l,ljj2),
     333             :      >            dp(:,:,lp,ljj1),dq_b(:,:,l,ljj2),
     334           0 :      >            lmaxd,lzb,uvd_sph(lh,l,lp,ljj2,ljj1),irank)
     335             : 
     336             :                 !if(irank.eq.0) write(*,*)'dvu',0.0
     337             :                 call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
     338             :      >            e,vr0(1,n),
     339             :      >            q(:,:,lp,ljj1),p_b(:,:,l,ljj2),
     340             :      >            dq(:,:,lp,ljj1),dp_b(:,:,l,ljj2),
     341           0 :      >            lmaxd,lzb,dvu_sph(lh,l,lp,ljj2,ljj1),irank)
     342             : 
     343             : c                ! TODO: check influence of symmetrization on OM and DMI
     344             : c                uvd_sph(lh,l,lp,ljj2,ljj1) = (uvd_sph(lh,l,lp,ljj2,ljj1)
     345             : c     >                                  +dvu_sph(lh,l,lp,ljj2,ljj1))/2.0
     346             : c                dvu_sph(lh,l,lp,ljj2,ljj1) = uvd_sph(lh,l,lp,ljj2,ljj1)
     347             : 
     348             :                 !if(irank.eq.0) write(*,*)'dvd',epar(lp,n,jspin)
     349             :                 call wann_uHu_radintsra5(jmtd,jri(n),rmsh(1,n),dx(n),
     350             :      >            e,vr0(1,n),
     351             :      >            q(:,:,lp,ljj1),q_b(:,:,l,ljj2),
     352             :      >            dq(:,:,lp,ljj1),dq_b(:,:,l,ljj2),
     353           0 :      >            lmaxd,lzb,dvd_sph(lh,l,lp,ljj2,ljj1),irank)
     354             :              ENDDO
     355             :           ENDDO
     356             :          ENDDO
     357             :         ENDDO
     358             :        ENDDO
     359             :        ENDIF
     360             : 
     361             : 
     362             :        ! nonspherical part
     363           0 :        IF(.not.l_nns) then
     364           0 :        uvu_non = 0.
     365           0 :        uvd_non = 0.
     366           0 :        dvu_non = 0.
     367           0 :        dvd_non = 0.
     368           0 :        DO ljj1 = 0, lwn
     369           0 :         DO ljj2 = 0, 0!lwn
     370           0 :          DO lp = 0, lwn
     371             : c          DO l = 0, lp
     372           0 :           DO l = 0, lwn
     373           0 :                DO lh = 1, nh
     374           0 :                      DO i = 1,jri(n)
     375             :                      x(i)=(p(i,1,lp,ljj1)*p_b(i,1,l,ljj2)
     376             :      >                    +p(i,2,lp,ljj1)*p_b(i,2,l,ljj2))
     377           0 :      >                    * vr(i,lh,n)
     378             :                      ENDDO
     379             :                      CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
     380           0 :      >                           uvu_non(lh,l,lp,ljj2,ljj1))
     381           0 :                      DO i = 1,jri(n)
     382             :                      x(i)=(q(i,1,lp,ljj1)*p_b(i,1,l,ljj2)
     383             :      >                    +q(i,2,lp,ljj1)*p_b(i,2,l,ljj2))
     384           0 :      >                   * vr(i,lh,n)
     385             :                      ENDDO
     386             :                      CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
     387           0 :      >                           dvu_non(lh,l,lp,ljj2,ljj1))
     388           0 :                      DO i = 1,jri(n)
     389             :                      x(i)=(p(i,1,lp,ljj1)*q_b(i,1,l,ljj2)
     390             :      >                    +p(i,2,lp,ljj1)*q_b(i,2,l,ljj2))
     391           0 :      >                   * vr(i,lh,n)
     392             :                      ENDDO
     393             :                      CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
     394           0 :      >                           uvd_non(lh,l,lp,ljj2,ljj1))
     395           0 :                      DO i = 1,jri(n)
     396             :                      x(i)=(q(i,1,lp,ljj1)*q_b(i,1,l,ljj2)
     397             :      >                    +q(i,2,lp,ljj1)*q_b(i,2,l,ljj2))
     398           0 :      >                   * vr(i,lh,n)
     399             :                      ENDDO
     400             :                      CALL intgr3(x(1:jri(n)),rmsh(1,n),dx(n),jri(n),
     401           0 :      >                           dvd_non(lh,l,lp,ljj2,ljj1))
     402             :                ENDDO
     403             :           ENDDO
     404             :          ENDDO
     405             :         ENDDO
     406             :        ENDDO
     407             :        ! symmetry can be used if f=f_b and g=g_b
     408             :        ! i.e. if jspin = jspin_b
     409             : c       DO lp=0,lwn
     410             : c        DO l=lp+1,lwn
     411             : c            uvu_non(:,l,lp,:,:) = uvu_non(:,lp,l,:,:)
     412             : c            dvd_non(:,l,lp,:,:) = dvd_non(:,lp,l,:,:)
     413             : c            uvd_non(:,l,lp,:,:) = dvu_non(:,lp,l,:,:)
     414             : c            dvu_non(:,l,lp,:,:) = uvd_non(:,lp,l,:,:)
     415             : c        ENDDO
     416             : c       ENDDO
     417             :        ENDIF
     418             : 
     419           0 :        call cpu_time(t0)
     420           0 :        t_radint = t_radint + t0-t1 
     421             : 
     422           0 :        tuu(:,:,n,indexx) = cmplx(0.,0.)
     423           0 :        tdu(:,:,n,indexx) = cmplx(0.,0.)
     424           0 :        tud(:,:,n,indexx) = cmplx(0.,0.)
     425           0 :        tdd(:,:,n,indexx) = cmplx(0.,0.)
     426             : 
     427           0 :        if(l_skip_sph) GOTO 444
     428             : c************** SPHERICAL CONTRIBUTION *******************c
     429             : c   compute product of the two Gaunt coefficients         c
     430             : c   with the radial integrals (+prefactors)               c
     431             : c*********************************************************c
     432             : c   We deal with two Gaunt coefficients:                  c
     433             : c  G1 = G( (ltil1,mtil1), (ljj1,mjj1) , ( lp, mp) )       c
     434             : c  G2 = G( (l , m)      , (ljj2, mjj2), (ltil1, mtil1) )  c
     435             : c*********************************************************c
     436             : c   use Gaunt conditions to reduce number of operations.  c
     437             : c   coefficient G(L1,L2,L3) only nonzero if               c
     438             : c       a)  l1 + l2 + l3 = even                           c
     439             : c       b)  |l2-l3| <= l1 <= l2+l3                        c
     440             : c       c)  m1 = m2 + m3                                  c
     441             : c*********************************************************c
     442           0 :        DO ljj1=0,lwn
     443           0 :         lljj1 = ljj1 * (ljj1 + 1)
     444           0 :         cciljj1 = cci**ljj1
     445           0 :         DO mjj1=-ljj1,ljj1
     446           0 :          lmjj1 = lljj1 + mjj1
     447           0 :          fac1 = cciljj1*conjg(yl1(lmjj1 + 1))
     448           0 :          if(fac1.eq.0.0) cycle
     449             : 
     450           0 :          DO lp=0,lwn
     451           0 :           llp = lp*(lp+1)
     452           0 :           lmini = abs(ljj1-lp)
     453           0 :           lmaxi = ljj1+lp
     454           0 :           lmx = min(lmaxi,lwn)
     455           0 :           DO mp=-lp,lp
     456           0 :            lmp = llp+mp
     457           0 :            mtil1=mjj1+mp
     458           0 :            IF(abs(mtil1).gt.lmx) cycle
     459             : 
     460           0 :            DO ltil1=lmini,lmx
     461             :             ! Gaunt conditions (G1):
     462           0 :             if((mod(lmaxi+ltil1,2).eq.1).or.
     463             :      >         (abs(mtil1).gt.ltil1)) cycle
     464             : 
     465             :              gc1 = gaunt1(ltil1, ljj1, lp, 
     466           0 :      >                    mtil1, mjj1, mp, lmaxd)
     467             : 
     468           0 :              cil1 =  (ImagUnit**lp) * fac1 * gc1
     469             : 
     470           0 :          DO ljj2=0,0!lwn
     471           0 :           lljj2 = ljj2 * (ljj2 + 1)
     472           0 :           lmini2 = abs(ltil1-ljj2)
     473           0 :           lmaxi2 = ltil1+ljj2
     474           0 :           lmx2 = min(lmaxi2,lwn)
     475           0 :           DO mjj2=-ljj2,ljj2
     476           0 :            lmjj2 = lljj2 + mjj2
     477           0 :            m=mjj2+mtil1
     478           0 :            IF(abs(m).gt.lmx2) cycle
     479             : 
     480           0 :            DO l=lmini2,lmx2
     481             :              ! Gaunt conditions (G2):
     482           0 :              if((mod(lmaxi2+l,2).eq.1).or.
     483             :      >          (abs(m).gt.l)) cycle
     484           0 :             ll = l*(l+1)
     485           0 :             lm = ll+m
     486             : 
     487             :              gc2 = gaunt1( l, ljj2, ltil1,
     488           0 :      >                     m, mjj2, mtil1, lmaxd)
     489             : 
     490             :              ! set up prefactor
     491             :              cil =  cil1* ( ImagUnit ** (ljj2 - l) ) 
     492           0 :      +              * gc2 * conjg( yl2(lmjj2 + 1) )
     493             : 
     494           0 :              if(cil.eq.0.0) cycle
     495             : 
     496             :              ! additional factor from symmetrization (below)
     497             : c             cil = 0.5 * cil
     498             : 
     499             :              ! compute T-coefficients
     500             :              ! using symmetrized uvu,uvd,dvu,dvd integrals
     501             :              tuu(lm, lmp, n, indexx)
     502             :      >         = tuu(lm, lmp, n, indexx)
     503           0 :      >         + cil * ( uvu_sph(ltil1, l, lp, ljj2, ljj1) )
     504             : 
     505             :              tdd(lm, lmp, n, indexx)
     506             :      >         = tdd(lm, lmp, n, indexx)
     507           0 :      >         + cil * ( dvd_sph(ltil1, l, lp, ljj2, ljj1) )
     508             : 
     509             :              tud(lm, lmp, n, indexx)
     510             :      >         = tud(lm, lmp, n, indexx)
     511           0 :      >         + cil * ( uvd_sph(ltil1, l, lp, ljj2, ljj1) )
     512             : 
     513             :              tdu(lm, lmp, n, indexx)
     514             :      >         = tdu(lm, lmp, n, indexx)
     515           0 :      >         + cil * ( dvu_sph(ltil1, l, lp, ljj2, ljj1) )
     516             :            ENDDO!ljj2
     517             :           ENDDO!m
     518             :          ENDDO!l
     519             :            ENDDO!ljj1
     520             :           ENDDO!mp
     521             :          ENDDO!lp
     522             :         ENDDO!mtil1
     523             :        ENDDO!ltil1
     524             : 
     525             :   444  CONTINUE
     526           0 :        call cpu_time(t1)
     527           0 :        t_sphint = t_sphint + t1-t0
     528             :  
     529           0 :        IF( l_nns ) GOTO 555
     530             : c************** NON-SPHERICAL CONTRIBUTION ***************c
     531             : c   compute product of the three Gaunt coefficients       c
     532             : c   with the radial integrals (+prefactors)               c
     533             : c*********************************************************c
     534             : c   We deal with three Gaunt coefficients:                c
     535             : c  G1 = G( (ltil1,mtil1) , (lamda,mu)  , (ltil2,mtil2) )  c
     536             : c  G2 = G( (ltil1,mtil1) , (ljj1,mjj1) , (lp,mp)       )  c
     537             : c  G3 = G( (l,m)         , (ljj2,mjj2) , (ltil2,mtil2) )  c
     538             : c*********************************************************c
     539             : c   use Gaunt conditions to reduce number of operations.  c
     540             : c   coefficient G(L1,L2,L3) only nonzero if               c
     541             : c       a)  l1 + l2 + l3 = even                           c
     542             : c       b)  |l2-l3| <= l1 <= l2+l3                        c
     543             : c       c)  m1 = m2 + m3                                  c
     544             : c*********************************************************c
     545           0 :        DO ltil2 = 0, lwn         
     546           0 :         DO mtil2 = -ltil2, ltil2
     547           0 :          DO lh = 1, nh
     548           0 :             lamda = llh(lh,nsym)
     549           0 :             lmini = abs(ltil2 - lamda)
     550           0 :             lmaxi = ltil2 + lamda
     551           0 :             lmx = min(lmaxi,lwn)
     552           0 :             mems = nmem(lh,nsym)
     553           0 :             DO mem = 1, mems
     554           0 :                mu = mlh(mem,lh,nsym)
     555           0 :                mtil1 = mtil2+mu
     556           0 :                IF(abs(mtil1).GT.lmx) CYCLE
     557             : 
     558           0 :                DO ltil1 = lmini, lmx
     559             : 
     560             :                ! Gaunt conditions (G1):
     561           0 :                if((mod(lmaxi+ltil1,2).eq.1).or.
     562             :      >            (abs(mtil1).gt.ltil1)) cycle
     563             : 
     564             :                   gc1 = gaunt1(ltil1, lamda, ltil2, 
     565           0 :      >                         mtil1,    mu, mtil2, lmaxd)
     566           0 :                   fac1 = conjg(clnu(mem,lh,nsym)) * gc1 * invsfct
     567           0 :                   if(fac1.eq.0.0) CYCLE
     568             : 
     569           0 :          DO lp = 0, lnonsph(n)
     570           0 :             llp = lp * (lp + 1)
     571           0 :             cilp = ImagUnit**lp
     572           0 :             DO mp = -lp, lp
     573           0 :                lmp = llp + mp
     574           0 :                mjj1 = mtil1 - mp
     575           0 :                IF(abs(mjj1).GT.lwn) CYCLE
     576             : 
     577           0 :                DO ljj1 = 0, lwn
     578           0 :                   lmini2 = abs(lp-ljj1)
     579           0 :                   lmaxi2 = lp+ljj1
     580             : 
     581             :                   ! Gaunt conditions (G2):
     582             :                   if((lmini2.gt.ltil1).or.(lmaxi2.lt.ltil1).or.
     583           0 :      >               (mod(lmaxi2+ltil1,2).eq.1).or.
     584             :      >               (abs(mjj1).gt.ljj1)) cycle
     585             : 
     586           0 :                   lljj1 = ljj1 * (ljj1 + 1)
     587           0 :                   lmjj1 = lljj1 + mjj1
     588             : 
     589             :                   gc2 = gaunt1(ltil1,  ljj1,  lp,
     590           0 :      >                         mtil1,  mjj1,  mp, lmaxd)
     591             :                   fac2 = fac1 * cilp * (cci**ljj1) * gc2
     592           0 :      >                 * conjg( yl1(lmjj1 + 1) )
     593           0 :                   IF(fac2.eq.0.0) CYCLE
     594             : 
     595           0 :          DO ljj2 = 0, 0!lwn
     596           0 :             lljj2 = ljj2 * (ljj2 + 1)
     597           0 :             ciljj2 = ImagUnit**ljj2
     598           0 :             lmini3 = abs(ltil2 - ljj2)
     599           0 :             lmaxi3 = ltil2 + ljj2
     600           0 :             lmx3 = min(lmaxi3,lnonsph(n))
     601           0 :             DO mjj2 = -ljj2,ljj2
     602           0 :                m = mjj2+mtil2
     603           0 :                IF(abs(m).GT.lmx3) CYCLE
     604           0 :                lmjj2 = lljj2 + mjj2
     605           0 :                fac3 = fac2 * ciljj2 * conjg(yl2(lmjj2 + 1))
     606           0 :                if(fac3.eq.0.0) CYCLE
     607             : 
     608           0 :             DO l = lmini3,lmx3
     609             : 
     610             :                   ! Gaunt conditions (G3):
     611           0 :                   if((mod(lmaxi3+l,2).eq.1).or.
     612             :      >               (abs(m).gt.l)) cycle
     613             : 
     614           0 :                   ll = l * (l + 1)
     615           0 :                   lm = ll + m
     616             : 
     617             :                   gc3 = gaunt1(l,  ljj2,  ltil2,
     618           0 :      >                         m,  mjj2,  mtil2, lmaxd)
     619             : 
     620             :                   ! set up prefactor
     621           0 :                    cil = fac3 * (cci**l) * gc3
     622             : 
     623             :                   ! compute T-coefficients
     624             :                   tuu(lm, lmp, n, indexx)
     625             :      >              = tuu(lm, lmp, n, indexx)
     626           0 :      >              + cil * uvu_non(lh, l, lp, ljj2, ljj1)
     627             : 
     628             :                   tdd(lm, lmp, n, indexx)
     629             :      >              = tdd(lm, lmp, n, indexx)
     630           0 :      >              + cil * dvd_non(lh, l, lp, ljj2, ljj1)
     631             : 
     632             :                   tud(lm, lmp, n, indexx)
     633             :      >              = tud(lm, lmp, n, indexx)
     634           0 :      >              + cil * uvd_non(lh, l, lp, ljj2, ljj1)
     635             : 
     636             :                   tdu(lm, lmp, n, indexx)
     637             :      >              = tdu(lm, lmp, n, indexx)
     638           0 :      >              + cil * dvu_non(lh, l, lp, ljj2, ljj1)
     639             :                ENDDO !l
     640             :             ENDDO !mjj2
     641             :          ENDDO !ljj2
     642             :                ENDDO !lp
     643             :             ENDDO !mjj1
     644             :          ENDDO !mjj1
     645             :                ENDDO !ltil1
     646             :             ENDDO !mem
     647             :          ENDDO !lh
     648             :         ENDDO !mtil2
     649             :        ENDDO !ltil2
     650             : 
     651             :   555  CONTINUE
     652           0 :        call cpu_time(t0)
     653           0 :        t_gggint = t_gggint + t0-t1
     654             : 
     655             : c************** LOCAL ORBITAL CONTRIBUTION ***************c
     656             : c   determine contributions to T-matrix due to loc. orb.  c
     657             : c*********************************************************c
     658           0 :        IF((nlo(n).GE.1).AND.(.NOT.l_skip_loc)) THEN
     659             :         !write(*,*)'TODO: include loc. orb. for DMI'
     660             :           call wann_uHu_tlo(
     661             :      >             memd,nlhd,ntypsd,jmtd,lmaxd,dx(n),rmsh(1,n),jri(n),
     662             :      >             lmax(n),ntypsy,natd,lnonsph(n),lmd,lmplmd,clnu,
     663             :      >             mlh,nmem,llh,nlh,irank,vr(1,0,n),nlod,llod,
     664             :      >             loplod,ello(1,n),
     665             :      >             llo(1,n),nlo(n),lo1l(0,n),l_dulo(1,n),ulo_der(1,n),
     666             :      >             f(n,1:,1:,0:),g(n,1:,1:,0:),flo(n,1:,1:,1:),
     667             :      >             f_b(n,1:,1:,0:),g_b(n,1:,1:,0:),flo_b(n,1:,1:,1:),
     668             :      >             vr0(1,n),epar(0,n,jspin),l_skip_sph,
     669             :      >             l_nns,rk1,rk2,invsfct,yl1,yl2,jj1,jj2,
     670             :      >             p,dp,p_b,dp_b,q,dq,q_b,dq_b,
     671             :      >             tuulo(0,1,-llod,n,indexx),
     672             :      >             tdulo(0,1,-llod,n,indexx),
     673             :      >             tulou(0,1,-llod,n,indexx),
     674             :      >             tulod(0,1,-llod,n,indexx),
     675           0 :      >             tuloulo(1,-llod,1,-llod,n,indexx))
     676           0 :           call cpu_time(t1)
     677           0 :           t_lo = t_lo + t1-t0
     678             :        ENDIF
     679             : 
     680           0 :        na = na + neq(n)
     681           0 :   210 ENDDO !loop over atoms
     682             : 
     683             : c        write(*,*)ikpt_b,ikpt_b2
     684             : c        write(*,*)'t_radint=',t_radint-tt1
     685             : c        write(*,*)'t_sphint=',t_sphint-tt2
     686             : c        write(*,*)'t_gggint=',t_gggint-tt3
     687             : c        tt1 = t_radint
     688             : c        tt2 = t_sphint
     689             : c        tt3 = t_gggint
     690             : 
     691             :        enddo !ikpt_b2
     692             :       enddo !ikpt_b
     693             : 
     694             : #if (defined(CPP_MPI) && !defined(CPP_T90))
     695           0 :       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
     696             : #endif
     697             : 
     698           0 :       deallocate( yl1, yl2 )
     699           0 :       deallocate( dvd_sph, dvu_sph, uvd_sph, uvu_sph ) 
     700           0 :       deallocate( dvd_non, dvu_non, uvd_non, uvu_non ) 
     701           0 :       deallocate( p, dp, p_b, dp_b, q, dq, q_b, dq_b )
     702           0 :       deallocate( x,vv )
     703           0 :       deallocate( vm )
     704           0 :       deallocate( jlpp, jj1, jj2 )
     705             : 
     706             : c      if(irank.eq.0) then
     707             : c       write(*,*)'t_sphbes=',t_sphbes
     708             : c       write(*,*)'t_ylm   =',t_ylm
     709             : c       write(*,*)'t_radint=',t_radint
     710             : c       write(*,*)'t_sphint=',t_sphint
     711             : c       write(*,*)'t_gggint=',t_gggint
     712             : c       write(*,*)'t_locorb=',t_lo
     713             : c      endif
     714             : 
     715           0 :       END SUBROUTINE wann_uHu_tlmplm2
     716             :       END MODULE m_wann_uHu_tlmplm2

Generated by: LCOV version 1.13