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
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,sym,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_sym),   INTENT(IN)      :: sym
28
    TYPE(t_enpara),  INTENT(IN)    :: enpara
29
    TYPE(t_moments), INTENT(INOUT) :: moments
30

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

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

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

    CALL timestart("cdnmt")

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

Daniel Wortmann's avatar
Daniel Wortmann committed
67
    !$OMP PARALLEL DEFAULT(none) &
68
    !$OMP SHARED(usdus,rho,moments,qmtl) &
69
    !$OMP SHARED(atoms,sym,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
70
    !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
71
    !$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
72 73 74 75 76 77 78 79 80 81
    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
82
    DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
83 84 85 86
       na = 1
       DO i = 1, itype - 1
          na = na + atoms%neq(i)
       ENDDO
87 88 89
       !--->    spherical component
       DO ispin = jsp_start,jsp_end
          DO l = 0,atoms%lmax(itype)
90 91
             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)
92
             DO j = 1,atoms%jri(itype)
93 94 95
                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) )
96
                rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
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

107
          CALL rhosphnlo(itype,atoms,sphhar,sym,&
108
               usdus%uloulopn(1,1,itype,ispin),usdus%dulon(1,itype,ispin),&
109
               usdus%uulon(1,itype,ispin),enpara%ello0(1,itype,ispin),&
110 111 112
               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),&
113
               f(1,1,0,ispin),g(1,1,0,ispin),&
114
               rho(:,0:,itype,ispin),qmtllo)
115 116 117


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

126 127 128
          !+soc
          !--->       spherical angular component
          IF (noco%l_soc) THEN
129
             CALL orbmom2(atoms,itype,ispin,usdus%ddn(0,itype,ispin),&
130
                          orb,usdus%uulon(1,itype,ispin),usdus%dulon(1,itype,ispin),&
131
                          usdus%uloulopn(1,1,itype,ispin),moments%clmom(1,itype,ispin))!keep
132 133 134
          ENDIF
          !-soc
          !--->       non-spherical components
135
          nd = sym%ntypsy(na)
136 137 138 139 140
          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)
141
                      s = denCoeffs%uunmt(llp,lh,itype,ispin)*( &
142
                           f(j,1,l,ispin)*f(j,1,lp,ispin)+ f(j,2,l,ispin)*f(j,2,lp,ispin) )&
143
                           + denCoeffs%ddnmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*g(j,1,lp,ispin)&
144
                           + g(j,2,l,ispin)*g(j,2,lp,ispin) )&
145
                           + denCoeffs%udnmt(llp,lh,itype,ispin)*(f(j,1,l,ispin)*g(j,1,lp,ispin)&
146
                           + f(j,2,l,ispin)*g(j,2,lp,ispin) )&
147
                           + denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
148
                           + g(j,2,l,ispin)*f(j,2,lp,ispin) )
149
                      rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
150 151 152 153 154 155 156 157 158 159
                   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)
160
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
161 162 163 164
                  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)
165 166
          ENDDO
          DO ilo = 1, atoms%nlo(itype)
167
             moments%qa21(itype) = moments%qa21(itype) + conjg(&
168 169 170 171
                  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) )/&
172 173
                  atoms%neq(itype)
             DO ilop = 1, atoms%nlo(itype)
174
                moments%qa21(itype) = moments%qa21(itype) + conjg(&
175 176
                     denCoeffsOffdiag%uloulop21(ilo,ilop,itype) *&
                     denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
177 178 179
             ENDDO
          ENDDO

180
          IF (denCoeffsOffdiag%l_fmpl) THEN
181 182 183 184 185 186 187
             !--->        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)
188 189 190 191
                   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) )
192
                   !rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
193
                   rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
194 195
                   rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
                   rho(j,0,itype,4)=rho(j,0,itype,4)+imag(rho21)
196 197 198 199
                ENDDO
             ENDDO

             !--->        non-spherical components
200
             nd = sym%ntypsy(na)
201 202 203 204 205
             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)
206 207 208 209
                         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)&
210
                              + g(j,2,lp,2)*g(j,2,l,1) )
211
                         !rho21(j,lh,itype)= rho21(j,lh,itype)+ CONJG(cs)/atoms%neq(itype)
212
                         rho21=CONJG(cs)/atoms%neq(itype)
213 214
                         rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
                         rho(j,lh,itype,4)=rho(j,lh,itype,4)+imag(rho21)
215 216 217 218 219
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO

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

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

228 229 230
    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
231

232 233 234 235 236 237 238 239 240 241 242 243
    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
244 245
                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
       ENDDO
246
    ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
247

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

251
 
252 253
  END SUBROUTINE cdnmt
END MODULE m_cdnmt