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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_wann_plot_from_lapw
       8             :       use m_juDFT
       9             : 
      10             : c****************************************************************
      11             : c     read in WF1.lapw (WF2.lapw) and plot wannierfunctions
      12             : c     Frank Freimuth, October2006
      13             : c****************************************************************
      14             :       CONTAINS
      15             : 
      16           0 :       SUBROUTINE wann_plot_from_lapw(nv2d,jspins,odi,ods,n3d,nmzxyd,n2d,
      17             :      >     ntypsd,
      18           0 :      >     ntype,lmaxd,jmtd,natd,nmzd,neq,nq3,nvac,
      19             :      >     nmz,nmzxy,nq2,nop,nop2,volint,film,slice,symor,
      20           0 :      >     invs,invs2,z1,delz,ngopr,ntypsy,jri,pos,zatom,
      21           0 :      >     lmax,mrot,tau,rmsh,invtab,amat,bmat,bbmat,nnne,kk,
      22           0 :      >     nlod,llod,lmd,omtil,nlo,llo)
      23             : 
      24             :       USE m_types, ONLY : od_inp, od_sym
      25             :       USE m_xsf_io
      26             :       USE m_wann_lapw_int_plot
      27             :       USE m_wann_lapw_sph_plot
      28             :  
      29             :       implicit none
      30             : !     .. Scalar Arguments ..
      31             :       INTEGER, INTENT (IN) :: n3d,nmzxyd,n2d,ntypsd,nv2d
      32             :       INTEGER, INTENT (IN) :: lmaxd,jmtd,natd,nmzd
      33             :       INTEGER, INTENT (IN) :: nq3,nvac,nmz,nmzxy,nq2,nop,nop2,ntype
      34             :       INTEGER, INTENT (IN) :: lmd,llod,nlod
      35             :       INTEGER, INTENT (IN) :: nnne,kk,jspins
      36             :       LOGICAL, INTENT (IN) :: symor,invs,slice,invs2,film
      37             :       REAL,    INTENT (IN) :: z1,delz,volint,omtil
      38             : !     ..
      39             : !     .. Array Arguments ..
      40             :       INTEGER, INTENT (IN) :: ngopr(natd),ntypsy(natd),lmax(ntype)
      41             :       INTEGER, INTENT (IN) :: jri(ntype),neq(ntype),mrot(3,3,nop)
      42             :       INTEGER, INTENT (IN) :: invtab(nop),nlo(ntype),llo(nlod,ntype)
      43             :       REAL,    INTENT (IN) :: zatom(:),amat(3,3),bmat(3,3),pos(3,natd)
      44             :       REAL,    INTENT (IN) :: rmsh(jmtd,ntype),tau(3,nop)
      45             :       REAL,    INTENT (IN) :: bbmat(3,3)
      46             :       TYPE (od_inp), INTENT (IN) :: odi
      47             :       TYPE (od_sym), INTENT (IN) :: ods
      48             :       integer nplo,nplot,nbn,nbmin,nbmax,ix,iy,iz,ii1,ii2,ii3,i1,i2,i3
      49             :       integer nt,na,nq,ivac,jvac
      50             :       real amat_old(3,3)
      51           0 :       real pos_old(3,natd)
      52             :       real s
      53             :       logical l_file
      54             :       integer jspin
      55             :       integer num_wann,cell1,cell2,cell3
      56             :       character(len=3) :: spin12(2)
      57             :       data spin12/'WF1' , 'WF2'/
      58           0 :       complex,allocatable::wann_acof(:,:,:,:,:,:)
      59           0 :       complex,allocatable::wann_bcof(:,:,:,:,:,:)
      60           0 :       complex,allocatable::wann_ccof(:,:,:,:,:,:,:)
      61           0 :       real,allocatable::ff(:,:,:,:)
      62           0 :       real,allocatable::gg(:,:,:,:)
      63           0 :       real,allocatable::flo(:,:,:,:)
      64             :       complex xdnout
      65             :       LOGICAL       :: twodim
      66             :       LOGICAL :: cartesian,xsf
      67             :       REAL    :: pt(3),vec1(3),vec2(3),vec3(3),zero(3),v(3),point(3)
      68             :       INTEGER :: grid(3)
      69             :       CHARACTER(len=30):: filename
      70             :       NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
      71             :       CHARACTER(len=20):: name1,name2,name3
      72             :       CHARACTER(len=10):: vandername
      73             :       CHARACTER*8      :: name(10)
      74             :       integer ntype2,jmtd2,lmaxd2,nlod2,natd2,lmd2,llod2
      75             :       integer unigrid(4)
      76           0 :       real tmp_zatom(1:ntype),tmp_dx(1:ntype),tmp_rmt(1:ntype)
      77           0 :       integer tmp_jri(1:ntype),tmp_neq(1:ntype),tmp_rmsh(jmtd,ntype)
      78           0 :       complex,allocatable::wannint(:,:,:,:)
      79             :       real delta1,delta2,plot_time_int,plot_read_data_time
      80             :       real plot_time_sph,time_sqrt
      81           0 :       plot_time_sph=0.0
      82           0 :       plot_time_int=0.0
      83           0 :       plot_read_data_time=0.0
      84           0 :       time_sqrt=0.0
      85             : 
      86             : 
      87           0 :       INQUIRE(file ="plot_inp",exist= twodim)
      88           0 :       IF (.NOT.twodim) THEN !no input file exists, create a template and
      89             :                             !exit
      90           0 :          OPEN(20,file ="plot_inp")
      91           0 :          WRITE(20,'(i2,a5,l1)') 1,",xsf=",.false.
      92             : c         WRITE(20,*) "&PLOT twodim=t,cartesian=t"
      93             : c         WRITE(20,*) "  vec1(1)=10.0 vec2(2)=10.0"
      94             : c         WRITE(20,*) "  filename='plot1' /"
      95           0 :          WRITE(20,*) "&PLOT twodim=f,cartesian=f"
      96           0 :          WRITE(20,*) "  vec1(1)=1.0 vec1(2)=0.0 vec1(3)=0.0 "
      97           0 :          WRITE(20,*) "  vec2(1)=0.0 vec2(2)=1.0 vec2(3)=0.0 "
      98           0 :          WRITE(20,*) "  vec3(1)=0.0 vec3(2)=0.0 vec3(3)=1.0 "
      99           0 :          WRITE(20,*) "  grid(1)=30  grid(2)=30  grid(3)=30  "
     100             : c         WRITE(20,*) "  zero(1)=0.0 zero(2)=0.0 zero(3)=0.5 "
     101           0 :          WRITE(20,*) "  filename ='plot2' /"
     102           0 :          CLOSE(20)
     103           0 :          WRITE(*,*) "No plot_inp file found. Created a template"
     104             :          CALL juDFT_error("Missing input for plot; modify plot_inp"
     105           0 :      +        ,calledby ="wann_plot_from_lapw")
     106             :       ENDIF
     107             : 
     108           0 :       l_file=.false.
     109           0 :       do jspin=1,2  !spin loop
     110           0 :          call cpu_time(delta1)
     111             : c*******************************************************
     112             : c          read in data
     113             : c*******************************************************
     114           0 :        inquire(file=spin12(jspin)//'.lapw',exist=l_file)
     115           0 :        IF(.NOT.l_file) CALL juDFT_error("WF..lapw not found",calledby
     116           0 :      +      ="wann_plot_from_lapw")
     117           0 :        print*,"open file WF..lapw"
     118           0 :        open(344,file=spin12(jspin)//'.lapw',form='unformatted')
     119           0 :        read(344)num_wann,cell1,cell2,cell3
     120           0 :        read(344)amat_old(3,3)
     121           0 :        read(344)ntype2,jmtd2,lmaxd2,nlod2,natd2,lmd2,llod2
     122           0 :        IF(ntype2/=ntype) CALL juDFT_error("ntype2",calledby
     123           0 :      +      ="wann_plot_from_lapw")
     124           0 :        IF(jmtd2/=jmtd) CALL juDFT_error("jmtd2",calledby
     125           0 :      +      ="wann_plot_from_lapw")
     126           0 :        IF(lmaxd2/=lmaxd) CALL juDFT_error("lmaxd2",calledby
     127           0 :      +      ="wann_plot_from_lapw")
     128           0 :        IF(nlod2/=nlod) CALL juDFT_error("nlod2",calledby
     129           0 :      +      ="wann_plot_from_lapw")
     130           0 :        IF(natd2/=natd) CALL juDFT_error("natd2",calledby
     131           0 :      +      ="wann_plot_from_lapw")
     132           0 :        IF(lmd2/=lmd) CALL juDFT_error("lmd2",calledby
     133           0 :      +      ="wann_plot_from_lapw")
     134           0 :        IF(nlod2/=nlod) CALL juDFT_error("nlod2",calledby
     135           0 :      +      ="wann_plot_from_lapw")
     136           0 :        read(344)pos_old(1:3,1:natd)
     137           0 :        read(344)tmp_neq(1:ntype)
     138           0 :        read(344)tmp_zatom(1:ntype)
     139           0 :        read(344)tmp_dx(1:ntype)
     140           0 :        read(344)tmp_rmt(1:ntype)
     141           0 :        read(344)tmp_jri(1:ntype)
     142           0 :        read(344)tmp_rmsh(1:jmtd,1:ntype)
     143             : c*******************************************************
     144             : c         read in radial wavefunctions
     145             : c*******************************************************
     146           0 :        allocate(ff(1:ntype,1:jmtd,1:2,0:lmaxd))
     147           0 :        allocate(gg(1:ntype,1:jmtd,1:2,0:lmaxd))
     148           0 :        allocate(flo(1:ntype,1:jmtd,1:2,1:nlod))
     149           0 :        read(344)ff(1:ntype,1:jmtd,1:2,0:lmaxd)
     150           0 :        read(344)gg(1:ntype,1:jmtd,1:2,0:lmaxd)
     151           0 :        read(344)flo(1:ntype,1:jmtd,1:2,1:nlod)
     152             : c*******************************************************
     153             : c         read in a-,b-, and c-coefficients
     154             : c*******************************************************
     155             : 
     156             :        allocate(wann_acof(num_wann,0:lmd,natd,-cell1:cell1,-cell2:cell2,
     157           0 :      &   -cell3:cell3))
     158             :        allocate(wann_bcof(num_wann,0:lmd,natd,-cell1:cell1,-cell2:cell2,
     159           0 :      &   -cell3:cell3))
     160             :        allocate(wann_ccof(num_wann,-llod:llod,nlod,natd,-cell1:cell1,
     161           0 :      &   -cell2:cell2,-cell3:cell3))
     162             : 
     163             :        read(344)wann_acof(1:num_wann,0:lmd,1:natd,-cell1:cell1,
     164           0 :      &                    -cell2:cell2,-cell3:cell3)
     165             :        read(344)wann_bcof(1:num_wann,0:lmd,1:natd,-cell1:cell1,
     166           0 :      &                    -cell2:cell2,-cell3:cell3)
     167             :        read(344)wann_ccof(1:num_wann,-llod:llod,1:nlod,1:natd,
     168           0 :      &                -cell1:cell1,-cell2:cell2,-cell3:cell3)
     169             : c*********************************************************
     170             : c        read in planewave expansion
     171             : c*********************************************************
     172           0 :        read(344)unigrid(1:4)
     173           0 :        allocate(wannint(-unigrid(4):unigrid(4),-unigrid(4):unigrid(4),
     174           0 :      ,           -unigrid(4):unigrid(4),num_wann))
     175           0 :        read(344)wannint(:,:,:,:)
     176             : 
     177           0 :        print*,"num_wann=",num_wann
     178           0 :        print*,"wannier supercell: ",cell1,cell2,cell3
     179           0 :        close(344)
     180           0 :        call cpu_time(delta2)
     181           0 :        plot_read_data_time=delta2-delta1
     182             : 
     183             : c**********************************************************************
     184             : c       plot the wannierfunction
     185             : c**********************************************************************
     186             :       !<-- Open the plot_inp file for input
     187           0 :       OPEN (18,file='plot_inp')
     188           0 :       READ(18,'(i2,5x,l1)') nplot,xsf
     189             :       ! If xsf is specified we create an input file for xcrysden
     190           0 :       IF (nplot.ge.2) 
     191             :      &     CALL juDFT_error
     192             :      +     ("plots one by one, please, this is not charge density"
     193           0 :      +     ,calledby="wann_plot_from_lapw")
     194             :       !<-- Loop over all plots
     195           0 :       DO nplo=1,nplot
     196             :          ! the defaults
     197           0 :          twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
     198           0 :          vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
     199           0 :          zero = (/0.,0.,0./);filename="default"
     200           0 :          READ(18,plot)
     201           0 :          IF (twodim.AND.ANY(grid(1:2)<1)) 
     202             :      +        CALL juDFT_error("Illegal grid size in plot",calledby
     203           0 :      +        ="wann_plot_from_lapw")
     204           0 :          IF (.NOT.twodim.AND.ANY(grid<1)) 
     205             :      +        CALL juDFT_error("Illegal grid size in plot",calledby
     206           0 :      +        ="wann_plot_from_lapw")
     207           0 :          IF (twodim) grid(3) = 1
     208             :          !calculate cartesian coordinates if needed
     209           0 :          IF (.NOT.cartesian) THEN
     210           0 :             vec1=matmul(amat,vec1)
     211           0 :             vec2=matmul(amat,vec2)
     212           0 :             vec3=matmul(amat,vec3)
     213           0 :             zero=matmul(amat,zero)
     214             :          ENDIF
     215             :          !Open the file
     216           0 :          IF (filename =="default") WRITE(filename,'(a,i2)') "plot",nplo
     217             : c..loop by the wannierfunctions
     218           0 :          nbmin=1
     219           0 :          nbmax=num_wann
     220           0 :          bands:DO nbn = nbmin,nbmax
     221             : 
     222           0 :          IF (xsf) THEN
     223           0 :             write (name1,22) nbn,jspin
     224             :    22       format (i3.3,'.real.',i1,'.xsf')
     225           0 :             write (name2,23) nbn,jspin
     226             :    23       format (i3.3,'.imag.',i1,'.xsf')
     227           0 :             write (name3,24) nbn,jspin
     228             :    24       format (i3.3,'.absv.',i1,'.xsf')
     229           0 :             OPEN(55,file=name1)
     230             : !#if 1==1
     231           0 :             call judft_error("NOT INPLEMENTED")
     232             : c$$$#else
     233             : c$$$            CALL xsf_WRITE_atoms(
     234             : c$$$     >                        55,film,odi%d1,amat,neq(:ntype),
     235             : c$$$     >                        zatom(:ntype),pos)
     236             : c$$$            OPEN(56,file=name2)
     237             : c$$$            CALL xsf_WRITE_atoms(
     238             : c$$$     >                        56,film,odi%d1,amat,neq(:ntype),
     239             : c$$$     >                        zatom(:ntype),pos)
     240             : c$$$            OPEN(57,file=name3)
     241             : c$$$            CALL xsf_WRITE_atoms(
     242             : c$$$     >                        57,film,odi%d1,amat,neq(:ntype),
     243             : c$$$     >                        zatom(:ntype),pos)
     244             : c$$$            CALL xsf_WRITE_header(55,twodim,filename,(vec1),
     245             : c$$$     &       (vec2),(vec3),zero
     246             : c$$$     $           ,grid)
     247             : c$$$            CALL xsf_WRITE_header(56,twodim,filename,(vec1),
     248             : c$$$     &       (vec2),(vec3),zero
     249             : c$$$     $           ,grid)
     250             : c$$$            CALL xsf_WRITE_header(57,twodim,filename,(vec1),
     251             : c$$$     &       (vec2),(vec3),zero
     252             : c$$$     $           ,grid)
     253             : c$$$#endif
     254             :          ELSE
     255           0 :                WRITE (vandername,201) nbn,jspin
     256             :   201          FORMAT (i5.5,'.',i1)            
     257           0 :                OPEN(55,file=vandername)
     258           0 :                WRITE (55,7) grid(1),grid(2),grid(3)
     259             :     7          FORMAT (5i4)
     260             :          ENDIF
     261             : 
     262           0 :          DO iz = 0,grid(3)-1
     263           0 :           DO iy = 0,grid(2)-1
     264           0 :            xloop:DO ix = 0,grid(1)-1
     265             :             point = zero+vec1*REAL(ix)/(grid(1)-1)+vec2*REAL(iy)
     266           0 :      $                 /(grid(2)-1)
     267           0 :             IF (.NOT.twodim) point = point+vec3*REAL(iz)/(grid(3)-1)
     268             : !Check if the point is in MT-sphere
     269             : 
     270           0 :              ii1 = cell1
     271           0 :              ii2 = cell2
     272           0 :              ii3 = cell3
     273           0 :              IF (film .AND. .NOT.odi%d1) ii3 = 0
     274           0 :              IF (odi%d1) THEN
     275           0 :                 ii1 = 0 ; ii2 = 0
     276             :              END IF
     277           0 :              DO  i1 = -ii1,ii1
     278           0 :               DO  i2 = -ii2,ii2
     279           0 :                DO  i3 = -ii3,ii3
     280           0 :                 pt = point+MATMUL(amat,(/i1,i2,i3/))
     281           0 :                 na = 0
     282           0 :                 DO nt = 1,ntype
     283           0 :                  DO nq = 1,neq(nt)
     284           0 :                   na   = na + 1
     285           0 :                   s  = dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt)
     286           0 :                   IF (s<(rmsh(jri(nt),nt))**2) THEN
     287           0 :                    pt(:)=pt(:)-pos(:,na)
     288             :                    call wann_lapw_sph_plot(ff(nt,:,1,:),gg(nt,:,1,:),
     289             :      >                   flo(nt,:,1,:),wann_acof(nbn,:,na,-i1,-i2,-i3),
     290             :      >                   wann_bcof(nbn,:,na,-i1,-i2,-i3),
     291             :      >                   wann_ccof(nbn,:,:,na,-i1,-i2,-i3),pt,nlo(nt),
     292             :      >                   jmtd,lmaxd,nlod,llod,lmd,rmsh(:,nt),lmax(nt),
     293             :      >                   llo(:,nt),jri(nt),
     294           0 :      <                   xdnout)
     295           0 :                    IF (xsf) THEN
     296           0 :                       WRITE(55,*) real(xdnout)
     297           0 :                       WRITE(56,*) aimag(xdnout)
     298           0 :                       WRITE(57,*) real(xdnout*conjg(xdnout))
     299             :                    ELSE
     300           0 :                       WRITE(55,8) real(xdnout),aimag(xdnout)
     301             :                    ENDIF
     302             :                    CYCLE xloop
     303             :                   ENDIF
     304             :                  ENDDO
     305             :                 ENDDO !nt
     306             :                ENDDO
     307             :               ENDDO
     308             :              ENDDO !i1
     309             : !Check for point in vacuum
     310           0 :              IF (film.AND..NOT.odi%d1.AND.ABS(point(3))>=z1) THEN
     311           0 :                 ivac=1
     312             :                 if (point(3).lt. 0.0)ivac=2
     313           0 :                 jvac=ivac
     314             :                 if(nvac==1)jvac=1
     315           0 :                 xdnout=0.0
     316             : c     call wann_plot_vac
     317             :                 if(real(xdnout).gt.9.0 .or.real(xdnout).lt.-9.0
     318           0 :      &        .or.aimag(xdnout).gt.9.0 .or. aimag(xdnout).lt.-9.0)then
     319             :                 xdnout=cmplx(0.0,0.0)
     320             :                 print*,"vac-problem at z=",point(3)
     321             :                 endif
     322             : 
     323           0 :               IF (xsf) THEN
     324           0 :                  WRITE(55,*) real(xdnout)
     325           0 :                  WRITE(56,*) aimag(xdnout)
     326           0 :                  WRITE(57,*) real(xdnout*conjg(xdnout))
     327             :               ELSE
     328           0 :                  WRITE(55,8) real(xdnout),aimag(xdnout)
     329             :               ENDIF
     330             :               CYCLE xloop
     331             :              END IF
     332             :             
     333           0 :              IF (odi%d1) THEN
     334           0 :               IF (SQRT((pt(1))**2 + (pt(2))**2)>=z1) THEN
     335             : c     call wann_real
     336           0 :                IF (xsf) THEN
     337           0 :                   WRITE(55,*) real(xdnout)
     338           0 :                   WRITE(56,*) aimag(xdnout)
     339           0 :                   WRITE(57,*) real(xdnout*conjg(xdnout))
     340             :                ELSE
     341           0 :                   WRITE (55,8) real(xdnout),aimag(xdnout)
     342             :                ENDIF
     343             :                CYCLE xloop
     344             :               END IF
     345             :              END IF
     346             :                call wann_lapw_int_plot(point,bmat,unigrid,
     347             :      >            wannint(:,:,:,nbn),
     348           0 :      <            xdnout)
     349             : c               xdnout=cmplx(0.0,0.0)
     350           0 :              IF (xsf) THEN
     351           0 :                 WRITE(55,*) real(xdnout)
     352           0 :                 WRITE(56,*) aimag(xdnout)
     353           0 :                 WRITE(57,*) real(xdnout*conjg(xdnout))
     354             :              ELSE
     355           0 :                 WRITE(55,8) real(xdnout),aimag(xdnout)
     356             :              ENDIF
     357             :             ENDDO xloop
     358             :            ENDDO
     359             :           ENDDO !z-loop
     360             : 
     361           0 :           IF (xsf) THEN
     362           0 :               CALL xsf_WRITE_endblock(55,twodim)
     363           0 :               CALL xsf_WRITE_endblock(56,twodim)
     364           0 :               CALL xsf_WRITE_endblock(57,twodim)
     365           0 :               CLOSE (55) ; CLOSE (56) ; CLOSE (57)
     366             :           ENDIF
     367             : c..end of the loop by the bands
     368             :           ENDDO bands   
     369             :       ENDDO   !nplot      
     370           0 :       CLOSE(18)
     371           0 :       IF (.not.xsf) CLOSE(55)
     372             :     8 FORMAT (2f16.12)
     373             : 
     374             : 
     375           0 :        deallocate(wann_acof,wann_bcof,wann_ccof)
     376             : 
     377           0 :        if(jspins.eq.1)exit
     378             :       enddo !spinloop
     379           0 :       print*,"plot_time_int=",plot_time_int
     380           0 :       print*,"plot_time_sph=",plot_time_sph
     381           0 :       print*,"plot_read_data_time=",plot_read_data_time
     382           0 :       print*,"time_sqrt=",time_sqrt
     383             : 
     384           0 :       END SUBROUTINE wann_plot_from_lapw
     385             :       END MODULE m_wann_plot_from_lapw

Generated by: LCOV version 1.13