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