LCOV - code coverage report
Current view: top level - math - d_wigner.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 100 106 94.3 %
Date: 2019-09-08 04:53:50 Functions: 2 3 66.7 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_dwigner
       8             :       use m_juDFT
       9             : 
      10             : ! Calculate the Wigner rotation matrices for complex spherical
      11             : ! harmonics for all space-group rotations and l=1,2,3. Needed 
      12             : ! for the calculation of the density matrix in nmat.
      13             : !                 gb 2002
      14             : 
      15             : c******************************************************
      16             : c     interface for real and integer rotation matrices
      17             : c       FF, Oct 2006
      18             : c*****************************************************
      19             :       PRIVATE
      20             :       INTERFACE d_wigner
      21             :       MODULE PROCEDURE real_wigner, integer_wigner, integer_wigner_1op
      22             :       END INTERFACE
      23             :       LOGICAL :: written=.false.
      24             :       PUBLIC :: d_wigner
      25             : 
      26             :                 
      27             :       CONTAINS
      28             : c***************************************************
      29             : c        private routine for integer rotation where only
      30             : c        one rotation is passed 
      31             : c***************************************************
      32           0 :       SUBROUTINE integer_wigner_1op(
      33             :      >                          mrot,bmat,lmax,
      34           0 :      <                          d_wgn)
      35             : 
      36             :       INTEGER, INTENT(IN)  :: lmax
      37             :       INTEGER, INTENT(IN)  :: mrot(3,3)
      38             :       REAL,    INTENT(IN)  :: bmat(3,3)
      39             :       COMPLEX, INTENT(OUT) :: d_wgn(-lmax:lmax,-lmax:lmax,lmax)
      40             :       REAL                    realmrot(3,3,1)
      41             :       
      42           0 :       realmrot(:,:,1) = mrot(:,:)
      43             :       CALL real_wigner(
      44             :      >                 1,realmrot,bmat,lmax,
      45           0 :      <                 d_wgn)
      46             : 
      47           0 :       END SUBROUTINE
      48             : 
      49             : c***************************************************
      50             : c        private routine for integer rotation
      51             : c***************************************************
      52           3 :       SUBROUTINE integer_wigner(
      53           3 :      >                          nop,mrot,bmat,lmax,
      54           3 :      <                          d_wgn)
      55             : 
      56             :       INTEGER, INTENT(IN)  :: nop,lmax
      57             :       INTEGER, INTENT(IN)  :: mrot(3,3,nop)
      58             :       REAL,    INTENT(IN)  :: bmat(3,3)
      59             :       COMPLEX, INTENT(OUT) :: d_wgn(-lmax:lmax,-lmax:lmax,lmax,nop)
      60           6 :       REAL                    realmrot(3,3,nop)
      61             :       
      62          21 :       realmrot(:,:,:) = mrot(:,:,:)
      63             :       CALL real_wigner(
      64             :      >                 nop,realmrot,bmat,lmax,
      65           3 :      <                 d_wgn)
      66             : 
      67           3 :       END SUBROUTINE
      68             : 
      69             : c**************************************************
      70             : c         private routine for real rotation
      71             : c**************************************************
      72           3 :       SUBROUTINE real_wigner(
      73           3 :      >                    nop,mrot,bmat,lmax,
      74           3 :      <                    d_wgn)
      75             : 
      76             :       USE m_constants, ONLY : pimach, ImagUnit
      77             :       USE m_inv3
      78             : 
      79             :       IMPLICIT NONE
      80             : 
      81             : ! .. arguments:
      82             :       INTEGER, INTENT(IN)  :: nop,lmax
      83             :       REAL,    INTENT(IN)  :: mrot(3,3,nop)
      84             :       REAL,    INTENT(IN)  :: bmat(3,3)
      85             :       COMPLEX, INTENT(OUT) :: d_wgn(-lmax:lmax,-lmax:lmax,lmax,nop)
      86             : 
      87             : ! .. local variables:
      88             :       INTEGER              :: ns,signum
      89             :       INTEGER              :: i,j,k,l,m,mp,x_lo,x_up,x,e_c,e_s
      90             :       REAL                 :: fac_l_m,fac_l_mp,fac_lmpx,fac_lmx,
      91             :      +                        fac_x,fac_xmpm
      92             :       REAL                 :: pi,co_bh,si_bh,zaehler,nenner,cp,sp
      93             :       REAL                 :: sina,sinb,sinc,cosa,cosb,cosc,determ,dt
      94             :       COMPLEX              :: phase_g,phase_a,bas,
      95           6 :      +                        d(-lmax:lmax,-lmax:lmax)
      96             : 
      97          12 :       REAL                 :: alpha(nop),beta(nop),gamma(nop)
      98           6 :       REAL                 :: dmat(3,3),dmati(3,3),det(nop),bmati(3,3)
      99             : 
     100             :       INTRINSIC sqrt,max,min
     101             : 
     102             : 
     103           3 :       pi = pimach()
     104             : c
     105             : c determine the eulerian angles of all the rotations
     106             : c
     107           3 :       CALL inv3(bmat,bmati,dt)
     108          21 :       DO ns = 1, nop
     109             : c
     110             : c first determine the determinant of the rotation
     111             : c +1 for a proper rotation
     112             : c -1 for a proper rotation times inversion
     113             : c
     114             :         determ = 0.00
     115         126 :         DO i = 1,3
     116         396 :           DO j = 1,3
     117         216 :             IF (i.NE.j) THEN
     118         108 :               k = 6 - i - j
     119         108 :               signum = 1
     120         108 :               IF ( (i.EQ.(j+1)).OR.(j.EQ.(k+1)) ) signum=-signum
     121             :               determ = determ + signum*
     122         108 :      +                 mrot(i,1,ns)*mrot(j,2,ns)*mrot(k,3,ns)
     123             :             ENDIF
     124             :           ENDDO
     125             :         ENDDO
     126          18 :         IF (abs(1.0-abs(determ)).GT.1.0e-5)
     127             :      +       CALL juDFT_error("d_wigner: determ/==+-1",calledby
     128           0 :      +       ="d_wigner")
     129          18 :         det(ns) = determ
     130             : c
     131             : c store the proper part of the rotation in dmati and convert to
     132             : c cartesian coordinates
     133             : c
     134          18 :         dmati = determ*matmul(bmati,matmul(mrot(:,:,ns),bmat))
     135             : c
     136             : c the eulerian angles are derived from the inverse of
     137             : c dmati, because we use the convention that we rotate functions
     138             : c
     139          18 :         CALL inv3(dmati,dmat,dt)
     140             : c
     141             : c beta follows directly from d33
     142             : c
     143          18 :         cosb = dmat(3,3)
     144          18 :         sinb = 1.00 - cosb*cosb
     145          18 :         sinb = max(sinb,0.00)
     146          18 :         sinb = sqrt(sinb)
     147             : c
     148             : c if beta = 0 or pi , only alpha+gamma or -gamma have a meaning:
     149             : c
     150          21 :         IF ( abs(sinb).LT.1.0e-5 ) THEN 
     151             : 
     152           6 :           beta(ns) = 0.0
     153           6 :           IF ( cosb.lt.0.0 ) beta(ns) = pi
     154           6 :           gamma(ns) = 0.0
     155           6 :           cosa = dmat(1,1)/cosb
     156           6 :           sina = dmat(1,2)/cosb
     157           6 :           IF ( abs(sina).LT.1.0e-5 ) THEN
     158           3 :             alpha(ns)=0.0
     159           3 :             IF ( cosa.LT.0.0 ) alpha(ns)=alpha(ns)+pi
     160             :           ELSE
     161           3 :             alpha(ns) = 0.5*pi - atan(cosa/sina)
     162           3 :             IF ( sina.LT.0.0 ) alpha(ns)=alpha(ns)+pi
     163             :           ENDIF
     164             : 
     165             :         ELSE
     166             : 
     167          12 :           beta(ns) = 0.5*pi - atan(cosb/sinb)
     168             : c
     169             : c determine alpha and gamma from d13 d23 d32 d31
     170             : c
     171          12 :           cosa = dmat(3,1)/sinb
     172          12 :           sina = dmat(3,2)/sinb
     173          12 :           cosc =-dmat(1,3)/sinb
     174          12 :           sinc = dmat(2,3)/sinb
     175          12 :           IF ( abs(sina).lt.1.0e-5 ) THEN
     176           6 :             alpha(ns)=0.0
     177           6 :             IF ( cosa.LT.0.0 ) alpha(ns)=alpha(ns)+pi
     178             :           ELSE
     179           6 :             alpha(ns) = 0.5*pi - atan(cosa/sina)
     180           6 :             IF ( sina.LT.0.0 ) alpha(ns)=alpha(ns)+pi
     181             :           ENDIF
     182          12 :           IF ( abs(sinc).lt.1.0e-5 ) THEN
     183           6 :             gamma(ns) = 0.0
     184           6 :             IF ( cosc.LT.0.0 ) gamma(ns)=gamma(ns)+pi
     185             :           ELSE
     186           6 :             gamma(ns) = 0.5*pi - atan(cosc/sinc)
     187           6 :             IF ( sinc.LT.0.0 ) gamma(ns)=gamma(ns)+pi
     188             :           ENDIF
     189             : 
     190             :         ENDIF
     191             : 
     192             :       ENDDO ! loop over nop
     193             : 
     194             : #ifndef CPP_MPI
     195             :       IF(.NOT.written) THEN
     196             :         WRITE (6,8000)
     197             :         DO ns = 1, nop
     198             :           WRITE (6,8010) ns,alpha(ns),beta(ns),gamma(ns),det(ns)
     199             :         ENDDO
     200             :         written=.true.
     201             :       ENDIF
     202             :  8000 FORMAT(//,'   eulerian angles for the rotations ',
     203             :      $  //,'   ns   alpha     beta      gamma    determ ')
     204             :  8010 FORMAT(i5,4f10.5)
     205             : #endif
     206          21 :       DO ns = 1, nop
     207             : 
     208          18 :         co_bh = cos(beta(ns)*0.5)
     209          18 :         si_bh = sin(beta(ns)*0.5)
     210             : 
     211          75 :         DO l = 1, lmax
     212             : 
     213         324 :           DO m = -l,l
     214         810 :             fac_l_m = fac(l+m) * fac(l-m)
     215         270 :             phase_g = exp( - ImagUnit * gamma(ns) * m ) 
     216             : 
     217        1818 :             DO mp = -l,l
     218        4482 :               fac_l_mp = fac(l+mp) * fac(l-mp)
     219             : 
     220        1494 :               zaehler = sqrt( real(fac_l_m * fac_l_mp) )
     221        1494 :               phase_a = exp( - ImagUnit * alpha(ns) * mp ) 
     222        1494 :               x_lo = max(0, m-mp)
     223        1494 :               x_up = min(l-mp, l+m)
     224             : 
     225        1494 :               bas = zaehler * phase_a * phase_g 
     226        1494 :               d(m,mp) = cmplx(0.0,0.0)
     227        4086 :               DO x = x_lo,x_up
     228        4644 :                 fac_lmpx = fac(l-mp-x)
     229        4644 :                 fac_lmx  = fac(l+m-x)
     230        2322 :                 fac_x    = fac(x)
     231        4644 :                 fac_xmpm = fac(x+mp-m)
     232        2322 :                 nenner = fac_lmpx * fac_lmx * fac_x * fac_xmpm
     233        2322 :                 e_c = 2*l + m - mp - 2*x 
     234        2322 :                 e_s = 2*x + mp - m
     235        2322 :                 IF (e_c.EQ.0) THEN
     236             :                   cp = 1.0
     237             :                 ELSE
     238        2052 :                   cp = co_bh ** e_c
     239             :                 ENDIF
     240        2322 :                 IF (e_s.EQ.0) THEN
     241             :                   sp = 1.0
     242             :                 ELSE
     243        2052 :                   sp = si_bh ** e_s
     244             :                 ENDIF
     245        3816 :                 d(m,mp) = d(m,mp) + bas * (-1)**x * cp * sp / nenner
     246             :               ENDDO
     247             : 
     248             :             ENDDO ! loop over mp
     249             :           ENDDO   ! loop over m
     250             : 
     251         594 :           DO m = -l,l
     252        3312 :             DO mp = -l,l
     253        1764 :              d( m,mp ) = d( m,mp ) * (-1)**(m-mp)
     254             :             ENDDO
     255             :           ENDDO
     256         612 :           DO m = -l,l
     257        3312 :             DO mp = -l,l
     258        1764 :               IF(abs(det(ns)+1).lt.1e-5) THEN
     259         747 :                 d_wgn(m,mp,l,ns) = d( m,mp) * (-1)**l ! adds inversion
     260             :               ELSE
     261         747 :                 d_wgn(m,mp,l,ns) = d( m,mp)
     262             :               ENDIF
     263             :             ENDDO
     264             :           ENDDO
     265             : 
     266             :         ENDDO
     267             :       ENDDO
     268             : 
     269           3 :       END SUBROUTINE real_wigner
     270             : 
     271             :       ELEMENTAL REAL FUNCTION  fac(n)
     272             : 
     273             :       INTEGER, INTENT (IN) :: n
     274             :       INTEGER :: i
     275             :  
     276       12816 :       fac = 0
     277       12816 :       IF (n.LT.0) RETURN
     278       12816 :       fac = 1
     279       12816 :       IF (n.EQ.0) RETURN
     280       32544 :       DO i = 2,n
     281       11952 :         fac = fac * i
     282             :       ENDDO
     283             : 
     284             :       END FUNCTION  fac
     285             :       
     286             :       END MODULE m_dwigner

Generated by: LCOV version 1.13