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
  SUBROUTINE cdnmt(mpi,jspd,atoms,sphhar,noco,jsp_start,jsp_end,enpara,&
14
                   vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
15
    use m_constants,only: sfp_const
16 17 18 19
    USE m_rhosphnlo
    USE m_radfun
    USE m_orbmom2
    USE m_types
20
    USE m_xmlOutput
21
    IMPLICIT NONE
22
    TYPE(t_mpi),     INTENT(IN)    :: mpi
23
    TYPE(t_usdus),   INTENT(INOUT) :: usdus !in fact only the lo part is intent(in)
24 25 26
    TYPE(t_noco),    INTENT(IN)    :: noco
    TYPE(t_sphhar),  INTENT(IN)    :: sphhar
    TYPE(t_atoms),   INTENT(IN)    :: atoms
27
    TYPE(t_enpara),  INTENT(IN)    :: enpara
28
    TYPE(t_moments), INTENT(INOUT) :: moments
29

30 31
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
32

33
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
34
    REAL, INTENT    (IN) :: vr(atoms%jmtd,atoms%ntype,jspd)
35
    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 ..
47
    REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
48
    CHARACTER(LEN=20) :: attributes(6)
49

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

    CALL timestart("cdnmt")

58
   IF (mpi%irank==0) THEN
59
    IF (noco%l_mperp) THEN
60
       IF (denCoeffsOffdiag%l_fmpl) THEN
61
          !ALLOCATE ( rho21(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) )
62
          rho(:,:,:,3:4) = CMPLX(0.0,0.0)
63 64 65
       ENDIF
    ENDIF

Daniel Wortmann's avatar
Daniel Wortmann committed
66
    !$OMP PARALLEL DEFAULT(none) &
67
    !$OMP SHARED(usdus,rho,moments,qmtl) &
68
    !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
69
    !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
70
    !$OMP PRIVATE(itype,na,ispin,l,rho21,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 90
             CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
                  f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
91
             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
                rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
96 97 98 99 100 101 102 103 104 105 106
             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
               f(1,1,0,ispin),g(1,1,0,ispin),&
113
               rho(:,0:,itype,ispin),qmtllo)
114 115 116


          !--->       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
             CALL orbmom2(atoms,itype,ispin,usdus%ddn(0,itype,ispin),&
129
                          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
                           + g(j,2,l,ispin)*f(j,2,lp,ispin) )
148
                      rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
149 150 151 152 153 154 155 156 157 158
                   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
             ENDDO
          ENDDO

179
          IF (denCoeffsOffdiag%l_fmpl) THEN
180 181 182 183 184 185 186
             !--->        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
                   !rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
192
                   rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
193 194
                   rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
                   rho(j,0,itype,4)=rho(j,0,itype,4)+imag(rho21)
195 196 197 198 199 200 201 202 203 204
                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)
205 206 207 208
                         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)&
209
                              + g(j,2,lp,2)*g(j,2,l,1) )
210
                         !rho21(j,lh,itype)= rho21(j,lh,itype)+ CONJG(cs)/atoms%neq(itype)
211
                         rho21=CONJG(cs)/atoms%neq(itype)
212 213
                         rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
                         rho(j,lh,itype,4)=rho(j,lh,itype,4)+imag(rho21)
214 215 216 217 218
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO

219
          ENDIF ! denCoeffsOffdiag%l_fmpl
220 221 222
       ENDIF ! noco%l_mperp

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

227 228 229
    WRITE (6,FMT=8000)
8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
         &     t33,'d',t42,'f',t51,'total')
Daniel Wortmann's avatar
Daniel Wortmann committed
230

231 232 233 234 235 236 237 238 239 240 241 242
    DO itype = 1,atoms%ntype
       DO ispin = jsp_start,jsp_end
          WRITE ( 6,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),moments%chmom(itype,ispin)
8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
          attributes = ''
          WRITE(attributes(1),'(i0)') itype
          WRITE(attributes(2),'(f12.7)') moments%chmom(itype,ispin)
          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,&
Daniel Wortmann's avatar
Daniel Wortmann committed
243 244
                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
       ENDDO
245
    ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
246

247
   ENDIF !(mpi%irank==0) THEN
248
    CALL timestop("cdnmt")
Daniel Wortmann's avatar
Daniel Wortmann committed
249

250
 
251 252
  END SUBROUTINE cdnmt
END MODULE m_cdnmt