target_classical.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_classical_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 ion_dynamics_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use messages_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 states_elec_oct_m
  use string_oct_m
  use target_oct_m
  use td_oct_m
 
  implicit none
 
  private
  public :: target_classical_t
 
  type, extends(target_t) :: target_classical_t
    private
    character(len=4096)            :: classical_input_string
    character(len=1024), allocatable :: mom_der_array(:,:)
    character(len=1024), allocatable :: pos_der_array(:,:)
  contains
    procedure :: init      => target_init_classical
    procedure :: cleanup   => target_end_classical
    procedure :: j1        => target_j1_classical
    procedure :: apply_chi => target_chi_classical
    procedure :: output    => target_output_classical
  end type target_classical_t
 
 
contains
 
  ! ----------------------------------------------------------------------
  subroutine target_init_classical(tg, gr, kpoints, namespace, space, ions, qcs, td, w0, oct, ep, restart)
    class(target_classical_t), intent(inout) :: tg
    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(oct_t),               intent(in)    :: oct
    type(epot_t),              intent(inout) :: ep
    type(restart_t),           intent(inout) :: restart
 
    integer             :: jj, ist, jst
    type(block_t)       :: blk
    character(len=1024) :: expression
    push_sub(target_init_classical)
 
    tg%move_ions = td%ions_dyn%ions_move()
    tg%dt = td%dt
    assert(tg%move_ions)
 
    !%Variable OCTClassicalTarget
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% If <tt>OCTTargetOperator = oct_tg_classical</tt>, the you must supply this block.
    !% It should contain a string (e.g. "(q[1,1]-q[1,2])*p[2,1]") with a mathematical
    !% expression in terms of two arrays, q and p, that represent the position and momenta
    !% of the classical variables. The first index runs through the various classical particles,
    !% and the second index runs through the spatial dimensions.
    !%
    !% In principle, the block only contains one entry (string). However, if the expression is very
    !% long, you can split it into various lines (one column each) that will be concatenated.
    !%
    !% The QOCT algorithm will attempt to maximize this expression, at the end of the propagation.
    !%End
    if (parse_block(namespace, 'OCTClassicalTarget', blk) == 0) then
      tg%classical_input_string = " "
      do jj = 0, parse_block_n(blk) - 1
        call parse_block_string(blk, jj, 0, expression)
        tg%classical_input_string = trim(tg%classical_input_string) // trim(expression)
      end do
      call parse_block_end(blk)
    else
      message(1) = 'If OCTTargetOperator = oct_tg_classical, then you must give the shape'
      message(2) = 'of this target in the block "OCTClassicalTarget".'
      call messages_fatal(2, namespace=namespace)
    end if
 
    !%Variable OCTMomentumDerivatives
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% This block should contain the derivatives of the expression given in
    !% <tt>OCTClassicalTarget</tt> with respect to the p array components.
    !% Each line corresponds to a different classical particle, whereas the
    !% columns correspond to each spatial dimension: the (i,j) block component
    !% corresponds with the derivative wrt p[i,j].
    !%End
    if (parse_block(namespace, 'OCTMomentumDerivatives', blk) == 0) then
      safe_allocate(tg%mom_der_array(1:ions%natoms,1:ions%space%dim))
      do ist = 0, ions%natoms - 1
        do jst = 0, ions%space%dim - 1
          call parse_block_string(blk, ist, jst, tg%mom_der_array(ist+1, jst+1))
        end do
      end do
      call parse_block_end(blk)
    else if (oct%algorithm == option__octscheme__oct_cg .or. oct%algorithm == option__octscheme__oct_bfgs) then
      message(1) = 'If "OCTTargetOperator = oct_classical" and "OCTScheme = oct_cg" or'
      message(2) = '"OCTScheme = oct_bfgs", then you must define the blocks "OCTClassicalTarget",'
      message(3) = '"OCTPositionDerivatives" AND "OCTMomentumDerivatives"'
      call messages_fatal(3, namespace=namespace)
    end if
 
    !%Variable OCTPositionDerivatives
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% This block should contain the derivatives of the expression given in
    !% <tt>OCTClassicalTarget</tt> with respect to the q array components.
    !% Each line corresponds to a different classical particle, whereas the
    !% columns correspond to each spatial dimension: the (i,j) block component
    !% corresponds with the derivative wrt q[i,j].
    !%End
    if (parse_block(namespace, 'OCTPositionDerivatives', blk) == 0) then
      safe_allocate(tg%pos_der_array(1:ions%natoms,1:ions%space%dim))
      do ist = 0, ions%natoms-1
        do jst = 0, ions%space%dim-1
          call parse_block_string(blk, ist, jst, tg%pos_der_array(ist+1, jst+1))
        end do
      end do
      call parse_block_end(blk)
    else if (oct%algorithm == option__octscheme__oct_cg .or. oct%algorithm == option__octscheme__oct_bfgs) then
      message(1) = 'If "OCTTargetOperator = oct_tg_classical" and "OCTScheme = oct_cg" or'
      message(2) = '"OCTScheme = oct_bfgs", then you must define the blocks "OCTClassicalTarget",'
      message(3) = '"OCTPositionDerivatives" AND "OCTMomentumDerivatives"'
      call messages_fatal(3, namespace=namespace)
    end if
 
    pop_sub(target_init_classical)
  end subroutine target_init_classical
  ! ----------------------------------------------------------------------
 
  ! ----------------------------------------------------------------------
  subroutine target_end_classical(tg, oct)
    class(target_classical_t), intent(inout) :: tg
    type(oct_t),               intent(in)    :: oct
 
    push_sub(target_end_classical)
 
    safe_deallocate_a(tg%pos_der_array)
    safe_deallocate_a(tg%mom_der_array)
 
    pop_sub(target_end_classical)
  end subroutine target_end_classical
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_output_classical(tg, namespace, space, gr, dir, ions, hm, outp)
    class(target_classical_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_classical)
    pop_sub(target_output_classical)
  end subroutine target_output_classical
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  real(real64) function target_j1_classical(tg, namespace, gr, kpoints, qcpsi, ions) result(j1)
    class(target_classical_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
 
    real(real64), pointer :: q(:, :), p(:, :)
    real(real64) :: dummy(3)
    character(len=4096) :: inp_string
 
    push_sub(target_j1_classical)
 
    q => opt_control_point_q(qcpsi)
    p => opt_control_point_p(qcpsi)
 
    inp_string = tg%classical_input_string
    call parse_array(inp_string, q, 'q')
    call parse_array(inp_string, p, 'p')
    call conv_to_c_string(inp_string)
    call parse_expression(j1, dummy(1), 1, dummy(1:3), dummy(1), dummy(1), inp_string)
 
    nullify(q)
    nullify(p)
    pop_sub(target_j1_classical)
  end function target_j1_classical
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  subroutine target_chi_classical(tg, namespace, gr, kpoints, qcpsi_in, qcchi_out, ions)
    class(target_classical_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
 
    integer :: ist, jst, ib, iqn
    character(len=1024) :: temp_string
    real(real64) :: df_dv, dummy(3)
    real(real64), pointer :: q(:, :), p(:, :), tq(:, :), tp(:, :)
    type(states_elec_t), pointer :: chi
    push_sub(target_chi_classical)
 
    tq => opt_control_point_q(qcchi_out)
    tp => opt_control_point_p(qcchi_out)
    q => opt_control_point_q(qcpsi_in)
    p => opt_control_point_p(qcpsi_in)
    tq = m_zero
    tp = m_zero
 
    do ist = 1, ions%natoms
      do jst=1, ions%space%dim
        temp_string = tg%mom_der_array(ist, jst)
        call parse_array(temp_string, p, 'p')
        call parse_array(temp_string, q, 'q')
        call conv_to_c_string(temp_string)
        call parse_expression(df_dv, dummy(1), 1, dummy(1:3), dummy(1), dummy(1), temp_string)
        tq(ist, jst) = -df_dv
      end do
    end do
 
    do ist = 1, ions%natoms
      do jst=1, ions%space%dim
        temp_string = tg%pos_der_array(ist, jst)
        call parse_array(temp_string, p, 'p')
        call parse_array(temp_string, q, 'q')
        call conv_to_c_string(temp_string)
        call parse_expression(df_dv, dummy(1), 1, dummy(1:3), dummy(1), dummy(1), temp_string)
        tp(ist, jst) = df_dv
      end do
    end do
 
    chi => opt_control_point_qs(qcchi_out)
 
    !We assume that there is no time-independent operator.
    do iqn = chi%d%kpt%start, chi%d%kpt%end
      do ib = chi%group%block_start, chi%group%block_end
        call batch_set_zero(chi%group%psib(ib, iqn))
      end do
    end do
 
    nullify(tq)
    nullify(tp)
    nullify(q)
    nullify(p)
    nullify(chi)
    pop_sub(target_chi_classical)
  end subroutine target_chi_classical
  ! ----------------------------------------------------------------------
 
end module target_classical_oct_m
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: