lasers.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 N. Tancogne-Dejean
!!
!! 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 lasers_oct_m
  use clock_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use lalg_basic_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lattice_vectors_oct_m
  use lorentz_force_oct_m
  use math_oct_m
  use mpi_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use space_oct_m
  use symmetries_oct_m
  use symm_op_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use tdfunction_oct_m
 
  implicit none
 
  private
  public ::                       &
    lasers_t,                     &
    lasers_parse_external_fields, &
    lasers_generate_potentials,   &
    load_lasers,                  &
    lasers_check_symmetries,      &
    laser_t,                      &
    laser_write_info,             &
    laser_field,                  &
    lasers_nondipole_laser_field_step,   &
    lasers_set_nondipole_parameters,     &
    lasers_with_nondipole_field,   &
    laser_electric_field,         &
    laser_potential,              &
    laser_vector_potential,       &
    laser_to_numerical,           &
    laser_to_numerical_all,       &
    laser_kind,                   &
    laser_polarization,           &
    laser_get_f,                  &
    laser_set_f,                  &
    laser_get_phi,                &
    laser_set_phi,                &
    laser_set_empty_phi,          &
    laser_set_f_value,            &
    laser_carrier_frequency,      &
    laser_set_frequency,          &
    laser_set_polarization
 
 
  ! TODO: (Micael, Alex) Issue 836. Remove the following paramaters and use the
  ! corresponding quantities defined in the multisystem (E_FIELD, B_FIELD, etc)
  integer, public, parameter ::     &
    E_FIELD_NONE             =  0,  &
    e_field_electric         =  1,  &
    e_field_magnetic         =  2,  &
    e_field_vector_potential =  3,  &
    e_field_scalar_potential =  4
 
  type laser_t
    private
    integer :: field      = e_field_none  
    complex(real64), allocatable :: pol(:)
    real(real64), allocatable :: prop(:)
    type(tdf_t) :: f
    type(tdf_t) :: phi
    real(real64) :: omega        = m_zero        
 
    real(real64), allocatable :: v(:)
    real(real64), allocatable :: a(:, :)
    character(len=200) :: scalar_pot_expression
    character(len=200) :: phase_expression
    character(len=200) :: envelope_expression
  end type laser_t
 
  type, extends(interaction_partner_t) :: lasers_t
    private
 
    integer, public                    :: no_lasers
    type(laser_t), allocatable, public :: lasers(:)
 
    real(real64), allocatable :: e(:)
    real(real64), allocatable :: b(:)
    real(real64), allocatable :: integrated_nondipole_afield(:)
    real(real64) :: nd_integration_time
    real(real64) :: nd_integration_step
  contains
    procedure :: init_interaction_as_partner => lasers_init_interaction_as_partner
    procedure :: update_quantity => lasers_update_quantity
    procedure :: copy_quantities_to_interaction => lasers_copy_quantities_to_interaction
    final :: lasers_finalize
  end type
 
  interface lasers_t
    module procedure lasers_constructor
  end interface lasers_t
 
 
contains
 
  function lasers_constructor(namespace) result(this)
    class(lasers_t), pointer :: this
    type(namespace_t), intent(in) :: namespace
 
    push_sub(lasers_constructor)
 
    safe_allocate(this)
 
    this%namespace = namespace_t("Lasers", parent=namespace)
 
    safe_allocate(this%e(1:3))
    safe_allocate(this%b(1:3))
    this%e = m_zero
    this%b = m_zero
 
    pop_sub(lasers_constructor)
  end function lasers_constructor
 
  subroutine lasers_parse_external_fields(this)
    class(lasers_t),      intent(inout) :: this
 
    type(block_t)     :: blk
    integer           :: il, jj, ierr, k
    real(real64) :: omega0
    complex(real64) :: cprop(3)
 
    push_sub(lasers_parse_external_fields)
 
    call messages_obsolete_variable(this%namespace, "TDLasers", "TDExternalFields")
    !%Variable TDExternalFields
    !%Type block
    !%Section Time-Dependent
    !%Description
    !% The block <tt>TDExternalFields</tt> describes the type and shape of time-dependent
    !% external perturbations that are applied to the system, in the form
    !% <math>f(x,y,z) \cos(\omega t + \phi (t)) g(t)</math>, where <math>f(x,y,z)</math> is defined by
    !% by a field type and polarization or a scalar potential, as below; <math>\omega</math>
    !% is defined by <tt>omega</tt>; <math>g(t)</math> is defined by
    !% <tt>envelope_function_name</tt>; and <math>\phi(t)</math> is the (time-dependent) phase from <tt>phase</tt>.
    !%
    !% These perturbations are only applied for time-dependent runs. If
    !% you want the value of the perturbation at time zero to be
    !% applied for time-independent runs, use <tt>TimeZero = yes</tt>.
    !%
    !% Each line of the block describes an external field; this way you can actually have more
    !% than one laser (<i>e.g.</i> a "pump" and a "probe").
    !%
    !% There are two ways to specify <math>f(x,y,z)</math> but both use the same <tt>omega | envelope_function_name [| phase]</tt>
    !% for the time-dependence.
    !% The float <tt>omega</tt> will be the carrier frequency of the
    !% pulse (in energy units). The envelope of the field is a time-dependent function whose definition
    !% must be given in a <tt>TDFunctions</tt> block. <tt>envelope_function_name</tt> is a string (and therefore
    !% it must be surrounded by quotation marks) that must match one of the function names
    !% given in the first column of the <tt>TDFunctions</tt> block.
    !% <tt>phase</tt> is optional and is taken to be zero if not provided, and is also a string specifying
    !% a time-dependent function.
    !%
    !% (A) type = <tt>electric field, magnetic field, vector_potential</tt>
    !%
    !% For these cases, the syntax is:
    !%
    !% <tt>%TDExternalFields
    !% <br>&nbsp;&nbsp; type | nx | ny | nz | omega | envelope_function_name | phase (| px | py | pz )
    !% <br>%</tt>
    !%
    !% The <tt>vector_potential</tt> option (constant in space) permits us to describe
    !% an electric perturbation in the velocity gauge.
    !% The three (possibly complex) numbers (<tt>nx</tt>, <tt>ny</tt>, <tt>nz</tt>) mark the polarization
    !% direction of the field.
    !% By default, (<tt>nx</tt>, <tt>ny</tt>, <tt>nz</tt>) are defined in Cartesian space.
    !% However, it is possible for solids to define them using the Miller indices.
    !% This can be achieved by defining the block <tt>MillerIndicesBasis</tt>. It is also possible to activate
    !% the nondipole SFA correction by inserting the laser propagation direction px | py | pz as prescribed in
    !% Phys. Rev. A 101, 043408 (2020).
    !%
    !% (B) type = <tt>scalar_potential</tt>
    !%
    !% <tt>%TDExternalFields
    !% <br>&nbsp;&nbsp; scalar_potential | "spatial_expression" | omega | envelope_function_name | phase
    !% <br>%</tt>
    !%
    !% The scalar potential is any expression of the spatial coordinates given by the string
    !% "spatial_expression", allowing a field beyond the dipole approximation.
    !%
    !% For DFTB runs, only fields of type type = <tt>electric field</tt> are allowed for the moment, and the
    !% <tt>type</tt> keyword is omitted.
    !%
    !% A NOTE ON UNITS:
    !%
    !% It is very common to describe the strength of a laser field by its intensity, rather
    !% than using the electric-field amplitude. In atomic units (or, more precisely, in any
    !% Gaussian system of units), the relationship between instantaneous electric field
    !% and intensity is:
    !% <math> I(t) = \frac{c}{8\pi} E^2(t) </math>.
    !%
    !% It is common to read intensities in W/cm<math>^2</math>. The dimensions of intensities are
    !% [W]/(L<math>^2</math>T), where [W] are the dimensions of energy. The relevant conversion factors
    !% are:
    !%
    !% Hartree / (<math>a_0^2</math> atomic_time) = <math>6.4364086 \times 10^{15} \mathrm{W/cm}^2</math>
    !%
    !% eV / ( &Aring;<math>^2 (\hbar</math>/eV) ) = <math>2.4341348 \times 10^{12} \mathrm{W/cm}^2</math>
    !%
    !% If, in atomic units, we set the electric-field amplitude to <math>E_0</math>,
    !% then the intensity is:
    !%
    !% <math> I_0 = 3.51 \times 10^{16} \mathrm{W/cm}^2 (E_0^2) </math>
    !%
    !% If, working with <tt>Units = ev_angstrom</tt>, we set <math>E_0</math>, then the intensity is:
    !%
    !% <math> I_0 = 1.327 \times 10^{13} (E_0^2) \mathrm{W/cm}^2 </math>
    !%
    !%Option electric_field 1
    !% The external field is an electric field, the usual case when we want to describe a
    !% laser in the length gauge.
    !%Option magnetic_field 2
    !% The external field is a (homogeneous) time-dependent magnetic field.
    !%Option vector_potential 3
    !% The external field is a time-dependent homogeneous vector potential, which may describe
    !% a laser field in the velocity gauge.
    !%Option scalar_potential 4
    !% The external field is an arbitrary scalar potential, which may describe an
    !% inhomogeneous electrical field.
    !%End
 
    this%no_lasers = 0
    if (parse_block(this%namespace, 'TDExternalFields', blk) == 0) then
      this%no_lasers = parse_block_n(blk)
      safe_allocate(this%lasers(1:this%no_lasers))
 
      do il = 1, this%no_lasers
        safe_allocate(this%lasers(il)%pol(1:3))
        this%lasers(il)%pol = m_z0
 
        call parse_block_integer(blk, il-1, 0, this%lasers(il)%field)
 
        select case (this%lasers(il)%field)
        case (e_field_scalar_potential)
          call parse_block_string(blk, il-1, 1, this%lasers(il)%scalar_pot_expression)
          jj = 1
          this%lasers(il)%pol = m_z1
        case default
          call parse_block_cmplx(blk, il-1, 1, this%lasers(il)%pol(1))
          call parse_block_cmplx(blk, il-1, 2, this%lasers(il)%pol(2))
          call parse_block_cmplx(blk, il-1, 3, this%lasers(il)%pol(3))
          jj = 3
        end select
 
        call parse_block_float(blk, il-1, jj+1, omega0)
 
        this%lasers(il)%omega = omega0
 
        call parse_block_string(blk, il-1, jj+2, this%lasers(il)%envelope_expression)
        call tdf_read(this%lasers(il)%f, this%namespace, trim(this%lasers(il)%envelope_expression), ierr)
 
        ! Check if there is a phase.
        if (parse_block_cols(blk, il-1) > jj+3) then
          call parse_block_string(blk, il-1, jj+3, this%lasers(il)%phase_expression)
          call tdf_read(this%lasers(il)%phi, this%namespace, trim(this%lasers(il)%phase_expression), ierr)
          if (ierr /= 0) then
            write(message(1),'(3A)') 'Error in the "', trim(this%lasers(il)%envelope_expression), &
              '" field defined in the TDExternalFields block:'
            write(message(2),'(3A)') 'Time-dependent phase function "', trim(this%lasers(il)%phase_expression), &
              '" not found.'
            call messages_warning(2, namespace=this%namespace)
          end if
        else
          call tdf_init(this%lasers(il)%phi)
        end if
 
        ! Check if there is nondipoole fields to activate the nondipole SFA formalism of Phys. Rev. A 101, 043408 (2020)
        if (parse_block_cols(blk, il-1) > jj+4) then
          safe_allocate(this%lasers(il)%prop(1:size(cprop)))
          do k = 1, size(cprop)
            call parse_block_cmplx(blk, il-1, jj+size(cprop)+k, cprop(k))
 
          end do
          if (any(abs(aimag(cprop)) > m_epsilon)) then
            write(message(1),'(3A)') 'Error in the "', trim(this%lasers(il)%envelope_expression), &
              '" field defined in the TDExternalFields block:'
            write(message(2),'(A)') 'Propagation direction cannot be complex.'
            call messages_fatal(2, namespace=this%namespace)
          end if
          this%lasers(il)%prop(:) = real(cprop(:), real64)
          if (.not.allocated(this%integrated_nondipole_afield)) then
            safe_allocate(this%integrated_nondipole_afield(1:size(cprop)))
            this%integrated_nondipole_afield(1:3)=m_zero
            this%nd_integration_time=m_zero
            !%Variable NDSFATimeIntegrationStep
            !%Type float
            !%Default 0.01
            !%Section Hamiltonian
            !%Description
            !% Timestep for the integration of the NDSFA first order
            !% response numerically as prescribed in Phys. Rev. A 101, 043408 (2020).
            !% Typically needs to be an order of magniude less than
            !% the TDtimestep. However, for writing "laser" the TDtimestep
            !% is used instead.
            !%End
            call parse_variable(this%namespace, 'NDSFATimeIntegrationStep', 0.01_real64, this%nd_integration_step)
 
            if (is_close(this%nd_integration_step, 0.01_real64)) then
              message(1) = "The default timestep of 0.01 is utilized for the nondipole SFA integration."
              message(2) = "Be aware that this should be at least an order of magniude less than the TDtimestep"
              call messages_info(2, namespace=this%namespace)
            end if
          end if
        end if
 
      end do
 
      call parse_block_end(blk)
    end if
 
    pop_sub(lasers_parse_external_fields)
  end subroutine lasers_parse_external_fields
 
  subroutine lasers_generate_potentials(this, mesh, space, latt)
    class(lasers_t),      intent(inout) :: this
    class(mesh_t),           intent(in) :: mesh
    class(space_t),          intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
 
    type(block_t)     :: blk2
    integer           :: il, ip, idir, idir2
    real(real64) :: rr, pot_re, pot_im, xx(3)
    real(real64) :: miller(space%dim,space%dim), miller_red(space%dim,space%dim)
 
    push_sub(lasers_generate_potentials)
 
    do il = 1, this%no_lasers
      ! For periodic systems, we might want to define the polarization direction in
      ! terms of Miller indices
      ! In this case, the use must provide the coordinates of the X, Y, and Z high symmetry points
      ! such that the code applying a laser along say the [100] direction converts it to
      ! the proper coordinates in Cartesian space.
      ! This is usefull for arbitrarily rotated crystals.
 
      !%Variable MillerIndicesBasis
      !%Type block
      !%Section Time-Dependent
      !%Description
      !% When this block is given, the polarisation of the TDExternalFields is
      !% understood to be defined in terms of Miller indices.
      !% This block define the corresponding basis, by defining the reduced coordinates
      !% of the X, Y, and Z high symmetry points, such that the code can do the corresponding
      !% transformation.
      !%
      !% For example, in an FCC crystal with the conventional primitive cell,
      !% the following input allows to define the polarization in terms
      !% message(1) = "The lasers break (at least) one of the symmetries used to reduce the k-points  ."
      !%
      !% <tt>%MillerIndicesBasis
      !% <br> 0.0 | 0.5 | 0.5
      !% <br> 0.5 | 0.0 | 0.5
      !% <br> 0.5 | 0.5 | 0.0
      !% <br>%</tt>
      !%
      !% Indeed, in this case, the reciprocal lattice vectors are (-1, 1, 1), (1, -1, 1),
      !% and (1, 1, -1) in units of 2*pi/a.
      !% This directly gives that the [100] direction correspond to the x direction, [111]
      !% gives the vector (1,1,1), etc.
      !%
      !%End
      if (parse_block(this%namespace, 'MillerIndicesBasis', blk2) == 0) then
        if(.not. space%is_periodic()) then
          write(message(1),'(a)') 'MillerIndicesBasis can only be used for periodic systems.'
          call messages_fatal(1, namespace=this%namespace)
        end if
 
        do idir = 1, space%dim
          do idir2 = 1, space%dim
            call parse_block_float(blk2, idir-1, idir2-1, miller_red(idir2, idir))
          end do
          call kpoints_to_absolute(latt, miller_red(:, idir), miller(:, idir))
        end do
 
 
        this%lasers(il)%pol = matmul(miller, this%lasers(il)%pol)
        call parse_block_end(blk2)
      end if
 
      this%lasers(il)%pol(:) = this%lasers(il)%pol(:)/sqrt(sum(abs(this%lasers(il)%pol(:))**2))
 
      select case (this%lasers(il)%field)
      case (e_field_scalar_potential)
        safe_allocate(this%lasers(il)%v(1:mesh%np_part))
        this%lasers(il)%v = m_zero
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, coords = xx(1:space%dim))
          xx(space%dim+1:3) = m_zero
          call parse_expression(pot_re, pot_im, 3, xx, rr, m_zero, trim(this%lasers(il)%scalar_pot_expression))
          this%lasers(il)%v(ip) = pot_re
        end do
 
      case (e_field_magnetic)
        ! \warning: note that for the moment we are ignoring the possibility of a complex
        ! polarizability vector for the td magnetic-field case.
        safe_allocate(this%lasers(il)%a(1:mesh%np_part, 1:3))
        this%lasers(il)%a = m_zero
        do ip = 1, mesh%np
          xx(1:space%dim) = mesh%x(:, ip)
          xx(space%dim+1:3) = m_zero
          ! Compute the vector potential from a uniform B field.
          ! The sign is determined by the relation $\vec{B} = \nabla \times \vec{A}$.
          ! This leads to  $\vec{A} = -\frac{1}{2}\vec{r}\times\vec{B}$.
          select case (space%dim)
          case (2)
            this%lasers(il)%a(ip, 1:2) = (/xx(2), -xx(1)/) * sign(m_one, real(this%lasers(il)%pol(3)))
          case (3)
            this%lasers(il)%a(ip, :) = (/ xx(2)*real(this%lasers(il)%pol(3)) - xx(3)*real(this%lasers(il)%pol(2)), &
              xx(3)*real(this%lasers(il)%pol(1)) - xx(1)*real(this%lasers(il)%pol(3)), &
              xx(1)*real(this%lasers(il)%pol(2)) - xx(2)*real(this%lasers(il)%pol(1))  /)
          case default
            message(1) = "Magnetic fields only allowed in 2 or 3D."
            call messages_fatal(1, namespace=this%namespace)
          end select
        end do
        this%lasers(il)%a = -m_half * this%lasers(il)%a
 
      end select
 
    end do
 
    pop_sub(lasers_generate_potentials)
  end subroutine lasers_generate_potentials
 
 
  ! ---------------------------------------------------------
  subroutine lasers_check_symmetries(this, kpoints)
    type(lasers_t),  intent(in) :: this
    type(kpoints_t), intent(in) :: kpoints
 
    integer :: iop, il
 
    push_sub(lasers_check_symmetries)
 
    if (kpoints%use_symmetries) then
      do iop = 1, symmetries_number(kpoints%symm)
        if (iop == symmetries_identity_index(kpoints%symm)) cycle
        do il = 1, this%no_lasers
          if (.not. symm_op_invariant_cart(kpoints%symm%ops(iop), this%lasers(il)%pol(:), symprec)) then
            message(1) = "The lasers break (at least) one of the symmetries used to reduce the k-points  ."
            message(2) = "Set SymmetryBreakDir accordingly to your laser fields."
            call messages_fatal(2, namespace=this%namespace)
          end if
        end do
      end do
    end if
 
    pop_sub(lasers_check_symmetries)
  end subroutine lasers_check_symmetries
 
  ! ---------------------------------------------------------
  subroutine lasers_finalize(this)
    type(lasers_t), intent(inout) :: this
 
    push_sub(lasers_finalize)
 
    call lasers_deallocate(this)
 
    pop_sub(lasers_finalize)
  end subroutine lasers_finalize
 
  ! ---------------------------------------------------------
  subroutine lasers_deallocate(this)
    class(lasers_t), intent(inout) :: this
 
    integer :: il
 
    push_sub(lasers_deallocate)
 
    safe_deallocate_a(this%e)
    safe_deallocate_a(this%b)
    safe_deallocate_a(this%integrated_nondipole_afield)
 
    do il = 1, this%no_lasers
      safe_deallocate_a(this%lasers(il)%pol)
      call tdf_end(this%lasers(il)%f)
      call tdf_end(this%lasers(il)%phi)
      select case (this%lasers(il)%field)
      case (e_field_scalar_potential)
        safe_deallocate_a(this%lasers(il)%v)
      case (e_field_magnetic)
        safe_deallocate_a(this%lasers(il)%a)
      end select
      safe_deallocate_a(this%lasers(il)%prop)
    end do
    safe_deallocate_a(this%lasers)
 
    pop_sub(lasers_deallocate)
  end subroutine lasers_deallocate
 
 
  ! ---------------------------------------------------------
  real(real64) function laser_carrier_frequency(laser) result(w0)
    type(laser_t), intent(in) :: laser
 
    push_sub(laser_carrier_frequency)
 
    w0 = laser%omega
 
    pop_sub(laser_carrier_frequency)
  end function laser_carrier_frequency
 
  ! ---------------------------------------------------------
  subroutine lasers_init_interaction_as_partner(partner, interaction)
    class(lasers_t),      intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(lasers_init_interaction_as_partner)
 
    select type (interaction)
    type is (lorentz_force_t)
      ! Nothing to be initialized for the Lorentz force.
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(lasers_init_interaction_as_partner)
  end subroutine lasers_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine lasers_update_quantity(this, label)
    class(lasers_t),  intent(inout) :: this
    character(len=*), intent(in)    :: label
 
    type(quantity_t), pointer :: quantity
    integer :: il
 
    push_sub(lasers_update_quantity)
 
    if (allocated(this%lasers)) then
      if (any(laser_kind(this%lasers) == e_field_vector_potential) .or. &
        any(laser_kind(this%lasers) == e_field_scalar_potential)) then
        call messages_not_implemented("Laser vector potentials and scalar potentials in multi-system framework", &
          namespace=this%namespace)
      end if
    end if
 
    quantity => this%quantities%get(label)
    select case (label)
    case ("E field")
      this%e = m_zero
      do il = 1, this%no_lasers
        if (laser_kind(this%lasers(il)) == e_field_electric) then
          call laser_field(this%lasers(il), this%e, quantity%iteration%value())
        end if
      end do
 
    case ("B field")
      this%b = m_zero
      do il = 1, this%no_lasers
        if (laser_kind(this%lasers(il)) == e_field_magnetic) then
          call laser_field(this%lasers(il), this%b, quantity%iteration%value())
        end if
      end do
 
    case default
      message(1) = "Incompatible quantity."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(lasers_update_quantity)
  end subroutine lasers_update_quantity
 
  ! ---------------------------------------------------------
  subroutine lasers_copy_quantities_to_interaction(partner, interaction)
    class(lasers_t),      intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    integer :: ip
 
    push_sub(lasers_copy_quantities_to_interaction)
 
    select type (interaction)
    type is (lorentz_force_t)
      do ip = 1, interaction%system_np
        interaction%partner_e_field(:, ip) = partner%e
        interaction%partner_b_field(:, ip) = partner%b
      end do
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(lasers_copy_quantities_to_interaction)
  end subroutine lasers_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  integer pure elemental function laser_kind(laser)
    type(laser_t), intent(in) :: laser
 
    ! no push_sub allowed in pure function
    laser_kind = laser%field
 
  end function laser_kind
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  function laser_polarization(laser) result(pol)
    type(laser_t), intent(in) :: laser
    complex(real64) :: pol(3)
 
    push_sub(laser_polarization)
 
    pol = laser%pol
 
    pop_sub(laser_polarization)
  end function laser_polarization
  ! ---------------------------------------------------------
 
  function lasers_with_nondipole_field(lasers) result(isnondipole)
    type(lasers_t),           intent(in)    :: lasers
    logical :: isnondipole
 
    push_sub(lasers_with_nondipole_field)
    isnondipole = .false.
    if(allocated(lasers%integrated_nondipole_afield)) isnondipole = .true.
    pop_sub(lasers_with_nondipole_field)
  end function lasers_with_nondipole_field
  ! ---------------------------------------------------------
 
 
  subroutine lasers_set_nondipole_parameters(this, ndfield, nd_integration_time)
    type(lasers_t), intent(inout) :: this
    real(real64),          intent(in)    :: ndfield(:)
    real(real64),          intent(in)    :: nd_integration_time
    integer :: dim
    dim = size(ndfield)
 
    push_sub(lasers_set_nondipole_parameters)
 
    this%integrated_nondipole_afield(1:dim) = ndfield(1:dim)
    this%nd_integration_time = nd_integration_time
 
    pop_sub(lasers_set_nondipole_parameters)
  end subroutine lasers_set_nondipole_parameters
 
 
  ! ---------------------------------------------------------
  subroutine laser_get_f(laser, ff)
    type(laser_t), intent(in)    :: laser
    type(tdf_t),   intent(inout) :: ff
 
    push_sub(laser_get_f)
    call tdf_copy(ff, laser%f)
 
    pop_sub(laser_get_f)
  end subroutine laser_get_f
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_set_f(laser, ff)
    type(laser_t), intent(inout) :: laser
    type(tdf_t),   intent(inout) :: ff
 
    push_sub(laser_set_f)
 
    call tdf_end(laser%f)
    call tdf_copy(laser%f, ff)
 
    pop_sub(laser_set_f)
  end subroutine laser_set_f
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_get_phi(laser, phi)
    type(laser_t), intent(in)    :: laser
    type(tdf_t),   intent(inout) :: phi
 
    push_sub(laser_get_phi)
    call tdf_copy(phi, laser%phi)
 
    pop_sub(laser_get_phi)
  end subroutine laser_get_phi
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_set_phi(laser, phi)
    type(laser_t), intent(inout) :: laser
    type(tdf_t),   intent(inout) :: phi
 
    push_sub(laser_set_phi)
 
    call tdf_end(laser%phi)
    call tdf_copy(laser%phi, phi)
 
    pop_sub(laser_set_phi)
  end subroutine laser_set_phi
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_set_empty_phi(laser)
    type(laser_t),          intent(inout) :: laser
 
    push_sub(laser_set_empty_phi)
 
    call tdf_init(laser%phi)
 
    pop_sub(laser_set_empty_phi)
  end subroutine laser_set_empty_phi
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine laser_set_f_value(laser, ii, xx)
    type(laser_t), intent(inout) :: laser
    integer,       intent(in)    :: ii
    real(real64),  intent(in)    :: xx
 
    push_sub(laser_set_f_value)
    call tdf_set_numerical(laser%f, ii, xx)
 
    pop_sub(laser_set_f_value)
  end subroutine laser_set_f_value
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_set_frequency(laser, omega)
    type(laser_t), intent(inout) :: laser
    real(real64),  intent(in)    :: omega
 
    push_sub(laser_set_frequency)
    laser%omega = omega
 
    pop_sub(laser_set_frequency)
  end subroutine laser_set_frequency
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_set_polarization(laser, pol)
    type(laser_t), intent(inout) :: laser
    complex(real64),         intent(in)    :: pol(:)
 
    push_sub(laser_set_polarization)
 
    laser%pol = pol
 
    pop_sub(laser_set_polarization)
  end subroutine laser_set_polarization
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine laser_to_numerical_all(laser, dt, max_iter, omegamax)
    type(laser_t),   intent(inout)  :: laser
    real(real64),    intent(in)     :: dt
    integer,         intent(in)     :: max_iter
    real(real64),    intent(in)     :: omegamax
 
    integer :: iter
    real(real64)   :: tt, fj, phi
 
    push_sub(lasers_to_numerical_all)
 
    call tdf_to_numerical(laser%f, max_iter, dt, omegamax)
    do iter = 1, max_iter + 1
      tt = (iter-1)*dt
      fj = tdf(laser%f, iter)
      phi = tdf(laser%phi, tt)
      call tdf_set_numerical(laser%f, iter, fj*cos(laser%omega*tt+phi))
    end do
    call tdf_end(laser%phi)
    call tdf_init_cw(laser%phi, m_zero, m_zero)
    laser%omega = m_zero
 
    pop_sub(lasers_to_numerical_all)
  end subroutine laser_to_numerical_all
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine laser_to_numerical(laser, dt, max_iter, omegamax)
    type(laser_t),   intent(inout) :: laser
    real(real64),    intent(in)    :: dt
    integer,         intent(in)    :: max_iter
    real(real64),    intent(in)    :: omegamax
 
    push_sub(lasers_to_numerical)
 
    call tdf_to_numerical(laser%f, max_iter, dt, omegamax)
    call tdf_to_numerical(laser%phi,  max_iter, dt, omegamax)
 
    pop_sub(lasers_to_numerical)
  end subroutine laser_to_numerical
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine laser_write_info(lasers, namespace, dt, max_iter, iunit)
    type(laser_t),               intent(in) :: lasers(:)
    type(namespace_t),           intent(in) :: namespace
    real(real64),      optional, intent(in) :: dt
    integer,           optional, intent(in) :: max_iter
    integer,           optional, intent(in) :: iunit
 
    real(real64) :: tt, fluence, max_intensity, intensity, dt_, field(3), up, maxfield,tmp
    integer :: il, iter, no_l, max_iter_
 
    push_sub(laser_write_info)
 
    no_l = size(lasers)
 
    do il = 1, no_l
 
      if (present(dt)) then
        dt_ = dt
      else
        dt_ = tdf_dt(lasers(il)%f)
      end if
      if (present(max_iter)) then
        max_iter_ = max_iter
      else
        max_iter_ = tdf_niter(lasers(il)%f)
      end if
 
      write(message(1),'(i2,a)') il, ':'
      select case (lasers(il)%field)
      case (e_field_electric)
        message(2) = '   Electric Field.'
      case (e_field_magnetic)
        message(2) = '   Magnetic Field.'
      case (e_field_vector_potential)
        message(2) = '   Vector Potential.'
      case (e_field_scalar_potential)
        message(2) = '   Scalar Potential.'
      end select
      call messages_info(2, iunit=iunit, namespace=namespace)
 
      if (lasers(il)%field /= e_field_scalar_potential) then
        write(message(1),'(3x,a,3(a1,f7.4,a1,f7.4,a1))') 'Polarization: ', &
          '(', real(lasers(il)%pol(1), real64), ',', aimag(lasers(il)%pol(1)), '), ', &
          '(', real(lasers(il)%pol(2), real64), ',', aimag(lasers(il)%pol(2)), '), ', &
          '(', real(lasers(il)%pol(3), real64), ',', aimag(lasers(il)%pol(3)), ')'
        call messages_info(1, iunit=iunit, namespace=namespace)
      end if
 
      write(message(1),'(3x,a,f14.8,3a)') 'Carrier frequency = ', &
        units_from_atomic(units_out%energy, lasers(il)%omega), &
        ' [', trim(units_abbrev(units_out%energy)), ']'
      message(2) = '   Envelope: '
      call messages_info(2, iunit=iunit, namespace=namespace)
      call tdf_write(lasers(il)%f, iunit)
 
      if (.not. tdf_is_empty(lasers(il)%phi)) then
        message(1) = '   Phase: '
        call messages_info(1, iunit=iunit, namespace=namespace)
        call tdf_write(lasers(il)%phi, iunit)
      end if
 
      ! 1 atomic unit of intensity = 3.5094448e+16 W / cm^2
      ! In a Gaussian system of units,
      ! I(t) = (1/(8\pi)) * c * E(t)^2
      ! (1/(8\pi)) * c = 5.4525289841210 a.u.
      if (lasers(il)%field == e_field_electric .or. lasers(il)%field == e_field_vector_potential) then
        fluence = m_zero
        max_intensity = m_zero
        maxfield= m_zero
        do iter = 1, max_iter_
          tt = iter * dt_
          call laser_electric_field(lasers(il), field, tt, dt_)
          intensity = 5.4525289841210_real64*sum(field**2)
          fluence = fluence + intensity
          if (intensity > max_intensity) max_intensity = intensity
 
          tmp = sum(field(:)**2)
          if (tmp > maxfield) maxfield = tmp
        end do
        fluence = fluence * dt_
 
        write(message(1),'(a,es17.6,3a)') '   Peak intensity       = ', max_intensity, ' [a.u]'
        write(message(2),'(a,es17.6,3a)') '                        = ', &
          max_intensity * 6.4364086e+15_real64, ' [W/cm^2]'
        write(message(3),'(a,es17.6,a)')  '   Int. intensity       = ', fluence, ' [a.u]'
        write(message(4),'(a,es17.6,a)')  '   Fluence              = ', &
          fluence / 5.4525289841210_real64 , ' [a.u]'
        call messages_info(4, iunit=iunit, namespace=namespace)
 
        if (abs(lasers(il)%omega) > m_epsilon) then
          ! Ponderomotive Energy is the cycle-averaged kinetic energy of
          ! a free electron quivering in the field
          ! Up = E^2/(4*\omega^2)
          !
          ! subroutine laser_to_numerical_all sets lasers%omega to zero
          !
          up = maxfield/(4*lasers(il)%omega**2)
 
          write(message(1),'(a,es17.6,3a)') '   Ponderomotive energy = ', &
            units_from_atomic(units_out%energy, up) ,&
            ' [', trim(units_abbrev(units_out%energy)), ']'
          call messages_info(1, iunit=iunit, namespace=namespace)
        end if
      end if
 
    end do
 
    pop_sub(laser_write_info)
  end subroutine laser_write_info
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_potential(laser, mesh, pot, time)
    type(laser_t),     intent(in)    :: laser
    class(mesh_t),     intent(in)    :: mesh
    real(real64),      intent(inout) :: pot(:)
    real(real64), optional,   intent(in)    :: time
 
    complex(real64) :: amp
    integer :: ip
    real(real64) :: field(3)
 
    push_sub(laser_potential)
 
    if (present(time)) then
      amp = tdf(laser%f, time) * exp(m_zi * (laser%omega * time + tdf(laser%phi, time)))
    else
      amp = m_z1
    end if
 
    select case (laser%field)
    case (e_field_scalar_potential)
      call lalg_axpy(mesh%np, real(amp, real64), laser%v, pot)
    case default
      field(:) = real(amp * laser%pol(:), real64)
      do ip = 1, mesh%np
        ! The -1 sign is missing here. Check epot.F90 for the explanation.
        pot(ip) = pot(ip) + sum(field(1:mesh%box%dim) * mesh%x(1:mesh%box%dim, ip))
      end do
    end select
 
    pop_sub(laser_potential)
  end subroutine laser_potential
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_vector_potential(laser, mesh, aa, time)
    type(laser_t),   intent(in)    :: laser
    type(mesh_t),    intent(in)    :: mesh
    real(real64),    intent(inout) :: aa(:, :)
    real(real64), optional, intent(in)    :: time
 
    real(real64)   :: amp
    integer :: ip, idir
 
    push_sub(laser_vector_potential)
 
    if (present(time)) then
      amp = tdf(laser%f, time)*cos((laser%omega*time + tdf(laser%phi, time)))
      do idir = 1, mesh%box%dim
        do ip = 1, mesh%np
          aa(ip, idir) = aa(ip, idir) + amp*laser%a(ip, idir)
        end do
      end do
    else
      do idir = 1, mesh%box%dim
        do ip = 1, mesh%np
          aa(ip, idir) = aa(ip, idir) + laser%a(ip, idir)
        end do
      end do
    end if
 
    pop_sub(laser_vector_potential)
  end subroutine laser_vector_potential
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine laser_field(laser, field, time)
    type(laser_t),     intent(in)    :: laser
    real(real64),      intent(inout) :: field(:)
    real(real64), optional,   intent(in)    :: time
 
    integer :: dim
    complex(real64) :: amp
 
    !no PUSH SUB, called too often
 
    dim = size(field)
 
    if (present(time)) then
      amp = tdf(laser%f, time) * exp(m_zi * (laser%omega * time + tdf(laser%phi, time)))
    else
      amp = m_z1
    end if
    if (laser%field == e_field_scalar_potential) then
      ! In this case we will just return the value of the time function. The "field", in fact,
      ! should be a function of the position in space (thus, a true "field"), given by the
      ! gradient of the scalar potential.
      field(1) = field(1) + real(amp, real64)
    else
      field(1:dim) = field(1:dim) + real(amp*laser%pol(1:dim), real64)
    end if
  end subroutine laser_field
 
! ---------------------------------------------------------
  subroutine lasers_nondipole_laser_field_step(this, field, time)
    type(lasers_t),     intent(in)    :: this
    real(real64),             intent(out) :: field(:)
    real(real64),             intent(in)    :: time
 
    real(real64) :: a0(3)
    real(real64) :: e0(3)
    integer :: ilaser
    integer :: jlaser
    integer :: dim
    integer :: iter
    dim = size(field)
 
    field(1:dim) = this%integrated_nondipole_afield(1:dim)
    if (time - this%nd_integration_time < 0) then
      return
    endif
    do iter = 1, nint((time-this%nd_integration_time)/this%nd_integration_step)
      do ilaser = 1, this%no_lasers
        if(laser_kind(this%lasers(ilaser)) /= e_field_vector_potential) cycle
        do jlaser = 1, this%no_lasers
          if(laser_kind(this%lasers(jlaser)) /= e_field_vector_potential) cycle
          if(allocated(this%lasers(jlaser)%prop)) then
            a0 = m_zero
            call laser_field(this%lasers(ilaser), a0, &
              iter*this%nd_integration_step+this%nd_integration_time)
            call laser_electric_field(this%lasers(jlaser), e0, &
              iter*this%nd_integration_step+this%nd_integration_time, this%nd_integration_step)
            ! We have front factor q/(mc) with q=-abs(e)  =-1 a.u. and m = 1 a.u.
            field(1:dim) = field(1:dim) - m_one/(p_c) * this%nd_integration_step &
              * dot_product(a0,e0) * this%lasers(jlaser)%prop(1:dim)
          end if
        end do
      end do
    end do
  end subroutine lasers_nondipole_laser_field_step
 
 
  ! ---------------------------------------------------------
  subroutine laser_electric_field(laser, field, time, dt)
    type(laser_t),     intent(in)    :: laser
    real(real64),      intent(out)   :: field(:)
    real(real64),      intent(in)    :: time
    real(real64),      intent(in)    :: dt
 
    integer :: dim
    real(real64), allocatable :: field1(:), field2(:)
 
    !no PUSH SUB, called too often
 
    dim = size(field)
 
    select case (laser%field)
    case (e_field_electric)
      field = m_zero
      call laser_field(laser, field(1:dim), time)
    case (e_field_vector_potential)
      safe_allocate(field1(1:dim))
      safe_allocate(field2(1:dim))
      field1 = m_zero
      field2 = m_zero
      call laser_field(laser, field1(1:dim), time - dt)
      call laser_field(laser, field2(1:dim), time + dt)
      field = - (field2 - field1) / (m_two * p_c * dt)
      safe_deallocate_a(field1)
      safe_deallocate_a(field2)
    case default
      field = m_zero
    end select
 
  end subroutine laser_electric_field
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  ! Loading of the lasers for the multisystem framework
  ! Ultimately this should be done in laser_init
  subroutine load_lasers(partners, namespace)
    class(partner_list_t), intent(inout)  :: partners
    type(namespace_t),    intent(in)     :: namespace
 
    class(lasers_t), pointer :: lasers
    integer :: il
 
    push_sub(load_lasers)
 
    lasers => lasers_t(namespace)
 
    call lasers_parse_external_fields(lasers)
 
    ! TODO: see how to do this in the multisystem framework
    if (parse_is_defined(namespace, 'MillerIndicesBasis')) then
      call messages_not_implemented("MillerIndicesBasis with load_lasers routine")
    end if
 
    ! This is done normally in lasers_generate_potential, but we cannot call this routine here,
    ! so we do it here
    do il = 1, lasers%no_lasers
      lasers%lasers(il)%pol(:) = lasers%lasers(il)%pol(:)/sqrt(sum(abs(lasers%lasers(il)%pol(:))**2))
    end do
 
    call lasers%quantities%add(quantity_t("E field", always_available = .true., updated_on_demand = .true., iteration = clock_t()))
    call lasers%quantities%add(quantity_t("B field", always_available = .true., updated_on_demand = .true., iteration = clock_t()))
 
    lasers%supported_interactions_as_partner = [lorentz_force]
 
    if (lasers%no_lasers > 0) then
      call partners%add(lasers)
    else
      safe_deallocate_p(lasers)
    end if
 
    pop_sub(load_lasers)
  end subroutine load_lasers
 
end module lasers_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: