LCOV - code coverage report
Current view: top level - cdn_mt - cdnmt.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 56 82 68.3 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       2             : ! This file is part of FLEUR and available as free software under the conditions
       3             : ! of the MIT license as expressed in the LICENSE file in more detail.
       4             : !--------------------------------------------------------------------------------
       5             : 
       6             : MODULE m_cdnmt
       7             :   !***********************************************************************
       8             :   !     This subroutine calculates the spherical and non-spherical charge-
       9             :   !     density and the orbital moment inside the muffin-tin spheres.
      10             :   !     Philipp Kurz 2000-02-03
      11             :   !***********************************************************************
      12             : CONTAINS
      13         608 :   SUBROUTINE cdnmt(mpi,jspd,atoms,sphhar,noco,jsp_start,jsp_end,enpara,&
      14         608 :                    vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
      15             :     use m_constants,only: sfp_const
      16             :     USE m_rhosphnlo
      17             :     USE m_radfun
      18             :     USE m_orbmom2
      19             :     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             :     TYPE(t_noco),    INTENT(IN)    :: noco
      25             :     TYPE(t_sphhar),  INTENT(IN)    :: sphhar
      26             :     TYPE(t_atoms),   INTENT(IN)    :: atoms
      27             :     TYPE(t_enpara),  INTENT(IN)    :: enpara
      28             :     TYPE(t_moments), INTENT(INOUT) :: moments
      29             : 
      30             :     !     .. Scalar Arguments ..
      31             :     INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
      32             : 
      33             :     !     .. Array Arguments ..
      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             :     TYPE (t_orb),              INTENT(IN) :: orb
      37             :     TYPE (t_denCoeffs),        INTENT(IN) :: denCoeffs
      38             :     TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
      39             :     !     ..
      40             :     !     .. Local Scalars ..
      41             :     INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
      42             :     INTEGER ilo,ilop,i
      43             :     REAL s,wronk,sumlm,qmtt
      44             :     COMPLEX cs
      45             :     !     ..
      46             :     !     .. Local Arrays ..
      47        1216 :     REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
      48             :     CHARACTER(LEN=20) :: attributes(6)
      49             : 
      50             :     !     ..
      51             :     !     .. Allocatable Arrays ..
      52         608 :     REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:)
      53             :     COMPLEX :: rho21
      54             :     !
      55             : 
      56         608 :     CALL timestart("cdnmt")
      57             : 
      58         608 :    IF (mpi%irank==0) THEN
      59         304 :     IF (noco%l_mperp) THEN
      60           0 :        IF (denCoeffsOffdiag%l_fmpl) THEN
      61             :           !ALLOCATE ( rho21(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) )
      62           0 :           rho(:,:,:,3:4) = CMPLX(0.0,0.0)
      63             :        ENDIF
      64             :     ENDIF
      65             : 
      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         608 :     !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,cs)
      71         304 :     IF (noco%l_mperp) THEN
      72           0 :        ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
      73             :     ELSE
      74         304 :        ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
      75         304 :        ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
      76             :     ENDIF
      77             : 
      78        1468 :     qmtl = 0
      79             :     
      80         304 :     !$OMP DO
      81             :     DO itype = 1,atoms%ntype
      82         582 :        na = 1
      83         910 :        DO i = 1, itype - 1
      84         910 :           na = na + atoms%neq(i)
      85             :        ENDDO
      86             :        !--->    spherical component
      87        1164 :        DO ispin = jsp_start,jsp_end
      88        5904 :           DO l = 0,atoms%lmax(itype)
      89             :              CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
      90        5322 :                   f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
      91     3224270 :              DO j = 1,atoms%jri(itype)
      92             :                 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) )&
      93             :                      +   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) )&
      94     3218366 :                      + 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     3223688 :                 rho(j,0,itype,ispin) = rho(j,0,itype,ispin)+ s/(atoms%neq(itype)*sfp_const)
      96             :              ENDDO
      97             :           ENDDO
      98             : 
      99             :           !--->       add the contribution of the local orbitals and flapw - lo
     100             :           !--->       cross-terms to rho, qmtl. the latter are stored in
     101             :           !--->       qmtllo. initialize qmtllo
     102        5968 :           DO l = 0,atoms%lmaxd
     103        5968 :              qmtllo(l) = 0.0
     104             :           END DO
     105             : 
     106             :           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             :                vr(1,itype,ispin),denCoeffs%aclo(1,itype,ispin),denCoeffs%bclo(1,itype,ispin),&
     110             :                denCoeffs%cclo(1,1,itype,ispin),denCoeffs%acnmt(0,1,1,itype,ispin),&
     111             :                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         582 :                rho(:,0:,itype,ispin),qmtllo)
     114             : 
     115             : 
     116             :           !--->       l-decomposed density for each atom type
     117         582 :           qmtt = 0.0
     118        5904 :           DO l = 0,atoms%lmax(itype)
     119             :              qmtl(l,ispin,itype) = ( denCoeffs%uu(l,itype,ispin)+denCoeffs%dd(l,itype,ispin)&
     120        5322 :                   &              *usdus%ddn(l,itype,ispin) )/atoms%neq(itype) + qmtllo(l)
     121        5904 :              qmtt = qmtt + qmtl(l,ispin,itype)
     122             :           END DO
     123         582 :           moments%chmom(itype,ispin) = qmtt
     124             : 
     125             :           !+soc
     126             :           !--->       spherical angular component
     127         582 :           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         184 :                           usdus%uloulopn(1,1,itype,ispin),moments%clmom(1,itype,ispin))!keep
     131             :           ENDIF
     132             :           !-soc
     133             :           !--->       non-spherical components
     134         582 :           nd = atoms%ntypsy(na)
     135        9966 :           DO lh = 1,sphhar%nlh(nd)
     136       90702 :              DO l = 0,atoms%lmax(itype)
     137      930192 :                 DO lp = 0,l
     138      420036 :                    llp = (l* (l+1))/2 + lp
     139   258950982 :                    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   258449628 :                            + g(j,2,l,ispin)*f(j,2,lp,ispin) )
     148   258869664 :                       rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ s/atoms%neq(itype)
     149             :                    ENDDO
     150             :                 ENDDO
     151             :              ENDDO
     152             :           ENDDO
     153             :        ENDDO ! end of spin loop (ispin = jsp_start,jsp_end)
     154             : 
     155         582 :        IF (noco%l_mperp) THEN
     156             : 
     157             :           !--->      calculate off-diagonal integrated density
     158           0 :           DO l = 0,atoms%lmax(itype)
     159             :              moments%qa21(itype) = moments%qa21(itype) + conjg(&
     160             :                   denCoeffsOffdiag%uu21(l,itype) * denCoeffsOffdiag%uu21n(l,itype) +&
     161             :                   denCoeffsOffdiag%ud21(l,itype) * denCoeffsOffdiag%ud21n(l,itype) +&
     162             :                   denCoeffsOffdiag%du21(l,itype) * denCoeffsOffdiag%du21n(l,itype) +&
     163           0 :                   denCoeffsOffdiag%dd21(l,itype) * denCoeffsOffdiag%dd21n(l,itype) )/atoms%neq(itype)
     164             :           ENDDO
     165           0 :           DO ilo = 1, atoms%nlo(itype)
     166             :              moments%qa21(itype) = moments%qa21(itype) + conjg(&
     167             :                   denCoeffsOffdiag%ulou21(ilo,itype) * denCoeffsOffdiag%ulou21n(ilo,itype) +&
     168             :                   denCoeffsOffdiag%ulod21(ilo,itype) * denCoeffsOffdiag%ulod21n(ilo,itype) +&
     169             :                   denCoeffsOffdiag%uulo21(ilo,itype) * denCoeffsOffdiag%uulo21n(ilo,itype) +&
     170             :                   denCoeffsOffdiag%dulo21(ilo,itype) * denCoeffsOffdiag%dulo21n(ilo,itype) )/&
     171           0 :                   atoms%neq(itype)
     172           0 :              DO ilop = 1, atoms%nlo(itype)
     173             :                 moments%qa21(itype) = moments%qa21(itype) + conjg(&
     174             :                      denCoeffsOffdiag%uloulop21(ilo,ilop,itype) *&
     175           0 :                      denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
     176             :              ENDDO
     177             :           ENDDO
     178             : 
     179           0 :           IF (denCoeffsOffdiag%l_fmpl) THEN
     180             :              !--->        the following part can be used to calculate the full magnet.
     181             :              !--->        density without the atomic sphere approximation for the
     182             :              !--->        magnet. density, e.g. for plotting.
     183             :              !--->        calculate off-diagonal part of the density matrix
     184             :              !--->        spherical component
     185           0 :              DO l = 0,atoms%lmax(itype)
     186           0 :                 DO j = 1,atoms%jri(itype)
     187             :                    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) )&
     188             :                         + denCoeffsOffdiag%ud21(l,itype)*( f(j,1,l,2)*g(j,1,l,1) +f(j,2,l,2)*g(j,2,l,1) )&
     189             :                         + denCoeffsOffdiag%du21(l,itype)*( g(j,1,l,2)*f(j,1,l,1) +g(j,2,l,2)*f(j,2,l,1) )&
     190           0 :                         + 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           0 :                    rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
     193           0 :                    rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
     194           0 :                    rho(j,0,itype,4)=rho(j,0,itype,4)+imag(rho21)
     195             :                 ENDDO
     196             :              ENDDO
     197             : 
     198             :              !--->        non-spherical components
     199           0 :              nd = atoms%ntypsy(na)
     200           0 :              DO lh = 1,sphhar%nlh(nd)
     201           0 :                 DO l = 0,atoms%lmax(itype)
     202           0 :                    DO lp = 0,atoms%lmax(itype)
     203           0 :                       llp = lp*(atoms%lmax(itype)+1)+l+1
     204           0 :                       DO j = 1,atoms%jri(itype)
     205             :                          cs = denCoeffsOffdiag%uunmt21(llp,lh,itype)*(f(j,1,lp,2)*f(j,1,l,1)&
     206             :                               + 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)&
     207             :                               + 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)&
     208             :                               + 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           0 :                               + 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           0 :                          rho21=CONJG(cs)/atoms%neq(itype)
     212           0 :                          rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
     213           0 :                          rho(j,lh,itype,4)=rho(j,lh,itype,4)+imag(rho21)
     214             :                       ENDDO
     215             :                    ENDDO
     216             :                 ENDDO
     217             :              ENDDO
     218             : 
     219             :           ENDIF ! denCoeffsOffdiag%l_fmpl
     220             :        ENDIF ! noco%l_mperp
     221             : 
     222             :     ENDDO ! end of loop over atom types
     223             :     !$OMP END DO
     224         304 :     DEALLOCATE ( f,g)
     225             :     !$OMP END PARALLEL
     226             : 
     227         304 :     WRITE (6,FMT=8000)
     228             : 8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
     229             :          &     t33,'d',t42,'f',t51,'total')
     230             : 
     231         886 :     DO itype = 1,atoms%ntype
     232        1468 :        DO ispin = jsp_start,jsp_end
     233         582 :           WRITE ( 6,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),moments%chmom(itype,ispin)
     234             : 8100      FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
     235        4074 :           attributes = ''
     236         582 :           WRITE(attributes(1),'(i0)') itype
     237         582 :           WRITE(attributes(2),'(f12.7)') moments%chmom(itype,ispin)
     238         582 :           WRITE(attributes(3),'(f12.7)') qmtl(0,ispin,itype)
     239         582 :           WRITE(attributes(4),'(f12.7)') qmtl(1,ispin,itype)
     240         582 :           WRITE(attributes(5),'(f12.7)') qmtl(2,ispin,itype)
     241         582 :           WRITE(attributes(6),'(f12.7)') qmtl(3,ispin,itype)
     242             :           CALL writeXMLElementForm('mtCharge',(/'atomType','total   ','s       ','p       ','d       ','f       '/),attributes,&
     243        1164 :                                    reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
     244             :        ENDDO
     245             :     ENDDO
     246             : 
     247             :    ENDIF !(mpi%irank==0) THEN
     248         608 :     CALL timestop("cdnmt")
     249             : 
     250             :  
     251         608 :   END SUBROUTINE cdnmt
     252             : END MODULE m_cdnmt

Generated by: LCOV version 1.13