brysh1.f90 7.28 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 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 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 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
MODULE m_brysh1
  USE m_juDFT
  !******************************************************
  !      shifts the charge density of the interstitial, m.t
  !      and vacuum part in one single vector
  !      in the spin polarized case the arrays consist of 
  !      spin up and spin down densities
  !******************************************************
CONTAINS
  SUBROUTINE brysh1(&
       &                  input,stars,atoms,sphhar,noco,vacuum,sym,oneD,&
       &                  intfac,vacfac,qpw,rho,rht,rhtxy,cdom,cdomvz,cdomvxy,n_mmp,&
       &                  nmap,nmaph,mapmt,mapvac,mapvac2,sout) 

    USE m_types
    IMPLICIT NONE
    TYPE(t_oneD),INTENT(IN)    :: oneD
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_noco),INTENT(IN)    :: noco
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_stars),INTENT(IN)   :: stars
    TYPE(t_sphhar),INTENT(IN)  :: sphhar
    TYPE(t_atoms),INTENT(IN)   :: atoms
    !     ..
    !     .. Scalar Arguments ..
    REAL,    INTENT (IN) :: intfac,vacfac
    INTEGER, INTENT (OUT):: mapmt,mapvac,mapvac2,nmap,nmaph
    !     ..
    !     .. Array Arguments ..
    !-odim
    !+odim
    COMPLEX, INTENT (IN) :: qpw(stars%n3d,input%jspins)
    COMPLEX, INTENT (IN) :: cdomvz(vacuum%nmz,2),cdomvxy(vacuum%nmzxy,oneD%odi%n2d-1,2)
    COMPLEX, INTENT (IN) :: cdom(stars%n3d),rhtxy(vacuum%nmzxy,oneD%odi%n2d-1,2,input%jspins)
    COMPLEX, INTENT (IN) :: n_mmp(-3:3,-3:3,atoms%n_u,input%jspins)
    REAL,    INTENT (IN) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,input%jspins)
    REAL,    INTENT (IN) :: rht(vacuum%nmz,2,input%jspins)
    REAL,    INTENT (OUT):: sout(:)
    !     ..
    !     .. Local Scalars ..
    INTEGER i,iv,j,js,k,l,n,nall,na,nvaccoeff,nvaccoeff2,mapmtd
    !
    !--->  put input into arrays sout 
    !      in the spin polarized case the arrays consist of 
    !      spin up and spin down densities

    j=0
    DO  js = 1,input%jspins
       DO i = 1,stars%ng3
          j = j + 1
          sout(j) = REAL(qpw(i,js))
       END DO
       IF (.NOT.sym%invs) THEN
          DO i = 1,stars%ng3
             j = j + 1
             sout(j) = AIMAG(qpw(i,js))
          END DO
       ENDIF
       mapmt=0
       na = 1
       DO n = 1,atoms%ntype
          DO l = 0,sphhar%nlh(atoms%ntypsy(na))
             DO i = 1,atoms%jri(n)
                mapmt = mapmt +1
                j = j + 1
                sout(j) = rho(i,l,n,js)
             END DO
          END DO
          na = na + atoms%neq(n)
       END DO
       mapvac=0
       IF (input%film) THEN
          DO iv = 1,vacuum%nvac
             DO k = 1,vacuum%nmz
                mapvac = mapvac + 1
                j = j + 1
                sout(j) = rht(k,iv,js)
             END DO
             DO k = 1,oneD%odi%nq2-1
                DO i = 1,vacuum%nmzxy
                   mapvac = mapvac + 1
                   j = j + 1
                   sout(j) =  REAL(rhtxy(i,k,iv,js))
                END DO
             END DO
             IF (.NOT.sym%invs2) THEN
                DO k = 1,oneD%odi%nq2-1
                   DO i = 1,vacuum%nmzxy
                      mapvac = mapvac + 1
                      j = j + 1
                      sout(j) =  AIMAG(rhtxy(i,k,iv,js))
                   END DO
                END DO
             END IF
          END DO
       END IF
       IF (js .EQ. 1) nmaph = j
    ENDDO

    mapvac2=0
    IF (noco%l_noco) THEN
       !--->    off-diagonal part of the density matrix
       DO i = 1,stars%ng3
          j = j + 1
          sout(j) = REAL(cdom(i))
       END DO
       DO i = 1,stars%ng3
          j = j + 1
          sout(j) = AIMAG(cdom(i))
       END DO
       IF (input%film) THEN
          DO iv = 1,vacuum%nvac
             DO k = 1,vacuum%nmz
                mapvac2 = mapvac2 + 1
                j = j + 1
                sout(j) = REAL(cdomvz(k,iv))
             END DO
             DO k = 1,oneD%odi%nq2-1
                DO i = 1,vacuum%nmzxy
                   mapvac2 = mapvac2 + 1
                   j = j + 1
                   sout(j) =  REAL(cdomvxy(i,k,iv))
                END DO
             END DO
          END DO
          DO iv = 1,vacuum%nvac
             DO k = 1,vacuum%nmz
                mapvac2 = mapvac2 + 1
                j = j + 1
                sout(j) = AIMAG(cdomvz(k,iv))
             END DO
             DO k = 1,oneD%odi%nq2-1
                DO i = 1,vacuum%nmzxy
                   mapvac2 = mapvac2 + 1
                   j = j + 1
                   sout(j) =  AIMAG(cdomvxy(i,k,iv))
                END DO
             END DO
          END DO
          nvaccoeff2 = 2*vacuum%nmzxy*(oneD%odi%nq2-1)*vacuum%nvac + 2*vacuum%nmz*vacuum%nvac
          IF (mapvac2 .NE. nvaccoeff2) THEN
             WRITE (6,*)'The number of vaccum coefficients off the'
             WRITE (6,*)'off-diagonal part of the density matrix is'
             WRITE (6,*)'inconsitent:'
             WRITE (6,8000) mapvac2,nvaccoeff2
8000         FORMAT ('mapvac2= ',i12,'nvaccoeff2= ',i12)
             CALL juDFT_error("brysh1:# of vacuum coeff. inconsistent"&
                  &              ,calledby ="brysh1")
          ENDIF
       END IF
    ENDIF ! noco

    IF (atoms%n_u > 0 ) THEN     ! lda+U
       DO js = 1,input%jspins
          DO n = 1, atoms%n_u
             DO k = -3, 3
                DO i = -3, 3
                   j = j + 1 
                   sout(j) = REAL(n_mmp(i,k,n,js))
                   j = j + 1 
                   sout(j) = AIMAG(n_mmp(i,k,n,js))
                ENDDO
             ENDDO
          ENDDO
       ENDDO
    ENDIF

    IF (input%film) THEN
       nvaccoeff = vacfac*vacuum%nmzxy*(oneD%odi%nq2-1)*vacuum%nvac + vacuum%nmz*vacuum%nvac
       IF (mapvac .NE. nvaccoeff) THEN
          WRITE(6,*)'The number of vaccum coefficients is'
          WRITE(6,*)'inconsitent:'
          WRITE (6,8010) mapvac,nvaccoeff
8010      FORMAT ('mapvac= ',i12,'nvaccoeff= ',i12)
          CALL juDFT_error("brysh1: # of vacuum coeff. inconsistent"&
               &           ,calledby ="brysh1")
       ENDIF
    ENDIF

    mapmtd = atoms%ntypd*(sphhar%nlhd+1)*atoms%jmtd
    IF (mapmt .GT. mapmtd) THEN
       WRITE(6,*)'The number of mt coefficients is larger than the'
       WRITE(6,*)'dimensions:'
       WRITE (6,8040) mapmt,mapmtd
8040   FORMAT ('mapmt= ',i12,' > mapmtd= ',i12)
       CALL juDFT_error("brysh1: mapmt > mapmtd (dimensions)",calledby&
            &        ="brysh1")
    ENDIF

    nmap = j
    nall = (intfac*stars%ng3 + mapmt + mapvac + 49*2*atoms%n_u )*input%jspins
    IF (noco%l_noco) nall = nall + 2*stars%ng3 + mapvac2
    IF (nall.NE.nmap) THEN
       WRITE(6,*)'The input%total number of charge density coefficients is'
       WRITE(6,*)'inconsitent:'
       WRITE (6,8020) nall,nmap
8020   FORMAT ('nall= ',i12,'not equal nmap= ',i12)
       WRITE (6,'(a,i5,a,i5)') 'nall = ',nall,' nmap = ',nmap
       CALL juDFT_error&
            &        ("brysh1: input # of charge density coeff. inconsistent"&
            &        ,calledby ="brysh1")
    ENDIF
    IF (nmap.GT.SIZE(sout)) THEN 
       WRITE(6,*)'The input%total number of charge density coefficients is'
       WRITE(6,*)'larger than the dimensions:'
       WRITE (6,8030) nmap,SIZE(sout)
8030   FORMAT ('nmap= ',i12,' > size(sout)= ',i12)
       CALL juDFT_error("brysh1: nmap > mmap (dimensions)",calledby&
            &        ="brysh1")
    ENDIF

  END SUBROUTINE brysh1
END MODULE m_brysh1