n_mat.f90 5.7 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
  !     ************************************************************
CONTAINS
  SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we, acof,bcof,ccof, n_mmp)
    !

    USE m_types
Gregor Michalicek's avatar
Gregor Michalicek committed
21
    USE m_constants
22 23 24 25 26 27 28 29 30 31 32 33
    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
34
    COMPLEX, INTENT (IN) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
Gregor Michalicek's avatar
Gregor Michalicek committed
35
    COMPLEX, INTENT (INOUT) :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
36 37 38
    !     ..
    !     .. Local Scalars ..
    COMPLEX c_0
Daniel Wortmann's avatar
Daniel Wortmann committed
39
    INTEGER i,j,k,l ,mp,n,it,is,isi,natom,natomTemp,n_ldau,lp,m,i_u
40 41 42 43 44 45 46 47 48 49 50
    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
51
    i_u = 1
52
    DO n = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
53 54 55 56 57
       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
58 59
          ll1 = (l+1)*l 
          DO nn = 1, atoms%neq(n)
Daniel Wortmann's avatar
Daniel Wortmann committed
60
             natomTemp = natomTemp + 1
61 62 63 64 65 66 67 68 69 70 71
             !
             !  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
72 73
                           conjg(bcof(i,lmp,natomTemp))*bcof(i,lm,natomTemp) +&
                           conjg(acof(i,lmp,natomTemp))*acof(i,lm,natomTemp) )
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
                   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
91 92
                                 conjg(acof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
                                 conjg(ccof(mp,i,ilo,natomTemp))*acof(i,lm,natomTemp) )&
93
                                 + usdus%dulon(ilo,n,jspin) * (&
Daniel Wortmann's avatar
Daniel Wortmann committed
94 95
                                 conjg(bcof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
                                 conjg(ccof(mp,i,ilo,natomTemp))*bcof(i,lm,natomTemp)))
96 97 98 99 100
                         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
101
                                       conjg(ccof(mp,i,ilop,natomTemp)) *ccof(m ,i,ilo ,natomTemp)
102 103 104 105 106 107 108 109 110 111 112 113
                               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
114
             DO it = 1, sym%invarind(natomTemp)
115

Daniel Wortmann's avatar
Daniel Wortmann committed
116 117
                fac = 1.0  /  ( sym%invarind(natomTemp) * atoms%neq(n) )
                is = sym%invarop(natomTemp,it)
118 119 120 121 122 123 124 125 126 127 128
                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
129
                      n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) * fac
130 131 132 133 134 135
                   ENDDO
                ENDDO

             ENDDO

          ENDDO ! sum  over equivalent atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
136 137 138
          i_u = i_u + 1
       END DO
       natom = natom + atoms%neq(n)
139 140 141 142 143 144 145 146 147
    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