n_mat.f90 5.78 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
  !     ************************************************************
CONTAINS
17
  SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we,eigVecCoeffs,n_mmp)
18 19 20
    !

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

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

             ENDDO

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