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

          Line data    Source code
       1             : c*****************************************c
       2             : c    Spin-orbit coupling cont. to uHu     c
       3             : c  < u_{k+b1} | V_{k}^{soc} | 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                J.-P. Hanke, Dec. 2015   c
      10             : c*****************************************c
      11             :       MODULE m_wann_uHu_soc_tlo
      12             :       CONTAINS
      13           0 :       SUBROUTINE wann_uHu_soc_tlo(
      14             :      >                  input,atoms,
      15             :      >                  ntypd,jmtd,lmaxd,jspd,
      16           0 :      >                  ntype,dx,rmsh,jri,lmax,natd,
      17             :      >                  lmd,irank,nlod,llod,
      18           0 :      >                  loplod,ello,llo,nlo,lo1l,l_dulo,ulo_der,
      19           0 :      >                  flo,flo_b,kdiff,kdiff2,nntot,nntot2,
      20           0 :      >                  vr,epar,jspin,jspin_b,jspins,l_spav,theta,phi,
      21           0 :      >                  yl1,yl2,jj1,jj2,p,p_b,q,q_b,ntyp,angso,
      22             :      >                  l_socscreen,socscreen_fac,
      23           0 :      <                  tuulo,tdulo,tulou,tulod,tuloulo)
      24             : 
      25             :       USE m_intgr, ONLY : intgr0
      26             :       USE m_sphbes
      27             :       USE m_ylm
      28             :       USE m_gaunt, ONLY: gaunt1
      29             :       USE m_sointg
      30             :       USE m_constants
      31             :       USE m_types
      32             :       
      33             : 
      34             :       IMPLICIT NONE
      35             : C     .. Intrinsic Functions ..
      36             :       INTRINSIC abs,cmplx,max,mod
      37             : 
      38             :       TYPE(t_input),INTENT(IN)   :: input
      39             :       TYPE(t_atoms),INTENT(IN)   :: atoms
      40             : 
      41             : C     ..
      42             : C     .. Scalar Arguments ..
      43             :       INTEGER, INTENT (IN) :: ntypd,jmtd,lmaxd,jspd,jspins
      44             :       INTEGER, INTENT (IN) :: lmd,ntype,irank,jspin,jspin_b
      45             :       INTEGER, INTENT (IN) :: nlod,llod,loplod,natd
      46             :       INTEGER, INTENT (IN) :: nntot,nntot2,ntyp
      47             :       REAL,    INTENT (IN) :: theta,phi,socscreen_fac
      48             :       LOGICAL, INTENT (IN) :: l_socscreen
      49             : C     ..
      50             : C     .. Array Arguments ..
      51             :       COMPLEX, INTENT (INOUT)::
      52             :      >        tdulo(0:lmd,nlod,-llod:llod),
      53             :      >        tuulo(0:lmd,nlod,-llod:llod),
      54             :      >        tulou(0:lmd,nlod,-llod:llod),
      55             :      >        tulod(0:lmd,nlod,-llod:llod),
      56             :      >        tuloulo(nlod,-llod:llod,nlod,-llod:llod)
      57             :       COMPLEX, INTENT (IN) :: angso(0:lmaxd,0:lmaxd)
      58             :       COMPLEX, INTENT (IN) :: yl1((lmaxd+1)**2),yl2((lmaxd+1)**2)
      59             :       INTEGER, INTENT (IN) :: llo(nlod,ntypd),nlo(ntypd)
      60             :       INTEGER, INTENT (IN) :: lo1l(0:llod,ntypd)
      61             :       INTEGER, INTENT (IN) :: jri(ntypd),lmax(ntypd)
      62             :       INTEGER, INTENT (IN) :: ulo_der(nlod,ntypd)
      63             :       REAL,    INTENT (IN) :: dx(ntypd),rmsh(jmtd,ntypd)
      64             :       REAL,    INTENT (IN) :: ello(nlod,ntypd,max(jspd,2))
      65             :       REAL,    INTENT (IN) :: epar(0:lmaxd,ntypd,max(jspd,2))
      66             :       REAL,    INTENT (IN) :: vr(jmtd,ntypd,jspd)
      67             :       REAL,    INTENT (IN) :: flo  (ntypd,jmtd,2,nlod),
      68             :      >                        flo_b(ntypd,jmtd,2,nlod)
      69             :       REAL,    INTENT (IN) :: p(jmtd,0:lmaxd,0:lmaxd),
      70             :      >                        q(jmtd,0:lmaxd,0:lmaxd),
      71             :      >                        p_b(jmtd,0:lmaxd,0:lmaxd),
      72             :      >                        q_b(jmtd,0:lmaxd,0:lmaxd)
      73             :       REAL,    INTENT (IN) :: kdiff (3,nntot)
      74             :       REAL,    INTENT (IN) :: kdiff2(3,nntot2)
      75             :       REAL,    INTENT (IN) :: jj1(0:lmaxd,jmtd),jj2(0:lmaxd,jmtd)
      76             :       LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd),l_spav
      77             : C     ..
      78             : C     .. Local Scalars ..
      79             :       REAL gc1,gc2,c,e,r0
      80             :       COMPLEX cil
      81             :       INTEGER i,l,l2,lh,lm,lmp,lp,lo,lop
      82             :       INTEGER lp1,lpl,m,mp,mu
      83             :       INTEGER ljj1,ljj2,mjj1,mjj2
      84             :       INTEGER ltil1,mtil1,mtil2
      85             :       INTEGER lmini,lmini2, lmaxi,lmaxi2
      86             :       INTEGER ll,llp,lljj1,lljj2,lmjj1,lmjj2
      87             :       INTEGER lltil1,lmtil1,lmtil2,sp1,sp2
      88             :       INTEGER mmin,mmax,lwn, ikpt_b,ikpt_b2
      89             : C     ..
      90             : C     .. Local Arrays ..
      91           0 :       REAL, ALLOCATABLE :: dvulo(:,:,:,:,:),ulovd(:,:,:,:,:)
      92           0 :       REAL, ALLOCATABLE :: uvulo(:,:,:,:,:),ulovu(:,:,:,:,:)
      93           0 :       REAL, ALLOCATABLE :: ulovulo(:,:,:,:,:)
      94           0 :       REAL, ALLOCATABLE :: plo(:,:,:),plo_b(:,:,:)
      95           0 :       REAL, ALLOCATABLE :: x(:)
      96           0 :       REAL, ALLOCATABLE :: v0(:),vso(:,:)
      97             :       INTEGER :: ispjsp(2)
      98             :       DATA ispjsp/1,-1/
      99             : 
     100             : #if (defined(CPP_MPI) && !defined(CPP_T90))
     101             :       INCLUDE 'mpif.h'
     102             :       INTEGER ierr(3)
     103             : #endif
     104             : 
     105           0 :       c = c_light(1.0)
     106           0 :       sp1 = ispjsp(jspin)
     107           0 :       sp2 = ispjsp(jspin_b)
     108             : 
     109           0 :       allocate( ulovulo(0:lmaxd,0:lmaxd,0:lmaxd,nlod,nlod) )
     110           0 :       allocate( ulovu(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
     111           0 :       allocate( ulovd(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
     112           0 :       allocate( uvulo(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
     113           0 :       allocate( dvulo(0:lmaxd,0:lmaxd,0:lmaxd,0:lmaxd,nlod) )
     114           0 :       allocate( x(jmtd) )
     115           0 :       allocate( v0(jmtd), vso(jmtd,2) )
     116           0 :       allocate( plo(jmtd,0:lmaxd,nlod), plo_b(jmtd,0:lmaxd,nlod) )
     117             : 
     118           0 :       lwn = lmax(ntyp)
     119           0 :       r0 = rmsh(1,ntyp)
     120             : 
     121           0 :        DO lo = 1, nlo(ntyp)
     122           0 :         DO ljj1 = 0, lwn
     123           0 :          DO i=1,jri(ntyp)
     124           0 :           plo(i,ljj1,lo) = flo(ntyp,i,1,lo)*jj1(ljj1,i)
     125           0 :           plo_b(i,ljj1,lo) = flo_b(ntyp,i,1,lo)*jj2(ljj1,i)
     126             :          ENDDO
     127             :         ENDDO
     128             :        ENDDO
     129             : 
     130             : c****************************************************c
     131             : c   compute radial integrals                         c
     132             : c     <u(l')jj1(ljj1) | vso(lh) | u(l)jj2(ljj2)>     c
     133             : c****************************************************c
     134             : 
     135             : ! lo-lo
     136           0 :        DO lop = 1, nlo(ntyp)
     137           0 :         lp = llo(lop,ntyp)
     138           0 :         DO lo = 1, nlo(ntyp)
     139           0 :          l = llo(lop,ntyp)
     140             : 
     141             :         ! construct vso
     142           0 :         v0 = 0.
     143           0 :         IF (jspins.EQ.1) THEN
     144           0 :            v0(1:jri(ntyp)) = vr(1:jri(ntyp),ntyp,1)
     145           0 :            e = ( ello(lo,ntyp,1) + ello(lop,ntyp,1) )/2.
     146             :         ELSE
     147           0 :            DO i = 1,jri(ntyp)
     148           0 :              v0(i) = (vr(i,ntyp,1)+vr(i,ntyp,jspins))/2.
     149             :            END DO
     150             :            e = ( ello(lo,ntyp,1) + ello(lo,ntyp,jspins)
     151           0 :      >          +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
     152             :         END IF
     153             : 
     154           0 :         CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)  
     155             : 
     156             :         ! spin-averaging
     157           0 :         if(l_spav)then
     158           0 :             DO i= 1,jmtd
     159           0 :              vso(i,1)= (vso(i,1)+vso(i,2))/2.
     160           0 :              vso(i,2)= vso(i,1)
     161             :            ENDDO        
     162             :         endif
     163             : 
     164           0 :         if(l_socscreen) vso(:,:) = vso(:,:)*socscreen_fac
     165             : 
     166             :         ! calculate radial integrals
     167           0 :          DO ljj1 = 0, lwn
     168           0 :           DO ljj2 = 0, lwn
     169           0 :            DO lh=0,lwn
     170           0 :               lmini = abs(lp-ljj1)
     171           0 :               lmini2= abs(lh-ljj2)
     172           0 :               lmaxi = lp+ljj1
     173           0 :               lmaxi2= lh+ljj2
     174             :               if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
     175             :      >           (lmini2.gt.l) .or. (lmaxi2.lt.l).or.
     176           0 :      >           (mod(lp+lh+ljj1,2).ne.0).or.
     177           0 :      >           (mod( l+lh+ljj2,2).ne.0) ) then
     178             : 
     179           0 :                ulovulo(lh,ljj2,ljj1,lo,lop) = 0.
     180             : 
     181             :               else
     182             : 
     183           0 :                DO i=1,jri(ntyp)
     184           0 :                 x(i) = plo(i,ljj1,lop)*plo_b(i,ljj2,lo)*vso(i,jspin)
     185             :                ENDDO
     186             :                call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
     187           0 :      >                     ulovulo(lh,ljj2,ljj1,lo,lop))
     188             : 
     189             :               endif
     190             :            ENDDO
     191             :           ENDDO
     192             :          ENDDO
     193             :         ENDDO
     194             :        ENDDO!lh
     195             : 
     196             : ! apw-lo
     197           0 :        DO lop = 1, nlo(ntyp)
     198           0 :         lp = llo(lop,ntyp)
     199           0 :         DO l = 0, lwn
     200             : 
     201             :         ! construct vso
     202           0 :         v0 = 0.
     203           0 :         IF (jspins.EQ.1) THEN
     204           0 :            v0(1:jri(ntyp)) = vr(1:jri(ntyp),ntyp,1)
     205           0 :            e = ( epar(l,ntyp,1) + ello(lop,ntyp,1) )/2.
     206             :         ELSE
     207           0 :            DO i = 1,jri(ntyp)
     208           0 :              v0(i) = (vr(i,ntyp,1)+vr(i,ntyp,jspins))/2.
     209             :            END DO
     210             :            e = ( epar(l,ntyp,1) + epar(l,ntyp,jspins)
     211           0 :      >          +ello(lop,ntyp,1) + ello(lop,ntyp,jspins) )/4.
     212             :         END IF
     213             : 
     214           0 :         CALL sointg(ntyp,e,vr(:,ntyp,:),v0,atoms,input,vso)  
     215             : 
     216             :         ! spin-averaging
     217           0 :         if(l_spav)then
     218           0 :             DO i= 1,jmtd
     219           0 :              vso(i,1)= (vso(i,1)+vso(i,2))/2.
     220           0 :              vso(i,2)= vso(i,1)
     221             :            ENDDO        
     222             :         endif
     223             : 
     224             :         ! calculate radial integrals
     225           0 :          DO ljj1 = 0, lwn
     226           0 :           DO ljj2 = 0, lwn
     227           0 :            DO lh=0,lwn
     228           0 :               lmini = abs(lp-ljj1)
     229           0 :               lmini2= abs(lh-ljj2)
     230           0 :               lmaxi = lp+ljj1
     231           0 :               lmaxi2= lh+ljj2
     232             :               if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
     233             :      >           (lmini2.gt.l) .or. (lmaxi2.lt.l).or.
     234           0 :      >           (mod(lp+lh+ljj1,2).ne.0).or.
     235             :      >           (mod( l+lh+ljj2,2).ne.0) ) then
     236             : 
     237           0 :                ulovu(lh,l,ljj2,ljj1,lop) = 0.
     238           0 :                ulovd(lh,l,ljj2,ljj1,lop) = 0.
     239             : 
     240             :               else
     241             : 
     242           0 :                DO i=1,jri(ntyp)
     243           0 :                 x(i) = plo(i,ljj1,lop)*p_b(i,l,ljj2)*vso(i,jspin)
     244             :                ENDDO
     245             :                call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
     246           0 :      >                     ulovu(lh,l,ljj2,ljj1,lop))
     247           0 :                DO i=1,jri(ntyp)
     248           0 :                 x(i) = plo(i,ljj1,lop)*q_b(i,l,ljj2)*vso(i,jspin)
     249             :                ENDDO
     250             :                call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
     251           0 :      >                     ulovd(lh,l,ljj2,ljj1,lop))
     252             : 
     253             :               endif
     254             : 
     255           0 :               lmini = abs(l-ljj1)
     256           0 :               lmini2= abs(lh-ljj2)
     257           0 :               lmaxi = l+ljj1
     258           0 :               lmaxi2= lh+ljj2
     259             :               if((lmini.gt.lh) .or. (lmaxi.lt.lh).or.
     260             :      >           (lmini2.gt.lp) .or. (lmaxi2.lt.lp).or.
     261           0 :      >           (mod( l+lh+ljj1,2).ne.0).or.
     262           0 :      >           (mod(lp+lh+ljj2,2).ne.0) ) then
     263             : 
     264           0 :                uvulo(lh,l,ljj2,ljj1,lop) = 0.
     265           0 :                dvulo(lh,l,ljj2,ljj1,lop) = 0.
     266             : 
     267             :               else
     268             : 
     269           0 :                DO i=1,jri(ntyp)
     270           0 :                 x(i) = p(i,l,ljj1)*plo_b(i,ljj2,lop)*vso(i,jspin)
     271             :                ENDDO
     272             :                call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
     273           0 :      >                     uvulo(lh,l,ljj2,ljj1,lop))
     274           0 :                DO i=1,jri(ntyp)
     275           0 :                 x(i) = q(i,l,ljj1)*plo_b(i,ljj2,lop)*vso(i,jspin)
     276             :                ENDDO
     277             :                call intgr0(x(1:jri(ntyp)),r0,dx(ntyp),jri(ntyp),
     278           0 :      >                     dvulo(lh,l,ljj2,ljj1,lop))
     279             :              
     280             :               endif
     281             :            ENDDO
     282             :           ENDDO
     283             :          ENDDO
     284             :         ENDDO
     285             :        ENDDO!lh
     286             : 
     287             : c***************** SOC CONTRIBUTION **********************c
     288             : c   compute product of the two Gaunt coefficients         c
     289             : c   with the radial integrals (+prefactors)               c
     290             : c*********************************************************c
     291             : c   We deal with two Gaunt coefficients:                  c
     292             : c  G1 = G( (ltil1, mtil1), (ljj1, mjj1), (lp, mp)       ) c
     293             : c  G2 = G( (l , m)       , (ljj2, mjj2), (ltil1, mtil2) ) c
     294             : c*********************************************************c
     295             : c   use Gaunt conditions to reduce number of operations.  c
     296             : c   coefficient G(L1,L2,L3) only nonzero if               c
     297             : c       a)  l1 + l2 + l3 = even                           c
     298             : c       b)  |l2-l3| <= l1 <= l2+l3                        c
     299             : c       c)  m1 = m2 + m3                                  c
     300             : c*********************************************************c
     301             : 
     302             : ! lo-lo
     303           0 :        DO lop=1,nlo(ntyp)
     304           0 :         lp = llo(lop,ntyp)
     305           0 :         DO mp=-lp,lp
     306           0 :          DO ltil1=0,lwn
     307           0 :           lltil1 = ltil1*(ltil1+1)
     308           0 :           DO mtil1=-ltil1,ltil1
     309           0 :            lmtil1 = lltil1+mtil1
     310           0 :            mmin = max(-ltil1,mtil1-1)
     311           0 :            mmax = min( ltil1,mtil1+1)
     312           0 :            mjj1=mtil1-mp
     313           0 :            DO ljj1=0,lwn
     314           0 :             lljj1 = ljj1 * (ljj1 + 1)
     315           0 :             lmjj1 = lljj1 + mjj1
     316           0 :             lmini = abs(lp-ljj1)
     317           0 :             lmaxi = lp+ljj1
     318             :             ! Gaunt conditions (G1):
     319             :             if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
     320           0 :      >         (mod(lp+ltil1+ljj1,2).ne.0).or.
     321             :      >         (abs(mjj1).gt.ljj1)) cycle
     322             : 
     323             :              gc1 = gaunt1(ltil1, ljj1, lp, 
     324           0 :      >                    mtil1, mjj1, mp, lmaxd)
     325             : 
     326           0 :         DO lo=1,nlo(ntyp)
     327           0 :          l = llo(lo,ntyp)
     328           0 :          DO m=-l,l
     329           0 :           DO mtil2=mmin,mmax     ! mtil1-1 <= mtil2 <= mtil1+1
     330           0 :            lmtil2 = lltil1+mtil2
     331           0 :            mjj2=m-mtil2
     332           0 :            DO ljj2=0,lwn
     333           0 :             lljj2 = ljj2 * (ljj2 + 1)
     334           0 :             lmjj2 = lljj2 + mjj2
     335           0 :             lmini2 = abs(ltil1-ljj2)
     336           0 :             lmaxi2 = ltil1+ljj2
     337             :              ! Gaunt conditions (G2):
     338             :              if((lmini2.gt.l).or.(lmaxi2.lt.l).or.
     339           0 :      >          (mod(l+ltil1+ljj2,2).ne.0).or.
     340             :      >          (abs(mjj2).gt.ljj2)) cycle
     341             : 
     342             :              gc2 = gaunt1( l, ljj2, ltil1,
     343           0 :      >                     m, mjj2, mtil2, lmaxd)
     344             : 
     345             :              ! set up prefactor
     346             :              cil =   ( ImagUnit ** (lp - l + ljj1 + ljj2) )
     347             :      +              *( (-1.0) ** ljj1 )
     348             :      +              * gc1 * gc2
     349             :      +              * conjg( yl1(lmjj1 + 1) )
     350             :      +              * conjg( yl2(lmjj2 + 1) )
     351           0 :      +              * conjg( angso(lmtil1,lmtil2) )
     352             : 
     353             :              ! compute T-coefficients
     354             :              tuloulo(lo, m, lop, mp) = tuloulo(lo, m, lop, mp)
     355           0 :      >         + cil * ulovulo(ltil1, ljj2, ljj1, lo, lop)
     356             : 
     357             :            ENDDO!ljj2
     358             :           ENDDO!m
     359             :          ENDDO!l
     360             :         ENDDO!mtil2
     361             :            ENDDO!ljj1
     362             :           ENDDO!mp
     363             :          ENDDO!lp
     364             :         ENDDO!mtil1
     365             :        ENDDO!ltil1
     366             : 
     367             : ! apw-lo
     368           0 :        DO lop=1,nlo(ntyp)
     369           0 :         lp = llo(lop,ntyp)
     370           0 :         DO mp=-lp,lp
     371           0 :          DO ltil1=0,lwn
     372           0 :           lltil1 = ltil1*(ltil1+1)
     373           0 :           DO mtil1=-ltil1,ltil1
     374           0 :            lmtil1 = lltil1+mtil1
     375           0 :            mmin = max(-ltil1,mtil1-1)
     376           0 :            mmax = min( ltil1,mtil1+1)
     377           0 :            mjj1=mtil1-mp
     378           0 :            DO ljj1=0,lwn
     379           0 :             lljj1 = ljj1 * (ljj1 + 1)
     380           0 :             lmjj1 = lljj1 + mjj1
     381           0 :             lmini = abs(lp-ljj1)
     382           0 :             lmaxi = lp+ljj1
     383             :             ! Gaunt conditions (G1):
     384             :             if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
     385           0 :      >         (mod(lp+ltil1+ljj1,2).ne.0).or.
     386             :      >         (abs(mjj1).gt.ljj1)) cycle
     387             : 
     388             :              gc1 = gaunt1(ltil1, ljj1, lp, 
     389           0 :      >                    mtil1, mjj1, mp, lmaxd)
     390             : 
     391           0 :         DO mtil2=mmin,mmax     ! mtil1-1 <= mtil2 <= mtil1+1
     392           0 :          lmtil2 = lltil1+mtil2
     393           0 :          DO l=0,lwn
     394           0 :           ll = l*(l+1)
     395           0 :           DO m=-l,l
     396           0 :            lm = ll+m
     397           0 :            mjj2=m-mtil2
     398           0 :            DO ljj2=0,lwn
     399           0 :             lljj2 = ljj2 * (ljj2 + 1)
     400           0 :             lmjj2 = lljj2 + mjj2
     401           0 :             lmini2 = abs(ltil1-ljj2)
     402           0 :             lmaxi2 = ltil1+ljj2
     403             :              ! Gaunt conditions (G2):
     404             :              if((lmini2.gt.l).or.(lmaxi2.lt.l).or.
     405           0 :      >          (mod(l+ltil1+ljj2,2).ne.0).or.
     406             :      >          (abs(mjj2).gt.ljj2)) cycle
     407             : 
     408             :              gc2 = gaunt1( l, ljj2, ltil1,
     409           0 :      >                     m, mjj2, mtil2, lmaxd)
     410             : 
     411             :              ! set up prefactor
     412             :              cil =   ( ImagUnit ** (lp - l + ljj1 + ljj2) )
     413             :      +              *( (-1.0) ** ljj1 )
     414             :      +              * gc1 * gc2
     415             :      +              * conjg( yl1(lmjj1 + 1) )
     416             :      +              * conjg( yl2(lmjj2 + 1) )
     417           0 :      +              * conjg( angso(lmtil1,lmtil2) )
     418             : 
     419             :              ! compute T-coefficients
     420             :              tulou(lm, lop, mp) = tulou(lm, lop, mp)
     421           0 :      >         + cil * ulovu(ltil1, l, ljj2, ljj1, lop)
     422             : 
     423             :              tulod(lm, lop, mp) = tulod(lm, lop, mp)
     424           0 :      >         + cil * ulovd(ltil1, l, ljj2, ljj1, lop)
     425             : 
     426             :            ENDDO!ljj2
     427             :           ENDDO!m
     428             :          ENDDO!l
     429             :         ENDDO!mtil2
     430             :            ENDDO!ljj1
     431             :           ENDDO!mp
     432             :          ENDDO!lp
     433             :         ENDDO!mtil1
     434             :        ENDDO!ltil1
     435             : 
     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 :          DO ltil1=0,lwn
     442           0 :           lltil1 = ltil1*(ltil1+1)
     443           0 :           DO mtil1=-ltil1,ltil1
     444           0 :            lmtil1 = lltil1+mtil1
     445           0 :            mmin = max(-ltil1,mtil1-1)
     446           0 :            mmax = min( ltil1,mtil1+1)
     447           0 :            mjj1=mtil1-m
     448           0 :            DO ljj1=0,lwn
     449           0 :             lljj1 = ljj1 * (ljj1 + 1)
     450           0 :             lmjj1 = lljj1 + mjj1
     451           0 :             lmini = abs(l-ljj1)
     452           0 :             lmaxi = l+ljj1
     453             :             ! Gaunt conditions (G1):
     454             :             if((lmini.gt.ltil1).or.(lmaxi.lt.ltil1).or.
     455           0 :      >         (mod(l+ltil1+ljj1,2).ne.0).or.
     456             :      >         (abs(mjj1).gt.ljj1)) cycle
     457             : 
     458             :              gc1 = gaunt1(ltil1, ljj1, l, 
     459           0 :      >                    mtil1, mjj1, m, lmaxd)
     460             : 
     461           0 :         DO mtil2=mmin,mmax     ! mtil1-1 <= mtil2 <= mtil1+1
     462           0 :          lmtil2 = lltil1+mtil2
     463           0 :          DO lop=1,nlo(ntyp)
     464           0 :           lp = llo(lop,ntyp)
     465           0 :           DO mp=-lp,lp
     466           0 :            mjj2=mp-mtil2
     467           0 :            DO ljj2=0,lwn
     468           0 :             lljj2 = ljj2 * (ljj2 + 1)
     469           0 :             lmjj2 = lljj2 + mjj2
     470           0 :             lmini2 = abs(ltil1-ljj2)
     471           0 :             lmaxi2 = ltil1+ljj2
     472             :              ! Gaunt conditions (G2):
     473             :              if((lmini2.gt.lp).or.(lmaxi2.lt.lp).or.
     474           0 :      >          (mod(lp+ltil1+ljj2,2).ne.0).or.
     475             :      >          (abs(mjj2).gt.ljj2)) cycle
     476             : 
     477             :              gc2 = gaunt1( lp, ljj2, ltil1,
     478           0 :      >                     mp, mjj2, mtil2, lmaxd)
     479             : 
     480             :              ! set up prefactor
     481             :              cil =   ( ImagUnit ** (l - lp + ljj1 + ljj2) )
     482             :      +              *( (-1.0) ** ljj1 )
     483             :      +              * gc1 * gc2
     484             :      +              * conjg( yl1(lmjj1 + 1) )
     485             :      +              * conjg( yl2(lmjj2 + 1) )
     486           0 :      +              * conjg( angso(lmtil1,lmtil2) )
     487             : 
     488             :              ! compute T-coefficients
     489             :              tuulo(lm, lop, mp) = tuulo(lm, lop, mp)
     490           0 :      >         + cil * uvulo(ltil1, l, ljj2, ljj1, lop)
     491             : 
     492             :              tdulo(lm, lop, mp) = tdulo(lm, lop, mp)
     493           0 :      >         + cil * dvulo(ltil1, l, ljj2, ljj1, lop)
     494             : 
     495             :            ENDDO!ljj2
     496             :           ENDDO!m
     497             :          ENDDO!l
     498             :         ENDDO!mtil2
     499             :            ENDDO!ljj1
     500             :           ENDDO!mp
     501             :          ENDDO!lp
     502             :         ENDDO!mtil1
     503             :        ENDDO!ltil1
     504             : 
     505             : #if (defined(CPP_MPI) && !defined(CPP_T90))
     506           0 :       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
     507             : #endif
     508             : 
     509           0 :       deallocate( ulovulo )
     510           0 :       deallocate( uvulo, dvulo, ulovu, ulovd )
     511           0 :       deallocate( plo,plo_b )
     512           0 :       deallocate( x )
     513           0 :       deallocate( v0, vso )
     514             : 
     515           0 :       END SUBROUTINE wann_uHu_soc_tlo
     516             :       END MODULE m_wann_uHu_soc_tlo

Generated by: LCOV version 1.13