pes_mask.F90

Go to the documentation of this file.

!! Copyright (C) 2006-2011 U. De Giovannini, M. Marques
!!
!! 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_mask_oct_m
  use absorbing_boundaries_oct_m
  use box_oct_m
  use box_sphere_oct_m
  use box_parallelepiped_oct_m
  use comm_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use debug_oct_m
  use density_oct_m
  use electron_space_oct_m
  use ext_partner_list_oct_m
  use fft_oct_m
  use global_oct_m
  use lalg_basic_oct_m
  use grid_oct_m
  use interaction_partner_oct_m
  use io_binary_oct_m
  use io_function_oct_m
  use io_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use lasers_oct_m
  use loct_oct_m
  use math_oct_m
  use mesh_cube_parallel_map_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
#if defined(HAVE_NETCDF)
  use netcdf
#endif
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use qshep_oct_m
  use restart_oct_m
  use space_oct_m
  use sort_oct_m
  use states_elec_oct_m
  use string_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use vtk_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                             &
    pes_mask_t,                         &
    pes_mask_init,                      &
    pes_mask_end,                       &
    pes_mask_calc,                      &
    pes_mask_mesh_to_cube,              &
    pes_mask_cube_to_mesh,              &
    pes_mask_generate_mask_function,    &
    pes_mask_x_to_k,                    &
    pes_mask_k_to_x,                    &
    pes_mask_read_info,                 &
    pes_mask_output_full_mapm,            &
    pes_mask_output_ar_spherical_cut_m,   &
    pes_mask_output_ar_plane_m,           &
    pes_mask_output_ar_polar_m,           &
    pes_mask_output_full_mapm_cut,        &
    pes_mask_output_power_totalm,         &
    pes_mask_load,                      &
    pes_mask_dump,                      &
    pes_mask_output,                    &
    pes_mask_pmesh,                     &
    pes_mask_map_from_states
 
 
  type pes_mask_t
    private
    complex(real64), allocatable, public :: k(:,:,:,:,:,:)
    !                                            !< mask%k(ll(1),ll(2),ll(3),st%d%dim, st%nst, st%nik)
 
    ! mesh- and cube-related stuff
    integer          :: np
    integer          :: ll(3)
    integer          :: fs_n_global(1:3)
    integer          :: fs_n(1:3)
    integer          :: fs_istart(1:3)
 
    real(real64)     :: spacing(3)
 
    type(mesh_t), pointer, public  :: mesh => null()   
    type(cube_t)                   :: cube
 
    real(real64), allocatable, public :: vec_pot(:,:)
 
    real(real64), allocatable, public :: Mk(:,:,:)
    type(cube_function_t)      :: cM
    real(real64), allocatable, public :: mask_R(:)
    integer                    :: shape
    real(real64), allocatable         :: ufn(:)
    logical                    :: user_def
 
    real(real64), allocatable, public :: Lk(:,:)
 
    real(real64)     :: enlarge(3)
    real(real64)     :: enlarge_2p(3)
 
    real(real64) :: start_time
    real(real64) :: energyMax
    real(real64) :: energyStep
 
    integer :: sw_evolve
    logical :: back_action
    logical :: add_psia
    logical :: filter_k
 
    integer :: mode
    integer :: pw_map_how
 
    type(fft_t)    :: fft
 
    type(mesh_cube_parallel_map_t) :: mesh_cube_map
 
 
  end type pes_mask_t
 
 
  integer, public, parameter ::        &
    PES_MASK_MODE_MASK         =   1,  &
    pes_mask_mode_backaction   =   2,  &
    pes_mask_mode_passive      =   3
 
  integer, parameter ::       &
    PW_MAP_FFT    =  3,  &    !< FFT - normally from fftw3
    pw_map_nfft        =  5,  &    
    pw_map_pfft        =  6,  &    
    pw_map_pnfft       =  7        
 
  integer, parameter ::       &
    M_SIN2            =  1,   &
    m_step            =  2,   &
    m_erf             =  3
 
  integer, parameter ::       &
    IN                =  1,   &
    out               =  2
 
  integer, public, 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_mask_init(mask, namespace, space, mesh, box, st, ext_partners, abs_boundaries, max_iter,dt)
    type(pes_mask_t),             intent(out) :: mask
    type(namespace_t),            intent(in)  :: namespace
    class(space_t),               intent(in)  :: space
    type(mesh_t), target,         intent(in)  :: mesh
    class(box_t),                 intent(in)  :: box
    type(states_elec_t),          intent(in)  :: st
    type(partner_list_t),         intent(in)  :: ext_partners
    type(absorbing_boundaries_t), intent(in)  :: abs_boundaries
    integer,                      intent(in)  :: max_iter
    real(real64),                 intent(in)  :: dt
 
    type(block_t) :: blk
 
    integer :: il, it, ll(3)
    real(real64) :: field(3)
    real(real64) :: deltae, maxe, pcutoff, tmp
    integer :: defaultmask,k1,k2,st1,st2
    integer :: cols_pesmask_block, idim, ip
 
    real(real64) :: xx(space%dim), r
    real(real64) :: ufn_re, ufn_im
    character(len=1024) :: user_def_expr
    type(lasers_t), pointer :: lasers
 
    push_sub(pes_mask_init)
 
    mask%mesh => mesh
 
    if (space%is_periodic()) then
      call messages_experimental("PES_mask with periodic dimensions", namespace=namespace)
    end if
 
 
    write(message(1),'(a,i1,a)') 'Info: Calculating PES using mask technique.'
    call messages_info(1, namespace=namespace)
 
 
    select type (box => box)
    type is (box_sphere_t)
    class default
      if (.not. space%is_periodic()) then
        message(1) = 'PhotoElectronSpectrum = pes_mask usually requires BoxShape = sphere.'
        message(2) = 'Unless you know what you are doing modify this parameter and rerun.'
        call messages_warning(2, namespace=namespace)
      end if
    end select
 
    if (abs_boundaries%abtype /= not_absorbing) then
      message(1) = 'PhotoElectronSpectrum = pes_mask already contains absorbing boundaries.'
      message(2) = 'Set AbsorbingBoundaries = no and rerun.'
      call messages_fatal(2, namespace=namespace)
    end if
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Calculation mode
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    !%Variable PESMaskMode
    !%Type integer
    !%Default mask_mode
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% PES calculation mode.
    !%Option mask_mode 1
    !% Mask method.
    !%Option fullmask_mode 2
    !% Full mask method. This includes a back action of the momentum-space states on the
    !% interaction region. This enables electrons to come back from the continuum.
    !%Option passive_mode 3
    !% Passive analysis of the wf. Simply analyze the plane-wave components of the
    !% wavefunctions on the region <i>r</i> > <i>R1</i>. This mode employs a step masking function by default.
    !%End
    call parse_variable(namespace, 'PESMaskMode', pes_mask_mode_mask, mask%mode)
    if (.not. varinfo_valid_option('PESMaskMode', mask%mode)) call messages_input_error(namespace, 'PESMaskMode')
    call messages_print_var_option("PESMaskMode", mask%mode, namespace=namespace)
 
    select case (mask%mode)
    case (pes_mask_mode_passive)
      defaultmask = m_step
      mask%back_action = .false.
 
    case (pes_mask_mode_backaction)
      defaultmask = m_sin2
      mask%back_action = .true.
      mask%mode = pes_mask_mode_mask
    case default
      defaultmask = m_sin2
      mask%back_action = .false.
    end select
 
 
    !%Variable PESMaskStartTime
    !%Type float
    !%Default -1.0
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% The time photoelectrons start to be recorded. In pump-probe simulations, this allows
    !% getting rid of an unwanted ionization signal coming from the pump.
    !% NOTE: This will enforce the mask boundary conditions for all times.
    !%End
    call parse_variable(namespace, 'PESMaskStartTime', -m_one, mask%start_time, unit = units_inp%time)
 
    !%Variable PESMaskPlaneWaveProjection
    !%Type integer
    !%Default fft_map
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% With the mask method, wavefunctions in the continuum are treated as plane waves.
    !% This variable sets how to calculate the plane-wave projection in the buffer
    !% region. We perform discrete Fourier transforms (DFT) in order to approximate
    !% a continuous Fourier transform. The major drawback of this approach is the built-in
    !% periodic boundary condition of DFT. Choosing an appropriate plane-wave projection
    !% for a given simulation in addition to <tt>PESMaskEnlargeFactor</tt> and
    !% <tt>PESMask2PEnlargeFactor</tt>will help to converge the results.
    !%
    !% NOTE: depending on the value of <tt>PESMaskMode</tt> <tt>PESMaskPlaneWaveProjection</tt>,
    !% may affect not only performance but also the time evolution of the density.
    !%Option fft_out 2
    !% FFT filtered in order to keep only outgoing waves. 1D only.
    !%Option fft_map 3
    !% FFT transform.
    !%Option nfft_map 5
    !% Non-equispaced FFT map.
    !%Option pfft_map 6
    !% Use PFFT library.
    !%Option pnfft_map 7
    !% Use PNFFT library.
    !%End
    call parse_variable(namespace, 'PESMaskPlaneWaveProjection', pw_map_fft, mask%pw_map_how)
 
    if (.not. varinfo_valid_option('PESMaskPlaneWaveProjection', mask%pw_map_how)) then
      call messages_input_error(namespace, 'PESMaskPlaneWaveProjection')
    end if
 
    call messages_print_var_option("PESMaskPlaneWaveProjection", mask%pw_map_how, namespace=namespace)
 
    if (mask%pw_map_how ==  pw_map_pfft .and. (.not. mask%mesh%parallel_in_domains)) then
      message(1)= "Trying to use PESMaskPlaneWaveProjection = pfft_map with no domain parallelization."
      message(2)= "Projection method changed to more efficient fft_map."
      call messages_warning(2, namespace=namespace)
      mask%pw_map_how = pw_map_fft
    end if
 
    if (mask%pw_map_how ==  pw_map_pnfft .and. (.not. mask%mesh%parallel_in_domains)) then
      message(1)= "Trying to use PESMaskPlaneWaveProjection = pnfft_map with no domain parallelization."
      message(2)= "Projection method changed to more efficient nfft_map."
      call messages_warning(2, namespace=namespace)
      mask%pw_map_how = pw_map_nfft
    end if
 
#if !defined(HAVE_NFFT)
    if (mask%pw_map_how ==  pw_map_nfft) then
      message(1) = "PESMaskPlaneWaveProjection = nfft_map requires FINUFFT but that library was not linked."
      call messages_fatal(1, namespace=namespace)
    end if
#endif
 
#if !defined(HAVE_PFFT)
    if (mask%pw_map_how ==  pw_map_pfft) then
      message(1) = "PESMaskPlaneWaveProjection = pfft_map requires PFFT but that library was not linked."
      call messages_fatal(1, namespace=namespace)
    end if
#endif
 
#if !defined(HAVE_PNFFT)
    if (mask%pw_map_how ==  pw_map_pnfft) then
      message(1) = "PESMaskPlaneWaveProjection = pnfft_map requires PNFFT but that library was not linked."
      call messages_fatal(1, namespace=namespace)
    end if
#endif
 
 
    !%Variable PESMaskEnlargeFactor
    !%Type float
    !%Default 1
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% Mask box enlargement level. Enlarges the mask bounding box by a <tt>PESMaskEnlargeFactor</tt>.
    !% This helps to avoid wavefunction wrapping at the boundaries.
    !%End
 
    mask%enlarge = m_one
    call parse_variable(namespace, 'PESMaskEnlargeFactor', m_one, mask%enlarge(1))
 
    if (.not. is_close(mask%enlarge(1), m_one)) then
 
      mask%enlarge(space%periodic_dim + 1:space%dim) = mask%enlarge(1)
      mask%enlarge(1:space%periodic_dim) = m_one
 
      if (space%is_periodic()) then
        call messages_print_var_value("PESMaskEnlargeFactor", mask%enlarge(1:space%dim), namespace=namespace)
      else
        call messages_print_var_value("PESMaskEnlargeFactor", mask%enlarge(1), namespace=namespace)
      end if
 
    end if
    if (mask%enlarge(1) < m_one) then
      message(1) = "PESMaskEnlargeFactor must be bigger than one."
      call messages_fatal(1, namespace=namespace)
    end if
 
    call messages_obsolete_variable(namespace, 'PESMaskEnlargeLev', 'PESMaskEnlargeFactor')
 
    !%Variable PESMask2PEnlargeFactor
    !%Type float
    !%Default 1.0
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% Mask two points enlargement factor. Enlarges the mask box by adding two
    !% points at the edges of the box in each direction (x,y,z) at a distance
    !% L=Lb*<tt>PESMask2PEnlargeFactor</tt> where <i>Lb</i> is the box size.
    !% This allows to run simulations with an additional void space at a price of
    !% adding few points. The Fourier space associated with the new box is restricted
    !% by the same factor.
    !%
    !% Note: needs <tt> PESMaskPlaneWaveProjection = nfft_map or pnfft_map </tt>.
    !%End
 
    mask%enlarge_2p = m_one
    call parse_variable(namespace, 'PESMask2PEnlargeFactor', m_one, mask%enlarge_2p(1))
 
 
    if (.not. is_close(mask%enlarge_2p(1), m_one)) then
 
      mask%enlarge_2p(space%periodic_dim + 1:space%dim) = mask%enlarge_2p(1)
      mask%enlarge_2p(1:space%periodic_dim) = m_one
 
      if (space%is_periodic()) then
        call messages_print_var_value("PESMask2PEnlargeFactor", mask%enlarge_2p(1:space%dim), namespace=namespace)
      else
        call messages_print_var_value("PESMask2PEnlargeFactor", mask%enlarge_2p(1), namespace=namespace)
      end if
 
      if (mask%pw_map_how /=  pw_map_nfft .and. mask%pw_map_how /=  pw_map_pnfft) then
        message(1) = "PESMask2PEnlargeFactor requires PESMaskPlaneWaveProjection = nfft_map"
        message(2) = "or pnfft_map in order to run properly."
        call messages_warning(2, namespace=namespace)
      end if
    end if
 
    if (mask%enlarge_2p(1) < m_one) then
      message(1) = "PESMask2PEnlargeFactor must be bigger than one."
      call messages_fatal(1, namespace=namespace)
    end if
 
    call messages_obsolete_variable(namespace, 'PESMaskNFFTEnlargeLev', 'PESMask2PEnlargeFactor')
 
 
    mask%ll = 1
    mask%spacing = -m_one
 
    mask%spacing(1:space%dim) = mesh%spacing(1:space%dim)
    mask%ll(1:space%dim) = mesh%idx%ll(1:space%dim)
 
    !Enlarge the cube region
    mask%ll(1:space%dim) = int(mask%ll(1:space%dim) * mask%enlarge(1:space%dim))
 
    mask%np = mesh%np_part ! the mask is local
 
    select case (mask%pw_map_how)
 
    case (pw_map_pfft)
      assert(mask%mesh%parallel_in_domains)
      call cube_init(mask%cube, mask%ll, namespace, space, &
        mesh%spacing, mesh%coord_system, &
        fft_type = fft_complex, fft_library = fftlib_pfft, nn_out = ll, &
        mpi_grp = mask%mesh%mpi_grp, need_partition=.true.)
      call cube_init_cube_map(mask%cube, mesh)
 
 
!       mask%ll(1) = mask%cube%fs_n(3)
!       mask%ll(2) = mask%cube%fs_n(1)
!       mask%ll(3) = mask%cube%fs_n(2)
!     Note: even if tempting, setting
!     mask%ll(1:3) = mask%cube%fs_n(1:3)
!     results in the wrong index mapping! (1->3, 2->1, 3->2)
!     Well.. very much against intuition it turns out to be
      mask%ll(1:3) = mask%cube%rs_n(1:3)
 
      mask%fft = mask%cube%fft
      if (mask%mesh%parallel_in_domains .and. mask%cube%parallel_in_domains) then
        call mesh_cube_parallel_map_init(mask%mesh_cube_map, mask%mesh, mask%cube)
      end if
 
    case (pw_map_fft)
      call cube_init(mask%cube, mask%ll, namespace, space, &
        mesh%spacing, mesh%coord_system, &
        fft_type = fft_complex, fft_library = fftlib_fftw, &
        nn_out = ll)
      call cube_init_cube_map(mask%cube, mesh)
      mask%ll = ll
      mask%fft = mask%cube%fft
 
 
    case (pw_map_nfft)
 
      !NFFT initialization
      ! we just add 2 points for the enlarged region
      if (.not. is_close(mask%enlarge_2p(1), m_one)) mask%ll(1:space%dim) = mask%ll(1:space%dim) + 2
 
      call cube_init(mask%cube, mask%ll, namespace, space, &
        mesh%spacing, mesh%coord_system, &
        fft_type = fft_complex, fft_library = fftlib_nfft, &
        nn_out = ll, tp_enlarge = mask%enlarge_2p, nthreads=1)
      call cube_init_cube_map(mask%cube, mesh)
 
      mask%ll = ll
      mask%fft = mask%cube%fft
 
 
    case (pw_map_pnfft)
 
      if (.not. is_close(mask%enlarge_2p(1), m_one))  mask%ll(1:space%dim) = mask%ll(1:space%dim) + 2
 
      call cube_init(mask%cube, mask%ll, namespace, space, &
        mask%mesh%spacing, mask%mesh%coord_system, &
        fft_type = fft_complex, fft_library = fftlib_pnfft, &
        nn_out = ll, tp_enlarge = mask%enlarge_2p, &
        mpi_grp = mask%mesh%mpi_grp, need_partition=.true.)
      call cube_init_cube_map(mask%cube, mesh)
 
      mask%ll(1:3) = mask%cube%fs_n(1:3)
 
      mask%fft = mask%cube%fft
      if (mask%mesh%parallel_in_domains .and. mask%cube%parallel_in_domains) then
        call mesh_cube_parallel_map_init(mask%mesh_cube_map, mask%mesh, mask%cube)
      end if
 
    case default
      !Program should die before coming here
      write(message(1),'(a)') "PESMaskPlaneWaveProjection unrecognized option."
      call messages_fatal(1, namespace=namespace)
 
    end select
 
    !Indices
 
    if (mask%pw_map_how == pw_map_pfft) then
      mask%fs_istart = mask%cube%rs_istart
      mask%fs_n = mask%cube%rs_n
      mask%fs_n_global = mask%cube%rs_n_global
    else
      mask%fs_istart = mask%cube%fs_istart
      mask%fs_n = mask%cube%fs_n
      mask%fs_n_global = mask%cube%fs_n_global
    end if
 
    !Allocations
    call zcube_function_alloc_rs(mask%cube, mask%cM, force_alloc = .true.)
 
    safe_allocate(mask%Lk(1:maxval(mask%fs_n_global(:)),1:3))
    mask%Lk(:,:) = m_zero
 
    st1 = st%st_start
    st2 = st%st_end
    k1 = st%d%kpt%start
    k2 = st%d%kpt%end
    safe_allocate(mask%k(1:mask%ll(1),1:mask%ll(2),1:mask%ll(3),1:st%d%dim,st1:st2,k1:k2))
    mask%k = m_z0
 
 
    ! generate the map between mesh and cube
    call pes_mask_generate_lk(mask) ! generate the physical momentum vector
 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Mask Function options
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    safe_allocate(mask%mask_R(1:2))
 
    !%Variable PESMaskShape
    !%Type integer
    !%Default m_sin2
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% The mask function shape.
    !%Option m_sin2 1
    !% sin2 mask.
    !%Option m_step 2
    !%step function.
    !%Option m_erf 3
    !%Error function. Not Implemented.
    !%End
    call parse_variable(namespace, 'PESMaskShape', defaultmask, mask%shape)
    if (.not. varinfo_valid_option('PESMaskShape', mask%shape)) call messages_input_error(namespace, 'PESMaskShape')
    call messages_print_var_option("PESMaskShape", mask%shape, namespace=namespace)
 
    !%Variable PESMaskSize
    !%Type block
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% Set the size of the mask function.
    !% Here you can set the inner (R1) and outer (R2) radius by setting
    !% the block as follows:
    !%
    !% <tt>%PESMaskSize
    !% <br>&nbsp;&nbsp; R1 | R2 | "user-defined"
    !% <br>%</tt>
    !%
    !% The optional 3rd column is a user-defined expression for the mask
    !% function. For example, <i>r</i> creates a spherical mask (which is the
    !% default for <tt>BoxShape = sphere</tt>). Note, values R2 larger than
    !% the box size may lead in this case to unexpected reflection
    !% behaviours.
    !%End
    cols_pesmask_block = 0
    if (parse_block(namespace, 'PESMaskSize', blk) == 0) then
      cols_pesmask_block = parse_block_cols(blk, 0)
    end if
 
    mask%user_def = .false.
 
    select case (cols_pesmask_block)
    case (0)
      select type (box => box)
      type is (box_sphere_t)
        mask%mask_R(1) = box%radius/m_two
        mask%mask_R(2) = box%radius
        message(1) = "Input: PESMaskSize R(1) and R(2) not specified. Using default values for spherical mask."
      type is (box_parallelepiped_t)
        mask%mask_R(1) = box%half_length(1)/m_two
        mask%mask_R(2) = box%half_length(1)
        message(1) = "Input: PESMaskSize R(1) and R(2) not specified. Using default values for cubic mask."
      end select
      call messages_info(1, namespace=namespace)
    case (1)
      call parse_block_float(blk, 0, 0, mask%mask_R(1), units_inp%length)
      select type (box => box)
      type is (box_sphere_t)
        mask%mask_R(2) = box%radius
        message(1) = "Input: PESMaskSize R(2) not specified. Using default value for spherical mask."
      type is (box_parallelepiped_t)
        mask%mask_R(2) = box%half_length(1)
        message(1) = "Input: PESMaskSize R(2) not specified. Using default value for cubic mask."
      end select
      call messages_info(1, namespace=namespace)
    case (2)
      call parse_block_float(blk, 0, 0, mask%mask_R(1), units_inp%length)
      call parse_block_float(blk, 0, 1, mask%mask_R(2), units_inp%length)
 
      select type (box => box)
      type is (box_sphere_t)
        if (mask%mask_R(2) > box%radius)  mask%mask_R(2) = box%radius
        message(1) = "Info: using spherical mask."
      type is (box_parallelepiped_t)
        if (mask%mask_R(2) > box%half_length(1)) mask%mask_R(2) = box%half_length(1)
        message(1) = "Info: using cubic mask."
      end select
      call messages_info(1, namespace=namespace)
 
    case (3)
      mask%user_def = .true.
      safe_allocate(mask%ufn(1:mask%np))
      mask%ufn = m_zero
      call parse_block_float(blk, 0, 0, mask%mask_R(1), units_inp%length)
      call parse_block_float(blk, 0, 1, mask%mask_R(2), units_inp%length)
      call parse_block_string(blk, 0, 2, user_def_expr)
      do ip = 1, mask%np
        xx(1:space%dim) = mesh%x(1:space%dim, ip)
        r = units_from_atomic(units_inp%length, norm2(xx(1:space%dim)))
        do idim = 1, space%dim
          xx(idim) = units_from_atomic(units_inp%length, xx(idim))
        end do
        call parse_expression(ufn_re, ufn_im, space%dim, xx, r, m_zero, user_def_expr)
        mask%ufn(ip) = ufn_re
      end do
      message(1) = "Input: using user-defined mask function from expression:"
      write(message(2),'(a,a)') '   R = ', trim(user_def_expr)
      call messages_info(2, namespace=namespace)
    end select
 
    call parse_block_end(blk)
 
    write(message(1),'(a,es10.3,3a)') &
      "  R1 = ", units_from_atomic(units_inp%length, mask%mask_R(1)),&
      ' [', trim(units_abbrev(units_inp%length)), ']'
    write(message(2),'(a,es10.3,3a)') &
      "  R2 = ", units_from_atomic(units_inp%length, mask%mask_R(2)),&
      ' [', trim(units_abbrev(units_inp%length)), ']'
    call messages_info(2, namespace=namespace)
 
 
    call pes_mask_generate_mask(mask, namespace, mesh)
 
    !%Variable PESMaskFilterCutOff
    !%Type float
    !%Default -1
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% In calculation with <tt>PESMaskMode = fullmask_mode</tt> and NFFT, spurious frequencies
    !% may lead to numerical instability of the algorithm. This option gives the possibility
    !% to filter out the unwanted components by setting an energy cut-off.
    !% If <tt>PESMaskFilterCutOff = -1</tt> no filter is applied.
    !%End
    call parse_variable(namespace, 'PESMaskFilterCutOff', -m_one, pcutoff, unit = units_inp%energy)
 
    mask%filter_k = .false.
 
    if (pcutoff > m_zero) then
      call messages_print_var_value("PESMaskFilterCutOff", pcutoff, unit = units_out%energy, namespace=namespace)
      mask%filter_k = .true.
 
      safe_allocate(mask%Mk(1:mask%ll(1),1:mask%ll(2),1:mask%ll(3)))
 
      call pes_mask_generate_filter(mask,pcutoff)
    end if
 
!!  Output
 
    !%Variable PESMaskIncludePsiA
    !%Type logical
    !%Default false
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% Add the contribution of <math>\Psi_A</math> in the mask region to the photo-electron spectrum.
    !% Literally adds the Fourier components of:
    !% <math>\Theta(r-R1) \Psi_A(r)</math>
    !% with <math>\Theta</math> being the Heaviside step function.
    !% With this option PES will contain all the contributions starting from the inner
    !% radius <math>R1</math>. Use this option to improve convergence with respect to the box size
    !% and total simulation time.
    !% Note: Carefully choose <math>R1</math> in order to avoid contributions from returning electrons.
    !%End
    call parse_variable(namespace, 'PESMaskIncludePsiA', .false., mask%add_psia)
    if (mask%add_psia) then
      message(1)= "Input: Include contribution from Psi_A."
      call messages_info(1, namespace=namespace)
    end if
 
 
    !%Variable PESMaskSpectEnergyMax
    !%Type float
    !%Default maxval(mask%Lk)<math>^2/2</math>
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% The maximum energy for the PES spectrum.
    !%End
    maxe = m_epsilon
    do idim = 1, space%dim
      tmp = maxval(mask%Lk(1:mask%ll(idim),1:space%dim))**2/m_two
      if (tmp > maxe) maxe = tmp
    end do
    call parse_variable(namespace, 'PESMaskSpectEnergyMax', maxe, mask%energyMax, unit = units_inp%energy)
    call messages_print_var_value("PESMaskSpectEnergyMax", mask%energyMax, unit = units_out%energy, namespace=namespace)
 
    !%Variable PESMaskSpectEnergyStep
    !%Type float
    !%Section Time-Dependent::PhotoElectronSpectrum
    !%Description
    !% The PES spectrum energy step.
    !%End
    deltae = minval(mask%Lk(2,1:space%dim) - mask%Lk(1,1:space%dim))**2/m_two
    call parse_variable(namespace, 'PESMaskSpectEnergyStep', deltae, mask%energyStep, unit = units_inp%energy)
    call messages_print_var_value("PESMaskSpectEnergyStep", mask%energyStep, unit = units_out%energy, namespace=namespace)
 
 
!!  Set external fields
 
    safe_allocate(mask%vec_pot(0:max_iter,1:3))
    mask%vec_pot=m_zero
 
    lasers => list_get_lasers(ext_partners)
    if(associated(lasers)) then
      do il = 1, lasers%no_lasers
        select case (laser_kind(lasers%lasers(il)))
        case (e_field_vector_potential)
          do it = 1, max_iter
            field=m_zero
            call laser_field(lasers%lasers(il), field, it*dt)
            ! We must sum with a -1 sign to account for the
            ! electron charge.
            mask%vec_pot(it,:)= mask%vec_pot(it,:) - field(:)
          end do
 
        case default
          write(message(1),'(a)') 'PESMask should work only with TDExternalFields = vector_potential.'
          write(message(2),'(a)') 'Unless PESMaskMode = passive_mode the results are likely to be wrong. '
          call messages_warning(2, namespace=namespace)
 
        end select
      end do
    end if
 
 
 
    ! Compensate the sign for the forward <-> backward FT inversion
    ! used with NFFT. See pes_mask_X_to_K comment for more info
    if (mask%pw_map_how == pw_map_pnfft .or. &
      mask%pw_map_how == pw_map_nfft) mask%Lk = - mask%Lk
 
 
    pop_sub(pes_mask_init)
  end subroutine pes_mask_init
 
 
  ! ---------------------------------------------------------
  subroutine pes_mask_end(mask)
    type(pes_mask_t), intent(inout) :: mask
 
    push_sub(pes_mask_end)
 
    safe_deallocate_a(mask%k)
 
    safe_deallocate_a(mask%vec_pot)
    safe_deallocate_a(mask%mask_R)
    safe_deallocate_a(mask%Lk)
 
 
    if (mask%filter_k) then
      safe_deallocate_a(mask%Mk)
    end if
 
    if (mask%mesh%parallel_in_domains .and. mask%cube%parallel_in_domains) then
      call mesh_cube_parallel_map_end(mask%mesh_cube_map)
    end if
 
    call zcube_function_free_rs(mask%cube, mask%cM)
    call cube_end(mask%cube)
 
    if (mask%user_def) then
      safe_deallocate_a(mask%ufn)
    end if
 
    pop_sub(pes_mask_end)
  end subroutine pes_mask_end
 
 
  ! --------------------------------------------------------
  subroutine pes_mask_generate_lk(mask)
    type(pes_mask_t),      intent(inout) :: mask
 
    integer :: ii, dim
 
    push_sub(pes_mask_generate_lk)
 
    dim = mask%mesh%box%dim
 
    do ii = 1, maxval(mask%ll(:))
      mask%Lk(ii,1:dim)= matmul(mask%cube%latt%klattice_primitive(1:dim,1:dim), mask%cube%Lfs(ii,1:dim))
    end do
 
    pop_sub(pes_mask_generate_lk)
  end subroutine pes_mask_generate_lk
 
  ! ======================================
  ! ======================================
  subroutine pes_mask_generate_filter(mask,cutOff)
    type(pes_mask_t), intent(inout) :: mask
    real(real64),     intent(in)    ::cutoff
 
 
    integer :: kx,ky,kz, power
    real(real64)   :: kk(3), ee, emax
 
    push_sub(pes_mask_generate_filter)
 
    mask%Mk = m_zero
 
    emax = maxval(mask%Lk(:,:))**2 / m_two
 
    power = 8
 
    do kx = 1, mask%ll(1)
      kk(1) = mask%Lk(kx,1)
      do ky = 1, mask%ll(2)
        kk(2) = mask%Lk(ky,2)
        do kz = 1, mask%ll(3)
          kk(3) = mask%Lk(kz,3)
 
          ee = sum(kk(1:mask%mesh%box%dim)**2) / m_two
 
          if (ee > cutoff .and. ee < emax) then
            mask%Mk(kx,ky,kz) = m_one * sin((emax-ee) * m_pi / (m_two * (cutoff)))**power
          else
            if (ee <= cutoff) mask%Mk(kx,ky,kz) = m_one
          end if
 
        end do
      end do
    end do
 
 
    pop_sub(pes_mask_generate_filter)
  end subroutine pes_mask_generate_filter
 
 
  ! --------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine pes_mask_generate_mask(mask, namespace, mesh)
    type(pes_mask_t),  intent(inout) :: mask
    type(namespace_t), intent(in)    :: namespace
    type(mesh_t),      intent(in)    :: mesh
 
    push_sub(pes_mask_generate_mask)
 
    call pes_mask_generate_mask_function(mask, namespace, mesh, mask%shape, mask%mask_R)
 
    pop_sub(pes_mask_generate_mask)
 
  end subroutine pes_mask_generate_mask
 
  ! --------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine pes_mask_generate_mask_function(mask, namespace, mesh, shape, R, mask_sq)
    type(pes_mask_t),  intent(inout) :: mask
    type(namespace_t), intent(in)    :: namespace
    type(mesh_t),      intent(in)    :: mesh
    integer,           intent(in)    :: shape
    real(real64),      intent(in)    :: r(2)
    real(real64), optional,   intent(out)   :: mask_sq(:,:,:)
 
    integer :: ip, dir
    real(real64)   :: width
    real(real64)   :: xx(1:mesh%box%dim), rr, dd, ddv(1:mesh%box%dim), tmp(1:mesh%box%dim)
    complex(real64), allocatable :: mask_fn(:)
 
    push_sub(pes_mask_generate_mask_function)
 
    assert(mask%cM%batch_size == 1)
 
    ! generate the mask function on the mesh
    safe_allocate(mask_fn(1:mask%np))
 
    mask_fn = m_zero
    width = r(2) - r(1)
    xx = m_zero
 
    select case (shape)
    case (m_sin2)
      do ip = 1, mask%np
        call mesh_r(mesh, ip, rr, coords=xx)
 
        if (mask%user_def) then
          dd = mask%ufn(ip) - r(1)
          if (dd > m_zero) then
            if (mask%ufn(ip) < r(2)) then
              mask_fn(ip) = m_one * sin(dd * m_pi / (m_two * (width)))**2
            else
              mask_fn(ip) = m_one
            end if
          end if
 
        else ! mask%user_def == .false.
 
          select type (box => mesh%box)
          type is (box_sphere_t)
            dd = rr -  r(1)
            if (dd > m_zero) then
              if (dd  <  width) then
                mask_fn(ip) = m_one * sin(dd * m_pi / (m_two * (width)))**2
              else
                mask_fn(ip) = m_one
              end if
            end if
 
          type is (box_parallelepiped_t)
 
            ! We are filling from the center opposite to the spherical case
            tmp = m_one
            mask_fn(ip) = m_one
            ddv = abs(xx) -  r(1)
            do dir=1, mesh%box%dim
              if (ddv(dir) > m_zero) then
                if (ddv(dir)  <  width) then
                  tmp(dir) = m_one - sin(ddv(dir) * m_pi / (m_two * width))**2
                else
                  tmp(dir) = m_zero
                end if
              end if
              mask_fn(ip) = mask_fn(ip) * tmp(dir)
            end do
            mask_fn(ip) = m_one - mask_fn(ip)
 
          end select
        end if
      end do
 
    case (m_step)
      do ip = 1, mask%np
        call mesh_r(mesh, ip, rr, coords=xx)
        dd = rr - r(1)
        if (dd > m_zero) then
          if (dd  <  width) then
            mask_fn(ip) = m_one
          else
            mask_fn(ip) = m_zero
          end if
        end if
      end do
 
    case (m_erf)
 
    case default
      message(1)="PhotoElectronSpectrum = pes_mask. Unrecognized mask type."
      call messages_fatal(1, namespace=namespace)
    end select
 
    mask_fn(:) = m_one - mask_fn(:)
    call pes_mask_mesh_to_cube(mask, mask_fn, mask%cM)
 
    if (present(mask_sq)) mask_sq = real(mask%cM%zRS(:,:,:,1), real64)
 
    safe_deallocate_a(mask_fn)
 
    pop_sub(pes_mask_generate_mask_function)
 
  end subroutine pes_mask_generate_mask_function
 
 
  ! --------------------------------------------------------
  subroutine pes_mask_apply_mask(mask, st)
    type(states_elec_t), intent(inout) :: st
    type(pes_mask_t),    intent(in)    :: mask
 
    integer :: ik, ist, idim, ip
    complex(real64), allocatable :: mmask(:), psi(:)
 
    push_sub(pes_mask_apply_mask)
    safe_allocate(mmask(1:mask%mesh%np))
    safe_allocate(psi(1:mask%mesh%np))
 
    call pes_mask_cube_to_mesh(mask, mask%cM, mmask)
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ist = st%st_start, st%st_end
        do idim = 1, st%d%dim
          call states_elec_get_state(st, mask%mesh, idim, ist, ik, psi)
          !$omp parallel do
          do ip = 1, mask%mesh%np
            psi(ip) = psi(ip) * mmask(ip)
          end do
          !$omp end parallel do
          call states_elec_set_state(st, mask%mesh, idim, ist, ik, psi)
        end do
      end do
    end do
 
    safe_deallocate_a(mmask)
 
    pop_sub(pes_mask_apply_mask)
  end subroutine pes_mask_apply_mask
 
 
 
 
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine pes_mask_volkov_time_evolution_wf(mask, space, mesh, kpoints, dt, iter, wf, ikpoint)
    type(pes_mask_t), intent(in)    :: mask
    class(space_t),   intent(in)    :: space
    type(mesh_t),     intent(in)    :: mesh
    type(kpoints_t),  intent(in)    :: kpoints
    real(real64),     intent(in)    :: dt
    integer,          intent(in)    :: iter
    complex(real64),  intent(inout) :: wf(:,:,:)
    integer,          intent(in)    :: ikpoint
 
    integer ::  ix, iy, iz
    real(real64) :: vec
    real(real64) :: kk(1:3), kpoint(1:3)
 
    push_sub(pes_mask_volkov_time_evolution_wf)
 
    kpoint = m_zero
    if (space%is_periodic()) then
      kpoint(1:mesh%box%dim) = kpoints%get_point(ikpoint)
    end if
 
    do ix = 1, mask%ll(1)
      kk(1) = mask%Lk(ix + mask%fs_istart(1) - 1, 1)
      do iy = 1, mask%ll(2)
        kk(2) = mask%Lk(iy + mask%fs_istart(2) - 1, 2)
        do iz = 1, mask%ll(3)
          kk(3) = mask%Lk(iz + mask%fs_istart(3) - 1, 3)
          ! The k-points have the same sign as the vector potential consistently
          ! with what is done to generate the phase (hm%phase) in hamiltonian_elec_update()
          vec = sum(( kk(1:mesh%box%dim) &
            - kpoint(1:mesh%box%dim) &
            - mask%vec_pot(iter,1:mesh%box%dim)/p_c)**2) / m_two
          wf(ix, iy, iz) = wf(ix, iy, iz) * exp(-m_zi * dt * vec)
        end do
      end do
    end do
 
 
!     do ix = 1, mask%ll(1)
!       KK(1) = mask%Lk(ix + mask%fs_istart(1) - 1)
!       do iy = 1, mask%ll(2)
!         KK(2) = mask%Lk(iy + mask%fs_istart(2) - 1)
!         do iz = 1, mask%ll(3)
!           KK(3) = mask%Lk(iz + mask%fs_istart(3) - 1)
!           vec = sum(( KK(1:mesh%box%dim) - mask%vec_pot(iter,1:mesh%box%dim)/P_C)**2) / M_TWO
!           wf(ix, iy, iz) = wf(ix, iy, iz) * exp(-M_zI * dt * vec)
!         end do
!       end do
!     end do
 
    pop_sub(pes_mask_volkov_time_evolution_wf)
  end subroutine pes_mask_volkov_time_evolution_wf
 
 
 
 
  !---------------------------------------------------------
  !---------------------------------------------------------
  subroutine pes_mask_x_to_k(mask, wfin, wfout)
    type(pes_mask_t), intent(in)    :: mask
    complex(real64),  intent(inout) :: wfin(:,:,:)
    complex(real64),  intent(out)   :: wfout(:,:,:)
 
    type(cube_function_t) :: cf_tmp
    real(real64)          :: norm
 
    call profiling_in("PESMASK_X_to_K")
 
    push_sub(pes_mask_x_to_k)
 
    assert(cf_tmp%batch_size == 1)
 
    wfout=m_z0
 
    select case (mask%pw_map_how)
 
    case (pw_map_fft)
      call zfft_forward(mask%cube%fft, wfin, wfout)
 
    case (pw_map_nfft)
      ! By definition the forward non-equispaced transform generates a function defined on an
      ! unstructured grid from its Fourier components (defined on a cubic equispaced grid).
      ! This is not what we want since for us it is the real-space grid the one that can be
      ! unstructured.
      ! In order to preserve the possibility to have an unstructured rs-grid we use the
      ! backward transform and flip the sign of all the momenta (Lk = -Lk).
      call zfft_backward(mask%cube%fft, wfin, wfout, norm)
      wfout = wfout * norm
 
    case (pw_map_pfft)
      call zcube_function_alloc_rs(mask%cube, cf_tmp)
      call cube_function_alloc_fs(mask%cube, cf_tmp)
      cf_tmp%zRs(:,:,:,1) = wfin
      call zfft_forward(mask%cube%fft, cf_tmp%zRs, cf_tmp%fs)
      wfout = cf_tmp%fs(:,:,:,1)
      call zcube_function_free_rs(mask%cube, cf_tmp)
      call cube_function_free_fs(mask%cube, cf_tmp)
 
    case (pw_map_pnfft)
      call zfft_backward(mask%cube%fft, wfin, wfout, norm)
      wfout = wfout * norm
 
    case default
 
    end select
 
    pop_sub(pes_mask_x_to_k)
    call profiling_out("PESMASK_X_to_K")
 
  end subroutine pes_mask_x_to_k
 
  ! ------------------------------------------------
  subroutine pes_mask_k_to_x(mask, wfin, wfout)
    type(pes_mask_t), intent(in)    :: mask
    complex(real64),  intent(inout) :: wfin(:,:,:)
    complex(real64),  intent(out)   :: wfout(:,:,:)
 
    type(cube_function_t) :: cf_tmp
    real(real64) :: norm
 
    call profiling_in("PESMASK_K_toX")
 
    push_sub(pes_mask_k_to_x)
 
    assert(cf_tmp%batch_size == 1)
 
    wfout=m_z0
 
    select case (mask%pw_map_how)
 
    case (pw_map_fft)
      call zfft_backward(mask%cube%fft, wfin, wfout)
 
    case (pw_map_nfft)
      call zfft_forward(mask%cube%fft, wfin, wfout, norm)
      wfout = wfout / norm
 
    case (pw_map_pfft)
      call zcube_function_alloc_rs(mask%cube, cf_tmp)
      call cube_function_alloc_fs(mask%cube, cf_tmp)
      cf_tmp%fs(:,:,:,1)  = wfin
      call zfft_backward(mask%cube%fft, cf_tmp%fs, cf_tmp%zRs)
      wfout = cf_tmp%zRs(:,:,:,1)
      call zcube_function_free_rs(mask%cube, cf_tmp)
      call cube_function_free_fs(mask%cube, cf_tmp)
 
    case (pw_map_pnfft)
      call zfft_forward(mask%cube%fft, wfin, wfout, norm)
      wfout = wfout / norm
 
    case default
 
    end select
 
    pop_sub(pes_mask_k_to_x)
 
    call profiling_out("PESMASK_K_toX")
 
  end subroutine pes_mask_k_to_x
 
  !---------------------------------------------------------
  subroutine pes_mask_mesh_to_cube(mask, mf, cf)
    type(pes_mask_t),            intent(in)    :: mask
    complex(real64), contiguous, intent(in)    :: mf(:)
    type(cube_function_t),       intent(inout) :: cf
 
    push_sub(pes_mask_mesh_to_cube)
 
    if (mask%cube%parallel_in_domains) then
      call zmesh_to_cube_parallel(mask%mesh, mf, mask%cube, cf, mask%mesh_cube_map)
    else
      call zmesh_to_cube(mask%mesh, mf, mask%cube, cf)
    end if
 
    pop_sub(pes_mask_mesh_to_cube)
  end subroutine pes_mask_mesh_to_cube
 
 
  !---------------------------------------------------------
  subroutine pes_mask_cube_to_mesh(mask, cf, mf)
    type(pes_mask_t),            intent(in) :: mask
    complex(real64), contiguous, intent(out):: mf(:)
    type(cube_function_t),       intent(in) :: cf
 
    push_sub(pes_mask_cube_to_mesh)
 
    if (mask%cube%parallel_in_domains) then
      call zcube_to_mesh_parallel(mask%cube, cf, mask%mesh, mf, mask%mesh_cube_map)
    else
      call zcube_to_mesh(mask%cube, cf, mask%mesh, mf)
    end if
 
    pop_sub(pes_mask_cube_to_mesh)
  end subroutine pes_mask_cube_to_mesh
 
 
  !---------------------------------------------------------
  !
  !            Performs all the dirty work
  !
  !---------------------------------------------------------
  subroutine pes_mask_calc(mask, namespace, space, mesh, st, kpoints, dt, iter)
    type(pes_mask_t),    intent(inout) :: mask
    type(namespace_t),   intent(in)    :: namespace
    class(space_t),      intent(in)    :: space
    type(mesh_t),        intent(in)    :: mesh
    type(states_elec_t), intent(inout) :: st
    type(kpoints_t),     intent(in)    :: kpoints
    real(real64),        intent(in)    :: dt
    integer,             intent(in)    :: iter
 
    integer :: idim, ist, ik
    type(cube_function_t):: cf1, cf2
    complex(real64), allocatable :: mf(:), psi(:)
 
    real(real64) :: time
 
 
    call profiling_in("PESMASK_calc")
 
    push_sub(pes_mask_calc)
 
    time = iter *dt
 
    if (time > mask%start_time) then ! record photoelectrons only after mask%start_time
 
      call zcube_function_alloc_rs(mask%cube, cf1, force_alloc = .true.)
      call cube_function_alloc_fs(mask%cube, cf1, force_alloc = .true.)
      call zcube_function_alloc_rs(mask%cube, cf2, force_alloc = .true.)
      call cube_function_alloc_fs(mask%cube, cf2, force_alloc = .true.)
 
      select case (mask%mode)
      case (pes_mask_mode_mask)
        if (mask%back_action .eqv. .true.) then
          safe_allocate(mf(1:mask%mesh%np_part))
        end if
      end select
 
      safe_allocate(psi(1:mask%mesh%np_part))
 
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = st%st_start, st%st_end
          do idim = 1, st%d%dim
 
            call states_elec_get_state(st, mask%mesh, idim, ist, ik, psi)
 
            cf1%zRs = m_z0
            cf2%zRS= m_z0
            cf1%Fs = m_z0
            cf2%Fs = m_z0
 
            call pes_mask_mesh_to_cube(mask, psi, cf1)
 
            select case (mask%mode)
              !-----------------------------------------
              ! Mask Method
              !----------------------------------------
            case (pes_mask_mode_mask)
 
              cf1%zRs = (m_one - mask%cM%zRs) * cf1%zRs                               ! cf1 =(1-M)*\Psi_A(x,t2)
              call pes_mask_x_to_k(mask, cf1%zRs(:,:,:,1), cf2%Fs(:,:,:,1))     ! cf2 = \tilde{\Psi}_A(k,t2)
 
              if (mask%filter_k) then ! apply a filter to the Fourier transform to remove unwanted energies
                assert(allocated(mask%Mk))
                cf2%Fs(:,:,:,1) = cf2%Fs(:,:,:,1) * mask%Mk
              end if
 
              cf1%Fs(:,:,:,1) = mask%k(:,:,:, idim, ist, ik)                            ! cf1 = \Psi_B(k,t1)
              mask%k(:,:,:, idim, ist, ik) =  cf2%Fs(:,:,:,1)                           ! mask%k = \tilde{\Psi}_A(k,t2)
              ! cf1 = \tilde{\Psi}_B(k,t2)
              call pes_mask_volkov_time_evolution_wf(mask, space, mesh, kpoints, dt,iter-1, cf1%Fs(:,:,:,1), &
                st%d%get_kpoint_index(ik))
 
              ! mask%k = \tilde{\Psi}_A(k,t2) + \tilde{\Psi}_B(k,t2)
              mask%k(:,:,:, idim, ist, ik) =  mask%k(:,:,:, idim, ist, ik) + cf1%Fs(:,:,:,1)
 
              if (mask%back_action .eqv. .true.) then
                ! Apply Back-action to wavefunction in A
                call pes_mask_k_to_x(mask, cf1%Fs(:,:,:,1), cf2%zRs(:,:,:,1))      ! cf2 = \Psi_B(x,t2)
                call pes_mask_cube_to_mesh(mask, cf2, mf)
                call lalg_axpy(mask%mesh%np, m_one, mf, psi)
                call states_elec_set_state(st, mask%mesh, idim, ist, ik, psi)
 
                ! Apply correction to wavefunction in B
                cf2%zRs= (mask%cM%zRs) * cf2%zRs                                     ! cf2 = M*\Psi_B(x,t1)
                call pes_mask_x_to_k(mask, cf2%zRs(:,:,:,1), cf1%Fs(:,:,:,1))
 
                mask%k(:,:,:, idim, ist, ik) = mask%k(:,:,:, idim, ist, ik) - cf1%Fs(:,:,:,1)
              end if
 
              !-----------------------------------------
              ! Passive Mask method
              !----------------------------------------
            case (pes_mask_mode_passive)
              cf1%zRs = (m_one-mask%cM%zRs) * cf1%zRs
              call pes_mask_x_to_k(mask, cf1%zRs(:,:,:,1), cf2%Fs(:,:,:,1))
              mask%k(:,:,:, idim, ist, ik) = cf2%Fs(:,:,:,1)
 
            case default
              !Program should die before coming here
              write(message(1),'(a)') "PhotoElectroSpectrum = pes_mask. Unrecognized calculation mode."
              call messages_fatal(1, namespace=namespace)
 
            end select
 
          end do
        end do
      end do
 
      safe_deallocate_a(psi)
 
      call zcube_function_free_rs(mask%cube, cf1)
      call cube_function_free_fs(mask%cube, cf1)
      call zcube_function_free_rs(mask%cube, cf2)
      call cube_function_free_fs(mask%cube, cf2)
 
      select case (mask%mode)
      case (pes_mask_mode_mask)
        if (mask%back_action .eqv. .true.) then
          safe_deallocate_a(mf)
        end if
      end select
 
    end if ! time > mask%start_time
 
    if (mask%mode == pes_mask_mode_mask) call pes_mask_apply_mask(mask, st)  !apply the mask to all the KS orbitals
 
    pop_sub(pes_mask_calc)
 
    call profiling_out("PESMASK_calc")
  end subroutine pes_mask_calc
 
#include "pes_mask_out_inc.F90"
 
end module pes_mask_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: