outcdn.f90 6.87 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 20 21

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

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 98
      IF (iflag.NE.1) THEN
         IF (iflag.NE.0) THEN
            ! Interstitial part:
            !CALL cotra1(p(1),rcc,cell%bmat)
            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
99
            END IF
100 101
         
         ! +odim
102
         ELSE
103 104 105 106 107 108
            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))
109
            END IF
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
            !CALL cotra1(p,rcc,cell%bmat)
            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
136
            END IF
137 138
         ! Vacuum part finished.
         ENDIF
139

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

206
      RETURN
207 208
   END SUBROUTINE outcdn
END MODULE m_outcdn