wann_nocoplot.F 10.6 KB
Newer Older
1
      MODULE m_wann_nocoplot
2
      USE m_juDFT
3 4 5 6 7 8 9
c     ++++++++++++++++++++++++++++++++++++++++++++++++
c     +  If noco=t, transform wave functions within  +
c     +  mu-th MT back to the global frame to plot   +
c     +  to plot WFs in a meaningful way.            +
c     +                                              +
c     ++++++++++++++++++++++++++++++++++++++++++++++++
      CONTAINS
10
      SUBROUTINE wann_nocoplot(atoms,slice,nnne,amat,bmat
11 12 13 14
     >                        ,nkpts,odi,film,natd,ntypd,jmtd
     >                        ,ntype,neq,pos,jri,rmsh,alph,beta
     >                        ,nqpts,qss,z1,zatom)
!    *****************************************************
15
      USE m_types
16 17
      USE m_xsf_io
      USE m_wann_gwf_tools, ONLY : get_index_q
18
      USE m_constants
19 20

      IMPLICIT NONE
21 22 23

      TYPE(t_atoms),INTENT(IN):: atoms

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
      LOGICAL, INTENT(IN) :: slice,film
      INTEGER, INTENT(IN) :: nnne,nkpts,ntypd,jmtd,natd
      INTEGER, INTENT(IN) :: ntype,neq(ntypd)
      INTEGER, INTENT(IN) :: jri(ntypd),nqpts
      REAL,    INTENT(IN) :: amat(3,3),bmat(3,3),qss(3,nqpts)
      REAL,    INTENT(IN) :: rmsh(jmtd,ntypd)
      REAL,    INTENT(IN) :: pos(3,natd),z1,zatom(:)
      REAL,    INTENT(IN) :: alph(ntypd),beta(ntypd)

      TYPE(od_inp), INTENT(in) :: odi

      LOGICAL :: twodim,xsf,cartesian
      INTEGER :: nbmin,nbmax,nplot,nplo,nbn
      INTEGER :: nslibd,nkqpts
      INTEGER :: ix,iy,iz,ikpt,jspin
      INTEGER :: ii1,ii2,ii3
      INTEGER :: i1,i2,i3,iintsp
      INTEGER :: gx,gy,gz,help_kpt
      INTEGER :: na,nt,nq,iqpt
      INTEGER :: grid(3),count_mt,count_int,count_vac
      REAL :: vec1(3),vec2(3),vec3(3),zero(3)
      REAL :: pt(3),point(3),rcc2(3)
      REAL :: s,arg,pi,u_r,u_i
      COMPLEX :: ci,phasfac
      COMPLEX :: U(2,2)
      COMPLEX :: wf_local(2),wf_global,xdnout
      CHARACTER(len=30) :: filename
      CHARACTER(len=20) :: vandername,name1,name2,name3

      NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
    
      intrinsic real,aimag,conjg,exp,cmplx

      pi = pimach()
      ci = cmplx(0.0,1.0)
      nkqpts = nkpts*nqpts

      INQUIRE(file ="plot_inp",exist=twodim)
      IF(.NOT.twodim) THEN
63
         CALL juDFT_error("Need the plot_inp from RNK generation!"
64 65 66 67 68 69 70 71
     >                  ,calledby="wann_nocoplot")
      ENDIF

      !<-- Open the plot_inp file for input
      OPEN (18,file='plot_inp')
      READ(18,'(i2,5x,l1)') nplot,xsf

      IF (nplot.ge.2) 
72
     &     CALL juDFT_error
73 74 75 76 77 78 79 80 81
     +     ("plots one by one, please, this is not charge density"
     +     ,calledby="wann_nocoplot")

         ! the defaults
         twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
         vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
         zero = (/0.,0.,0./);filename="default"
         READ(18,plot)
         IF (twodim.AND.ANY(grid(1:2)<1)) 
82
     +        CALL juDFT_error("Illegal grid size in plot",calledby
83 84
     +        ="wann_nocoplot")
         IF (.NOT.twodim.AND.ANY(grid<1)) 
85
     +        CALL juDFT_error("Illegal grid size in plot",calledby
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
     +        ="wann_nocoplot")
         IF (twodim) grid(3) = 1
         !calculate cartesian coordinates if needed
         IF (.NOT.cartesian) THEN
            vec1=matmul(amat,vec1)
            vec2=matmul(amat,vec2)
            vec3=matmul(amat,vec3)
            zero=matmul(amat,zero)
         ENDIF
         IF (filename =="default") WRITE(filename,'(a,i2)') "plot",1
      CLOSE(18)

         ! loop over k-points
         DO ikpt=1,nkqpts
            count_mt=0; count_int=0; count_vac=0
            iqpt = get_index_q(ikpt,nkpts)
            write(*,*)'kq',ikpt,' q',iqpt,' qz',qss(3,iqpt)

            ! open the old RNK.ikpt.jspin files
            DO jspin=1,2  ! local frame axis (inside MT)
               WRITE(vandername,202) ikpt,jspin
               OPEN(400+jspin,file=vandername)
               READ(400+jspin,7)gx,gy,gz,help_kpt,nslibd
            ENDDO 

            IF(.not.xsf) THEN
            ! open the new UNK.ikpt.iintsp files
            DO iintsp=1,2  ! global frame axis (INT and vacuum)              
               WRITE(vandername,201) ikpt,iintsp
               OPEN(500+iintsp,file=vandername,status='unknown')
               WRITE(500+iintsp,7)grid(1),grid(2),grid(3),ikpt,nslibd
            ENDDO
            ENDIF

         ! loop over all bands
               nbmin=1
               nbmax=nslibd
         bands:DO nbn = nbmin,nbmax

         IF (xsf) THEN
           do jspin=1,2
            write (name1,22) ikpt,nbn,jspin
   22       format (i5.5,'.',i3.3,'.real.',i1,'.xsf')
            write (name2,23) ikpt,nbn,jspin
   23       format (i5.5,'.',i3.3,'.imag.',i1,'.xsf')
            write (name3,24) ikpt,nbn,jspin
   24       format (i5.5,'.',i3.3,'.absv.',i1,'.xsf')
            OPEN(600+jspin,file=name1)
134
            CALL xsf_WRITE_atoms(600+jspin,atoms,film,odi%d1,amat)
135
            OPEN(602+jspin,file=name2)
136
            CALL xsf_WRITE_atoms(602+jspin,atoms,film,odi%d1,amat)
137
            OPEN(604+jspin,file=name3)
138
            CALL xsf_WRITE_atoms(604+jspin,atoms,film,odi%d1,amat)
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
            CALL xsf_WRITE_header(600+jspin,twodim,filename,(vec1),
     &       (vec2),(vec3),zero
     $           ,grid)
            CALL xsf_WRITE_header(602+jspin,twodim,filename,(vec1),
     &       (vec2),(vec3),zero
     $           ,grid)
            CALL xsf_WRITE_header(604+jspin,twodim,filename,(vec1),
     &       (vec2),(vec3),zero
     $           ,grid)
            enddo
         ENDIF

         ! loop over real space grid points
         DO iz = 0,grid(3)-1
          DO iy = 0,grid(2)-1
           xloop:DO ix = 0,grid(1)-1
            point = zero+vec1*REAL(ix)/grid(1)+vec2*REAL(iy)
     $                 /grid(2)
            IF (.NOT.twodim) point = point+vec3*REAL(iz)/grid(3)

            ! read old RNK information at grid point
            DO jspin=1,2
               READ(400+jspin,8)u_r,u_i
               wf_local(jspin) = cmplx(u_r,u_i)
            ENDDO

             ! is point in MT?
             ii1 = 3
             ii2 = 3
             ii3 = 3
             IF (film .AND. .NOT.odi%d1) ii3 = 0
             IF (odi%d1) THEN
                ii1 = 0 ; ii2 = 0
             END IF
             DO  i1 = -ii1,ii1
              DO  i2 = -ii2,ii2
               DO  i3 = -ii3,ii3
                pt = point+MATMUL(amat,(/i1,i2,i3/))
                na = 0
                DO nt = 1,ntype
                 DO nq = 1,neq(nt)
                  na   = na + 1
                  s  = SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))
                  IF (s<rmsh(jri(nt),nt)) THEN
                    count_mt=count_mt+1
                    ! we are inside the MT with alph(nt),beta(nt)
                    ! set up transformation local -> global
                    U(1,1) =  exp(-ci*alph(nt)/2.)*cos(beta(nt)/2.) 
                    U(1,2) = -exp(-ci*alph(nt)/2.)*sin(beta(nt)/2.)
                    U(2,1) =  exp( ci*alph(nt)/2.)*sin(beta(nt)/2.)
                    U(2,2) =  exp( ci*alph(nt)/2.)*cos(beta(nt)/2.)
                    
                    ! transform wfs to global frame
                    DO iintsp=1,2
193 194
                       pt = matmul(bmat,rcc2)/tpi_const
!                       CALL cotra1(pt(1),rcc2,bmat)
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
                       arg = -pi*real(2*iintsp-3)*( qss(1,iqpt)*rcc2(1)
     >                                             +qss(2,iqpt)*rcc2(2)
     >                                             +qss(3,iqpt)*rcc2(3))
                       phasfac = cmplx(cos(arg),sin(arg))

                       wf_global= phasfac*(  U(iintsp,1)*wf_local(1)
     >                                     + U(iintsp,2)*wf_local(2) )

                       if(xsf) THEN
                          xdnout=wf_global
                          WRITE(600+iintsp,*) real(xdnout)
                          WRITE(602+iintsp,*) aimag(xdnout)
                          WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))
                       ELSE
                          WRITE(500+iintsp,8)real(wf_global),
     >                                       aimag(wf_global)
                       ENDIF
                    ENDDO

                   CYCLE xloop
                  ENDIF
                 ENDDO
                ENDDO !nt
               ENDDO
              ENDDO
             ENDDO !i1

             ! VACUUM region
             IF (SQRT((pt(1))**2+(pt(2))**2)>=z1)THEN
             count_vac=count_vac+1
             DO iintsp=1,2
                phasfac=cmplx(1.,0.)
                wf_global = phasfac*wf_local(iintsp)

                if(xsf) THEN
                   xdnout=wf_global
                   WRITE(600+iintsp,*) real(xdnout)
                   WRITE(602+iintsp,*) aimag(xdnout)
                   WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))      
                ELSE
                   WRITE(500+iintsp,8)real(wf_global),
     >                  aimag(wf_global)
                ENDIF
             ENDDO
 
             CYCLE xloop
             ENDIF

             ! if we are here, point is in INTERSTITIAL
             ! therefore just copy wfs
             count_int = count_int+1
             DO iintsp=1,2
                phasfac=cmplx(1.,0.)
                wf_global = phasfac*wf_local(iintsp)

                if(xsf) THEN
                   xdnout=wf_global
                   WRITE(600+iintsp,*) real(xdnout)
                   WRITE(602+iintsp,*) aimag(xdnout)
                   WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))      
                ELSE
                   WRITE(500+iintsp,8)real(wf_global),
     >                  aimag(wf_global)
                ENDIF
             ENDDO

            ENDDO xloop
           ENDDO
          ENDDO !z-loop

          IF (xsf) THEN    
             DO iintsp=1,2
              CALL xsf_WRITE_endblock(600+iintsp,twodim)
              CALL xsf_WRITE_endblock(602+iintsp,twodim)
              CALL xsf_WRITE_endblock(604+iintsp,twodim)
              CLOSE(600+iintsp); CLOSE(602+iintsp); CLOSE(604+iintsp)
             ENDDO
          ENDIF

               ENDDO bands

               IF(.not.xsf) THEN
               ! close new UNK.ikpt.iintsp files
               DO iintsp=1,2
                  CLOSE(500+iintsp)
               ENDDO
               ENDIF

               ! delete RNK.ikpt.jspin files
               DO jspin=1,2   
                  IF(xsf) THEN
                     CLOSE(400+jspin)
                  ELSE
                     CLOSE(400+jspin,status='delete')
                  ENDIF
               ENDDO

               write(*,*)count_mt,count_int,count_vac

         ENDDO  ! end ikpt loop

  201 FORMAT ('UNK',i5.5,'.',i1)
  202 FORMAT ('RNK',i5.5,'.',i1)
    7 FORMAT (5i4)
    8 FORMAT (f20.12,1x,f20.12)

      END SUBROUTINE wann_nocoplot
!------------------------------------------
      
      END MODULE m_wann_nocoplot