cdntot.f90 5.53 KB
Newer Older
1 2 3 4 5 6
      MODULE m_cdntot
!     ********************************************************
!     calculate the total charge density in the interstial.,
!     vacuum, and mt regions      c.l.fu
!     ********************************************************
      CONTAINS
7 8
      SUBROUTINE cdntot(stars,atoms,sym,vacuum,input,cell,oneD,&
                        den,l_printData,qtot,qistot)
9 10 11 12 13 14

      USE m_intgr, ONLY : intgr3
      USE m_constants
      USE m_qsf
      USE m_pwint
      USE m_types
15
      USE m_juDFT
16
      USE m_convol
17
      USE m_xmlOutput
18
      IMPLICIT NONE
19

20
!     .. Scalar Arguments ..
21 22 23 24 25 26 27 28 29 30 31
      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

32 33
!     .. Local Scalars ..
      COMPLEX x(stars%ng3)
34
      REAL q,qis,w,mtCharge
35 36 37
      INTEGER i,ivac,j,jspin,n,nz
!     ..
!     .. Local Arrays ..
38
      REAL qmt(atoms%ntype),qvac(2),q2(vacuum%nmz),rht1(vacuum%nmzd,2,input%jspins)
39 40
      INTEGER, ALLOCATABLE :: lengths(:,:)
      CHARACTER(LEN=20) :: attributes(6), names(6)
41 42 43 44 45
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC real
!     ..
!
46 47 48 49 50
      IF (input%film) THEN
         ALLOCATE(lengths(4+vacuum%nvac,2))
      ELSE
         ALLOCATE(lengths(4,2))
      END IF
51 52 53
      CALL timestart("cdntot")
      qtot = 0.e0
      qistot = 0.e0
54
      DO jspin = 1,input%jspins
55 56 57 58
         q = 0.e0
!     -----mt charge
         CALL timestart("MT")
         DO 10 n = 1,atoms%ntype
59
            CALL intgr3(den%mt(:,0,n,jspin),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),w)
60 61 62 63 64 65 66 67 68 69
            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)*&
70
     &                    den%vacz(nz,ivac,jspin)
71
                  ELSE
72
                     rht1(nz,ivac,jspin) =  den%vacz(nz,ivac,jspin)
73 74 75 76 77 78 79 80 81 82 83 84
                  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
85
         IF (.not.judft_was_Argument("-oldfix")) THEN
86
            CALL convol(stars,x,den%pw(:,jspin),stars%ufft)
87 88 89
            qis = x(1)*cell%omtil
         ELSE
          qis = 0.
90 91 92 93 94 95 96
!         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)
97
!            qis = qis + den%pw(j,jspin)*x*nstr(j)
98 99 100 101 102 103
!   30    CONTINUE
         CALL pwint_all(&
     &                 stars,atoms,sym,oneD,&
     &                 cell,&
     &                 x)
         DO j = 1,stars%ng3
104
             qis = qis + den%pw(j,jspin)*x(j)*stars%nstr(j)
105
         ENDDO
106
         endif
107 108 109 110 111 112
         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)
113 114 115 116 117
         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
118 119 120 121
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
122
                  WRITE(attributes(4+i),'(f14.7)') qvac(i)
123
                  lengths(4+i,1)=7
124
                  lengths(4+i,2)=14
125 126 127 128 129 130
               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
131
         END IF
132
         qtot = qtot + q
133 134
      END DO ! loop over spins
      DEALLOCATE (lengths)
135 136
      WRITE (6,FMT=8020) qtot
      WRITE (16,FMT=8020) qtot
137 138 139
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
140 141 142 143 144 145 146 147 148
 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