n_mat.f90 5.64 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13
MODULE m_nmat
  !     ************************************************************
  !     This subroutine calculates the density matrix n^{s}_{m,m'}
  !     for a given atom 'n' and l-quantum number 'l'. The l's for
  !     all atoms are stored in lda_u(), if lda_u()<0, no +U is used.
  !     For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999)
  !     Part of the LDA+U package                   G.B., Oct. 2000
Daniel Wortmann's avatar
Daniel Wortmann committed
14
  !     Extension to multiple U per atom type by G.M. 2017
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
  !     ************************************************************
CONTAINS
  SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we, acof,bcof,ccof, n_mmp)
    !

    USE m_types
    IMPLICIT NONE
    TYPE(t_usdus),INTENT(IN)   :: usdus
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_atoms),INTENT(IN)   :: atoms
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ne,jspin 
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: we(:)!(dimension%neigd)
    COMPLEX, INTENT (IN) :: acof(:,0:,:)!(nobd,0:atoms%lmaxd*(lmaxd+2) ,natd)
    COMPLEX, INTENT (IN) :: bcof(:,0:,:)!(nobd,0:atoms%lmaxd*(lmaxd+2) ,natd)
Daniel Wortmann's avatar
Daniel Wortmann committed
33
    COMPLEX, INTENT (IN) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
34 35 36 37
    COMPLEX, INTENT (INOUT) :: n_mmp(-3:3,-3:3,atoms%n_u)
    !     ..
    !     .. Local Scalars ..
    COMPLEX c_0
Daniel Wortmann's avatar
Daniel Wortmann committed
38
    INTEGER i,j,k,l ,mp,n,it,is,isi,natom,natomTemp,n_ldau,lp,m,i_u
39 40 41 42 43 44 45 46 47 48 49
    INTEGER ilo,ilop,ll1,nn,lmp,lm
    REAL fac
    !     ..
    !     .. Local Arrays ..
    COMPLEX n_tmp(-3:3,-3:3),nr_tmp(-3:3,-3:3),d_tmp(-3:3,-3:3)
    COMPLEX n1_tmp(-3:3,-3:3)
    !     ..
    !
    ! calculate n_mat:
    !
    natom = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
50
    i_u = 1
51
    DO n = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
52 53 54 55 56
       DO WHILE (i_u.LE.atoms%n_u)
          IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
          natomTemp = natom
          n_tmp(:,:) = cmplx(0.0,0.0)
          l = atoms%lda_u(i_u)%l
57 58
          ll1 = (l+1)*l 
          DO nn = 1, atoms%neq(n)
Daniel Wortmann's avatar
Daniel Wortmann committed
59
             natomTemp = natomTemp + 1
60 61 62 63 64 65 66 67 68 69 70
             !
             !  prepare n_mat in local frame (in noco-calculations this depends 
             !                                also on alpha(n) and beta(n) )
             !
             DO m = -l,l
                lm = ll1+m
                DO mp = -l,l
                   lmp = ll1+mp
                   c_0 = cmplx(0.0,0.0)
                   DO i = 1,ne
                      c_0 = c_0 +  we(i) * ( usdus%ddn(l,n,jspin) *&
Daniel Wortmann's avatar
Daniel Wortmann committed
71 72
                           conjg(bcof(i,lmp,natomTemp))*bcof(i,lm,natomTemp) +&
                           conjg(acof(i,lmp,natomTemp))*acof(i,lm,natomTemp) )
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
                   ENDDO
                   n_tmp(m,mp) = c_0 
                ENDDO
             ENDDO
             !
             !  add local orbrbital contribution (if there is one) (untested so far)
             !
             DO ilo = 1, atoms%nlo(n)
                IF (atoms%llo(ilo,n).EQ.l) THEN

                   DO m = -l,l
                      lm = ll1+m
                      DO mp = -l,l
                         lmp = ll1+mp
                         c_0 = cmplx(0.0,0.0)
                         DO i = 1,ne
                            c_0 = c_0 +  we(i) * (  usdus%uulon(ilo,n,jspin) * (&
Daniel Wortmann's avatar
Daniel Wortmann committed
90 91
                                 conjg(acof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
                                 conjg(ccof(mp,i,ilo,natomTemp))*acof(i,lm,natomTemp) )&
92
                                 + usdus%dulon(ilo,n,jspin) * (&
Daniel Wortmann's avatar
Daniel Wortmann committed
93 94
                                 conjg(bcof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
                                 conjg(ccof(mp,i,ilo,natomTemp))*bcof(i,lm,natomTemp)))
95 96 97 98 99
                         ENDDO
                         DO ilop = 1, atoms%nlo(n)
                            IF (atoms%llo(ilop,n).EQ.l) THEN
                               DO i = 1,ne
                                  c_0 = c_0 +  we(i) * usdus%uloulopn(ilo,ilop,n,jspin) *&
Daniel Wortmann's avatar
Daniel Wortmann committed
100
                                       conjg(ccof(mp,i,ilop,natomTemp)) *ccof(m ,i,ilo ,natomTemp)
101 102 103 104 105 106 107 108 109 110 111 112
                               ENDDO
                            ENDIF
                         ENDDO
                         n_tmp(m,mp) = n_tmp(m,mp) + c_0
                      ENDDO
                   ENDDO

                ENDIF
             ENDDO
             !
             !  n_mmp should be rotated by D_mm' ; compare force_a21
             !
Daniel Wortmann's avatar
Daniel Wortmann committed
113
             DO it = 1, sym%invarind(natomTemp)
114

Daniel Wortmann's avatar
Daniel Wortmann committed
115 116
                fac = 1.0  /  ( sym%invarind(natomTemp) * atoms%neq(n) )
                is = sym%invarop(natomTemp,it)
117 118 119 120 121 122 123 124 125 126 127
                isi = sym%invtab(is)
                d_tmp(:,:) = cmplx(0.0,0.0)
                DO m = -l,l
                   DO mp = -l,l
                      d_tmp(m,mp) = sym%d_wgn(m,mp,l,isi)
                   ENDDO
                ENDDO
                nr_tmp = matmul( transpose( conjg(d_tmp) ) , n_tmp)
                n1_tmp =  matmul( nr_tmp, d_tmp )
                DO m = -l,l
                   DO mp = -l,l
Daniel Wortmann's avatar
Daniel Wortmann committed
128
                      n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) * fac
129 130 131 132 133 134
                   ENDDO
                ENDDO

             ENDDO

          ENDDO ! sum  over equivalent atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
135 136 137
          i_u = i_u + 1
       END DO
       natom = natom + atoms%neq(n)
138 139 140 141 142 143 144 145 146
    ENDDO     ! loop over atom types

    !     do m=-l,l
    !      write(*,'(14f12.6)') (n_mmp(m,mp),mp=-l,l)
    !     enddo
    !
    RETURN
  END SUBROUTINE n_mat
END MODULE m_nmat