cdnmt.f90 11.5 KB
Newer Older
1 2 3 4 5
! 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.
!--------------------------------------------------------------------------------

6 7 8 9 10 11 12
MODULE m_cdnmt
  !***********************************************************************
  !     This subroutine calculates the spherical and non-spherical charge-
  !     density and the orbital moment inside the muffin-tin spheres.
  !     Philipp Kurz 2000-02-03
  !***********************************************************************
CONTAINS
13 14
  SUBROUTINE cdnmt(jspd,atoms,sphhar,noco,l_fmpl,jsp_start,jsp_end,enpara,&
                   vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
15 16 17 18 19
    use m_constants,only: sfp_const
    USE m_rhosphnlo
    USE m_radfun
    USE m_orbmom2
    USE m_types
20
    USE m_xmlOutput
21
    IMPLICIT NONE
22 23 24 25
    TYPE(t_usdus),   INTENT(INOUT) :: usdus !in fact only the lo part is intent(in)
    TYPE(t_noco),    INTENT(IN)    :: noco
    TYPE(t_sphhar),  INTENT(IN)    :: sphhar
    TYPE(t_atoms),   INTENT(IN)    :: atoms
26
    TYPE(t_enpara),  INTENT(IN)    :: enpara
27
    TYPE(t_moments), INTENT(INOUT) :: moments
28 29 30 31 32 33

    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
    LOGICAL, INTENT (IN) :: l_fmpl
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
34 35
    REAL, INTENT    (IN) :: vr(atoms%jmtd,atoms%ntype,jspd)
    REAL, INTENT (INOUT) :: rho(:,0:,:,:)!(toms%jmtd,0:sphhar%nlhd,atoms%ntype,jspd)
36 37 38
    TYPE (t_orb),              INTENT(IN) :: orb
    TYPE (t_denCoeffs),        INTENT(IN) :: denCoeffs
    TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
39 40 41
    !     ..
    !     .. Local Scalars ..
    INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
Daniel Wortmann's avatar
Daniel Wortmann committed
42
    INTEGER ilo,ilop,i
43 44 45 46
    REAL s,wronk,sumlm,qmtt
    COMPLEX cs
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
47
    REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
48
    CHARACTER(LEN=20) :: attributes(6)
49 50 51 52 53 54

    !     ..
    !     .. Allocatable Arrays ..
    REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:)
    COMPLEX, ALLOCATABLE :: rho21(:,:,:)
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
55 56 57

    CALL timestart("cdnmt")

58 59
    IF (noco%l_mperp) THEN
       IF (l_fmpl) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
60
          ALLOCATE ( rho21(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) )
61 62 63 64
          rho21(:,:,:) = cmplx(0.0,0.0)
       ENDIF
    ENDIF

Daniel Wortmann's avatar
Daniel Wortmann committed
65
    !$OMP PARALLEL DEFAULT(none) &
66

67
    !$OMP SHARED(usdus,rho,moments,rho21,qmtl) &
68
    !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
69
    !$OMP SHARED(orb,noco,l_fmpl,denCoeffsOffdiag,jspd)&
70
    !$OMP PRIVATE(itype,na,ispin,l,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,cs)
Daniel Wortmann's avatar
Daniel Wortmann committed
71 72 73 74 75 76 77 78 79 80
    IF (noco%l_mperp) THEN
       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
    ELSE
       ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
       ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
    ENDIF

    qmtl = 0
    
    !$OMP DO
81
    DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
82 83 84 85
       na = 1
       DO i = 1, itype - 1
          na = na + atoms%neq(i)
       ENDDO
86 87 88
       !--->    spherical component
       DO ispin = jsp_start,jsp_end
          DO l = 0,atoms%lmax(itype)
89
             CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
90 91
                  f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
             DO j = 1,atoms%jri(itype)
92 93 94
                s = denCoeffs%uu(l,itype,ispin)*( f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin) )&
                     +   denCoeffs%dd(l,itype,ispin)*( g(j,1,l,ispin)*g(j,1,l,ispin)+g(j,2,l,ispin)*g(j,2,l,ispin) )&
                     + 2*denCoeffs%du(l,itype,ispin)*( f(j,1,l,ispin)*g(j,1,l,ispin)+f(j,2,l,ispin)*g(j,2,l,ispin) )
95 96 97 98 99 100 101 102 103 104 105 106
                rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
             ENDDO
          ENDDO

          !--->       add the contribution of the local orbitals and flapw - lo
          !--->       cross-terms to rho, qmtl. the latter are stored in
          !--->       qmtllo. initialize qmtllo
          DO l = 0,atoms%lmaxd
             qmtllo(l) = 0.0
          END DO

          CALL rhosphnlo(itype,atoms,sphhar,&
107
               usdus%uloulopn(1,1,itype,ispin),usdus%dulon(1,itype,ispin),&
108
               usdus%uulon(1,itype,ispin),enpara%ello0(1,itype,ispin),&
109 110 111
               vr(1,itype,ispin),denCoeffs%aclo(1,itype,ispin),denCoeffs%bclo(1,itype,ispin),&
               denCoeffs%cclo(1,1,itype,ispin),denCoeffs%acnmt(0,1,1,itype,ispin),&
               denCoeffs%bcnmt(0,1,1,itype,ispin),denCoeffs%ccnmt(1,1,1,itype,ispin),&
112 113 114 115 116
               f(1,1,0,ispin),g(1,1,0,ispin),&
               rho(:,0:,itype,ispin),qmtllo)


          !--->       l-decomposed density for each atom type
Daniel Wortmann's avatar
Daniel Wortmann committed
117
          qmtt = 0.0
118
          DO l = 0,atoms%lmax(itype)
119
             qmtl(l,ispin,itype) = ( denCoeffs%uu(l,itype,ispin)+denCoeffs%dd(l,itype,ispin)&
120
                  &              *usdus%ddn(l,itype,ispin) )/atoms%neq(itype) + qmtllo(l)
Daniel Wortmann's avatar
Daniel Wortmann committed
121
             qmtt = qmtt + qmtl(l,ispin,itype)
122
          END DO
123
          moments%chmom(itype,ispin) = qmtt
124

125 126 127
          !+soc
          !--->       spherical angular component
          IF (noco%l_soc) THEN
128 129
             CALL orbmom2(atoms,itype,ispin,usdus%ddn(0,itype,ispin),&
                          orb,usdus%uulon(1,itype,ispin),usdus%dulon(1,itype,ispin),&
130
                          usdus%uloulopn(1,1,itype,ispin),moments%clmom(1,itype,ispin))!keep
131 132 133 134 135 136 137 138 139
          ENDIF
          !-soc
          !--->       non-spherical components
          nd = atoms%ntypsy(na)
          DO lh = 1,sphhar%nlh(nd)
             DO l = 0,atoms%lmax(itype)
                DO lp = 0,l
                   llp = (l* (l+1))/2 + lp
                   DO j = 1,atoms%jri(itype)
140
                      s = denCoeffs%uunmt(llp,lh,itype,ispin)*( &
141
                           f(j,1,l,ispin)*f(j,1,lp,ispin)+ f(j,2,l,ispin)*f(j,2,lp,ispin) )&
142
                           + denCoeffs%ddnmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*g(j,1,lp,ispin)&
143
                           + g(j,2,l,ispin)*g(j,2,lp,ispin) )&
144
                           + denCoeffs%udnmt(llp,lh,itype,ispin)*(f(j,1,l,ispin)*g(j,1,lp,ispin)&
145
                           + f(j,2,l,ispin)*g(j,2,lp,ispin) )&
146
                           + denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
147 148 149 150 151 152 153 154 155 156 157 158
                           + g(j,2,l,ispin)*f(j,2,lp,ispin) )
                      rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
                   ENDDO
                ENDDO
             ENDDO
          ENDDO
       ENDDO ! end of spin loop (ispin = jsp_start,jsp_end)

       IF (noco%l_mperp) THEN

          !--->      calculate off-diagonal integrated density
          DO l = 0,atoms%lmax(itype)
159
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
160 161 162 163
                  denCoeffsOffdiag%uu21(l,itype) * denCoeffsOffdiag%uu21n(l,itype) +&
                  denCoeffsOffdiag%ud21(l,itype) * denCoeffsOffdiag%ud21n(l,itype) +&
                  denCoeffsOffdiag%du21(l,itype) * denCoeffsOffdiag%du21n(l,itype) +&
                  denCoeffsOffdiag%dd21(l,itype) * denCoeffsOffdiag%dd21n(l,itype) )/atoms%neq(itype)
164 165
          ENDDO
          DO ilo = 1, atoms%nlo(itype)
166
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
167 168 169 170
                  denCoeffsOffdiag%ulou21(ilo,itype) * denCoeffsOffdiag%ulou21n(ilo,itype) +&
                  denCoeffsOffdiag%ulod21(ilo,itype) * denCoeffsOffdiag%ulod21n(ilo,itype) +&
                  denCoeffsOffdiag%uulo21(ilo,itype) * denCoeffsOffdiag%uulo21n(ilo,itype) +&
                  denCoeffsOffdiag%dulo21(ilo,itype) * denCoeffsOffdiag%dulo21n(ilo,itype) )/&
171 172
                  atoms%neq(itype)
             DO ilop = 1, atoms%nlo(itype)
173
                moments%qa21(itype) = moments%qa21(itype) + conjg(&
174 175
                     denCoeffsOffdiag%uloulop21(ilo,ilop,itype) *&
                     denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
176 177 178 179 180 181 182 183 184 185 186
             ENDDO
          ENDDO

          IF (l_fmpl) THEN
             !--->        the following part can be used to calculate the full magnet.
             !--->        density without the atomic sphere approximation for the
             !--->        magnet. density, e.g. for plotting.
             !--->        calculate off-diagonal part of the density matrix
             !--->        spherical component
             DO l = 0,atoms%lmax(itype)
                DO j = 1,atoms%jri(itype)
187 188 189 190
                   cs = denCoeffsOffdiag%uu21(l,itype)*( f(j,1,l,2)*f(j,1,l,1) +f(j,2,l,2)*f(j,2,l,1) )&
                        + denCoeffsOffdiag%ud21(l,itype)*( f(j,1,l,2)*g(j,1,l,1) +f(j,2,l,2)*g(j,2,l,1) )&
                        + denCoeffsOffdiag%du21(l,itype)*( g(j,1,l,2)*f(j,1,l,1) +g(j,2,l,2)*f(j,2,l,1) )&
                        + denCoeffsOffdiag%dd21(l,itype)*( g(j,1,l,2)*g(j,1,l,1) +g(j,2,l,2)*g(j,2,l,1) )
191 192 193 194 195 196 197 198 199 200 201
                   rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
                ENDDO
             ENDDO

             !--->        non-spherical components
             nd = atoms%ntypsy(na)
             DO lh = 1,sphhar%nlh(nd)
                DO l = 0,atoms%lmax(itype)
                   DO lp = 0,atoms%lmax(itype)
                      llp = lp*(atoms%lmax(itype)+1)+l+1
                      DO j = 1,atoms%jri(itype)
202 203 204 205
                         cs = denCoeffsOffdiag%uunmt21(llp,lh,itype)*(f(j,1,lp,2)*f(j,1,l,1)&
                              + f(j,2,lp,2)*f(j,2,l,1) )+ denCoeffsOffdiag%udnmt21(llp,lh,itype)*(f(j,1,lp,2)*g(j,1,l,1)&
                              + f(j,2,lp,2)*g(j,2,l,1) )+ denCoeffsOffdiag%dunmt21(llp,lh,itype)*(g(j,1,lp,2)*f(j,1,l,1)&
                              + g(j,2,lp,2)*f(j,2,l,1) )+ denCoeffsOffdiag%ddnmt21(llp,lh,itype)*(g(j,1,lp,2)*g(j,1,l,1)&
206 207 208 209 210 211 212 213 214 215 216
                              + g(j,2,lp,2)*g(j,2,l,1) )
                         rho21(j,lh,itype)= rho21(j,lh,itype)+ conjg(cs)/atoms%neq(itype)
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO

          ENDIF ! l_fmpl
       ENDIF ! noco%l_mperp

    ENDDO ! end of loop over atom types
Daniel Wortmann's avatar
Daniel Wortmann committed
217 218 219 220 221 222 223 224 225 226 227
    !$OMP END DO
    DEALLOCATE ( f,g)
    !$OMP END PARALLEL

    WRITE (6,FMT=8000)
    WRITE (16,FMT=8000)
8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
         &     t33,'d',t42,'f',t51,'total')

    DO itype = 1,atoms%ntype
       DO ispin = jsp_start,jsp_end
228 229
          WRITE ( 6,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),moments%chmom(itype,ispin)
          WRITE (16,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),moments%chmom(itype,ispin)
Daniel Wortmann's avatar
Daniel Wortmann committed
230 231 232
8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
          attributes = ''
          WRITE(attributes(1),'(i0)') itype
233
          WRITE(attributes(2),'(f12.7)') moments%chmom(itype,ispin)
Daniel Wortmann's avatar
Daniel Wortmann committed
234 235 236 237 238 239 240 241 242
          WRITE(attributes(3),'(f12.7)') qmtl(0,ispin,itype)
          WRITE(attributes(4),'(f12.7)') qmtl(1,ispin,itype)
          WRITE(attributes(5),'(f12.7)') qmtl(2,ispin,itype)
          WRITE(attributes(6),'(f12.7)') qmtl(3,ispin,itype)
          CALL writeXMLElementForm('mtCharge',(/'atomType','total   ','s       ','p       ','d       ','f       '/),attributes,&
                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
       ENDDO
    ENDDO

243
    CALL timestop("cdnmt")
Daniel Wortmann's avatar
Daniel Wortmann committed
244

245 246 247 248 249 250 251 252 253 254 255 256
    !---> for testing: to plot the offdiag. part of the density matrix it
    !---> is written to the file rhomt21. This file can read in pldngen.
    IF (l_fmpl) THEN
       OPEN (26,file='rhomt21',form='unformatted',status='unknown')
       WRITE (26) rho21
       CLOSE (26)
       DEALLOCATE ( rho21 )
    ENDIF
    !---> end of test output

  END SUBROUTINE cdnmt
END MODULE m_cdnmt