cdntot.f90 7.55 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 11 12
   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
      USE m_types
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
13
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
14 15 16 17 18 19 20 21 22 23
      CLASS(t_xcpot),INTENT(IN)     :: xcpot
      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
24 25 26 27 28 29 30 31 32
      REAL, intent(in)   :: is_inte(:,:), mt_inte(:,:)
      REAL, INTENT(out)  :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
                            qvac(2,input%jspins), qtot, qistot

      TYPE(t_potden)     :: integrand
      
      INTEGER            :: n

      !allocate potden type
Matthias Redies's avatar
Matthias Redies committed
33
      call integrand%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
Matthias Redies's avatar
Matthias Redies committed
34 35 36 37 38

      !put is in potden-basis
      call pw_from_grid(xcpot, stars,.True., is_inte, integrand%pw, integrand%pw_w)
      !put mt in potden-basis
      do n = 1,atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
39
         call mt_from_grid(atoms,sphhar,n,input%jspins,mt_inte,integrand%mt(:,0:,n,:))
Matthias Redies's avatar
Matthias Redies committed
40 41 42 43 44 45 46 47
      enddo

      ! integrate my integrand
      call integrate_cdn(stars, atoms, sym, vacuum, input, cell, oneD, integrand,&
                        q, qis, qmt, qvac, qtot, qistot)
   END SUBROUTINE integrate_grid

   SUBROUTINE integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, integrand, &
48
                                   q, qis, qmt, qvac, qtot, qistot)
49 50 51 52 53
      USE m_intgr, ONLY : intgr3
      USE m_constants
      USE m_qsf
      USE m_pwint
      USE m_types
54
      USE m_juDFT
Daniel Wortmann's avatar
Daniel Wortmann committed
55
      USE m_convol
56
      IMPLICIT NONE
57 58 59 60 61 62
      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
63 64 65 66 67 68 69 70
      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)
      
71
      CALL timestart("cdntot")
72 73 74 75
      qtot = 0.0
      qistot = 0.0
      DO jsp = 1,input%jspins
         q(jsp) = 0.0
76 77
!     -----mt charge
         CALL timestart("MT")
78
         DO n = 1,atoms%ntype
79 80 81
            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)
82
         ENDDO
83 84 85
         CALL timestop("MT")
!     -----vacuum region
         IF (input%film) THEN
86
            DO ivac = 1,vacuum%nvac
87 88
               DO nz = 1,vacuum%nmz
                  IF (oneD%odi%d1) THEN
89 90
                     rht1(nz,ivac,jsp) = (cell%z1+(nz-1)*vacuum%delz)*&
                                           integrand%vacz(nz,ivac,jsp)
91
                  ELSE
92
                     rht1(nz,ivac,jsp) =  integrand%vacz(nz,ivac,jsp)
93 94
                  END IF
               END DO
95 96
               CALL qsf(vacuum%delz,rht1(1,ivac,jsp),q2,vacuum%nmz,0)
               qvac(ivac,jsp) = q2(1)*cell%area
97
               IF (.NOT.oneD%odi%d1) THEN
98
                  q(jsp) = q(jsp) + qvac(ivac,jsp)*2./real(vacuum%nvac)
99
               ELSE
100
                  q(jsp) = q(jsp) + cell%area*q2(1)
101
               END IF
102
            ENDDO
103 104
         END IF
!     -----is region
105
         qis(jsp) = 0.
Matthias Redies's avatar
Matthias Redies committed
106 107 108

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

112 113 114 115
         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
116
   END SUBROUTINE integrate_cdn
117 118 119 120 121 122 123 124 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

   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
153
      call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, den, &
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
                                   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
169
         names(4) = 'mtSpheres'    ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9  ; lengths(4,2)=14
170 171 172 173
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
174
                  WRITE(attributes(4+i),'(f14.7)') qvac(i,jsp)
175
                  lengths(4+i,1)=7
176
                  lengths(4+i,2)=14
177 178 179 180 181 182
               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
183 184
         END IF
      END DO ! loop over spins
185
      WRITE (6,FMT=8020) qtot
186 187 188
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
189 190 191 192 193
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)
194 195

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