cdntot.f90 7.91 KB
Newer Older
1
MODULE m_cdntot
2 3 4 5
!     ********************************************************
!     calculate the total charge density in the interstial.,
!     vacuum, and mt regions      c.l.fu
!     ********************************************************
6
CONTAINS
Matthias Redies's avatar
Matthias Redies committed
7 8 9 10
   SUBROUTINE integrate_grid(xcpot, stars, atoms, sym, vacuum, input, cell, oneD, sphhar,noco,&
                             is_inte, mt_inte, &
                             q, qis, qmt, qvac, qtot, qistot)
      USE m_pw_tofrom_grid
Matthias Redies's avatar
Matthias Redies committed
11
      USE m_mt_tofrom_grid
Matthias Redies's avatar
Matthias Redies committed
12 13
      USE m_types
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
14
      !USE m_types_xcpot
Matthias Redies's avatar
Matthias Redies committed
15
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
16
      CLASS(t_xcpot),INTENT(IN) :: xcpot
Matthias Redies's avatar
Matthias Redies committed
17 18 19 20 21 22 23 24 25
      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_cell),INTENT(IN)   :: cell
      TYPE(t_oneD),INTENT(IN)   :: oneD
      TYPE(t_sphhar), INTENT(IN):: sphhar
      TYPE(t_noco), INTENT(INOUT)   :: noco
Matthias Redies's avatar
Matthias Redies committed
26
      TYPE(t_grid), INTENT(in)  :: is_inte, mt_inte(:)
Matthias Redies's avatar
Matthias Redies committed
27 28 29 30 31
      REAL, INTENT(out)  :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
                            qvac(2,input%jspins), qtot, qistot

      TYPE(t_potden)     :: integrand
      
Matthias Redies's avatar
Matthias Redies committed
32
      TYPE(t_grid)       :: is_inte_mut
Matthias Redies's avatar
Matthias Redies committed
33
      INTEGER            :: n
Matthias Redies's avatar
Matthias Redies committed
34 35
      
      call init_pw_grid(xcpot, stars, sym, cell)
Matthias Redies's avatar
Matthias Redies committed
36 37
      call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)

Matthias Redies's avatar
Matthias Redies committed
38
      is_inte_mut = is_inte
Matthias Redies's avatar
Matthias Redies committed
39 40

      !allocate potden type
Matthias Redies's avatar
Matthias Redies committed
41
      call integrand%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
Matthias Redies's avatar
Matthias Redies committed
42
      allocate(integrand%pw_w, mold=integrand%pw)
Matthias Redies's avatar
Matthias Redies committed
43 44

      !put is in potden-basis
Matthias Redies's avatar
Matthias Redies committed
45 46
      call pw_from_grid(xcpot, stars,.True., is_inte_mut%grid, integrand%pw, integrand%pw_w)

Matthias Redies's avatar
Matthias Redies committed
47 48
      !put mt in potden-basis
      do n = 1,atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
49
         call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte(n)%grid,integrand%mt(:,0:,n,:))
Matthias Redies's avatar
Matthias Redies committed
50 51 52 53 54
      enddo

      ! integrate my integrand
      call integrate_cdn(stars, atoms, sym, vacuum, input, cell, oneD, integrand,&
                        q, qis, qmt, qvac, qtot, qistot)
Matthias Redies's avatar
Matthias Redies committed
55
      call finish_pw_grid()
Matthias Redies's avatar
Matthias Redies committed
56
      call finish_mt_grid()
Matthias Redies's avatar
Matthias Redies committed
57 58 59
   END SUBROUTINE integrate_grid

   SUBROUTINE integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, integrand, &
60
                                   q, qis, qmt, qvac, qtot, qistot)
61 62 63 64 65
      USE m_intgr, ONLY : intgr3
      USE m_constants
      USE m_qsf
      USE m_pwint
      USE m_types
66
      USE m_juDFT
Daniel Wortmann's avatar
Daniel Wortmann committed
67
      USE m_convol
68
      IMPLICIT NONE
69 70 71 72 73 74
      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_cell),INTENT(IN)   :: cell
75 76 77 78 79 80 81 82 83 84 85 86
      TYPE(t_oneD),INTENT(IN)   :: oneD
      TYPE(t_potden),INTENT(IN) :: integrand
      REAL, INTENT(out)         :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
                                   qvac(2,input%jspins), qtot, qistot
      INTEGER                   :: jsp, j, ivac, nz, n
      REAL                      :: q2(vacuum%nmz), w, rht1(vacuum%nmzd,2,input%jspins)
      COMPLEX                   :: x(stars%ng3)
      
      qtot = 0.0
      qistot = 0.0
      DO jsp = 1,input%jspins
         q(jsp) = 0.0
87 88
!     -----mt charge
         CALL timestart("MT")
89
         DO n = 1,atoms%ntype
90 91 92
            CALL intgr3(integrand%mt(:,0,n,jsp),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),w)
            qmt(n, jsp) = w*sfp_const
            q(jsp) = q(jsp) + atoms%neq(n)*qmt(n,jsp)
93
         ENDDO
94 95 96
         CALL timestop("MT")
!     -----vacuum region
         IF (input%film) THEN
97
            DO ivac = 1,vacuum%nvac
98 99
               DO nz = 1,vacuum%nmz
                  IF (oneD%odi%d1) THEN
100 101
                     rht1(nz,ivac,jsp) = (cell%z1+(nz-1)*vacuum%delz)*&
                                           integrand%vacz(nz,ivac,jsp)
102
                  ELSE
103
                     rht1(nz,ivac,jsp) =  integrand%vacz(nz,ivac,jsp)
104 105
                  END IF
               END DO
106 107
               CALL qsf(vacuum%delz,rht1(1,ivac,jsp),q2,vacuum%nmz,0)
               qvac(ivac,jsp) = q2(1)*cell%area
108
               IF (.NOT.oneD%odi%d1) THEN
109
                  q(jsp) = q(jsp) + qvac(ivac,jsp)*2./real(vacuum%nvac)
110
               ELSE
111
                  q(jsp) = q(jsp) + cell%area*q2(1)
112
               END IF
113
            ENDDO
114 115
         END IF
!     -----is region
116
         qis(jsp) = 0.
Matthias Redies's avatar
Matthias Redies committed
117 118 119

         CALL pwint_all(stars,atoms,sym,oneD,cell,1,stars%ng3,x)
         DO j = 1,stars%ng3
120
            qis(jsp) = qis(jsp) + integrand%pw(j,jsp)*x(j)*stars%nstr(j)
Matthias Redies's avatar
Matthias Redies committed
121
         ENDDO
122

123 124 125 126
         qistot = qistot + qis(jsp)
         q(jsp) = q(jsp) + qis(jsp)
         qtot = qtot + q(jsp)
      END DO ! loop over spins
Matthias Redies's avatar
Matthias Redies committed
127
   END SUBROUTINE integrate_cdn
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163

   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(input%jspins),qis(input%jspins),w,mtCharge
      INTEGER i,ivac,j,jsp,n,nz
!     ..
!     .. Local Arrays ..
      REAL qmt(atoms%ntype,input%jspins),qvac(2,input%jspins)
      INTEGER, ALLOCATABLE :: lengths(:,:)
      CHARACTER(LEN=20) :: attributes(6), names(6)
      
Matthias Redies's avatar
Matthias Redies committed
164
      CALL timestart("cdntot")
Matthias Redies's avatar
Matthias Redies committed
165
      call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, den, &
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
                                   q, qis, qmt, qvac, qtot, qistot)
 
      IF (input%film) THEN
         ALLOCATE(lengths(4+vacuum%nvac,2))
      ELSE
         ALLOCATE(lengths(4,2))
      END IF

      DO jsp = 1,input%jspins
         WRITE (6,FMT=8000) jsp,q(jsp),qis, (qmt(n,jsp),n=1,atoms%ntype)
         IF (input%film) WRITE (6,FMT=8010) (i,qvac(i,jsp),i=1,vacuum%nvac)
         mtCharge = SUM(qmt(1:atoms%ntype,jsp) * atoms%neq(1:atoms%ntype))
         names(1) = 'spin'         ; WRITE(attributes(1),'(i0)')    jsp      ; lengths(1,1)=4  ; lengths(1,2)=1
         names(2) = 'total'        ; WRITE(attributes(2),'(f14.7)') q(jsp)   ; lengths(2,1)=5  ; lengths(2,2)=14
         names(3) = 'interstitial' ; WRITE(attributes(3),'(f14.7)') qis(jsp) ; lengths(3,1)=12 ; lengths(3,2)=14
181
         names(4) = 'mtSpheres'    ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9  ; lengths(4,2)=14
182 183 184 185
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
186
                  WRITE(attributes(4+i),'(f14.7)') qvac(i,jsp)
187
                  lengths(4+i,1)=7
188
                  lengths(4+i,2)=14
189 190 191 192 193 194
               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
195 196
         END IF
      END DO ! loop over spins
197
      WRITE (6,FMT=8020) qtot
198 199 200
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
201 202 203 204 205
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)
206 207

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