outcdn.f90 6.71 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_outcdn

CONTAINS

   SUBROUTINE outcdn(p, n, na, iv, iflag, jsp, l_potential, stars, vacuum, &
                     sphhar, atoms, sym, cell, oneD, potDen, xdnout)
12
      USE m_types
13
      USE m_constants
14 15 16
      USE m_angle
      USE m_starf, ONLY : starf2,starf3
      USE m_ylm
17 18

      !--------------------------------------------------------------------------
19
      ! Calculates the charge density at a given point p(i=1,3).
20 21
      !--------------------------------------------------------------------------

22
      IMPLICIT NONE
23

24 25 26 27 28 29 30
      TYPE(t_stars),INTENT(IN)     :: stars
      TYPE(t_vacuum),INTENT(IN)    :: vacuum
      TYPE(t_sphhar),INTENT(IN)    :: sphhar
      TYPE(t_atoms),INTENT(IN)     :: atoms
      TYPE(t_sym),INTENT(IN)       :: sym
      TYPE(t_cell),INTENT(IN)      :: cell
      TYPE(t_oneD),INTENT(IN)      :: oneD
31
      TYPE(t_potden),INTENT(IN)    :: potDen
32
  
33

34
      ! Scalar Arguments
35 36
      INTEGER, INTENT (IN) :: iflag,jsp,n,na,iv
      REAL,    INTENT (OUT) :: xdnout
37 38
      ! -odim
      ! +odim
39

40
      ! Array Arguments
41
      REAL,    INTENT (INOUT) :: p(3)
42 43 44 45 46

      ! Logical argument
      LOGICAL, INTENT (IN) :: l_potential 

      ! Local scalars 
47 48
      REAL delta,s,sx,xd1,xd2,xx1,xx2,rrr,phi
      INTEGER i,j,jp3,jr,k,lh,mem,nd,nopa,ivac,ll1,lm ,gzi,m
49

50
      ! Local arrays
51
      COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm((atoms%lmaxd+1)**2)
52
      REAL rcc(3),x(3)
53
      
54 55
      ivac=iv
     
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
      IF (iflag.NE.1) THEN
         IF (iflag.NE.0) THEN
            ! Interstitial part:
            rcc=matmul(cell%bmat,p)/tpi_const
            CALL starf3(sym%nop, stars%ng3, sym%symor, stars%kv3, sym%mrot, &
                        sym%tau, rcc, sym%invtab, sf3)

            xdnout=dot_product(real(potDen%pw(:,jsp)*sf3(:)),stars%nstr)
            RETURN

         ENDIF

         ! Vacuum part:
         xdnout = 0.

         ! -odim
         IF (oneD%odi%d1) THEN
            rrr = sqrt( p(1)**2 + p(2)**2 )
            phi = angle(p(1),p(2))
            jp3 = (rrr-cell%z1)/vacuum%delz
            delta = (rrr-cell%z1)/vacuum%delz - jp3
            ! We count 0 as point 1.
            jp3 = jp3 + 1
            IF (jp3.LT.vacuum%nmz) THEN
               xdnout = potDen%vacz(jp3,ivac,jsp) + delta*(potDen%vacz(jp3+1,ivac,jsp)-potDen%vacz(jp3,ivac,jsp))
               IF (jp3.LT.vacuum%nmzxy) THEN
                  xx1 = 0.
                  xx2 = 0.
                  DO  k = 2,oneD%odi%nq2
                     m = oneD%odi%kv(2,k)
                     gzi = oneD%odi%kv(1,k)
                     xx1 = xx1 + REAL(potDen%vacxy(jp3,k-1,ivac,jsp)*EXP( &
                           ImagUnit*m*phi)*EXP(ImagUnit*gzi*cell%bmat(3,3)* &
                           p(3)))*oneD%odi%nst2(k)
                     xx2 = xx2 + REAL(potDen%vacxy(jp3+1,k-1,ivac,jsp)*EXP( &
                           ImagUnit*m*phi)*EXP(ImagUnit*gzi*cell%bmat(3,3)* &
                           p(3)))*oneD%odi%nst2(k)
                  END DO
                  xdnout = xdnout + xx1 + delta* (xx2-xx1)
               END IF
            ELSE
               xdnout = 0.0
98
            END IF
99 100
         
         ! +odim
101
         ELSE
102 103 104 105 106 107
            IF (p(3).LT.0.0) THEN
               ivac = vacuum%nvac
               IF (sym%invs) THEN
                  p(1:2) = -p(1:2)
               END IF
               p(3) = abs(p(3))
108
            END IF
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
            rcc=matmul(cell%bmat,p)/tpi_const
            CALL starf2(sym%nop2, stars%ng2, stars%kv3, sym%mrot, sym%symor, &
                        sym%tau,rcc,sym%invtab,sf2)

            jp3 = (p(3)-cell%z1)/vacuum%delz
            delta = (p(3)-cell%z1)/vacuum%delz - jp3
            ! We count 0 as point 1.
            jp3 = jp3 + 1
            IF (jp3.LT.vacuum%nmz) THEN
               xdnout = potDen%vacz(jp3,ivac,jsp) + &
                        delta*(potDen%vacz(jp3+1,ivac,jsp) - &
                        potDen%vacz(jp3,ivac,jsp))
               IF (jp3.LT.vacuum%nmzxy) THEN
                  xx1 = 0.
                  xx2 = 0.
                  DO  k = 2,stars%ng2
                     xx1 = xx1 + REAL(potDen%vacxy(jp3,k-1,ivac,jsp)*sf2(k))* &
                           stars%nstr2(k)
                     xx2 = xx2 + REAL(potDen%vacxy(jp3+1,k-1,ivac,jsp)*sf2(k))* &
                           stars%nstr2(k)
                  ENDDO
                  xdnout = xdnout + xx1 + delta* (xx2-xx1)
               END IF
            ELSE
               xdnout = 0.0
134
            END IF
135 136
         ! Vacuum part finished.
         ENDIF
137

138
         RETURN
139
      ENDIF
140
      ! MT part:
141
      
142 143
      nd = sym%ntypsy(na)
      nopa = sym%ngopr(na)
144 145 146 147 148
      IF (oneD%odi%d1) nopa = oneD%ods%ngopr(na)
      sx = 0.0
      DO  i = 1,3
         x(i) = p(i) - atoms%pos(i,na)
         sx = sx + x(i)*x(i)
149
      END DO
150 151
      sx = sqrt(sx)
      IF (nopa.NE.1) THEN
152
         ! Switch to internal units.
153
         rcc=matmul(cell%bmat,x)/tpi_const
154
         ! Rotate into representative.
155 156 157
         DO  i = 1,3
            p(i) = 0.
            DO  j = 1,3
158 159 160 161 162 163 164 165
               IF (.NOT.oneD%odi%d1) THEN
                  p(i) = p(i) + sym%mrot(i,j,nopa)*rcc(j)
               ELSE
                  p(i) = p(i) + oneD%ods%mrot(i,j,nopa)*rcc(j)
               END IF
            END DO
         END DO
         ! Switch back to cartesian units.
166 167 168 169 170 171
         x=matmul(cell%amat,p)
      END IF
      DO j = atoms%jri(n),2,-1
         IF (sx.GE.atoms%rmsh(j,n)) EXIT
      ENDDO
      jr = j
172
      CALL ylm4(atoms%lmax(n),x,ylm)
173 174 175 176 177 178
      xd1 = 0.0
      xd2 = 0.0
      DO  lh = 0, sphhar%nlh(nd)
         ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
         s = 0.0
         DO mem = 1,sphhar%nmem(lh,nd)
179 180
            lm = ll1 + sphhar%mlh(mem,lh,nd)
            s = s + real( sphhar%clnu(mem,lh,nd)*ylm(lm) )
181
         ENDDO
182
         IF (l_potential) THEN
183
            xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s
184
         ELSE
185
            xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s/(atoms%rmsh(jr,n)**2)
186 187
         END IF
         IF (jr.EQ.atoms%jri(n)) CYCLE
188
         IF (l_potential) THEN
189
            xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s
190
         ELSE
191
            xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s/(atoms%rmsh(jr+1,n)**2)
192
         END IF
193
      ENDDO
194 195 196
      IF (jr.EQ.atoms%jri(n)) THEN
         xdnout = xd1
      ELSE
197 198
         xdnout = xd1 + (xd2-xd1) * (sx-atoms%rmsh(jr,n))/ &
                  (atoms%rmsh(jr+1,n)-atoms%rmsh(jr,n))
199
      END IF
200 201
      8000 FORMAT (2f10.6)

202
      RETURN
203 204
   END SUBROUTINE outcdn
END MODULE m_outcdn