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

Generated by: LCOV version 1.13