PROGRAM diel_vis
! diel_vis reads in the dielectric matrix, stored in
! G-space and transforms it in to real space. Also with frequency
! dependence we can look at plasmon modes.
  USE environment,ONLY : environment_start
  USE io_global,  ONLY : stdout, ionode, ionode_id
  USE mp_global,  ONLY : nproc_pool
  USE mp,         ONLY : mp_bcast
  USE parameters, ONLY : ntypx
  USE constants,  ONLY :  pi, fpi, tpi
  USE cell_base
  USE ions_base,         ONLY : nat, ityp, atm, ntyp => nsp, tau, zv
  USE lsda_mod,          ONLY: nspin
  USE gvect
  USE gsmooth
  USE wavefunctions_module,  ONLY: psic
  USE io_files,              ONLY: nd_nmbr, prefix, outdir
  USE fft_base,              ONLY : grid_scatter
  USE printout_base,         ONLY : title
  USE control_flags,         ONLY : gamma_only
  USE scf,                   ONLY : rho, vltot, v

IMPLICIT NONE

integer, parameter :: DIRECT_IO_FACTOR = 8
real(DP), parameter :: ryd2ev = 13.6058
complex(DP), allocatable :: a(:)
REAL(DP), allocatable    ::  wim(:), rhor(:)
REAL(DP)                 :: ar, ai
INTEGER :: nwim, iq, ios, iw
INTEGER :: iuncoul, iunsex, iunsigma
INTEGER :: ig, igp, iwim, recl
INTEGER*8 :: unf_recl
REAL (DP) :: eta, w_ryd, RYTOEV
integer :: glist(3), gplist(3), ngcount, ngpcount
integer :: ounit, style, output_format
complex(DP), allocatable :: invepsz(:)
complex(DP) :: arg
integer :: nx, ny, ngmpol, ipol, idum, na
integer :: ngmsex, nwsigma
complex(DP), allocatable :: rhog (:), fiu(:)
integer  :: plot_num
CHARACTER(len=256) :: filplot
COMPLEX(DP), ALLOCATABLE :: scrcoul_g(:,:,:), sigma(:,:,:), sigma_ex(:,:)
COMPLEX(DP), ALLOCATABLE :: fxyz(:,:,:)
INTEGER                  :: x, y, z
LOGICAL                  :: sigorcoul
LOGICAL                  :: plot_re
LOGICAL                  :: plot_rho

!e1,e2, e3 ar vector that
  real(DP) :: e1(3), e2(3), e3(3), x0 (3), radius, m1, m2, m3

!ounit just standard out.
  ounit = 6
  CALL environment_start ( 'POST-PROC' )

  CALL input_from_file ( )

  prefix = 'MoS2'
  outdir = './tmp'

  !    This subroutine reads the data for the output file produced by pw.x
  !    extracts and calculates the desired quantity (rho, V, ...)
  !    writes it to a file for further processing or plotting
  call extract (filplot, plot_num) 


!!!!!!!!! SHOULD FORMULATE INPUT FILE !!!!!!!!!!!!!!!!!
!ngmpol: number of G-vectors used to describe dielectric matrix.
!nwim: number of imaginary frequencies
!for Si:
  ngmpol    = 137
  ngmsex    = 411
  nwsigma   = 42

!sigorcoul
  sigorcoul =.true.
  plot_rho  =.false.

!output format
  output_format = 0

!vectors for 
 !!! x0, e1, e2 are in alat units !!!
     e1(:) = (/1.0, 1.0, 0.0/)
     e2(:) = (/0.0, 0.0, 1.0/)
   !!!orthorhombic crystals the bond mid point should
   !!!just be mid point of lattice vectors...
     e3(:) = (/0.0, 0.0, 0.0/)
   !!!origin of plane!!!
     x0(:) = (/0.0, 0.0, 0.0/)

!number of points in the plane
!     rho(i,j) = rho( x0 + e1 * (i-1)/(nx-1)
!              + e2 * (j-1)/(ny-1) ), i=1,nx ; j=1,ny

     nx    = 56
     ny    = 40

     if (e1(1)**2 + e1(2)**2 + e1(3)**2 <  1d-6 .or. &
         e2(1)**2 + e2(2)**2 + e2(3)**2 <  1d-6)     &
         call errore ('chdens', 'missing e1/e2 vectors', 1)

     if (abs(e1(1)*e2(1) + e1(2)*e2(2) + e1(3)*e2(3)) > 1d-6) &
         call errore ('chdens', 'e1 and e2 are not orthogonal', 1)

     if (nx <= 0 .or. ny <= 0 )   call errore ('chdens', 'wrong nx/ny', 2)

  !
  !    At this point we start the calculations, first we normalize the 
  !    vectors defining the plotting region. 
  !    If these vectors have 0 length, replace them with crystal axis
  !

  m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
  if (abs(m1) < 1.d-6) then
     e1 (:) = at(:,1)
     m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
  end if
  e1 (:) = e1 (:) / m1
  !
  m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
  if (abs(m2) < 1.d-6) then
     e2 (:) = at(:,2)
     m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
  end if
  e2 (:) = e2 (:) / m2
  !
  m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
  if (abs(m3) < 1.d-6) then
     e3 (:) = at(:,3)
     m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
  end if
  e3 (:) = e3 (:) / m3

!READ in dielectric matrix...
     open (45, file = "./imfreq.dat", form = 'formatted', status = 'unknown')
           read (45, *)
           read (45, *) nwim

     allocate(fiu(nwim))

     if (sigorcoul) then 
         ALLOCATE (sigma_ex (ngmsex, ngmsex))
         ALLOCATE (sigma(ngmpol,ngmpol,nwsigma))
     else
         ALLOCATE ( scrcoul_g (ngmpol, ngmpol, nwim) )
     endif

     do iw = 1, nwim
        read (45, *) ar, ai
        fiu(iw) = dcmplx(ar, ai)
        write(6,*) fiu(iw)
     enddo
     close(45)

     if(sigorcoul) then
        iunsex = 33
        recl = 2 * ngmsex * ngmsex
        unf_recl = 8 * int(recl, kind=kind(unf_recl))
        open ( iunsex, file = "./_gw0si.sigma_ex1", iostat = ios, form = 'unformatted', &
               status = 'unknown', access = 'direct', recl = unf_recl)

        iunsigma = 32
        recl = 2 * ngmpol * ngmpol * nwsigma
        unf_recl = 8 * int(recl, kind=kind(unf_recl))
        open ( iunsigma, file = "./_gw0si.sigma1", iostat = ios, form = 'unformatted', &
               status = 'unknown', access = 'direct', recl = unf_recl)

     else
       iuncoul = 28
       recl = 2 * ngmpol * ngmpol * nwim
       unf_recl = 8 * int(recl, kind=kind(unf_recl))
       open ( iuncoul, file = "./_gw0si.coul1", iostat = ios, form = 'unformatted', &
              status = 'unknown', access = 'direct', recl = unf_recl)
     endif

     write(6,*) ios

     iq = 1
     iw = 1

     if(sigorcoul) then
        READ ( iunsigma, rec = iq, iostat = ios) sigma
        READ ( iunsex,   rec = iq, iostat = ios) sigma_ex
     else
        READ ( iuncoul,  rec = iq, iostat = ios) scrcoul_g
     endif

     write(6,*) ios


     if(.not.sigorcoul) then
    !add one back since we store: \epsilon^{-1} - 1.
         do ig = 1, ngmpol
            scrcoul_g(ig,ig,iw) = scrcoul_g(ig,ig,iw) + 1.0d0
         enddo
     endif

ALLOCATE  (rhog( ngm))
ALLOCATE  (rhor(nrx1*nrx2*nrx3))

rhog(:) = (0.0d0,0.0d0)

style=4
if(style.eq.1) then
! response in cubefile format to a point test charge at a particular point in real space.
       do igp   = 1, ngmpol
          do ig = 1, ngmpol
            !arg = tpi*dcmplx(0.0d0, 1.0d0)*(e3(1)*(g(1,ig) + e(2)*g(2,ig) + e(3)*g(3,ig))
             arg = -tpi*(e3(1)*g(1,ig) + e3(2)*g(2,ig) + e3(3)*g(3,ig))
             !rhog(igp) = rhog(igp) + scrcoul_g(ig,igp,iw)*exp(dcmplx(0.0d0, 1.0d0)*arg)
             !to plot sigma_{\r}(\r').
             rhog(igp) = rhog(igp) + sigma(ig,igp,6)*exp(dcmplx(0.0d0, 1.0d0)*arg)
          enddo
       enddo

       rhog(:) = (1.0d0/sqrt(omega))*rhog(:)

       ALLOCATE  (invepsz(nrx1*nrx2*nrx3))

       invepsz(nl(:)) = rhog(:)

       call cft3 (invepsz, nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)

       rhor(:) = (1.0d0/sqrt(omega))*abs(real(invepsz(:)))

       call write_cubefile (alat, at, bg, nat, tau, atm, ityp, rhor, &
             nr1, nr2, nr3, nrx1, nrx2, nrx3, ounit)

else if (style.eq.2) then

!could also just FFT this whole thing in r but should try to do things
!with at least a modicum of class.
!GNUPLOT FORMAT OUTPUT FOR function_{\r}(\r').

     do igp   = 1, ngmpol
        do ig = 1, ngmpol
          !arg = tpi*dcmplx(0.0d0, 1.0d0)*(e3(1)*(g(1,ig) + e(2)*g(2,ig) + e(3)*g(3,ig))
           arg = -tpi*(e3(1)*g(1,ig) + e3(2)*g(2,ig) + e3(3)*g(3,ig))
          !rhog(igp) = rhog(igp) + scrcoul_g(ig,igp,iw)*exp(dcmplx(0.0d0, 1.0d0)*arg)
          !to plot sigma_{\r}(\r').
           rhog(igp) = rhog(igp) + sigma(ig,igp, 11)*exp(dcmplx(0.0d0, 1.0d0)*arg)
        enddo
     enddo

     rhog(:) = (1.0d0/sqrt(omega))*rhog(:)
     ALLOCATE  (invepsz(nrx1*nrx2*nrx3))
     invepsz(nl(:)) = rhog(:)
     call cft3 (invepsz, nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)
     rhor(:) = (1.0d0/sqrt(omega))*(real(invepsz(:)))

     allocate(fxyz(nr1,nr2,nr3))
  
     DO x = 1, nr1
          DO y = 1, nr2
             DO z = 1, nr3
                fxyz(x,y,z) = fxyz(x,y,z) + rhor(x + (y-1)*nr2 + (z-1)*nr3*nr2)
             ENDDO
         ENDDO
     ENDDO

!writing potential in xyz format should be middle of cell:
     write(6,*) int(nr1/2)
     !z = int(nr3/4)
     z = 10

     do x = int(nr1/2)+1, nr1
        do y = int(nr2/2)+1, nr2
              write (450, '(2i4, 1f12.7)') x-nr1, y-nr2, real(fxyz(x,y,z))
        enddo
        do y = 1, int(nr2/2)
              write (450, '(2i4, 1f12.7)') x-nr1, y, real(fxyz(x,y,z))
        enddo
           write (450,*)
     enddo

     do x = 1, int(nr1/2)
        do y = int(nr2/2)+1, nr2
              write (450, '(2i4, 1f12.7)') x, y-nr2, real(fxyz(x,y,z))
        enddo
        do y = 1, int(nr2/2)
              write (450, '(2i4, 1f12.7)') x, y, real(fxyz(x,y,z))
        enddo
            write (450,*)
     enddo

else if (style.eq.3) then
     !\inveps_{\G}(G').
     !ig = 9
     !  do igp = 1, ngmpol
     !      rhog(igp) = rhog(igp) + scrcoul_g(ig,igp,iw)
     !  enddo
       ALLOCATE  (invepsz(nrx1*nrx2*nrx3))
       invepsz(nl(:)) = rhog(:)
       call cft3 (invepsz, nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)
       rhor = (1.0d0/sqrt(omega))*real(invepsz(:))
       call write_cubefile (alat, at, bg, nat, tau, atm, ityp, rhor, &
             nr1, nr2, nr3, nrx1, nrx2, nrx3, ounit)

else if (style.eq.4) then
!planar fxns.
! If one wanted to plot charge density they would do this.
  if(plot_rho) then
      ALLOCATE  (invepsz(nrx1*nrx2*nrx3))
       invepsz = rho%of_r(:,1)
       call cft3 (invepsz, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
       rhog (:) = invepsz(nl(:))
  else
! If one wanted to plot Sigma they would do...
     !  do igp   = 1, ngmpol
     !     do ig = 1, ngmpol
     !        arg = -tpi*(x0(1)*g(1,ig) + x0(2)*g(2,ig) + x0(3)*g(3,ig))
     !        rhog(igp) = rhog(igp) + sigma(ig,igp,11)*exp(dcmplx(0.0d0, 1.0d0)*arg)
     !     enddo
     !  enddo
     !  rhog(:) = (1.0d0/sqrt(omega))*rhog(:)
!or sigma_x
       do igp   = 1, ngmsex
          do ig = 1, ngmsex
             arg = -tpi*(x0(1)*g(1,ig) + x0(2)*g(2,ig) + x0(3)*g(3,ig))
             rhog(igp) = rhog(igp) + sigma_ex(ig,igp)*exp(dcmplx(0.0d0, 1.0d0)*arg)
          enddo
       enddo
       rhog(:) = (1.0d0/sqrt(omega))*rhog(:)
  endif!plotrho

       call plot_2d (nx, ny, m1, m2, x0, e1, e2, ngm, g, rhog, alat, &
                     at, nat, tau, atm, ityp, output_format, ounit)
endif!style

        write (ounit, '(i4)') nat
        write (ounit, '(3f8.4,i3)') ( (tau(ipol,na), ipol=1,3), 1, na=1,nat)
        write (ounit, '(f10.6)') celldm (1)
        write (ounit, '(3(3f12.6/))') at
END PROGRAM diel_vis


subroutine plot_2d (nx, ny, m1, m2, x0, e1, e2, ngm, g, rhog, alat, &
     at, nat, tau, atm, ityp, output_format, ounit)
  !-----------------------------------------------------------------------
  !
  USE kinds, only : DP
  use constants, only : pi
  use io_global, only : stdout, ionode
  USE mp_global,  ONLY : intra_pool_comm
  USE mp,         ONLY : mp_sum
  implicit none
  integer :: nx, ny, ngm, nat, ityp (nat), output_format, ounit
  ! number of points along x
  ! number of points along y
  ! number of G vectors
  ! number of atoms
  ! types of atoms
  ! output unit
  ! output format
  character(len=3) :: atm(*) ! atomic symbols
  real(DP) :: e1(3), e2(3), x0(3), m1, m2, g(3,ngm), alat, &
       tau(3,nat), at(3,3)
  ! vectors e1, e2 defining the plane
  ! origin
  ! modulus of e1
  ! modulus of e2
  ! G-vectors

  complex(DP) :: rhog (ngm)
  ! rho or polarization in G space
  integer :: i, j, ig

  real(DP) :: rhomin, rhomax, rhoim, deltax, deltay
  ! minimum value of the charge
  ! maximum value of the charge
  ! integrated imaginary charge
  ! steps along e1
  ! steps along e2
  complex(DP), allocatable :: eigx (:), eigy (:), carica(:,:)

  allocate (eigx(  nx))    
  allocate (eigy(  ny))    
  allocate (carica( nx , ny))    

  deltax = m1 / (nx - 1)
  deltay = m2 / (ny - 1)

  carica(:,:) = (0.d0,0.d0)
  do ig = 1, ngm
     !
     ! eigx=exp(iG*e1+iGx0), eigy=(iG*e2)
     ! These factors are calculated and stored in order to save CPU time
     !
     do i = 1, nx
        eigx (i) = exp ( (0.d0, 1.d0) * 2.d0 * pi * ( (i - 1) * deltax * &
             (e1(1) * g(1,ig) + e1(2) * g(2,ig) + e1(3) * g(3,ig) ) + &
             (x0 (1) * g(1,ig) + x0 (2) * g(2,ig) + x0 (3) * g(3,ig) ) ) )
     enddo
     do j = 1, ny
        eigy (j) = exp ( (0.d0, 1.d0) * 2.d0 * pi * (j - 1) * deltay * &
             (e2(1) * g(1,ig) + e2(2) * g(2,ig) + e2(3) * g(3,ig) ) )
     enddo
     do j = 1, ny
        do i = 1, nx
           carica (i, j) = carica (i, j) + rhog (ig) * eigx (i) * eigy (j)
        enddo
     enddo
  enddo
  !HL
  !call mp_sum( carica, intra_pool_comm ) 
  !
  !    Here we check the value of the resulting charge
  !
  rhomin =  1.0d10
  rhomax = -1.0d10

  rhoim = 0.d0
  do i = 1, nx
     do j = 1, ny
        rhomin = min (rhomin,  DBLE (carica (i, j) ) )
        rhomax = max (rhomax,  DBLE (carica (i, j) ) )
        rhoim = rhoim + abs (AIMAG (carica (i, j) ) )
     enddo

  enddo

  rhoim = rhoim / nx / ny
  write(stdout, '(5x,"Min, Max, imaginary charge: ",3f12.6)') &
                 rhomin, rhomax, rhoim

  !
  !     and we print the charge on output
  !
  if (ionode) then
     if (output_format == 0) then
        !
        !     gnuplot format
        !
        !         write(ounit,'(2i6)') nx,ny
        !HL
        !original
        !do i = 1, nx
        !   write (ounit, '(e25.14)') (  DBLE(carica(i,j)), j = 1, ny )
        !   write (ounit, * )
        !enddo
        !with folding along y
        do i = 1, nx
           write (ounit, '(e25.14)') (  DBLE(carica(i,j)), j = int(ny/2)+1, ny )
           write (ounit, '(e25.14)') (  DBLE(carica(i,j)), j = 1, (ny/2) )
           write (ounit, * )
        enddo

     elseif (output_format == 1) then
        !
        !     contour.x format
        !
        write (ounit, '(3i5,2e25.14)') nx, ny, 1, deltax, deltay
        write (ounit, '(4e25.14)') ( (  DBLE(carica(i,j)), j = 1, ny ), i = 1, nx )
     elseif (output_format == 2) then
        !
        !     plotrho format
        !
        write (ounit, '(2i4)') nx - 1, ny - 1
        write (ounit, '(8f8.4)') (deltax * (i - 1) , i = 1, nx)
        write (ounit, '(8f8.4)') (deltay * (j - 1) , j = 1, ny)
        write (ounit, '(6e12.4)') ( (  DBLE(carica(i,j)), i = 1, nx ), j = 1, ny )
        write (ounit, '(3f8.4)') x0
        write (ounit, '(3f8.4)') (m1 * e1 (i) , i = 1, 3)
        write (ounit, '(3f8.4)') (m2 * e2 (i) , i = 1, 3)

    ! elseif (output_format == 3) then
        !
        ! XCRYSDEN's XSF format
        !
    !    call xsf_struct (alat, at, nat, tau, atm, ityp, ounit)
    !    call xsf_datagrid_2d (carica, nx, ny, m1, m2, x0, e1, e2, alat, ounit)
     else
        call errore('plot_2d', 'wrong output_format', 1)
     endif
  endif

  deallocate (carica)
  deallocate (eigy)
  deallocate (eigx)
  return
end subroutine plot_2d