LCOV - code coverage report
Current view: top level - wannier - wann_nocoplot.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 137 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       MODULE m_wann_nocoplot
       2             :       USE m_juDFT
       3             : c     ++++++++++++++++++++++++++++++++++++++++++++++++
       4             : c     +  If noco=t, transform wave functions within  +
       5             : c     +  mu-th MT back to the global frame to plot   +
       6             : c     +  to plot WFs in a meaningful way.            +
       7             : c     +                                              +
       8             : c     ++++++++++++++++++++++++++++++++++++++++++++++++
       9             :       CONTAINS
      10           0 :       SUBROUTINE wann_nocoplot(atoms,slice,nnne,amat,bmat
      11             :      >                        ,nkpts,odi,film,natd,ntypd,jmtd
      12           0 :      >                        ,ntype,neq,pos,jri,rmsh,alph,beta
      13           0 :      >                        ,nqpts,qss,z1,zatom)
      14             : !    *****************************************************
      15             :       USE m_types
      16             :       USE m_xsf_io
      17             :       USE m_wann_gwf_tools, ONLY : get_index_q
      18             :       USE m_constants
      19             : 
      20             :       IMPLICIT NONE
      21             : 
      22             :       TYPE(t_atoms),INTENT(IN):: atoms
      23             : 
      24             :       LOGICAL, INTENT(IN) :: slice,film
      25             :       INTEGER, INTENT(IN) :: nnne,nkpts,ntypd,jmtd,natd
      26             :       INTEGER, INTENT(IN) :: ntype,neq(ntypd)
      27             :       INTEGER, INTENT(IN) :: jri(ntypd),nqpts
      28             :       REAL,    INTENT(IN) :: amat(3,3),bmat(3,3),qss(3,nqpts)
      29             :       REAL,    INTENT(IN) :: rmsh(jmtd,ntypd)
      30             :       REAL,    INTENT(IN) :: pos(3,natd),z1,zatom(:)
      31             :       REAL,    INTENT(IN) :: alph(ntypd),beta(ntypd)
      32             : 
      33             :       TYPE(od_inp), INTENT(in) :: odi
      34             : 
      35             :       LOGICAL :: twodim,xsf,cartesian
      36             :       INTEGER :: nbmin,nbmax,nplot,nplo,nbn
      37             :       INTEGER :: nslibd,nkqpts
      38             :       INTEGER :: ix,iy,iz,ikpt,jspin
      39             :       INTEGER :: ii1,ii2,ii3
      40             :       INTEGER :: i1,i2,i3,iintsp
      41             :       INTEGER :: gx,gy,gz,help_kpt
      42             :       INTEGER :: na,nt,nq,iqpt
      43             :       INTEGER :: grid(3),count_mt,count_int,count_vac
      44             :       REAL :: vec1(3),vec2(3),vec3(3),zero(3)
      45             :       REAL :: pt(3),point(3),rcc2(3)
      46             :       REAL :: s,arg,pi,u_r,u_i
      47             :       COMPLEX :: phasfac
      48             :       COMPLEX :: U(2,2)
      49             :       COMPLEX :: wf_local(2),wf_global,xdnout
      50             :       CHARACTER(len=30) :: filename
      51             :       CHARACTER(len=20) :: vandername,name1,name2,name3
      52             : 
      53             :       NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
      54             :     
      55             :       intrinsic real,aimag,conjg,exp,cmplx
      56             : 
      57           0 :       pi = pimach()
      58           0 :       nkqpts = nkpts*nqpts
      59             : 
      60           0 :       INQUIRE(file ="plot_inp",exist=twodim)
      61           0 :       IF(.NOT.twodim) THEN
      62             :          CALL juDFT_error("Need the plot_inp from RNK generation!"
      63           0 :      >                  ,calledby="wann_nocoplot")
      64             :       ENDIF
      65             : 
      66             :       !<-- Open the plot_inp file for input
      67           0 :       OPEN (18,file='plot_inp')
      68           0 :       READ(18,'(i2,5x,l1)') nplot,xsf
      69             : 
      70           0 :       IF (nplot.ge.2) 
      71             :      &     CALL juDFT_error
      72             :      +     ("plots one by one, please, this is not charge density"
      73           0 :      +     ,calledby="wann_nocoplot")
      74             : 
      75             :          ! the defaults
      76           0 :          twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
      77           0 :          vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
      78           0 :          zero = (/0.,0.,0./);filename="default"
      79           0 :          READ(18,plot)
      80           0 :          IF (twodim.AND.ANY(grid(1:2)<1)) 
      81             :      +        CALL juDFT_error("Illegal grid size in plot",calledby
      82           0 :      +        ="wann_nocoplot")
      83           0 :          IF (.NOT.twodim.AND.ANY(grid<1)) 
      84             :      +        CALL juDFT_error("Illegal grid size in plot",calledby
      85           0 :      +        ="wann_nocoplot")
      86           0 :          IF (twodim) grid(3) = 1
      87             :          !calculate cartesian coordinates if needed
      88           0 :          IF (.NOT.cartesian) THEN
      89           0 :             vec1=matmul(amat,vec1)
      90           0 :             vec2=matmul(amat,vec2)
      91           0 :             vec3=matmul(amat,vec3)
      92           0 :             zero=matmul(amat,zero)
      93             :          ENDIF
      94           0 :          IF (filename =="default") WRITE(filename,'(a,i2)') "plot",1
      95           0 :       CLOSE(18)
      96             : 
      97             :          ! loop over k-points
      98           0 :          DO ikpt=1,nkqpts
      99           0 :             count_mt=0; count_int=0; count_vac=0
     100           0 :             iqpt = get_index_q(ikpt,nkpts)
     101           0 :             write(*,*)'kq',ikpt,' q',iqpt,' qz',qss(3,iqpt)
     102             : 
     103             :             ! open the old RNK.ikpt.jspin files
     104           0 :             DO jspin=1,2  ! local frame axis (inside MT)
     105           0 :                WRITE(vandername,202) ikpt,jspin
     106           0 :                OPEN(400+jspin,file=vandername)
     107           0 :                READ(400+jspin,7)gx,gy,gz,help_kpt,nslibd
     108             :             ENDDO 
     109             : 
     110           0 :             IF(.not.xsf) THEN
     111             :             ! open the new UNK.ikpt.iintsp files
     112           0 :             DO iintsp=1,2  ! global frame axis (INT and vacuum)              
     113           0 :                WRITE(vandername,201) ikpt,iintsp
     114           0 :                OPEN(500+iintsp,file=vandername,status='unknown')
     115           0 :                WRITE(500+iintsp,7)grid(1),grid(2),grid(3),ikpt,nslibd
     116             :             ENDDO
     117             :             ENDIF
     118             : 
     119             :          ! loop over all bands
     120           0 :                nbmin=1
     121           0 :                nbmax=nslibd
     122           0 :          bands:DO nbn = nbmin,nbmax
     123             : 
     124           0 :          IF (xsf) THEN
     125           0 :            do jspin=1,2
     126           0 :             write (name1,22) ikpt,nbn,jspin
     127             :    22       format (i5.5,'.',i3.3,'.real.',i1,'.xsf')
     128           0 :             write (name2,23) ikpt,nbn,jspin
     129             :    23       format (i5.5,'.',i3.3,'.imag.',i1,'.xsf')
     130           0 :             write (name3,24) ikpt,nbn,jspin
     131             :    24       format (i5.5,'.',i3.3,'.absv.',i1,'.xsf')
     132           0 :             OPEN(600+jspin,file=name1)
     133           0 :             CALL xsf_WRITE_atoms(600+jspin,atoms,film,odi%d1,amat)
     134           0 :             OPEN(602+jspin,file=name2)
     135           0 :             CALL xsf_WRITE_atoms(602+jspin,atoms,film,odi%d1,amat)
     136           0 :             OPEN(604+jspin,file=name3)
     137           0 :             CALL xsf_WRITE_atoms(604+jspin,atoms,film,odi%d1,amat)
     138             :             CALL xsf_WRITE_header(600+jspin,twodim,filename,(vec1),
     139             :      &       (vec2),(vec3),zero
     140           0 :      $           ,grid)
     141             :             CALL xsf_WRITE_header(602+jspin,twodim,filename,(vec1),
     142             :      &       (vec2),(vec3),zero
     143           0 :      $           ,grid)
     144             :             CALL xsf_WRITE_header(604+jspin,twodim,filename,(vec1),
     145             :      &       (vec2),(vec3),zero
     146           0 :      $           ,grid)
     147             :             enddo
     148             :          ENDIF
     149             : 
     150             :          ! loop over real space grid points
     151           0 :          DO iz = 0,grid(3)-1
     152           0 :           DO iy = 0,grid(2)-1
     153           0 :            xloop:DO ix = 0,grid(1)-1
     154             :             point = zero+vec1*REAL(ix)/grid(1)+vec2*REAL(iy)
     155           0 :      $                 /grid(2)
     156           0 :             IF (.NOT.twodim) point = point+vec3*REAL(iz)/grid(3)
     157             : 
     158             :             ! read old RNK information at grid point
     159           0 :             DO jspin=1,2
     160           0 :                READ(400+jspin,8)u_r,u_i
     161           0 :                wf_local(jspin) = cmplx(u_r,u_i)
     162             :             ENDDO
     163             : 
     164             :              ! is point in MT?
     165           0 :              ii1 = 3
     166           0 :              ii2 = 3
     167           0 :              ii3 = 3
     168           0 :              IF (film .AND. .NOT.odi%d1) ii3 = 0
     169           0 :              IF (odi%d1) THEN
     170           0 :                 ii1 = 0 ; ii2 = 0
     171             :              END IF
     172           0 :              DO  i1 = -ii1,ii1
     173           0 :               DO  i2 = -ii2,ii2
     174           0 :                DO  i3 = -ii3,ii3
     175           0 :                 pt = point+MATMUL(amat,(/i1,i2,i3/))
     176           0 :                 na = 0
     177           0 :                 DO nt = 1,ntype
     178           0 :                  DO nq = 1,neq(nt)
     179           0 :                   na   = na + 1
     180           0 :                   s  = SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))
     181           0 :                   IF (s<rmsh(jri(nt),nt)) THEN
     182           0 :                     count_mt=count_mt+1
     183             :                     ! we are inside the MT with alph(nt),beta(nt)
     184             :                     ! set up transformation local -> global
     185           0 :                    U(1,1) =  exp(-ImagUnit*alph(nt)/2.)*cos(beta(nt)/2.) 
     186           0 :                    U(1,2) = -exp(-ImagUnit*alph(nt)/2.)*sin(beta(nt)/2.)
     187           0 :                    U(2,1) =  exp( ImagUnit*alph(nt)/2.)*sin(beta(nt)/2.)
     188           0 :                    U(2,2) =  exp( ImagUnit*alph(nt)/2.)*cos(beta(nt)/2.)
     189             :                     
     190             :                     ! transform wfs to global frame
     191           0 :                     DO iintsp=1,2
     192           0 :                        pt = matmul(bmat,rcc2)/tpi_const
     193             : !                       CALL cotra1(pt(1),rcc2,bmat)
     194             :                        arg = -pi*real(2*iintsp-3)*( qss(1,iqpt)*rcc2(1)
     195             :      >                                             +qss(2,iqpt)*rcc2(2)
     196           0 :      >                                             +qss(3,iqpt)*rcc2(3))
     197           0 :                        phasfac = cmplx(cos(arg),sin(arg))
     198             : 
     199             :                        wf_global= phasfac*(  U(iintsp,1)*wf_local(1)
     200           0 :      >                                     + U(iintsp,2)*wf_local(2) )
     201             : 
     202           0 :                        if(xsf) THEN
     203           0 :                           xdnout=wf_global
     204           0 :                           WRITE(600+iintsp,*) real(xdnout)
     205           0 :                           WRITE(602+iintsp,*) aimag(xdnout)
     206           0 :                           WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))
     207             :                        ELSE
     208           0 :                           WRITE(500+iintsp,8)real(wf_global),
     209           0 :      >                                       aimag(wf_global)
     210             :                        ENDIF
     211             :                     ENDDO
     212             : 
     213             :                    CYCLE xloop
     214             :                   ENDIF
     215             :                  ENDDO
     216             :                 ENDDO !nt
     217             :                ENDDO
     218             :               ENDDO
     219             :              ENDDO !i1
     220             : 
     221             :              ! VACUUM region
     222           0 :              IF (SQRT((pt(1))**2+(pt(2))**2)>=z1)THEN
     223           0 :              count_vac=count_vac+1
     224           0 :              DO iintsp=1,2
     225           0 :                 phasfac=cmplx(1.,0.)
     226           0 :                 wf_global = phasfac*wf_local(iintsp)
     227             : 
     228           0 :                 if(xsf) THEN
     229           0 :                    xdnout=wf_global
     230           0 :                    WRITE(600+iintsp,*) real(xdnout)
     231           0 :                    WRITE(602+iintsp,*) aimag(xdnout)
     232           0 :                    WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))      
     233             :                 ELSE
     234           0 :                    WRITE(500+iintsp,8)real(wf_global),
     235           0 :      >                  aimag(wf_global)
     236             :                 ENDIF
     237             :              ENDDO
     238             :  
     239             :              CYCLE xloop
     240             :              ENDIF
     241             : 
     242             :              ! if we are here, point is in INTERSTITIAL
     243             :              ! therefore just copy wfs
     244           0 :              count_int = count_int+1
     245           0 :              DO iintsp=1,2
     246           0 :                 phasfac=cmplx(1.,0.)
     247           0 :                 wf_global = phasfac*wf_local(iintsp)
     248             : 
     249           0 :                 if(xsf) THEN
     250           0 :                    xdnout=wf_global
     251           0 :                    WRITE(600+iintsp,*) real(xdnout)
     252           0 :                    WRITE(602+iintsp,*) aimag(xdnout)
     253           0 :                    WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))      
     254             :                 ELSE
     255           0 :                    WRITE(500+iintsp,8)real(wf_global),
     256           0 :      >                  aimag(wf_global)
     257             :                 ENDIF
     258             :              ENDDO
     259             : 
     260             :             ENDDO xloop
     261             :            ENDDO
     262             :           ENDDO !z-loop
     263             : 
     264           0 :           IF (xsf) THEN    
     265           0 :              DO iintsp=1,2
     266           0 :               CALL xsf_WRITE_endblock(600+iintsp,twodim)
     267           0 :               CALL xsf_WRITE_endblock(602+iintsp,twodim)
     268           0 :               CALL xsf_WRITE_endblock(604+iintsp,twodim)
     269           0 :               CLOSE(600+iintsp); CLOSE(602+iintsp); CLOSE(604+iintsp)
     270             :              ENDDO
     271             :           ENDIF
     272             : 
     273             :                ENDDO bands
     274             : 
     275           0 :                IF(.not.xsf) THEN
     276             :                ! close new UNK.ikpt.iintsp files
     277           0 :                DO iintsp=1,2
     278           0 :                   CLOSE(500+iintsp)
     279             :                ENDDO
     280             :                ENDIF
     281             : 
     282             :                ! delete RNK.ikpt.jspin files
     283           0 :                DO jspin=1,2   
     284           0 :                   IF(xsf) THEN
     285           0 :                      CLOSE(400+jspin)
     286             :                   ELSE
     287           0 :                      CLOSE(400+jspin,status='delete')
     288             :                   ENDIF
     289             :                ENDDO
     290             : 
     291           0 :                write(*,*)count_mt,count_int,count_vac
     292             : 
     293             :          ENDDO  ! end ikpt loop
     294             : 
     295             :   201 FORMAT ('UNK',i5.5,'.',i1)
     296             :   202 FORMAT ('RNK',i5.5,'.',i1)
     297             :     7 FORMAT (5i4)
     298             :     8 FORMAT (f20.12,1x,f20.12)
     299             : 
     300           0 :       END SUBROUTINE wann_nocoplot
     301             : !------------------------------------------
     302             :       
     303             :       END MODULE m_wann_nocoplot

Generated by: LCOV version 1.13