cdntot.f90 8.79 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
   SUBROUTINE integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, integrand, &
8
                                   q, qis, qmt, qvac, qtot, qistot)
9 10 11 12 13
      USE m_intgr, ONLY : intgr3
      USE m_constants
      USE m_qsf
      USE m_pwint
      USE m_types
14
      USE m_juDFT
Daniel Wortmann's avatar
Daniel Wortmann committed
15
      USE m_convol
16
      IMPLICIT NONE
17 18 19 20 21 22
      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
23 24 25 26 27 28 29 30 31 32 33 34
      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
35 36
!     -----mt charge
         CALL timestart("MT")
37
         DO n = 1,atoms%ntype
38 39 40
            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)
41
         ENDDO
42 43 44
         CALL timestop("MT")
!     -----vacuum region
         IF (input%film) THEN
45
            DO ivac = 1,vacuum%nvac
46 47
               DO nz = 1,vacuum%nmz
                  IF (oneD%odi%d1) THEN
48 49
                     rht1(nz,ivac,jsp) = (cell%z1+(nz-1)*vacuum%delz)*&
                                           integrand%vacz(nz,ivac,jsp)
50
                  ELSE
51
                     rht1(nz,ivac,jsp) =  integrand%vacz(nz,ivac,jsp)
52 53
                  END IF
               END DO
54 55
               CALL qsf(vacuum%delz,rht1(1,ivac,jsp),q2,vacuum%nmz,0)
               qvac(ivac,jsp) = q2(1)*cell%area
56
               IF (.NOT.oneD%odi%d1) THEN
57
                  q(jsp) = q(jsp) + qvac(ivac,jsp)*2./real(vacuum%nvac)
58
               ELSE
59
                  q(jsp) = q(jsp) + cell%area*q2(1)
60
               END IF
61
            ENDDO
62 63
         END IF
!     -----is region
64
         qis(jsp) = 0.
Matthias Redies's avatar
Matthias Redies committed
65 66 67

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

71 72 73 74
         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
75
   END SUBROUTINE integrate_cdn
76

Matthias Redies's avatar
Matthias Redies committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
   SUBROUTINE integrate_realspace(xcpot, atoms, sym, sphhar, input, &
                                  stars, cell, oneD, vacuum, noco, mt, is, hint)
      use m_types
      use m_mt_tofrom_grid
      use m_pw_tofrom_grid
      use m_constants
      implicit none
      CLASS(t_xcpot), INTENT(inout)   :: xcpot
      TYPE(t_atoms),INTENT(IN)      :: atoms
      TYPE(t_sym), INTENT(in)       :: sym
      TYPE(t_sphhar), INTENT(IN)    :: sphhar
      TYPE(t_input), INTENT(IN)     :: input
      TYPE(t_stars), INTENT(IN)     :: stars
      TYPE(t_cell), INTENT(IN)      :: cell
      TYPE(t_oneD), INTENT(in)      :: oneD
      TYPE(t_vacuum), INTENT(in)    :: vacuum
      TYPE(t_noco), INTENT(in)      :: noco
      real, intent(inout)           :: mt(:,:,:), is(:,:)
      character(len=*), intent(in), optional :: hint
      integer                       :: n_atm, i

      TYPE(t_potden)                :: tmp_potden
      REAL                          :: q(input%jspins), qis(input%jspins), &
                                       qmt(atoms%ntype,input%jspins), qvac(2,input%jspins),&
                                       qtot, qistot

      call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
      call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
      do n_atm =1,atoms%ntype
         call mt_from_grid(atoms, sphhar, n_atm, input%jspins, mt(:,:,n_atm), &
                           tmp_potden%mt(:,0:,n_atm,:))

         do i=1,atoms%jri(n_atm)
            tmp_potden%mt(i,:,n_atm,:) = tmp_potden%mt(i,:,n_atm,:) * atoms%rmsh(i,n_atm)**2
         enddo
      enddo
      call finish_mt_grid()

      call init_pw_grid(xcpot, stars, sym, cell)
      call pw_from_grid(xcpot, stars, .False., is, tmp_potden%pw)
      call finish_pw_grid()

      call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, tmp_potden, &
                                   q, qis, qmt, qvac, qtot, qistot)

      call print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
   END SUBROUTINE integrate_realspace

125 126 127 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
   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
160
      CALL timestart("cdntot")
Matthias Redies's avatar
Matthias Redies committed
161
      call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, den, &
162 163 164 165 166 167 168 169 170
                                   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
171
         WRITE (6,FMT=8000) jsp,q(jsp),qis(jsp), (qmt(n,jsp),n=1,atoms%ntype)
172 173 174 175 176
         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
177
         names(4) = 'mtSpheres'    ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9  ; lengths(4,2)=14
178 179 180 181
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
182
                  WRITE(attributes(4+i),'(f14.7)') qvac(i,jsp)
183
                  lengths(4+i,1)=7
184
                  lengths(4+i,2)=14
185 186 187 188 189 190
               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
191 192
         END IF
      END DO ! loop over spins
193
      WRITE (6,FMT=8020) qtot
194 195 196
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
197 198 199 200 201
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)
202 203

      CALL timestop("cdntot")
204
   END SUBROUTINE cdntot
Matthias Redies's avatar
Matthias Redies committed
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227

   SUBROUTINE print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
      use  ieee_arithmetic
      implicit none
      REAL, INTENT(in)                       :: q(:), qis(:), qmt(:,:), qvac(:,:), qtot, qistot
      character(len=*), intent(in), optional :: hint
      integer                                :: n_mt
      

      if(present(hint)) write (*,*) "DEN of ", hint
      write (*,*) "q   = ", q
      write (*,*) "qis = ", qis

      write (*,*) "qmt"
      do n_mt = 1,size(qmt, dim=1)
         write (*,*) "mt = ", n_mt, qmt(n_mt,:)
      enddo

      if(.not. any(ieee_is_nan(qvac))) then
         write (*, *) "qvac",    qvac
      endif
      write (*, *) "qtot",    qtot
      write (*, *) "qis_tot", qistot
Matthias Redies's avatar
Matthias Redies committed
228
      write (*, *) "-------------------------"
Matthias Redies's avatar
Matthias Redies committed
229
   END SUBROUTINE print_cdn_inte
230
END MODULE m_cdntot