cdntot.f90 6.01 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
7
   SUBROUTINE cdntot_integrate(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
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
      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)
      
31
      CALL timestart("cdntot")
32 33 34 35
      qtot = 0.0
      qistot = 0.0
      DO jsp = 1,input%jspins
         q(jsp) = 0.0
36 37
!     -----mt charge
         CALL timestart("MT")
38
         DO n = 1,atoms%ntype
39 40 41
            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)
42
         ENDDO
43 44 45
         CALL timestop("MT")
!     -----vacuum region
         IF (input%film) THEN
46
            DO ivac = 1,vacuum%nvac
47 48
               DO nz = 1,vacuum%nmz
                  IF (oneD%odi%d1) THEN
49 50
                     rht1(nz,ivac,jsp) = (cell%z1+(nz-1)*vacuum%delz)*&
                                           integrand%vacz(nz,ivac,jsp)
51
                  ELSE
52
                     rht1(nz,ivac,jsp) =  integrand%vacz(nz,ivac,jsp)
53 54
                  END IF
               END DO
55 56
               CALL qsf(vacuum%delz,rht1(1,ivac,jsp),q2,vacuum%nmz,0)
               qvac(ivac,jsp) = q2(1)*cell%area
57
               IF (.NOT.oneD%odi%d1) THEN
58
                  q(jsp) = q(jsp) + qvac(ivac,jsp)*2./real(vacuum%nvac)
59
               ELSE
60
                  q(jsp) = q(jsp) + cell%area*q2(1)
61
               END IF
62
            ENDDO
63 64
         END IF
!     -----is region
65
         qis(jsp) = 0.
Matthias Redies's avatar
Matthias Redies committed
66 67 68

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

72 73 74 75
         qistot = qistot + qis(jsp)
         q(jsp) = q(jsp) + qis(jsp)
         qtot = qtot + q(jsp)
      END DO ! loop over spins
76
   END SUBROUTINE cdntot_integrate
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

   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)
      
113
      call cdntot_integrate(stars,atoms,sym,vacuum,input,cell,oneD, den, &
114 115 116 117 118 119 120 121 122
                                   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
123
         WRITE (6,FMT=8000) jsp,q(jsp),qis(jsp), (qmt(n,jsp),n=1,atoms%ntype)
124 125 126 127 128
         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
129
         names(4) = 'mtSpheres'    ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9  ; lengths(4,2)=14
130 131 132 133
         IF(l_printData) THEN
            IF(input%film) THEN
               DO i = 1, vacuum%nvac
                  WRITE(names(4+i),'(a6,i0)') 'vacuum', i
134
                  WRITE(attributes(4+i),'(f14.7)') qvac(i,jsp)
135
                  lengths(4+i,1)=7
136
                  lengths(4+i,2)=14
137 138 139 140 141 142
               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
143 144
         END IF
      END DO ! loop over spins
145
      WRITE (6,FMT=8020) qtot
146 147 148
      IF(l_printData) THEN
         CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
      END IF
149 150 151 152 153
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)
154 155

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