cdnmt.f90 11.6 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 13
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
  SUBROUTINE cdnmt(jspd,atoms,sphhar,llpd, noco,l_fmpl,jsp_start,jsp_end, epar,&
14
                   ello,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 26
    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
    TYPE(t_moments), INTENT(INOUT) :: moments
27 28 29 30 31 32 33

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

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

    CALL timestart("cdnmt")

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

Daniel Wortmann's avatar
Daniel Wortmann committed
67
    !$OMP PARALLEL DEFAULT(none) &
68

69
    !$OMP SHARED(usdus,rho,moments,rho21,qmtl) &
70 71
    !$OMP SHARED(atoms,jsp_start,jsp_end,epar,vr,denCoeffs,sphhar,ello)&
    !$OMP SHARED(orb,noco,l_fmpl,denCoeffsOffdiag,jspd)&
72
    !$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
73 74 75 76 77 78 79 80 81 82
    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
83
    DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
84 85 86 87
       na = 1
       DO i = 1, itype - 1
          na = na + atoms%neq(i)
       ENDDO
88 89 90 91 92 93
       !--->    spherical component
       DO ispin = jsp_start,jsp_end
          DO l = 0,atoms%lmax(itype)
             CALL radfun(l,itype,ispin,epar(l,itype,ispin),vr(1,itype,ispin),atoms,&
                  f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
             DO j = 1,atoms%jri(itype)
94 95 96
                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) )
97 98 99 100 101 102 103 104 105 106 107 108
                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,&
109
               usdus%uloulopn(1,1,itype,ispin),usdus%dulon(1,itype,ispin),&
110
               usdus%uulon(1,itype,ispin),ello(1,itype,ispin),&
111 112 113
               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),&
114 115 116 117 118
               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
119
          qmtt = 0.0
120
          DO l = 0,atoms%lmax(itype)
121
             qmtl(l,ispin,itype) = ( denCoeffs%uu(l,itype,ispin)+denCoeffs%dd(l,itype,ispin)&
122
                  &              *usdus%ddn(l,itype,ispin) )/atoms%neq(itype) + qmtllo(l)
Daniel Wortmann's avatar
Daniel Wortmann committed
123
             qmtt = qmtt + qmtl(l,ispin,itype)
124
          END DO
125
          moments%chmom(itype,ispin) = qmtt
126

127 128 129
          !+soc
          !--->       spherical angular component
          IF (noco%l_soc) THEN
130 131
             CALL orbmom2(atoms,itype,ispin,usdus%ddn(0,itype,ispin),&
                          orb,usdus%uulon(1,itype,ispin),usdus%dulon(1,itype,ispin),&
132
                          usdus%uloulopn(1,1,itype,ispin),moments%clmom(1,itype,ispin))!keep
133 134 135 136 137 138 139 140 141
          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)
142
                      s = denCoeffs%uunmt(llp,lh,itype,ispin)*( &
143
                           f(j,1,l,ispin)*f(j,1,lp,ispin)+ f(j,2,l,ispin)*f(j,2,lp,ispin) )&
144
                           + denCoeffs%ddnmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*g(j,1,lp,ispin)&
145
                           + g(j,2,l,ispin)*g(j,2,lp,ispin) )&
146
                           + denCoeffs%udnmt(llp,lh,itype,ispin)*(f(j,1,l,ispin)*g(j,1,lp,ispin)&
147
                           + f(j,2,l,ispin)*g(j,2,lp,ispin) )&
148
                           + denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
149 150 151 152 153 154 155 156 157 158 159 160
                           + 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)
161
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
162 163 164 165
                  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)
166 167
          ENDDO
          DO ilo = 1, atoms%nlo(itype)
168
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
169 170 171 172
                  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) )/&
173 174
                  atoms%neq(itype)
             DO ilop = 1, atoms%nlo(itype)
175
                moments%qa21(itype) = moments%qa21(itype) + conjg(&
176 177
                     denCoeffsOffdiag%uloulop21(ilo,ilop,itype) *&
                     denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
178 179 180 181 182 183 184 185 186 187 188
             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)
189 190 191 192
                   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) )
193 194 195 196 197 198 199 200 201 202 203
                   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)
204 205 206 207
                         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)&
208 209 210 211 212 213 214 215 216 217 218
                              + 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
219 220 221 222 223 224 225 226 227 228 229
    !$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
230 231
          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
232 233 234
8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
          attributes = ''
          WRITE(attributes(1),'(i0)') itype
235
          WRITE(attributes(2),'(f12.7)') moments%chmom(itype,ispin)
Daniel Wortmann's avatar
Daniel Wortmann committed
236 237 238 239 240 241 242 243 244
          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

245
    CALL timestop("cdnmt")
Daniel Wortmann's avatar
Daniel Wortmann committed
246

247 248 249 250 251 252 253 254 255 256 257 258
    !---> 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