wann_real.F 8.15 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14
      MODULE m_wann_real
c     ********************************************************
c     calculates the value of the periodic part of the
c     wavefunction at the given real-grid point p(:)
c                          Y.Mokrousov 16.8.6
c     ********************************************************
      CONTAINS
      SUBROUTINE wann_real(
15
     >                  p,n,na,iv,iflag,bkpt,iband,
16 17 18 19 20
     >                  n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
     >                  natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
     >                  nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
     >                  lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
     >                  omtil,amat,bmat,odi,ods,nlod,llod,nlo,llo,
21
     >                  ff,gg,flo,acof,bcof,ccof,zMat,
22
     >                  nv,k1,k2,k3,lmd,nbasfcn,l_ss,qss,jspin,addnoco,
23
     <                  xdnout)
24 25

      USE m_types
26
      USE m_ylm
27
      USE m_constants
28
      IMPLICIT NONE
29

30
      TYPE(t_mat),INTENT(IN)        :: zMat
31

32
C     .. Scalar Arguments ..
33
      INTEGER, INTENT (IN) :: n3d,nmzxyd,n2d,ntypsd,llod,nlod,iband
34 35 36
      INTEGER, INTENT (IN) :: lmaxd,jmtd,ntypd,natd,nmzd
      INTEGER, INTENT (IN) :: iflag,n,na,iv,lmd,nv,nvd,nbasfcn
      INTEGER, INTENT (IN) :: nq3,nvac,nmz,nmzxy,nq2,nop,nop2
37 38 39
      LOGICAL, INTENT (IN) :: invs,l_ss
      REAL,    INTENT (IN) :: z1,delz,omtil,bkpt(3),qss(3)
      INTEGER, INTENT (IN) :: jspin,addnoco
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
      COMPLEX, INTENT (OUT):: xdnout
c-odim
      TYPE (od_inp), INTENT (IN) :: odi
      TYPE (od_sym), INTENT (IN) :: ods
c+odim
C     ..
C     .. Array Arguments ..
      INTEGER, INTENT (IN) :: ngopr(natd),ntypsy(natd),jri(ntypd)
      INTEGER, INTENT (IN) :: lmax(ntypd),mrot(3,3,nop),invtab(nop)
      INTEGER, INTENT (IN) :: nlo(ntypd),llo(nlod,ntypd)
      REAL,    INTENT (IN) :: amat(3,3),bmat(3,3),pos(3,natd)
      REAL,    INTENT (IN) :: rmsh(jmtd,ntypd),tau(3,nop)
      INTEGER, INTENT (IN) :: k1(nvd),k2(nvd),k3(nvd) 
      COMPLEX, INTENT (IN) :: ccof(-llod:llod,nlod,natd)
      COMPLEX, INTENT (IN) :: acof(0:lmd,natd)
      COMPLEX, INTENT (IN) :: bcof(0:lmd,natd)
      REAL,    INTENT (IN) :: ff(ntypd,jmtd,2,0:lmaxd)
      REAL,    INTENT (IN) :: gg(ntypd,jmtd,2,0:lmaxd)
      REAL,    INTENT (IN) :: flo(ntypd,jmtd,2,nlod)
      REAL,    INTENT (INOUT) :: p(3)
60

61 62 63 64 65
C     ..
C     .. Local Scalars ..
      REAL delta,sx,xx1,xx2,rrr,phi,const,arg,tpi,arg1
      INTEGER i,j,jp3,jr,k,l,nd,nopa,ivac,lm,m,gzi,kk
      INTEGER kk1,kk2,kk3
66
      COMPLEX const2,s,xd1,xd2,const3
67 68 69 70
C     ..
C     .. Local Arrays ..
      COMPLEX sf2(n2d),sf3(n3d),ylm((lmaxd+1)**2)
      REAL rcc(3),x(3),rcc2(3)
71
      REAL bqpt(3)
72

73

74 75 76 77
      tpi   = 2 * pimach()
      const = 1./(sqrt(omtil))

c..define the factor e^{-ikr}
78
      rcc2=matmul(bmat,p)/tpi_const
79 80 81 82 83 84 85 86 87 88 89 90

      bqpt = 0.0
!      if(l_ss.and.jspin.eq.1) then 
!         bqpt = -qss/2.0
!      elseif(l_ss.and.jspin.eq.2) then
!         bqpt = +qss/2.0
!      endif

      arg = -tpi*(   (bkpt(1)+bqpt(1))*rcc2(1) 
     >             + (bkpt(2)+bqpt(2))*rcc2(2)
     >             + (bkpt(3)+bqpt(3))*rcc2(3)  )

91 92 93 94 95 96 97 98 99 100
      arg1 = tpi*( bkpt(1)*rcc2(1) + bkpt(2)*rcc2(2) + bkpt(3)*rcc2(3) )
      const2 = cmplx(cos(arg),sin(arg))
      const3 = cmplx(cos(arg1),sin(arg1))
c     write (6,*) 'bkpt,const2,const3=',bkpt(:),const2,const3

      ivac=iv

      IF (iflag.EQ.0) GO TO 20
      IF (iflag.EQ.1) GO TO 40
c     ---> interstitial part
101
      rcc=matmul(bmat,p)/tpi_const
102 103
      xdnout = cmplx(0.,0.)
c     write (6,*) 'nv,nvd=',nv,nvd
104 105 106 107 108
      IF (zMat%l_real) THEN
         DO k = 1,nv
c           write (6,*) 'k1,k2,k3=',k1(k),k2(k),k3(k)
c           write (6,*) 'z(k,iband)=', z(k,iband)
            arg = tpi * ((k1(k))*rcc(1)+(k2(k))*rcc(2)+(k3(k))*rcc(3))
109
            xdnout = xdnout + zMat%data_r(k+addnoco,iband)*
110 111 112 113 114 115 116 117 118 119 120 121 122 123
     +                        cmplx(cos(arg),sin(arg))*const
            IF (((abs(p(1)-2.2).le.0.0001).and.(abs(p(2)).le.0.0001))
     &    .or.((abs(p(2)-2.2).le.0.0001).and.(abs(p(1)).le.0.0001)))then
c              write (6,*) 'p(i)=',p(1:2)
c              write (6,*) 'G=',k1(k),k2(k),k3(k)
c              write (6,*) 'z(k,iband)=',z(k,iband)
c              write (6,*) 'val=',z(k,iband)*cmplx(cos(arg),sin(arg))
            ENDIF
         END DO
      ELSE
         DO k = 1,nv
c           write (6,*) 'k1,k2,k3=',k1(k),k2(k),k3(k)
c           write (6,*) 'z(k,iband)=', z(k,iband)
            arg = tpi * ((k1(k))*rcc(1)+(k2(k))*rcc(2)+(k3(k))*rcc(3))
124
            xdnout = xdnout + zMat%data_c(k+addnoco,iband)*
125 126 127 128 129 130 131 132 133 134
     +                        cmplx(cos(arg),sin(arg))*const
            IF (((abs(p(1)-2.2).le.0.0001).and.(abs(p(2)).le.0.0001))
     &    .or.((abs(p(2)-2.2).le.0.0001).and.(abs(p(1)).le.0.0001)))then
c              write (6,*) 'p(i)=',p(1:2)
c              write (6,*) 'G=',k1(k),k2(k),k3(k)
c              write (6,*) 'z(k,iband)=',z(k,iband)
c              write (6,*) 'val=',z(k,iband)*cmplx(cos(arg),sin(arg))
            ENDIF
         END DO
      END IF
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
c     write (6,*) 'ir:p(i)',p(:)
      RETURN
c     ---> vacuum part
   20 CONTINUE
      xdnout = cmplx(0.,0.)
      return
c     ----> m.t. part
   40 CONTINUE
      

      nd = ntypsy(na)
      nopa = ngopr(na)
      nopa=1
     
      IF (odi%d1) nopa = ods%ngopr(na)
      sx = 0.0
      DO 50 i = 1,3
         x(i) = p(i) - pos(i,na)
         sx = sx + x(i)*x(i)
   50 CONTINUE
      sx = sqrt(sx)
      IF (nopa.NE.1) THEN
c... switch to internal units
158
         rcc=matmul(bmat,p)/tpi_const
159 160 161 162 163 164 165 166 167 168 169 170
c... rotate into representative
         DO 70 i = 1,3
            p(i) = 0.
            DO 60 j = 1,3
              IF (.NOT.odi%d1) THEN
               p(i) = p(i) + mrot(i,j,nopa)*rcc(j)
              ELSE
               p(i) = p(i) + ods%mrot(i,j,nopa)*rcc(j)
              END IF
   60       CONTINUE
   70    CONTINUE
c... switch back to cartesian units
171
         x=matmul(amat,p)/tpi_const
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
      END IF
      DO 80 j = jri(n),2,-1
         IF (sx.GE.rmsh(j,n)) GO TO 90
   80 CONTINUE
   90 jr = j
      CALL ylm4(
     >          lmax(n),x,
     <          ylm)
      xd1 = cmplx(0.,0.)
      xd2 = cmplx(0.,0.)
      DO l = 0,lmax(n)
c        if (p(1).eq.0. .and. p(2).eq.0. .and. p(3).eq.0)then
c               write (6,*) 'ff(l,300)=',ff(1,300,1,l)
c               write (6,*) 'ff(l,300)=',ff(1,300,2,l)
c               write (6,*) 'gg(l,300)=',gg(1,300,1,l)
c               write (6,*) 'gg(l,300)=',gg(1,300,2,l)
c        endif
       DO 110 m = -l,l
        lm = l*(l+1)+m
191
        s = ylm(lm+1)*(ImagUnit)**l
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
c       if (p(1).eq.0. .and. p(2).eq.0. .and. p(3).eq.0)then
c              write (6,*) 'acof=',acof(lm,1)
c              write (6,*) 'bcof=',bcof(lm,1)
c       endif
        xd1 = xd1 + (acof(lm,na)*cmplx(ff(n,jr,1,l),0.)+
     +               bcof(lm,na)*cmplx(gg(n,jr,1,l),0.))*s/
     /               (rmsh(jr,n)) 
c    /               (rmsh(jr,n)*rmsh(jr,n))
        IF (jr.EQ.1) GO TO 110
        xd2 = xd2 + (acof(lm,na)*cmplx(ff(n,jr+1,1,l),0.)+
     +               bcof(lm,na)*cmplx(gg(n,jr+1,1,l),0.))*s/  
     /               (rmsh(jr+1,n))
c    /               (rmsh(jr+1,n)*rmsh(jr+1,n))
  110  CONTINUE
      ENDDO
c..contributions from the local orbitals
      IF (nlo(n).GE.1) THEN
       DO l = 1,nlo(n)
        DO 111 m = -llo(l,n),llo(l,n)
         lm = llo(l,n)*(llo(l,n)+1)+m
212
         s = ylm(lm+1)*(ImagUnit)**llo(l,n) 
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
         xd1 = xd1 + ccof(m,l,na)*flo(n,jr,1,l)*s/
     /               (rmsh(jr,n))         
         IF (jr.EQ.1) GO TO 111
         xd2 = xd2 + ccof(m,l,na)*flo(n,jr+1,1,l)*s/
     /               (rmsh(jr+1,n))         
  111   CONTINUE
       ENDDO
      ENDIF    
      IF (jr.EQ.1) THEN
         xdnout = xd1
      ELSE
         xdnout = xd1 + (xd2-xd1) *
     +                  (sx-rmsh(jr,n)) / (rmsh(jr+1,n)-rmsh(jr,n))
         
      END IF
      xdnout = xdnout*const2
c     write (6,*) 'mt:p(i)',p(:)
 8000 FORMAT (2f10.6)
c
      RETURN
      END SUBROUTINE wann_real
      END MODULE m_wann_real