MODULE m_cdntot
!     ********************************************************
!     calculate the total charge density in the interstial.,
!     vacuum, and mt regions      c.l.fu
!     ********************************************************
      CONTAINS
      SUBROUTINE cdntot(stars,atoms,sym,vacuum,input,cell,oneD,&
                        den,l_printData,qtot,qistot)

      USE m_intgr, ONLY : intgr3
      USE m_constants
      USE m_qsf
      USE m_pwint
      USE m_types
      USE m_juDFT
      USE m_convol
      USE m_xmlOutput
      IMPLICIT NONE

!     .. Scalar Arguments ..
      TYPE(t_stars),INTENT(IN)  :: stars
      TYPE(t_atoms),INTENT(IN)  :: atoms
      TYPE(t_sym),INTENT(IN)    :: sym
      TYPE(t_vacuum),INTENT(IN) :: vacuum
      TYPE(t_input),INTENT(IN)  :: input
      TYPE(t_oneD),INTENT(IN)   :: oneD
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_potden),INTENT(IN) :: den
      LOGICAL,INTENT(IN)        :: l_printData
      REAL,INTENT(OUT)          :: qtot,qistot

!     .. Local Scalars ..
      COMPLEX x(stars%ng3)
      REAL q,qis,w,mtCharge
      INTEGER i,ivac,j,jspin,n,nz
!     ..
!     .. Local Arrays ..
      REAL qmt(atoms%ntype),qvac(2),q2(vacuum%nmz),rht1(vacuum%nmzd,2,input%jspins)
      INTEGER, ALLOCATABLE :: lengths(:,:)
      CHARACTER(LEN=20) :: attributes(6), names(6)
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC real
!     ..
!
      IF (input%film) THEN
         ALLOCATE(lengths(4+vacuum%nvac,2))
      ELSE
         ALLOCATE(lengths(4,2))
      END IF
      CALL timestart("cdntot")
      qtot = 0.e0
      qistot = 0.e0
      DO jspin = 1,input%jspins
         q = 0.e0
!     -----mt charge
         CALL timestart("MT")
         DO 10 n = 1,atoms%ntype
            CALL intgr3(den%mt(:,0,n,jspin),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),w)
            qmt(n) = w*sfp_const
            q = q + atoms%neq(n)*qmt(n)
   10    CONTINUE
         CALL timestop("MT")
!     -----vacuum region
         IF (input%film) THEN
            DO 20 ivac = 1,vacuum%nvac
               DO nz = 1,vacuum%nmz
                  IF (oneD%odi%d1) THEN
                     rht1(nz,ivac,jspin) = (cell%z1+(nz-1)*vacuum%delz)*&
     &                    den%vacz(nz,ivac,jspin)
                  ELSE
                     rht1(nz,ivac,jspin) =  den%vacz(nz,ivac,jspin)
                  END IF
               END DO
               CALL qsf(vacuum%delz,rht1(1,ivac,jspin),q2,vacuum%nmz,0)
               qvac(ivac) = q2(1)*cell%area
               IF (.NOT.oneD%odi%d1) THEN
                  q = q + qvac(ivac)*2./real(vacuum%nvac)
               ELSE
                  q = q + cell%area*q2(1)
               END IF
   20       CONTINUE
         END IF
!     -----is region
         IF (.not.judft_was_Argument("-oldfix")) THEN
            CALL convol(stars,x,den%pw(:,jspin),stars%ufft)
            qis = x(1)*cell%omtil
         ELSE
          qis = 0.
!         DO 30 j = 1,nq3
!            CALL pwint(
!     >                 k1d,k2d,k3d,n3d,ntypd,natd,nop,invtab,odi,
!     >                 ntype,neq,volmts,taual,z1,vol,volint,
!     >                 symor,tau,mrot,rmt,sk3,bmat,ig2,ig,
!     >                 kv3(1,j),
!     <                 x)
!            qis = qis + den%pw(j,jspin)*x*nstr(j)
!   30    CONTINUE
         CALL pwint_all(&
     &                 stars,atoms,sym,oneD,&
     &                 cell,&
     &                 x)
         DO j = 1,stars%ng3
             qis = qis + den%pw(j,jspin)*x(j)*stars%nstr(j)
         ENDDO
         endif
         qistot = qistot + qis
         q = q + qis
         WRITE (6,FMT=8000) jspin,q,qis, (qmt(n),n=1,atoms%ntype)
         IF (input%film) WRITE (6,FMT=8010) (i,qvac(i),i=1,vacuum%nvac)
         WRITE (16,FMT=8000) jspin,q,qis, (qmt(n),n=1,atoms%ntype)
         IF (input%film) WRITE (16,FMT=8010) (i,qvac(i),i=1,vacuum%nvac)
         mtCharge = SUM(qmt(1:atoms%ntype) * atoms%neq(1:atoms%ntype))
         names(1) = 'spin'         ; WRITE(attributes(1),'(i0)') jspin       ; lengths(1,1)=4  ; lengths(1,2)=1
         names(2) = 'total'        ; WRITE(attributes(2),'(f14.7)') q        ; lengths(2,1)=5  ; lengths(2,2)=14
         names(3) = 'interstitial' ; WRITE(attributes(3),'(f14.7)') qis      ; lengths(3,1)=12 ; lengths(3,2)=14
         names(4) = 'mtSpheres'    ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9  ; lengths(4,2)=14
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
                  WRITE(attributes(4+i),'(f14.7)') qvac(i)
                  lengths(4+i,1)=7
                  lengths(4+i,2)=14
               END DO
               CALL writeXMLElementFormPoly('spinDependentCharge',names(1:4+vacuum%nvac),&
                                            attributes(1:4+vacuum%nvac),lengths)
            ELSE
               CALL writeXMLElementFormPoly('spinDependentCharge',names(1:4),attributes(1:4),lengths)
            END IF
         END IF
         qtot = qtot + q
      END DO ! loop over spins
      DEALLOCATE (lengths)
      WRITE (6,FMT=8020) qtot
      WRITE (16,FMT=8020) qtot
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
 8000 FORMAT (/,10x,'total charge for spin',i3,'=',f12.6,/,10x,&
     &       'interst. charge =   ',f12.6,/,&
     &       (10x,'mt charge=          ',4f12.6,/))
 8010 FORMAT (10x,'vacuum ',i2,'  charge=  ',f12.6)
 8020 FORMAT (/,10x,'total charge  =',f12.6)

      CALL timestop("cdntot")
      END SUBROUTINE cdntot
      END MODULE m_cdntot