target_hhg.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 target_hhg_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use epot_oct_m
  use fft_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use ions_oct_m
  use ion_dynamics_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use opt_control_global_oct_m
  use output_oct_m
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use spectrum_oct_m
  use states_elec_oct_m
  use target_low_oct_m
  use td_oct_m
  use td_calc_oct_m
 
 
  implicit none
 
  private
  public ::               &
    target_init_hhg,      &
    target_init_hhgnew,   &
    target_end_hhg,       &
    target_end_hhgnew,    &
    target_j1_hhg,        &
    target_j1_hhgnew,     &
    target_tdcalc_hhgnew, &
    target_tdcalc_hhg,    &
    target_output_hhg,    &
    target_chi_hhg
 
 
contains
 
 
 
  ! ----------------------------------------------------------------------
  subroutine target_init_hhg(tg, namespace, td, w0)
    type(target_t),    intent(inout) :: tg
    type(namespace_t), intent(in)    :: namespace
    type(td_t),        intent(in)    :: td
    real(real64),      intent(in)    :: w0
 
    integer       :: jj
    type(block_t) :: blk
    push_sub(target_init_hhg)
 
    tg%move_ions = ion_dynamics_ions_move(td%ions_dyn)
 
    !%Variable OCTOptimizeHarmonicSpectrum
    !%Type block
    !%Default no
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% (Experimental)
    !% If <tt>OCTTargetOperator = oct_tg_hhg</tt>, the target is the harmonic emission spectrum.
    !% In that case, you must supply an <tt>OCTOptimizeHarmonicSpectrum</tt> block in the <tt>inp</tt>
    !% file. The target is given, in general, by:
    !%
    !% <math>J_1 = \int_0^\infty d\omega \alpha(\omega) H(\omega)</math>,
    !%
    !% where <math>H(\omega)</math> is the harmonic spectrum generated by the system, and
    !% <math>\alpha(\omega)</math> is some function that determines what exactly we want
    !% to optimize. The role of the <tt>OCTOptimizeHarmonicSpectrum</tt> block is to determine
    !% this <math>\alpha(\omega)</math> function. Currently, this function is defined as:
    !%
    !% <math>\alpha(\omega) = \sum_{L=1}^{M} \frac{\alpha_L}{a_L} \sqcap( (\omega - L\omega_0)/a_L )</math>,
    !%
    !% where <math>\omega_0</math> is the carrier frequency. <math>M</math> is
    !% the number of columns in the <tt>OCTOptimizeHarmonicSpectrum</tt> block. The values of <i>L</i> will be listed
    !% in the first row of this block; <math>\alpha_L</math> in the second row, and <math>a_L</math> in
    !% the third.
    !%
    !% Example:
    !%
    !% <tt>%OCTOptimizeHarmonicSpectrum
    !% <br>&nbsp;&nbsp;  7    |  9    | 11
    !% <br>&nbsp;&nbsp; -1    |  1    | -1
    !% <br>&nbsp;&nbsp;  0.01 |  0.01 |  0.01
    !% <br>%</tt>
    !%
    !%End
    if (parse_is_defined(namespace, 'OCTOptimizeHarmonicSpectrum')) then
      if (parse_block(namespace, 'OCTOptimizeHarmonicSpectrum', blk) == 0) then
        tg%hhg_nks = parse_block_cols(blk, 0)
        safe_allocate(    tg%hhg_k(1:tg%hhg_nks))
        safe_allocate(tg%hhg_alpha(1:tg%hhg_nks))
        safe_allocate(    tg%hhg_a(1:tg%hhg_nks))
        do jj = 1, tg%hhg_nks
          call parse_block_integer(blk, 0, jj - 1, tg%hhg_k(jj))
          call parse_block_float(blk, 1, jj - 1, tg%hhg_alpha(jj))
          call parse_block_float(blk, 2, jj - 1, tg%hhg_a(jj))
        end do
        call parse_block_end(blk)
      else
        message(1) = '"OCTOptimizeHarmonicSpectrum" has to be specified as a block.'
        call messages_info(1, namespace=namespace)
        call messages_input_error(namespace, 'OCTOptimizeHarmonicSpectrum')
      end if
    else
      write(message(1), '(a)') 'If "OCTTargetMode = oct_targetmode_hhg", you must supply an'
      write(message(2), '(a)') '"OCTOptimizeHarmonicSpectrum" block.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    tg%hhg_w0 = w0
    tg%dt     = td%dt
    safe_allocate(tg%td_fitness(0:td%max_iter))
    tg%td_fitness = m_zero
 
    pop_sub(target_init_hhg)
  end subroutine target_init_hhg
 
 
  ! ----------------------------------------------------------------------
  subroutine target_init_hhgnew(gr, namespace, tg, td, ions, ep)
    type(grid_t),      intent(in)    :: gr
    type(namespace_t), intent(in)    :: namespace
    type(target_t),    intent(inout) :: tg
    type(td_t),        intent(in)    :: td
    type(ions_t),      intent(in)    :: ions
    type(epot_t),      intent(in)    :: ep
 
    integer :: ist, jst, iunit, jj, nn(3), optimize_parity(3)
    logical :: optimize(3)
    real(real64) :: dw, psi_re, psi_im, ww
    real(real64), allocatable  :: vl(:), vl_grad(:,:)
    push_sub(target_init_hhgnew)
 
    tg%move_ions = ion_dynamics_ions_move(td%ions_dyn)
 
    ! We allocate many things that are perhaps not necessary if we use a direct optimization scheme.
    safe_allocate(tg%vel(1:td%max_iter+1, 1:gr%box%dim))
    safe_allocate(tg%acc(1:td%max_iter+1, 1:gr%box%dim))
    safe_allocate(tg%gvec(1:td%max_iter+1, 1:gr%box%dim))
    safe_allocate(tg%alpha(1:td%max_iter))
 
    ! The following is a temporary hack, that assumes only one atom at the origin of coordinates.
    if (ions%natoms > 1) then
      message(1) = 'If "OCTTargetOperator = oct_tg_hhgnew", then you can only have one atom.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! The following is a temporary hack, that assumes only one atom at the origin of coordinates.
    if (ions%natoms > 1) then
      message(1) = 'If "OCTTargetOperator = oct_tg_hhgnew", then you can only have one atom.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(tg%grad_local_pot(1:ions%natoms, 1:gr%np, 1:gr%box%dim))
    safe_allocate(vl(1:gr%np_part))
    safe_allocate(vl_grad(1:gr%np, 1:gr%box%dim))
    safe_allocate(tg%rho(1:gr%np))
 
    vl(:) = m_zero
    vl_grad(:,:) = m_zero
    call epot_local_potential(ep, namespace, ions%space, ions%latt, gr, ions%atom(1)%species, &
      ions%pos(:, 1), 1, vl)
    call dderivatives_grad(gr%der, vl, vl_grad)
    do jst=1, gr%box%dim
      do ist = 1, gr%np
        tg%grad_local_pot(1, ist, jst) = vl_grad(ist, jst)
      end do
    end do
    ! Note that the calculation of the gradient of the potential
    ! is wrong at the borders of the box, since it assumes zero boundary
    ! conditions. The best way to solve this problems is to define the
    ! target making use of the definition of the forces based on the gradient
    ! of the density, rather than on the gradient of the potential.
 
 
    !%Variable OCTHarmonicWeight
    !%Type string
    !%Default "1"
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% (Experimental) If <tt>OCTTargetOperator = oct_tg_plateau</tt>, then the function to optimize is the integral of the
    !% harmonic spectrum <math>H(\omega)</math>, weighted with a function <math>f(\omega)</math>
    !% that is defined as a string here. For example, if
    !% you set <tt>OCTHarmonicWeight  = "step(w-1)"</tt>, the function to optimize is
    !% the integral of <math>step(\omega-1)*H(\omega)</math>, <i>i.e.</i>
    !% <math>\int_1^{\infty} H \left( \omega \right) d\omega</math>.
    !% In practice, it is better if you also set an upper limit, <i>e.g.</i>
    !% <math>f(\omega) = step(\omega-1) step(2-\omega)</math>.
    !%End
    call parse_variable(namespace, 'OCTHarmonicWeight', '1', tg%plateau_string)
    tg%dt = td%dt
    safe_allocate(tg%td_fitness(0:td%max_iter))
    tg%td_fitness = m_zero
 
    iunit = io_open('.alpha', namespace, action='write')
    dw = (m_two * m_pi) / (td%max_iter * tg%dt)
    do jj = 0, td%max_iter - 1
      ww = jj * dw
      call parse_expression(psi_re, psi_im, "w", ww, tg%plateau_string)
      tg%alpha(jj+1) = psi_re
      write(iunit, *) ww, psi_re
    end do
    call io_close(iunit)
 
    nn(1:3) = (/ td%max_iter, 1, 1 /)
    optimize(1:3) = .false.
    optimize_parity(1:3) = -1
    call fft_init(tg%fft_handler, nn(1:3), 1, fft_complex, fftlib_fftw, optimize, optimize_parity)
 
    pop_sub(target_init_hhgnew)
  end subroutine target_init_hhgnew
 
 
  ! ----------------------------------------------------------------------
  subroutine target_end_hhg(tg)
    type(target_t),   intent(inout) :: tg
 
    push_sub(target_end_hhg)
 
    safe_deallocate_a(tg%hhg_k)
    safe_deallocate_a(tg%hhg_alpha)
    safe_deallocate_a(tg%hhg_a)
    safe_deallocate_a(tg%td_fitness)
 
    pop_sub(target_end_hhg)
  end subroutine target_end_hhg
 
 
  ! ----------------------------------------------------------------------
  subroutine target_output_hhg(tg, namespace, space, gr, dir, ions, hm, outp)
    type(target_t),      intent(in) :: tg
    type(namespace_t),   intent(in) :: namespace
    class(space_t),      intent(in)    :: space
    type(grid_t),        intent(in) :: gr
    character(len=*),    intent(in) :: dir
    type(ions_t),        intent(in) :: ions
    type(hamiltonian_elec_t), intent(in) :: hm
    type(output_t),      intent(in) :: outp
 
    push_sub(target_output_hhg)
 
    call io_mkdir(trim(dir), namespace)
    call output_states(outp, namespace, space, trim(dir), tg%st, gr, ions, hm, -1)
 
    pop_sub(target_output_hhg)
  end subroutine target_output_hhg
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_end_hhgnew(tg, oct)
    type(target_t), intent(inout) :: tg
    type(oct_t),    intent(in)    :: oct
 
    push_sub(target_end_hhgnew)
 
    if ((oct%algorithm == option__octscheme__oct_cg) .or. (oct%algorithm == option__octscheme__oct_bfgs)) then
      safe_deallocate_a(tg%grad_local_pot)
      safe_deallocate_a(tg%rho)
      safe_deallocate_a(tg%vel)
      safe_deallocate_a(tg%acc)
      safe_deallocate_a(tg%gvec)
      safe_deallocate_a(tg%alpha)
      safe_deallocate_a(tg%td_fitness)
      call fft_end(tg%fft_handler)
    end if
 
    pop_sub(target_end_hhgnew)
  end subroutine target_end_hhgnew
 
 
  ! ----------------------------------------------------------------------
  real(real64) function target_j1_hhg(tg, namespace) result(j1)
    type(target_t),    intent(in) :: tg
    type(namespace_t), intent(in) :: namespace
 
    integer :: maxiter, jj
    real(real64) :: aa, ww, maxhh, omega
    complex(real64), allocatable :: ddipole(:)
    push_sub(target_j1_hhg)
 
    maxiter = size(tg%td_fitness) - 1
    safe_allocate(ddipole(0:maxiter))
    ddipole = m_z0
    ddipole = tg%td_fitness
 
    j1 = m_zero
 
    call spectrum_hsfunction_init(tg%dt, 0, maxiter, maxiter, ddipole)
    do jj = 1, tg%hhg_nks
      aa = tg%hhg_a(jj) * tg%hhg_w0
      ww = tg%hhg_k(jj) * tg%hhg_w0
      call spectrum_hsfunction_min(namespace, ww - aa, ww + aa, omega, maxhh)
      j1 = j1 + tg%hhg_alpha(jj) * log(-maxhh)
    end do
    call spectrum_hsfunction_end()
 
    safe_deallocate_a(ddipole)
    pop_sub(target_j1_hhg)
  end function target_j1_hhg
 
 
  ! ----------------------------------------------------------------------
  real(real64) function target_j1_hhgnew(gr, tg) result(j1)
    type(grid_t),   intent(in) :: gr
    type(target_t), intent(in) :: tg
 
    integer :: maxiter, i
    real(real64) :: dw, ww
    push_sub(target_j1_hhgnew)
 
    maxiter = size(tg%td_fitness) - 1
    dw = (m_two * m_pi) / (maxiter * tg%dt)
    j1 = m_zero
    do i = 0, maxiter - 1
      ww = i * dw
      j1 = j1 + dw * tg%alpha(i+1) * sum(abs(tg%vel(i+1, 1:gr%box%dim))**2)
    end do
 
    pop_sub(target_j1_hhgnew)
  end function target_j1_hhgnew
 
 
  ! ----------------------------------------------------------------------
  subroutine target_chi_hhg(chi_out)
    type(states_elec_t), intent(inout) :: chi_out
 
    integer :: ik, ib
    push_sub(target_chi_hhg)
 
    !we have a time-dependent target --> Chi(T)=0
    do ik = chi_out%d%kpt%start, chi_out%d%kpt%end
      do ib = chi_out%group%block_start, chi_out%group%block_end
        call batch_set_zero(chi_out%group%psib(ib, ik))
      end do
    end do
 
    pop_sub(target_chi_hhg)
  end subroutine target_chi_hhg
 
 
  ! ---------------------------------------------------------
  subroutine target_tdcalc_hhgnew(tg, gr, psi, time, max_time)
    type(target_t),      intent(inout) :: tg
    type(grid_t),        intent(in)    :: gr
    type(states_elec_t), intent(in)    :: psi
    integer,             intent(in)    :: time
    integer,             intent(in)    :: max_time
 
    complex(real64), allocatable :: opsi(:, :), zpsi(:, :)
    integer :: iw, ia, ist, idim, ik
    real(real64) :: acc(gr%box%dim), dt, dw
 
    push_sub(target_tdcalc_hhgnew)
 
    tg%td_fitness(time) = m_zero
 
    ! If the ions move, the tg is computed in the propagation routine.
    if (.not. target_move_ions(tg)) then
 
      safe_allocate(opsi(1:gr%np_part, 1))
      safe_allocate(zpsi(1:gr%np_part, 1))
 
      opsi = m_z0
      ! WARNING This does not work for spinors.
      ! The following is a temporary hack. It assumes only one atom at the origin.
      do ik = 1, psi%nik
        do ist = 1, psi%nst
          call states_elec_get_state(psi, gr, ist, ik, zpsi)
          do idim = 1, gr%box%dim
            opsi(1:gr%np, 1) = tg%grad_local_pot(1, 1:gr%np, idim)*zpsi(1:gr%np, 1)
            acc(idim) = acc(idim) + real( psi%occ(ist, ik) * zmf_dotp(gr, psi%d%dim, opsi, zpsi), real64)
            tg%acc(time+1, idim) = tg%acc(time+1, idim) + psi%occ(ist, ik)*zmf_dotp(gr, psi%d%dim, opsi, zpsi)
          end do
        end do
      end do
 
      safe_deallocate_a(opsi)
      safe_deallocate_a(zpsi)
 
    end if
 
    dt = tg%dt
    dw = (m_two * m_pi/(max_time * tg%dt))
    if (time == max_time) then
      tg%acc(1, 1:gr%box%dim) = m_half * (tg%acc(1, 1:gr%box%dim) + tg%acc(max_time+1, 1:gr%box%dim))
      do ia = 1, gr%box%dim
        call zfft_forward(tg%fft_handler, tg%acc(1:max_time, ia), tg%vel(1:max_time, ia))
      end do
      tg%vel = tg%vel * tg%dt
      do iw = 1, max_time
        ! We add the one-half dt term because when doing the propagation we want the value at interpolated times.
 
        tg%acc(iw, 1:gr%box%dim) = tg%vel(iw, 1:gr%box%dim) * tg%alpha(iw) * exp(m_zi * (iw-1) * dw * m_half * dt)
      end do
      do ia = 1, gr%box%dim
        call zfft_backward(tg%fft_handler, tg%acc(1:max_time, ia), tg%gvec(1:max_time, ia))
      end do
      tg%gvec(max_time + 1, 1:gr%box%dim) = tg%gvec(1, 1:gr%box%dim)
      tg%gvec = tg%gvec * (m_two * m_pi/ tg%dt)
    end if
 
    pop_sub(target_tdcalc_hhgnew)
  end subroutine target_tdcalc_hhgnew
  ! ----------------------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine target_tdcalc_hhg(tg, namespace, space, hm, gr, ions, ext_partners, psi, time)
    type(target_t),           intent(inout) :: tg
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(in)    :: psi
    integer,                  intent(in)    :: time
 
    real(real64) :: acc(space%dim)
 
    push_sub(target_tdcalc_hhg)
 
    call td_calc_tacc(namespace, space, gr, ions, ext_partners, psi, hm, acc, time*tg%dt)
    tg%td_fitness(time) = acc(1)
 
    pop_sub(target_tdcalc_hhg)
  end subroutine target_tdcalc_hhg
  ! ----------------------------------------------------------------------
 
end module target_hhg_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: