cdnmt.F90 13.9 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
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
11 12 13 14
  !
  !     Parallel MPI version.
  !     Jan. 2019          U.Alekseeva
  !
15 16
  !***********************************************************************
CONTAINS
17
  SUBROUTINE cdnmt(mpi,jspd,atoms,sphhar,noco,jsp_start,jsp_end,enpara,&
18
                   vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
19 20 21
#include "cpp_double.h"

    USE m_constants,ONLY: sfp_const
22 23 24 25
    USE m_rhosphnlo
    USE m_radfun
    USE m_orbmom2
    USE m_types
26
    USE m_xmlOutput
27

28
    IMPLICIT NONE
29

30
    TYPE(t_mpi),     INTENT(IN)    :: mpi
31
    TYPE(t_usdus),   INTENT(INOUT) :: usdus 
32 33 34
    TYPE(t_noco),    INTENT(IN)    :: noco
    TYPE(t_sphhar),  INTENT(IN)    :: sphhar
    TYPE(t_atoms),   INTENT(IN)    :: atoms
35
    TYPE(t_enpara),  INTENT(IN)    :: enpara
36
    TYPE(t_moments), INTENT(INOUT) :: moments
37
    !     ..
38 39
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
40
    !     ..
41
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
42
    REAL, INTENT    (IN) :: vr(atoms%jmtd,atoms%ntype,jspd)
43
    REAL, INTENT (INOUT) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,jspd)
44 45 46
    TYPE (t_orb),              INTENT(IN) :: orb
    TYPE (t_denCoeffs),        INTENT(IN) :: denCoeffs
    TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
47 48 49
    !     ..
    !     .. Local Scalars ..
    INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
Daniel Wortmann's avatar
Daniel Wortmann committed
50
    INTEGER ilo,ilop,i
51
    REAL s,wronk,sumlm,qmtt
52
    COMPLEX rho21
53 54 55
    COMPLEX cs
    !     ..
    !     .. Local Arrays ..
56
    REAL              :: qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
57
    CHARACTER(LEN=20) :: attributes(6)
58 59 60
    !     ..
    !     .. Allocatable Arrays ..
    REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:)
61 62 63 64 65 66 67 68 69 70
    REAL,    ALLOCATABLE :: rho_loc(:,:,:,:)
    REAL,    ALLOCATABLE :: clmom(:,:,:),chmom(:,:)
    COMPLEX, ALLOCATABLE :: qa21(:)

#ifdef CPP_MPI
    INCLUDE 'mpif.h'
    REAL,    ALLOCATABLE :: buf_r(:)
    COMPLEX, ALLOCATABLE :: buf_c(:)
    INTEGER n_buf, ierr
#endif
71
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
72 73 74

    CALL timestart("cdnmt")

75
    IF (mpi%irank==0) THEN
76 77
    ALLOCATE(rho_loc(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,jspd))
    ALLOCATE(chmom(atoms%ntype,jsp_start:jsp_end)) 
78
    IF (noco%l_mperp) THEN
79
       ALLOCATE(qa21(atoms%ntype))
80
       IF (denCoeffsOffdiag%l_fmpl) THEN
81
          rho(:,:,:,3:4) = CMPLX(0.0,0.0)
82 83
       ENDIF
    ENDIF
84 85 86
    IF (noco%l_soc) THEN
       ALLOCATE(clmom(3,atoms%ntype,jsp_start:jsp_end)) 
    ENDIF
87

Daniel Wortmann's avatar
Daniel Wortmann committed
88
    !$OMP PARALLEL DEFAULT(none) &
89
    !$OMP SHARED(mpi,usdus,rho_loc,qmtl) &
90
    !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
91
    !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd,qa21,chmom,clmom)&
92
    !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,cs)
93 94
    rho_loc = 0.0
    chmom = 0.0
Daniel Wortmann's avatar
Daniel Wortmann committed
95
    IF (noco%l_mperp) THEN
96
       qa21 = CMPLX(0.0,0.0)
Daniel Wortmann's avatar
Daniel Wortmann committed
97 98 99 100 101
       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
102 103 104
    IF (noco%l_soc) THEN
       clmom = 0.0
    ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
105 106 107 108

    qmtl = 0
    
    !$OMP DO
109
    DO itype = 1 + mpi%irank, atoms%ntype!, mpi%isize
Daniel Wortmann's avatar
Daniel Wortmann committed
110 111 112 113
       na = 1
       DO i = 1, itype - 1
          na = na + atoms%neq(i)
       ENDDO
114 115 116
       !--->    spherical component
       DO ispin = jsp_start,jsp_end
          DO l = 0,atoms%lmax(itype)
117 118 119 120 121 122
            ! The subroutine radfun is called here only because of f and g.
            ! Variables nodeu,noded,wronk are calculated but never used.
            ! Arrays usdus%[us,dus,uds,duds,ddn] are calculated again, although they
            ! are already there.
            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)  
123
             DO j = 1,atoms%jri(itype)
124 125 126
                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) )
127
                rho_loc(j,0,itype,ispin) = rho_loc(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
128 129 130 131 132 133 134 135 136 137 138
             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,&
139
               usdus%uloulopn(1,1,itype,ispin),usdus%dulon(1,itype,ispin),&   
140
               usdus%uulon(1,itype,ispin),enpara%ello0(1,itype,ispin),&
141 142 143
               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),&
144
               f(1,1,0,ispin),g(1,1,0,ispin),&
145
               rho_loc(:,0:,itype,ispin),qmtllo) !<>rho
146 147 148


          !--->       l-decomposed density for each atom type
Daniel Wortmann's avatar
Daniel Wortmann committed
149
          qmtt = 0.0
150
          DO l = 0,atoms%lmax(itype)
151
             qmtl(l,ispin,itype) = ( denCoeffs%uu(l,itype,ispin)+denCoeffs%dd(l,itype,ispin)&
152
                  &              *usdus%ddn(l,itype,ispin) )/atoms%neq(itype) + qmtllo(l)
Daniel Wortmann's avatar
Daniel Wortmann committed
153
             qmtt = qmtt + qmtl(l,ispin,itype)
154
          END DO
155
          chmom(itype,ispin) = qmtt
156

157 158 159
          !+soc
          !--->       spherical angular component
          IF (noco%l_soc) THEN
160
             CALL orbmom2(atoms,itype,ispin,usdus%ddn(0,itype,ispin),&  
161
                          orb,usdus%uulon(1,itype,ispin),usdus%dulon(1,itype,ispin),&
162
                          usdus%uloulopn(1,1,itype,ispin),clmom(1,itype,ispin))!keep  
163 164 165 166 167 168 169 170 171
          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)
172
                      s = denCoeffs%uunmt(llp,lh,itype,ispin)*( &
173
                           f(j,1,l,ispin)*f(j,1,lp,ispin)+ f(j,2,l,ispin)*f(j,2,lp,ispin) )&
174
                           + denCoeffs%ddnmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*g(j,1,lp,ispin)&
175
                           + g(j,2,l,ispin)*g(j,2,lp,ispin) )&
176
                           + denCoeffs%udnmt(llp,lh,itype,ispin)*(f(j,1,l,ispin)*g(j,1,lp,ispin)&
177
                           + f(j,2,l,ispin)*g(j,2,lp,ispin) )&
178
                           + denCoeffs%dunmt(llp,lh,itype,ispin)*(g(j,1,l,ispin)*f(j,1,lp,ispin)&
179
                           + g(j,2,l,ispin)*f(j,2,lp,ispin) )
180
                      rho_loc(j,lh,itype,ispin) = rho_loc(j,lh,itype,ispin)+ s/atoms%neq(itype)
181 182 183 184 185 186 187 188 189 190
                   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)
191
             qa21(itype) = qa21(itype) + conjg(&
192 193 194 195
                  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)
196 197
          ENDDO
          DO ilo = 1, atoms%nlo(itype)
198
             qa21(itype) = qa21(itype) + conjg(&
199 200 201 202
                  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) )/&
203 204
                  atoms%neq(itype)
             DO ilop = 1, atoms%nlo(itype)
205
                qa21(itype) = qa21(itype) + conjg(&
206 207
                     denCoeffsOffdiag%uloulop21(ilo,ilop,itype) *&
                     denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
208 209 210
             ENDDO
          ENDDO

211
          IF (denCoeffsOffdiag%l_fmpl) THEN
212 213 214 215 216 217 218
             !--->        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)
219 220 221 222
                   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) )
223
                   rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
224 225
                   rho_loc(j,0,itype,3)=rho_loc(j,0,itype,3)+REAL(rho21)
                   rho_loc(j,0,itype,4)=rho_loc(j,0,itype,4)+imag(rho21)
226 227 228 229 230 231 232 233 234 235
                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)
236 237 238 239
                         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)&
240
                              + g(j,2,lp,2)*g(j,2,l,1) )
241
                         rho21=CONJG(cs)/atoms%neq(itype)
242 243
                         rho_loc(j,lh,itype,3)=rho_loc(j,lh,itype,3)+REAL(rho21)
                         rho_loc(j,lh,itype,4)=rho_loc(j,lh,itype,4)+imag(rho21)
244 245 246 247 248
                      ENDDO
                   ENDDO
                ENDDO
             ENDDO

249
          ENDIF ! denCoeffsOffdiag%l_fmpl
250 251 252
       ENDIF ! noco%l_mperp

    ENDDO ! end of loop over atom types
Daniel Wortmann's avatar
Daniel Wortmann committed
253
    !$OMP END DO
254
    DEALLOCATE (f,g)
Daniel Wortmann's avatar
Daniel Wortmann committed
255 256
    !$OMP END PARALLEL

257
#ifdef CPP_MPI
258
    IF (.false.) THEN
259 260 261 262 263
    n_buf=atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*jspd
    ALLOCATE(buf_r(n_buf))
    CALL MPI_REDUCE(rho_loc(1,0,1,1),buf_r,n_buf,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
    IF (mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n_buf,buf_r,1,rho_loc(1,0,1,1),1)
    DEALLOCATE (buf_r)
Daniel Wortmann's avatar
Daniel Wortmann committed
264

265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
    n_buf=atoms%ntype*(jsp_end-jsp_start+1)
    ALLOCATE(buf_r(n_buf))
    CALL MPI_REDUCE(chmom(1,jsp_start),buf_r,n_buf,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
    IF (mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n_buf,buf_r,1,chmom(1,jsp_start),1)
    DEALLOCATE (buf_r)

    IF (noco%l_soc) THEN
       n_buf=3*atoms%ntype*(jsp_end-jsp_start+1)
       ALLOCATE(buf_r(n_buf))
       CALL MPI_REDUCE(clmom(1,1,jsp_start),buf_r,n_buf,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       IF (mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n_buf,buf_r,1,clmom(1,1,jsp_start),1)
       DEALLOCATE (buf_r)
    ENDIF

    IF (noco%l_mperp) THEN
       n_buf= atoms%ntype
       ALLOCATE(buf_c(n_buf))
       CALL MPI_REDUCE(qa21(1),buf_c,n_buf,CPP_MPI_COMPLEX,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       IF (mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n_buf,buf_c,1,qa21(1),1)
       DEALLOCATE (buf_c)
    ENDIF
286
    ENDIF
287 288 289 290 291 292
#endif 
    rho = rho + rho_loc
    moments%chmom(:,jsp_start:jsp_end) = chmom(:,jsp_start:jsp_end)
    IF (noco%l_mperp) moments%qa21 = moments%qa21 + qa21
    IF (noco%l_soc)   moments%clmom(:,:,jsp_start:jsp_end) = clmom(:,:,jsp_start:jsp_end)

293
    !IF (mpi%irank==0) THEN
294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
       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
             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)
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
312
                                   reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
313
          ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
314 315
       ENDDO

316
    ENDIF !(mpi%irank==0)
317
    CALL timestop("cdnmt")
Daniel Wortmann's avatar
Daniel Wortmann committed
318

319
 
320 321
  END SUBROUTINE cdnmt
END MODULE m_cdnmt