rhonmt.f90 4.38 KB
Newer Older
1 2
MODULE m_rhonmt
CONTAINS
3
  SUBROUTINE rhonmt(atoms,sphhar,we,ne,sym,eigVecCoeffs,denCoeffs,ispin)
4 5 6 7 8 9
    !     *************************************************************
    !     subroutine sets up the coefficients of non-sphereical
    !     muffin-tin density                          c.l.fu
    !     *************************************************************
    USE m_gaunt,ONLY:gaunt1
    USE m_types
10
    use m_constants
11
    IMPLICIT NONE
12 13 14 15 16
    TYPE(t_sym),          INTENT(IN)    :: sym
    TYPE(t_sphhar),       INTENT(IN)    :: sphhar
    TYPE(t_atoms),        INTENT(IN)    :: atoms
    TYPE(t_eigVecCoeffs), INTENT(IN)    :: eigVecCoeffs
    TYPE(t_denCoeffs),    INTENT(INOUT) :: denCoeffs
17 18 19

    INTEGER,           INTENT(IN)    :: ne, ispin

20
    REAL,    INTENT(IN) :: we(:)!(nobd)
21

22 23
    !     ..
    !     .. Local Scalars ..
24
    COMPLEX cconst,cil,cmv
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
    REAL coef
    INTEGER  :: jmem,l,lcond,lh,llp,llpmax,lm,lmp,lp,lphi,lplow,lplow0,lv
    INTEGER  :: mp,mv,na,natom,nb,nn,ns,nt,m

    !Initialize private variables in gaunt module before parallel region
    !$      coef = gaunt1(0,0,0,0,0,0,atoms%lmaxd)

    DO ns = 1,sym%nsymt
       !$OMP PARALLEL DO DEFAULT(SHARED) &
       !$OMP&  PRIVATE(lv,jmem,mv,cmv,l,lm,mp,m,llpmax,nt,na,nb,lplow0)&
       !$OMP&  PRIVATE(lphi,lplow,lcond,lp,cil,lmp,llp,coef,cconst&
       !$OMP& ,natom,nn) 
       DO  lh = 1,sphhar%nlh(ns)
          lv = sphhar%llh(lh,ns)
          DO  jmem = 1,sphhar%nmem(lh,ns)
             mv = sphhar%mlh(jmem,lh,ns)
             cmv = conjg(sphhar%clnu(jmem,lh,ns))
             DO  l = 0,atoms%lmaxd
                m_loop: DO m = -l,l
                   lm = l* (l+1) + m
                   mp = m - mv
                   !     -----> set up the lower and upper limit of lp in such a way that
                   !     -----> lp+l+lv is even, lp<=l, and (lp,l,lv) satisfies the
                   !     -----> triangular relation
                   lplow0 = iabs(l-lv)
                   lphi = l - mod(lv,2)
                   lplow = max(lplow0,iabs(mp))
                   lcond = iabs(lphi-lplow)
                   lplow = lplow + mod(lcond,2)
                   IF (lplow.GT.lphi) CYCLE m_loop
                   DO  lp = lplow,lphi,2
56
                      cil = ImagUnit** (l-lp)
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
                      lmp = lp* (lp+1) + mp
                      IF (lmp.GT.lm) CYCLE m_loop
                      llp = (l* (l+1))/2 + lp
                      !     -----> gaunt's coefficient
                      coef = 2.*gaunt1(l,lv,lp,m,mv,mp,atoms%lmaxd)
                      IF (lmp.EQ.lm) coef = coef/2.
                      cconst = coef* (cil*cmv)
                      natom = 0
                      DO  nn = 1,atoms%ntype
                         llpmax = (atoms%lmax(nn)* (atoms%lmax(nn)+3))/2
                         IF (llp.LE.llpmax) THEN 
                            nt = natom
                            DO  na = 1,atoms%neq(nn)
                               nt = nt + 1
                               IF (atoms%ntypsy(nt).EQ.ns) THEN
72 73
                                  DO nb = 1,ne
                                     denCoeffs%uunmt(llp,lh,nn,ispin) = denCoeffs%uunmt(llp,lh,nn,ispin)&
74
                                          +we(nb)*real(cconst*eigVecCoeffs%acof(nb,lm,nt,ispin)*conjg(eigVecCoeffs%acof(nb,lmp,nt,ispin)))
75
                                     denCoeffs%ddnmt(llp,lh,nn,ispin) = denCoeffs%ddnmt(llp,lh,nn,ispin) +&
76
                                          we(nb)*real(cconst*eigVecCoeffs%bcof(nb,lm,nt,ispin)*conjg(eigVecCoeffs%bcof(nb,lmp,nt,ispin)))
77
                                     denCoeffs%udnmt(llp,lh,nn,ispin) = denCoeffs%udnmt(llp,lh,nn,ispin) +&
78
                                          we(nb)*real(cconst*eigVecCoeffs%acof(nb,lm,nt,ispin)*conjg(eigVecCoeffs%bcof(nb,lmp,nt,ispin)))
79
                                     denCoeffs%dunmt(llp,lh,nn,ispin) = denCoeffs%dunmt(llp,lh,nn,ispin) +&
80
                                          we(nb)*real(cconst*eigVecCoeffs%bcof(nb,lm,nt,ispin)*conjg(eigVecCoeffs%acof(nb,lmp,nt,ispin)))
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
                                  ENDDO
                               ENDIF
                            ENDDO
                         ENDIF
                         natom = natom + atoms%neq(nn)
                      ENDDO
                   ENDDO
                ENDDO m_loop
             ENDDO
          ENDDO
       ENDDO
       !$OMP END PARALLEL DO
    ENDDO

  END SUBROUTINE rhonmt
END MODULE m_rhonmt