pes_out.F90

Go to the documentation of this file.

!! Copyright (C) 2016 U. De Giovannini
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module pes_out_oct_m
  use coordinate_system_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use io_function_oct_m
  use loct_oct_m
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
#if defined(HAVE_NETCDF)
  use netcdf
#endif
  use profiling_oct_m
  use qshep_oct_m
  use space_oct_m
  use sort_oct_m
  use string_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use vtk_oct_m
 
  implicit none
 
  private
 
  public ::                             &
    pes_out_velocity_map,               &
    pes_out_velocity_map_cut,           &
    pes_out_arpes_cut
 
  integer, parameter ::   &
    INTEGRATE_NONE    = -1,       &
    integrate_phi     =  1,       &
    integrate_theta   =  2,       &
    integrate_r       =  3,       &
    integrate_kx      =  4,       &
    integrate_ky      =  5,       &
    integrate_kz      =  6
 
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine pes_out_velocity_map(pesK, file, namespace, space, Lk, ll, how, spacing, coord_system, pmesh)
    real(real64),      intent(in) :: pesK(:,:,:)
    character(len=*),  intent(in) :: file
    type(namespace_t), intent(in) :: namespace
    class(space_t),    intent(in) :: space
    real(real64),      intent(in) :: Lk(:,:)
    integer,           intent(in) :: ll(:)
    integer(int64),       intent(in) :: how
    real(real64),      intent(in) :: spacing(3)
    class(coordinate_system_t), intent(in) :: coord_system
    real(real64), optional,   intent(in) :: pmesh(:,:,:,:)
 
    integer :: ierr
    character(len=512) :: filename
    type(cube_t) :: cube
    type(cube_function_t) :: cf
    real(real64) :: dk(3)
 
    push_sub(pes_out_velocity_map)
 
    call cube_init(cube, ll, namespace, space, spacing, coord_system)
    call dcube_function_alloc_rs(cube, cf, force_alloc = .true.)
    cf%dRS = pesk
 
 
    if (.not. present(pmesh)) then
      ! Ignore Lk and use pmesh
      dk(:) = m_zero
      dk(1:space%dim) = abs(lk(2, 1:space%dim) - lk(1,1:space%dim))
    end if
 
#if defined(HAVE_NETCDF)
 
    if (bitand(how, option__outputformat__netcdf) /= 0) then
      filename = trim(file)//".ncdf"
      write(message(1), '(a)') 'Writing netcdf format file: '
      call messages_info(1, namespace=namespace)
 
      call dout_cf_netcdf(filename, ierr, cf, cube, space, dk(:), .false., unit_one/units_out%length, namespace)
 
    end if
 
#endif
 
    if (bitand(how, option__outputformat__vtk) /= 0) then
      filename = trim(file)//".vtk"
      write(message(1), '(a)') 'Writing vtk format file: '
      call messages_info(1, namespace=namespace)
 
      if (present(pmesh)) then
        call dvtk_out_cf_structured(filename, namespace, 'PES_vel_map', ierr, cf, cube,&
          unit_one/units_out%length, pmesh, ascii = .false.)
      else
        call dvtk_out_cf(filename, namespace, 'PES_vel_map', ierr, cf, cube, dk(:),&
          unit_one/units_out%length)
      end if
    end if
 
    call cube_end(cube)
    call dcube_function_free_rs(cube, cf)
 
    pop_sub(pes_out_velocity_map)
 
  end subroutine pes_out_velocity_map
 
 
  ! ---------------------------------------------------------
  subroutine pes_out_arpes_cut(namespace, arpes, file, dim, ll, pmesh, Ekin)
    type(namespace_t), intent(in) :: namespace
    real(real64),      intent(in) :: arpes(:,:,:)
    character(len=*),  intent(in) :: file
    integer,           intent(in) :: dim
    integer,           intent(in) :: ll(:)
    real(real64),      intent(in) :: pmesh(:,:,:,:)
    real(real64),      intent(in) :: Ekin(:,:,:)
 
    integer :: iunit, ip, ie
    real(real64)   :: dp, length, pp(1:2), pdiff(1:2)
 
    push_sub(pes_out_arpes_cut)
 
    iunit = io_open(file, namespace, action='write')
    write(iunit, '(a)') '##################################################'
    write(iunit, '(a1,a18,2x,a18,2x,a18,2x,a18,2x, a18,2x,a18)') '#', &
      str_center("Ppath", 18), &
      str_center("px", 18), str_center("py", 18), &
      str_center("E", 18),  str_center("P[Ppath,E]", 18)
 
    write(iunit, '(a1,a18,2x,a18,2x,a18,2x,a18,2x, a18,2x,a18)') '#', &
      str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
      str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
      str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
      str_center('['//trim(units_abbrev(units_out%energy)) // ']', 18), &
      str_center('[1/' //trim(units_abbrev(units_out%energy))//']', 18)
    write(iunit, '(a)') '##################################################'
 
 
    length = m_zero
    pdiff(:) = m_zero
 
    do ip = 1, ll(1)
      pp(1:2) = pmesh(ip,1,1,1:2)
      if (ip > 1) pdiff = pp(1:2) - pmesh(ip-1,1,1,1:2)
      dp = norm2(pdiff(1:2))
      length = length + dp
 
      select case (dim)
      case (2)
        do ie = 1, ll(2)
          write(iunit, '(es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12)')   &
            units_from_atomic(unit_one/units_out%length, length),&
            units_from_atomic(unit_one/units_out%length, pp(1)), &
            units_from_atomic(unit_one/units_out%length, pp(2)), &
            units_from_atomic(units_out%energy, ekin(ip,ie,1)), &
            arpes(ip,ie,1)
        end do
 
      case (3)
        do ie = 1, ll(3)
          write(iunit, '(es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12)')   &
            units_from_atomic(unit_one/units_out%length, length),&
            units_from_atomic(unit_one/units_out%length, pp(1)), &
            units_from_atomic(unit_one/units_out%length, pp(2)), &
            units_from_atomic(units_out%energy, ekin(ip,1,ie)), &
            arpes(ip,1,ie)
        end do
 
      end select
 
 
      write(iunit, *)
 
 
    end do
 
 
    call io_close(iunit)
 
    pop_sub(pes_out_arpes_cut)
  end subroutine pes_out_arpes_cut
 
 
  ! ---------------------------------------------------------
  subroutine pes_out_velocity_map_cut(namespace, pesK, file, ll, dim, pol, dir, integrate, pos, Lk, pmesh)
    type(namespace_t), intent(in) :: namespace
    real(real64),      intent(in) :: pesk(:,:,:)
    character(len=*),  intent(in) :: file
    integer,           intent(in) :: ll(:)
    integer,           intent(in) :: dim
    real(real64),      intent(in) :: pol(3)
    integer,           intent(in) :: dir
    integer,           intent(in) :: integrate
    integer, optional, intent(in) :: pos(3)
    real(real64), optional,   intent(in) :: lk(:,:)
    real(real64), optional,   intent(in) :: pmesh(:,:,:,:)
 
    integer              :: ii, ix, iy, iunit,ldir(2), icut(3)
    real(real64)         :: kk(3),temp
    integer, allocatable :: idx(:,:)
    real(real64), allocatable   :: lk_(:,:)
    real(real64)         :: rotation(1:dim,1:dim)
    logical              :: aligned_axis
    ! integration
    real(real64)         :: k, kkk(3), theta, phi, dphi, dk(3)
    integer              :: iph, nphi
    ! progress
    integer              :: idone, ntodo
 
    real(real64), allocatable :: cube_f(:)
    type(qshep_t) :: interp
 
    push_sub(pes_out_velocity_map_cut)
 
    iunit = io_open(file, namespace, action='write')
 
 
 
    assert(size(pesk, 1) == ll(1))
 
    if (.not. present(pmesh)) then
      assert(present(lk))
 
      safe_allocate(idx(1:maxval(ll(:)), 1:3))
      safe_allocate(lk_(1:maxval(ll(:)), 1:3))
 
      dk(:) = m_zero
      dk(1:dim) = abs(lk(2,1:dim)-lk(1,1:dim))
 
      do ii = 1, 3
        lk_(:,ii) = lk(:,ii)
        call sort(lk_(1:ll(ii), ii), idx(1:ll(ii), ii)) !We need to sort the k-vectors in order to dump in gnuplot format
      end do
    end if
 
    aligned_axis = sum((pol-(/0 ,0 ,1/))**2)  <= m_epsilon  .or. &
      sum((pol-(/0 ,1 ,0/))**2)  <= m_epsilon  .or. &
      sum((pol-(/1 ,0 ,0/))**2)  <= m_epsilon
 
    if (present(pos)) then
      icut(1:3) = pos(1:3)
    else
      icut(1:3) = ll(1:3)/2 + 1
    end if
 
    if (aligned_axis .and. integrate == integrate_none) then !no need to rotate and interpolate
 
      select case (dir)
      case (1)
        ldir(:) =(/2,3/)
      case (2)
        ldir(:) =(/1,3/)
      case (3)
        ldir(:) =(/1,2/)
 
      end select
 
 
      if (present(pmesh)) then
        do ix = 1, ll(ldir(1))
          do iy = 1, ll(ldir(2))
 
            select case (dir)
            case (1)
              temp = pesk(icut(dir), ix, iy)
              kk(1) = pmesh(icut(dir), ix, iy, ldir(1))
              kk(2) = pmesh(icut(dir), ix, iy, ldir(2))
            case (2)
              temp = pesk(ix, icut(dir), iy)
              kk(1) = pmesh(ix, icut(dir), iy, ldir(1))
              kk(2) = pmesh(ix, icut(dir), iy, ldir(2))
            case (3)
              temp = pesk(ix, iy, icut(dir))
              kk(1) = pmesh(ix, iy, icut(dir), ldir(1))
              kk(2) = pmesh(ix, iy, icut(dir), ldir(2))
            end select
 
            write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
              units_from_atomic(sqrt(units_out%energy), kk(1)),&
              units_from_atomic(sqrt(units_out%energy), kk(2)),&
              temp
          end do
          write(iunit, *)
        end do
 
      else
        do ix = 1, ll(ldir(1))
          kk(1) = lk_(ix, ldir(1))
          do iy = 1, ll(ldir(2))
            kk(2) = lk_(iy, ldir(2))
 
            select case (dir)
            case (1)
              temp = pesk(idx(icut(dir), 1), idx(ix, 2), idx(iy, 3))
            case (2)
              temp = pesk(idx(ix, 1), idx(icut(dir), 2), idx(iy, 3))
            case (3)
              temp = pesk(idx(ix, 1), idx(iy, 2), idx(icut(dir), 3))
            end select
 
            write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
              units_from_atomic(sqrt(units_out%energy), kk(1)),&
              units_from_atomic(sqrt(units_out%energy), kk(2)),&
              temp
          end do
          write(iunit, *)
        end do
      end if
 
    else
      ! We set the z-axis along the pol vector
      call generate_rotation_matrix(rotation, (/ m_zero, m_zero, m_one /), pol)
 
      if (debug%info) then
        print *,"Rotate z-axis over the zenith axis"
        print *,rotation(1,:)
        print *,rotation(2,:)
        print *,rotation(3,:)
      end if
 
      call pes_out_interpolator_init(namespace, pesk, lk, ll, dim, cube_f, interp, pmesh)
 
      ntodo = product(ll(1:2))
      idone = 0
      call loct_progress_bar(-1, ntodo)
 
      do ix = 1, ll(1)
        do iy = 1, ll(2)
 
          !cut
          select case (dir)
          case (1)
            kk(1) = m_zero
            if (present(pmesh)) then
              kk(2) = pmesh(ix, 1, 1, 1)
              kk(3) = pmesh(1, iy, 1, 2)
            else
              kk(2) = lk_(ix, 1)
              kk(3) = lk_(iy, 2)
            end if
          case (2)
            kk(2) = m_zero
            if (present(pmesh)) then
              kk(1) = pmesh(ix,  1, 1, 1)
              kk(3) = pmesh(1,  iy, 1, 2)
            else
              kk(1) = lk_(ix, 1)
              kk(3) = lk_(iy, 2)
            end if
 
          case (3)
            kk(3) = m_zero
            if (present(pmesh)) then
              kk(1) = pmesh(ix,  1, 1, 1)
              kk(2) = pmesh(1,  iy, 1, 2)
            else
              kk(1) = lk_(ix, 1)
              kk(2) = lk_(iy, 2)
            end if
 
          end select
 
          temp = qshep_interpolate(interp, cube_f, matmul(rotation,kk(1:3)))
 
          select case (integrate)
          case (integrate_phi)
            temp = m_zero
            k = norm2(kk(1:3))
 
            nphi = 360
            dphi = m_two * m_pi / nphi
 
            do iph = 0, nphi
              phi = iph * dphi
              theta = atan2(norm2(kk(1:2)), kk(3))
 
              kkk(1) = k *sin(theta) * cos(phi)
              kkk(2) = k *sin(theta) * sin(phi)
              kkk(3) = k *cos(theta)
 
              temp = temp + &
                abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
            end do
            temp = temp * dphi
 
          case (integrate_kx)
            temp = m_zero
            do ii =1, ll(1)
              kkk(:) = kk(:) + (/lk(ii,1), m_zero, m_zero/)
              temp = temp + &
                abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
            end do
            temp = temp * dk(1)
 
          case (integrate_ky)
            temp = m_zero
            do ii =1, ll(2)
              kkk(:) = kk(:) + (/m_zero, lk(ii,2), m_zero/)
              temp = temp + &
                abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
            end do
            temp = temp * dk(2)
 
          case (integrate_kz)
            temp = m_zero
            do ii =1, ll(3)
              kkk(:) = kk(:) + (/m_zero, m_zero, lk(ii,3)/)
              temp = temp + &
                abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
            end do
            temp = temp * dk(3)
 
          end select
 
          write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
            units_from_atomic(sqrt(units_out%energy), kk(1)),&
            units_from_atomic(sqrt(units_out%energy), kk( 2)),&
            temp
 
 
          idone = idone +1
          call loct_progress_bar(idone, ntodo)
 
        end do
        write(iunit, *)
      end do
 
      write(stdout, '(1x)')
 
 
    end if
 
    call io_close(iunit)
 
    call pes_out_interpolator_end(cube_f, interp)
 
    safe_deallocate_a(idx)
    safe_deallocate_a(lk_)
 
    pop_sub(pes_out_velocity_map_cut)
  end subroutine pes_out_velocity_map_cut
 
  ! --------------------------------------------------------
  !
  !
  ! ---------------------------------------------------------
  subroutine pes_out_interpolator_init(namespace, pesK, Lk, ll, dim, cube_f, interp, pmesh)
    type(namespace_t),  intent(in)    :: namespace
    real(real64),       intent(in)    :: pesk(:,:,:)
    real(real64),       intent(in)    :: lk(:,:)
    integer,            intent(in)    :: ll(:)
    integer,            intent(in)    :: dim
    real(real64), allocatable, intent(out)   :: cube_f(:)
    type(qshep_t),      intent(out)   :: interp
    real(real64), optional,    intent(in)    :: pmesh(:,:,:,:)
 
    integer :: np, ii, ix, iy, iz
    real(real64)   :: kk(3)
    real(real64), allocatable ::  kx(:),ky(:),kz(:)
 
    push_sub(pes_out_interpolator_init)
 
 
    call messages_write("Initializing Qshep interpolator. Be patient it may take a while... ")
    call messages_info(namespace=namespace)
 
    np = ll(1)*ll(2)*ll(3)
 
    !check dim
    if (dim  <  2 .or. dim > 3) then
      message(1) = "This interpolator works only for 2 <= dim <= 3."
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(cube_f(1:np))
 
    safe_allocate(kx(1:np))
    safe_allocate(ky(1:np))
    safe_allocate(kz(1:np))
 
    cube_f = m_zero
    kx = m_zero
    ky = m_zero
    kz = m_zero
 
 
    ii=1
    if (present(pmesh)) then
      do ix = 1, ll(1)
        do iy = 1, ll(2)
          do iz = 1, ll(3)
 
            cube_f(ii) =  pesk(ix,iy,iz)
 
            kx(ii) = pmesh(ix,iy,iz,1)
            ky(ii) = pmesh(ix,iy,iz,2)
            kz(ii) = pmesh(ix,iy,iz,3)
 
            ii = ii +1
          end do
        end do
      end do
 
    else
 
      do ix = 1, ll(1)
        kk(1) = lk(ix, 1)
        do iy = 1, ll(2)
          kk(2) = lk(iy, 2)
          do iz = 1, ll(3)
            kk(3) = lk(iz, 3)
 
            cube_f(ii) =  pesk(ix,iy,iz)
 
            kx(ii) = kk(1)
            ky(ii) = kk(2)
            kz(ii) = kk(3)
 
            ii = ii +1
          end do
        end do
      end do
    end if
 
 
    select case (dim)
    case (2)
      call qshep_init(interp, i4_to_i8(np), cube_f, kx, ky)
    case (3)
      call qshep_init(interp, i4_to_i8(np), cube_f, kx, ky, kz)
    end select
 
    safe_deallocate_a(kx)
    safe_deallocate_a(ky)
    safe_deallocate_a(kz)
 
    call messages_write("done")
    call messages_new_line()
    call messages_info(namespace=namespace)
 
 
    pop_sub(pes_out_interpolator_init)
  end subroutine pes_out_interpolator_init
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine pes_out_interpolator_end(cube_f, interp)
    real(real64), allocatable, intent(inout) :: cube_f(:)
    type(qshep_t),      intent(inout) :: interp
 
    push_sub(pes_out_interpolator_end)
 
    call qshep_end(interp)
 
    safe_deallocate_a(cube_f)
 
    pop_sub(pes_out_interpolator_end)
  end subroutine pes_out_interpolator_end
 
end module pes_out_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: