nmat_rot.f 3.42 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
      MODULE m_nmat_rot

! Calculate the Wigner rotation matrices for complex spherical
! harmonics for all space-group rotations and l=1,2,3. Needed 
! for the calculation of the density matrix in nmat.
!
! also allows to use rotated "n_mmp_mat" file by specifying a 
! "n_mmp_rot" file (see its use in u_mix or u_setup) 
!                                                       gb10
      CONTAINS
      SUBROUTINE nmat_rot(
     >                    alpha,beta,gamma,l_in,n_u,jspins,lty,
     X                    n_mmp)

      USE m_inv3

      IMPLICIT NONE

! .. arguments:
      INTEGER, INTENT(IN)  :: l_in,n_u,jspins,lty(n_u)
      REAL,    INTENT(IN)  :: alpha(n_u),beta(n_u),gamma(n_u)
      COMPLEX, INTENT(INOUT) :: n_mmp(-3:3,-3:3,n_u,jspins)

! .. local variables:
      INTEGER ns,signum,ispin,n
      INTEGER i,j,k,l,m,mp,x_lo,x_up,x,e_c,e_s
      REAL fac_l_m,fac_l_mp,fac_lmpx,fac_lmx,fac_x,fac_xmpm
      REAL co_bh,si_bh,zaehler,nenner,cp,sp
      REAL sina,sinb,sinc,cosa,cosb,cosc,determ,dt
      COMPLEX ci,phase_g,phase_a,bas,d(-l_in:l_in,-l_in:l_in)
      COMPLEX d_wig(-l_in:l_in,-l_in:l_in,l_in,n_u)
      COMPLEX n_tmp(-l_in:l_in,-l_in:l_in)
      COMPLEX nr_tmp(-l_in:l_in,-l_in:l_in)
      LOGICAL, SAVE :: written = .false.

      REAL dmat(3,3),dmati(3,3)

      INTRINSIC sqrt,max,min

      ci = cmplx(0.0,1.0)
      
      DO n = 1, n_u

      co_bh = cos(beta(n)*0.5)
      si_bh = sin(beta(n)*0.5)

      DO l = 1, lty(n)
        d = (0.0,0.0)

        DO m = -l,l
          fac_l_m = fac(l+m) * fac(l-m)
          phase_g = exp( - ci * gamma(n) * m )

          DO mp = -l,l
            fac_l_mp = fac(l+mp) * fac(l-mp)

            zaehler = sqrt( real(fac_l_m * fac_l_mp) )
            phase_a = exp( - ci * alpha(n) * mp ) 
            x_lo = max(0, m-mp)
            x_up = min(l-mp, l+m)

            bas = zaehler * phase_a * phase_g 
            d(m,mp) = cmplx(0.0,0.0)
            DO x = x_lo,x_up
              fac_lmpx = fac(l-mp-x)
              fac_lmx  = fac(l+m-x)
              fac_x    = fac(x)
              fac_xmpm = fac(x+mp-m)
              nenner = fac_lmpx * fac_lmx * fac_x * fac_xmpm
              e_c = 2*l + m - mp - 2*x 
              e_s = 2*x + mp - m
              IF (e_c.EQ.0) THEN
                cp = 1.0
              ELSE
                cp = co_bh ** e_c
              ENDIF
              IF (e_s.EQ.0) THEN
                sp = 1.0
              ELSE
                sp = si_bh ** e_s
              ENDIF
              d(m,mp) = d(m,mp) + bas * (-1)**x * cp * sp / nenner
            ENDDO

          ENDDO ! loop over mp
        ENDDO   ! loop over m
        DO m = -l,l
          DO mp = -l,l
            d( m,mp ) = d( m,mp ) * (-1)**(m-mp)
          ENDDO
        ENDDO
        d_wig(:,:,l,n) = d(:,:)

      ENDDO ! l
      ENDDO ! n 

      DO ispin = 1, jspins
        DO n = 1, n_u
           n_tmp(:,:) = n_mmp(-l_in:l_in,-l_in:l_in,n,ispin)
           d(:,:) = d_wig(:,:,lty(n),n)

           nr_tmp = matmul( transpose( conjg(d) ) , n_tmp)
           n_tmp =  matmul( nr_tmp, d )

           n_mmp(-l_in:l_in,-l_in:l_in,n,ispin) = n_tmp(:,:)
         ENDDO
      ENDDO
      !write(*,'(14f8.4)') n_mmp

      END SUBROUTINE nmat_rot

      ELEMENTAL REAL FUNCTION  fac(n)

      INTEGER, INTENT (IN) :: n
      INTEGER :: i
 
      fac = 0
      IF (n.LT.0) RETURN
      fac = 1
      IF (n.EQ.0) RETURN
      DO i = 2,n
        fac = fac * i
      ENDDO

      END FUNCTION  fac
      
      END MODULE m_nmat_rot