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

Generated by: LCOV version 1.13