target.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_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use ions_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lasers_oct_m
  use mesh_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use opt_control_global_oct_m
  use opt_control_state_oct_m
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use species_oct_m
  use states_elec_oct_m
  use target_classical_oct_m
  use target_density_oct_m
  use target_excited_oct_m
  use target_exclude_oct_m
  use target_ground_state_oct_m
  use target_gstransformation_oct_m
  use target_hhg_oct_m
  use target_spin_oct_m
  use target_local_oct_m
  use target_tdlocal_oct_m
  use target_userdefined_oct_m
  use target_velocity_oct_m
  use target_low_oct_m
  use td_calc_oct_m
  use td_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                  &
    target_t,                &
    target_get_state,        &
    target_init,             &
    target_end,              &
    target_output,           &
    target_tdcalc,           &
    target_inh,              &
    target_mode,             &
    target_type,             &
    target_j1,               &
    target_chi,              &
    target_move_ions,        &
    target_curr_functional,  &
    target_init_propagation
 
contains
 
  ! ---------------------------------------------------------
  subroutine target_init_propagation(tg)
    type(target_t), intent(inout)    :: tg
    push_sub(target_init_propagation)
 
    select case (tg%type)
    case (oct_tg_hhgnew)
      tg%vel = m_z0
      tg%gvec = m_z0
      tg%acc = m_z0
    end select
 
    pop_sub(target_init_propagation)
  end subroutine target_init_propagation
 
 
  ! ----------------------------------------------------------------------
  subroutine target_get_state(tg, st)
    type(target_t),      intent(in)    :: tg
    type(states_elec_t), intent(inout) :: st
 
    push_sub(target_get_state)
    call states_elec_copy(st, tg%st)
 
    pop_sub(target_get_state)
  end subroutine target_get_state
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_init(gr, kpoints, namespace, space, ions, qcs, td, w0, tg, oct, ep, mc)
    type(grid_t),                intent(in)    :: gr
    type(kpoints_t),             intent(in)    :: kpoints
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    type(ions_t),                intent(in)    :: ions
    type(opt_control_state_t),   intent(inout) :: qcs
    type(td_t),                  intent(in)    :: td
    real(real64),                intent(in)    :: w0
    type(target_t),              intent(inout) :: tg
    type(oct_t),                 intent(in)    :: oct
    type(epot_t),                intent(inout) :: ep
    type(multicomm_t),           intent(in)    :: mc
 
    integer :: ierr
    type(states_elec_t), pointer :: stin
    type(restart_t) :: restart
 
    push_sub(target_init)
 
    stin => opt_control_point_qs(qcs)
 
    !%Variable OCTTargetOperator
    !%Type integer
    !%Section Calculation Modes::Optimal Control
    !%Default oct_tg_gstransformation
    !%Description
    !% The variable <tt>OCTTargetOperator</tt> prescribes which kind of target functional is
    !% to be used.
    !%Option oct_tg_groundstate 1
    !% The target operator is a projection operator on the ground state, <i>i.e.</i> the
    !% objective is to populate the ground state as much as possible.
    !%Option oct_tg_excited 2
    !% (Experimental) The target operator is an "excited state". This means that the target operator
    !% is a linear combination of Slater determinants, each one formed by replacing
    !% in the ground-state Slater determinant one occupied state with one excited
    !% state (<i>i.e.</i> "single excitations"). The description of which excitations are
    !% used, and with which weights, should be given in a file called
    !% <tt>oct-excited-state-target</tt>.
    !% See the documentation of subroutine <tt>excited_states_elec_init</tt> in the source
    !% code in order to use this feature.
    !%Option oct_tg_gstransformation 3
    !% The target operator is a projection operator on a transformation of the ground-state
    !% orbitals defined by the block <tt>OCTTargetTransformStates</tt>.
    !%Option oct_tg_userdefined 4
    !% (Experimental) Allows to define target state by using <tt>OCTTargetUserdefined</tt>.
    !%Option oct_tg_jdensity 5
    !% (Experimental)
    !%Option oct_tg_local 6
    !% (Experimental) The target operator is a local operator.
    !%Option oct_tg_td_local 7
    !% (Experimental) The target operator is a time-dependent local operator.
    !%Option oct_tg_exclude_state 8
    !% (Experimental) Target operator is the projection onto the complement of a given state, given by the
    !% block <tt>OCTTargetTransformStates</tt>. This means that the target operator is the unity
    !% operator minus the projector onto that state.
    !%Option oct_tg_hhg 9
    !% (Experimental) The target is the optimization of the HHG yield. You must supply the <tt>OCTOptimizeHarmonicSpectrum</tt>
    !% block, and it attempts to optimize the maximum of the spectrum around each harmonic peak. You may
    !% use only one of the gradient-less optimization schemes.
    !%Option oct_tg_velocity 10
    !% (Experimental) The target is a function of the velocities of the nuclei at the end of the influence of
    !% the external field, defined by <tt>OCTVelocityTarget</tt>
    !%Option oct_tg_hhgnew 12
    !% (Experimental) The target is the optimization of the HHG yield. You must supply the
    !% <tt>OCTHarmonicWeight</tt> string. It attempts to optimize the integral of the harmonic spectrum multiplied
    !% by some user-defined weight function.
    !%Option oct_tg_classical 13
    !% (Experimental)
    !%Option oct_tg_spin 14
    !% (Experimental)
    !%End
    call parse_variable(namespace, 'OCTTargetOperator', oct_tg_gstransformation, tg%type)
    if (tg%type == oct_tg_excited) call messages_experimental('OCTTargetOperator = oct_tg_excited', namespace=namespace)
    if (tg%type == oct_tg_userdefined) call messages_experimental('OCTTargetOperator = oct_tg_userdefined')
    if (tg%type == oct_tg_jdensity) call messages_experimental('OCTTargetOperator = oct_tg_jdensity', namespace=namespace)
    if (tg%type == oct_tg_local) call messages_experimental('OCTTargetOperator = oct_tg_local', namespace=namespace)
    if (tg%type == oct_tg_td_local) call messages_experimental('OCTTargetOperator = oct_tg_td_local', namespace=namespace)
    if (tg%type == oct_tg_exclude_state) call messages_experimental('OCTTargetOperator = oct_tg_exclude_state', namespace=namespace)
    if (tg%type == oct_tg_hhg) call messages_experimental('OCTTargetOperator = oct_tg_hhg', namespace=namespace)
    if (tg%type == oct_tg_velocity) call messages_experimental('OCTTargetOperator = oct_tg_velocity', namespace=namespace)
    if (tg%type == oct_tg_hhgnew) call messages_experimental('OCTTargetOperator = oct_tg_hhgnew', namespace=namespace)
    if (tg%type == oct_tg_classical) call messages_experimental('OCTTargetOperator = oct_tg_classical', namespace=namespace)
    if (tg%type == oct_tg_spin) call messages_experimental('OCTTargetOperator = oct_tg_spin', namespace=namespace)
 
 
    if (.not. varinfo_valid_option('OCTTargetOperator', tg%type)) then
      call messages_input_error(namespace, 'OCTTargetOperator')
    end if
 
    call states_elec_copy(tg%st, stin)
    call states_elec_deallocate_wfns(tg%st)
    call states_elec_allocate_wfns(tg%st, gr, type_cmplx)
    call restart_init(restart, namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
    if (ierr /= 0) then
      message(1) = "Could not read gs for OCTTargetOperator."
      call messages_fatal(1, namespace=namespace)
    end if
 
    tg%curr_functional = oct_no_curr
 
    select case (tg%type)
    case (oct_tg_groundstate)
      call target_init_groundstate(gr, namespace, space, tg, td, restart, kpoints)
    case (oct_tg_excited)
      call messages_experimental('OCTTargetOperator = oct_tg_excited', namespace=namespace)
      call target_init_excited(gr, namespace, space, tg, td, restart, kpoints)
    case (oct_tg_exclude_state)
      call target_init_exclude(gr, namespace, space, tg, td, restart, kpoints)
    case (oct_tg_gstransformation)
      call target_init_gstransformation(gr, namespace, space, tg, td, restart, kpoints)
    case (oct_tg_userdefined)
      call target_init_userdefined(gr, namespace, tg, td)
    case (oct_tg_jdensity)
      call target_init_density(gr, kpoints, namespace, space, tg, stin, td, restart)
    case (oct_tg_local)
      call target_init_local(gr, namespace, tg, td)
    case (oct_tg_td_local)
      call target_init_tdlocal(gr, namespace, tg, td)
    case (oct_tg_hhg)
      call target_init_hhg(tg, namespace, td, w0)
    case (oct_tg_hhgnew)
      call messages_experimental('OCTTargetOperator = oct_tg_hhgnew', namespace=namespace)
      call target_init_hhgnew(gr, namespace, tg, td, ions, ep)
    case (oct_tg_velocity)
      call target_init_velocity(gr, namespace, ions, tg, oct, td, ep)
    case (oct_tg_classical)
      call messages_experimental('OCTTargetOperator = oct_tg_classical', namespace=namespace)
      call target_init_classical(ions, namespace, tg, td, oct)
    case (oct_tg_spin)
      call messages_experimental('OCTTargetOperator = oct_tg_spin', namespace=namespace)
      call target_init_spin(tg, namespace)
    case default
      write(message(1),'(a)') "Target Operator not properly defined."
      call messages_fatal(1, namespace=namespace)
    end select
 
    call restart_end(restart)
 
    nullify(stin)
    pop_sub(target_init)
  end subroutine target_init
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_end(tg, oct)
    type(target_t), intent(inout) :: tg
    type(oct_t), intent(in)       :: oct
 
    push_sub(target_end)
 
    call states_elec_end(tg%st)
 
    select case (tg%type)
    case (oct_tg_exclude_state)
      call target_end_exclude()
    case (oct_tg_jdensity)
      call target_end_density(tg)
    case (oct_tg_local)
      call target_end_local(tg)
    case (oct_tg_td_local)
      call target_end_tdlocal(tg)
    case (oct_tg_hhg)
      call target_end_hhg(tg)
    case (oct_tg_hhgnew)
      call target_end_hhgnew(tg, oct)
    case (oct_tg_velocity)
      call target_end_velocity(tg, oct)
    case (oct_tg_classical)
      call target_end_classical(tg)
    end select
 
    pop_sub(target_end)
  end subroutine target_end
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_output(tg, namespace, space, gr, dir, ions, hm, outp)
    type(target_t),           intent(inout) :: 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)
 
    select case (tg%type)
    case (oct_tg_groundstate)
      call target_output_groundstate(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_excited)
      call target_output_excited(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_exclude_state)
      call target_output_exclude(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_gstransformation)
      call target_output_gstransformation(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_userdefined)
      call target_output_userdefined(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_jdensity)
      call target_output_density(tg, namespace, space, gr, dir, ions, outp)
    case (oct_tg_local)
      call target_output_local(tg, namespace, space, gr, dir, ions, outp)
    case (oct_tg_td_local)
      call target_output_tdlocal(tg, namespace, space, gr, dir, ions, outp)
    case (oct_tg_hhg)
      call target_output_hhg(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_hhgnew)
      call target_output_hhg(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_velocity)
      call target_output_velocity(tg, namespace, space, gr, dir, ions, hm, outp)
    case (oct_tg_classical)
      call target_output_classical()
    end select
 
    pop_sub(target_output)
  end subroutine target_output
  ! ----------------------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine target_tdcalc(tg, namespace, space, hm, gr, ions, ext_partners, psi, time, max_time)
    type(target_t),           intent(inout) :: tg
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(hamiltonian_elec_t), intent(inout) :: 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(inout) :: psi
    integer,                  intent(in)    :: time
    integer,                  intent(in)    :: max_time
 
    if (target_mode(tg)  /= oct_targetmode_td) return
 
    push_sub(target_tdcalc)
 
    tg%td_fitness(time) = m_zero
 
    select case (tg%type)
    case (oct_tg_hhgnew)
      call target_tdcalc_hhgnew(tg, gr, psi, time, max_time)
    case (oct_tg_velocity)
      call target_tdcalc_velocity(tg, hm, gr, ions, psi, time, max_time)
    case (oct_tg_td_local)
      call target_tdcalc_tdlocal(tg, gr, psi, time)
    case (oct_tg_hhg)
      call target_tdcalc_hhg(tg, namespace, space, hm, gr, ions, ext_partners, psi, time)
    case (oct_tg_jdensity)
      call target_tdcalc_density(tg, gr, hm%kpoints, psi, time)
    case default
      message(1) = 'Error in target.target_tdcalc: default.'
      call messages_fatal(1, namespace=namespace)
    end select
 
    pop_sub(target_tdcalc)
  end subroutine target_tdcalc
  ! ----------------------------------------------------------------------
 
 
 
  ! ---------------------------------------------------------------
  subroutine target_inh(psi, gr, kpoints, tg, time, inh, iter)
    type(states_elec_t), intent(inout)     :: psi
    type(grid_t),        intent(in)        :: gr
    type(kpoints_t),     intent(in)        :: kpoints
    type(target_t),      intent(inout)     :: tg
    real(real64),        intent(in)        :: time
    type(states_elec_t), intent(inout)     :: inh
    integer,             intent(in)        :: iter
 
    integer :: ik, ist, ip, idim, ib
    complex(real64), allocatable :: zpsi(:)
    complex(real64) :: gvec(gr%box%dim)
 
    push_sub(target_inh)
 
    safe_allocate(zpsi(1:gr%np))
 
    select case (tg%type)
    case (oct_tg_td_local)
 
      call target_build_tdlocal(tg, gr, time)
 
      do ik = inh%d%kpt%start, inh%d%kpt%end
        do ist = inh%st_start, inh%st_end
          do idim = 1, inh%d%dim
            call states_elec_get_state(psi, gr, idim, ist, ik, zpsi)
            zpsi(1:gr%np) = -psi%occ(ist, ik)*tg%rho(1:gr%np)*zpsi(1:gr%np)
            call states_elec_set_state(inh, gr, idim, ist, ik, zpsi)
          end do
        end do
      end do
 
    case (oct_tg_hhgnew)
      gvec(:) = real(tg%gvec(iter + 1, :), real64)
 
      do ik = inh%d%kpt%start, inh%d%kpt%end
        do ist = inh%st_start, inh%st_end
          do idim = 1, inh%d%dim
            call states_elec_get_state(psi, gr, idim, ist, ik, zpsi)
            do ip = 1, gr%np
              zpsi(ip) = -psi%occ(ist, ik)*m_two*sum(tg%grad_local_pot(1, ip, 1:gr%box%dim)*gvec(:))*zpsi(ip)
            end do
            call states_elec_set_state(inh, gr, idim, ist, ik, zpsi)
          end do
        end do
      end do
 
    case (oct_tg_velocity)
 
      do ik = inh%d%kpt%start, inh%d%kpt%end
        do ist = inh%st_start, inh%st_end
          do idim = 1, inh%d%dim
            call states_elec_get_state(psi, gr, idim, ist, ik, zpsi)
            do ip = 1, gr%np
              zpsi(ip) = -psi%occ(ist, ik)*tg%rho(ip)*zpsi(ip)
            end do
            call states_elec_set_state(inh, gr, idim, ist, ik, zpsi)
          end do
        end do
      end do
 
    case (oct_tg_jdensity)
 
      do ik = inh%d%kpt%start, inh%d%kpt%end
        do ib = inh%group%block_start, inh%group%block_end
          call batch_set_zero(inh%group%psib(ib, ik))
        end do
      end do
 
      if (abs(nint(time/tg%dt)) >= tg%strt_iter_curr_tg) then
        call chi_current(tg, gr, kpoints, -1.0_real64, psi, inh)
      end if
 
    case default
      write(message(1),'(a)') 'Internal error in target_inh'
      call messages_fatal(1)
 
    end select
 
    safe_deallocate_a(zpsi)
 
    pop_sub(target_inh)
  end subroutine target_inh
  !----------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  real(real64) function target_j1(tg, namespace, gr, kpoints, qcpsi, ions) result(j1)
    type(target_t),             intent(inout)   :: tg
    type(namespace_t),          intent(in)      :: namespace
    type(grid_t),               intent(in)      :: gr
    type(kpoints_t),            intent(in)      :: kpoints
    type(opt_control_state_t),  intent(inout)   :: qcpsi
    type(ions_t),     optional, intent(in)      :: ions
 
    type(states_elec_t), pointer :: psi
 
    psi => opt_control_point_qs(qcpsi)
 
    push_sub(target_j1)
 
    j1 = m_zero
    select case (tg%type)
    case (oct_tg_groundstate)
      j1 = target_j1_groundstate(tg, gr, psi)
    case (oct_tg_excited)
      j1 = target_j1_excited(tg, namespace, gr, psi)
    case (oct_tg_gstransformation)
      j1 = target_j1_gstransformation(tg, gr, psi)
    case (oct_tg_userdefined)
      j1 = target_j1_userdefined(tg, gr, psi)
    case (oct_tg_jdensity)
      j1 = target_j1_density(gr, kpoints, tg, psi)
    case (oct_tg_local)
      j1 = target_j1_local(gr, tg, psi)
    case (oct_tg_td_local)
      j1 = target_j1_tdlocal(tg)
    case (oct_tg_exclude_state)
      j1 = target_j1_exclude(gr, tg, psi)
    case (oct_tg_hhg)
      j1 = target_j1_hhg(tg, namespace)
    case (oct_tg_hhgnew)
      j1 = target_j1_hhgnew(gr, tg)
    case (oct_tg_velocity)
      j1 = target_j1_velocity(tg, ions)
    case (oct_tg_classical)
      j1 = target_j1_classical(tg, qcpsi)
    case (oct_tg_spin)
      j1 = target_j1_spin(tg, gr, psi)
    end select
 
    nullify(psi)
    pop_sub(target_j1)
  end function target_j1
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine target_chi(tg, namespace, gr, kpoints, qcpsi_in, qcchi_out, ions)
    type(target_t),                    intent(inout) :: tg
    type(namespace_t),                 intent(in)    :: namespace
    type(grid_t),                      intent(in)    :: gr
    type(kpoints_t),                   intent(in)    :: kpoints
    type(opt_control_state_t), target, intent(inout) :: qcpsi_in
    type(opt_control_state_t), target, intent(inout) :: qcchi_out
    type(ions_t),                      intent(in)    :: ions
 
    real(real64), pointer :: q(:, :), p(:, :)
    type(states_elec_t), pointer :: psi_in, chi_out
    push_sub(target_chi)
 
    psi_in => opt_control_point_qs(qcpsi_in)
    chi_out => opt_control_point_qs(qcchi_out)
 
    select case (tg%type)
    case (oct_tg_groundstate)
 
      call target_chi_groundstate(tg, gr, psi_in, chi_out)
    case (oct_tg_excited)
      call target_chi_excited(tg, namespace, gr, psi_in, chi_out)
    case (oct_tg_gstransformation)
      call target_chi_gstransformation(tg, gr, psi_in, chi_out)
    case (oct_tg_userdefined)
      call target_chi_userdefined(tg, gr, psi_in, chi_out)
    case (oct_tg_jdensity)
      call target_chi_density(tg, gr, kpoints, psi_in, chi_out)
    case (oct_tg_local)
      call target_chi_local(tg, gr, psi_in, chi_out)
    case (oct_tg_td_local)
      call target_chi_tdlocal(chi_out)
    case (oct_tg_exclude_state)
      call target_chi_exclude(tg, gr, psi_in, chi_out)
    case (oct_tg_hhg)
      call target_chi_hhg(chi_out)
    case (oct_tg_hhgnew)
      call target_chi_hhg(chi_out)
    case (oct_tg_velocity)
      call target_chi_velocity(gr, tg, chi_out, ions)
    case (oct_tg_classical)
      call target_chi_classical(tg, qcpsi_in, qcchi_out, ions)
    case (oct_tg_spin)
      call target_chi_spin(tg, gr, psi_in, chi_out)
    end select
 
    ! Unless the target is "classical", the co-state classical variables are zero at time t=T.
    if (tg%type .ne. oct_tg_classical) then
      q => opt_control_point_q(qcchi_out)
      p => opt_control_point_p(qcchi_out)
      q = m_zero
      p = m_zero
      nullify(q)
      nullify(p)
    end if
 
    nullify(psi_in)
    nullify(chi_out)
    pop_sub(target_chi)
  end subroutine target_chi
 
 
end module target_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: