target_velocity.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_velocity_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use io_oct_m
  use ions_oct_m
  use ion_dynamics_oct_m
  use kpoints_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use opt_control_global_oct_m
  use opt_control_state_oct_m
  use output_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_velocity_t
 
  type, extends(target_t) :: target_velocity_t
    private
    character(len=4096)            :: vel_input_string
    character(len=1024), allocatable :: vel_der_array(:,:)
    real(real64), allocatable      :: grad_local_pot(:,:,:)
    real(real64), allocatable      :: rho(:)
    real(real64), allocatable      :: td_fitness(:)
  contains
    procedure :: init      => target_init_velocity
    procedure :: cleanup   => target_end_velocity
    procedure :: j1        => target_j1_velocity
    procedure :: apply_chi => target_chi_velocity
    procedure :: output    => target_output_velocity
    procedure :: tdcalc    => target_tdcalc_velocity
    procedure :: inh       => target_inh_velocity
  end type target_velocity_t
 
 
contains
 
  ! ----------------------------------------------------------------------
  subroutine target_init_velocity(tg, gr, kpoints, namespace, space, ions, qcs, td, w0, oct, ep, restart)
    class(target_velocity_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             :: iatom, ist, jst, jj
    real(real64), allocatable  :: vl(:), vl_grad(:,:)
    type(block_t)       :: blk
    character(len=1024) :: expression
 
    push_sub(target_init_velocity)
 
    !%Variable OCTVelocityTarget
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% If <tt>OCTTargetOperator = oct_tg_velocity</tt>, then one must supply the
    !% target to optimize in terms of the ionic velocities. This is done by
    !% supplying a string through the block <tt>OCTVelocityTarget</tt>.
    !% Each velocity component is supplied by <tt>"v[n_atom,vec_comp]"</tt>,
    !% where <tt>n_atom</tt> is the atom number, corresponding to the
    !% <tt>Coordinates</tt> block, and <tt>vec_comp</tt> is the corresponding
    !% vector component of the velocity. The target string can be
    !% supplied by using several lines in this block.
    !% As an example, the following target can be used to maximize the
    !% velocity difference between atom 1 and 2 (in a 3D system):
    !%
    !% <tt>%OCTVelocityTarget
    !% <br> "(v[1,1]-v[2,1])^2 + (v[1,2]-v[2,2])^2 + "
    !% <br> "(v[1,3]-v[2,3])^2"
    !% <br>%</tt>
    !%
    !%End
 
    !%Variable OCTVelocityDerivatives
    !%Type block
    !%Section Calculation Modes::Optimal Control
    !%Description
    !% If <tt>OCTTargetOperator = oct_tg_velocity</tt>, and
    !% <tt>OCTScheme = oct_cg</tt> or <tt>OCTScheme = oct_bfgs</tt>
    !% then you must supply the target in terms of the ionic velocities AND
    !% the derivatives of the target with respect to the ionic velocity components.
    !% The derivatives are supplied via strings through the block
    !% <tt>OCTVelocityDerivatives</tt>.
    !% Each velocity component is supplied by <tt>"v[n_atom,vec_comp]"</tt>,
    !% while <tt>n_atom</tt> is the atom number, corresponding to the
    !% <tt>Coordinates</tt> block, and <tt>vec_comp</tt> is the corresponding
    !% vector component of the velocity. The first line of the
    !% <tt>OCTVelocityDerivatives</tt> block contains the derivatives
    !% with respect to <tt>v[1,*]</tt>, the second with respect to <tt>v[2,*]</tt> and so
    !% on. The first column contains all derivatives with respect <tt>v[*,1]</tt>,
    !% the second with respect to <tt>v[*,2]</tt> and the third w.r.t. <tt>v[*,3]</tt>.
    !% As an example, we show the <tt>OCTVelocityDerivatives</tt> block
    !% corresponding to the target shown in the <tt>OCTVelocityTarget</tt>
    !% help section:
    !%
    !% <tt>%OCTVelocityDerivatives
    !% <br> " 2*(v[1,1]-v[2,1])" | " 2*(v[1,2]-v[2,2])" | " 2*(v[1,3]-v[2,3])"
    !% <br> "-2*(v[1,1]-v[2,1])" | "-2*(v[1,2]-v[2,2])" | "-2*(v[1,3]-v[2,3])"
    !% <br>%</tt>
    !%
    !%End
 
    if (parse_block(namespace, 'OCTVelocityTarget', blk) == 0) then
      tg%vel_input_string = " "
      do jj=0, parse_block_n(blk)-1
        call parse_block_string(blk, jj, 0, expression)
        tg%vel_input_string = trim(tg%vel_input_string) // trim(expression)
      end do
      call parse_block_end(blk)
    else
      message(1) = 'If OCTTargetOperator = oct_tg_velocity, then you must give the shape'
      message(2) = 'of this target in the block "OCTVelocityTarget".'
      call messages_fatal(2, namespace=namespace)
    end if
 
    tg%move_ions = td%ions_dyn%ions_move()
    if (tg%move_ions) then
      message(1) = 'If OCTTargetOperator = oct_tg_velocity, then you must not allow the ions'
      message(2) = 'to move. If you want to move the ions, then you can get the same functionality'
      message(3) = 'with OCTTargetOperator = oct_tg_classical.'
      call messages_fatal(3, namespace=namespace)
    end if
 
    if (oct%algorithm == option__octscheme__oct_cg .or. oct%algorithm == option__octscheme__oct_bfgs) then
      if (parse_block(namespace, 'OCTVelocityDerivatives', blk) == 0) then
        safe_allocate(tg%vel_der_array(1:ions%natoms,1:gr%box%dim))
        do ist=0, ions%natoms-1
          do jst=0, gr%box%dim-1
            call parse_block_string(blk, ist, jst, tg%vel_der_array(ist+1, jst+1))
          end do
        end do
        call parse_block_end(blk)
      else
        message(1) = 'If OCTTargetOperator = oct_tg_velocity, and'
        message(2) = 'OCTScheme = oct_cg, or OCTScheme = oct_bfgs then you must define the'
        message(3) = 'blocks "OCTVelocityTarget" AND "OCTVelocityDerivatives"'
        call messages_fatal(3, namespace=namespace)
      end if
    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))
 
    ! calculate gradient of each species potential
    do iatom = 1, ions%natoms
      vl(:) = m_zero
      vl_grad(:,:) = m_zero
      call epot_local_potential(ep, namespace, ions%space, ions%latt, gr, ions%atom(iatom)%species, &
        ions%pos(:, iatom), iatom, vl)
      call dderivatives_grad(gr%der, vl, vl_grad)
      do ist = 1, gr%np
        do jst=1, gr%box%dim
          tg%grad_local_pot(iatom, ist, jst) = vl_grad(ist, jst)
        end do
      end do
    end do
    safe_deallocate_a(vl)
    safe_deallocate_a(vl_grad)
 
    ! 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.
 
    tg%dt = td%dt
    safe_allocate(tg%td_fitness(0:td%max_iter))
    tg%td_fitness = m_zero
 
    pop_sub(target_init_velocity)
  end subroutine target_init_velocity
 
 
  ! ----------------------------------------------------------------------
  subroutine target_end_velocity(tg, oct)
    class(target_velocity_t), intent(inout) :: tg
    type(oct_t),    intent(in)    :: oct
 
    push_sub(target_end_velocity)
 
    if (oct%algorithm == option__octscheme__oct_cg .or. oct%algorithm == option__octscheme__oct_bfgs) then
      safe_deallocate_a(tg%vel_der_array)
      safe_deallocate_a(tg%grad_local_pot)
      safe_deallocate_a(tg%rho)
    end if
    safe_deallocate_a(tg%td_fitness)
 
    pop_sub(target_end_velocity)
  end subroutine target_end_velocity
 
 
  ! ----------------------------------------------------------------------
  subroutine target_output_velocity(tg, namespace, space, gr, dir, ions, hm, outp)
    class(target_velocity_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_velocity)
 
    call io_mkdir(trim(dir), namespace)
    call output_states(outp, namespace, space, trim(dir), tg%st, gr, ions, hm, -1)
 
    pop_sub(target_output_velocity)
  end subroutine target_output_velocity
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  real(real64) function target_j1_velocity(tg, namespace, gr, kpoints, qcpsi, ions) result(j1)
    class(target_velocity_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
 
    integer :: i
    real(real64) :: f_re, dummy(3)
    real(real64), allocatable :: x(:, :)
    character(len=4096) :: inp_string
    push_sub(target_j1_velocity)
 
    safe_allocate(x(1:ions%natoms, 1:ions%space%dim))
    do i = 1, ions%natoms
      x(i, :) = ions%vel(:, i)
    end do
 
    f_re = m_zero
    dummy(:) = m_zero
    inp_string = tg%vel_input_string
    call parse_array(inp_string, x, 'v')
    call conv_to_c_string(inp_string)
    call parse_expression(f_re, dummy(1), 1, dummy(1:3), dummy(1), dummy(1), inp_string)
    j1 = f_re
 
    safe_deallocate_a(x)
    pop_sub(target_j1_velocity)
  end function target_j1_velocity
 
 
  ! ----------------------------------------------------------------------
  subroutine target_chi_velocity(tg, namespace, gr, kpoints, qcpsi_in, qcchi_out, ions)
    class(target_velocity_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 :: ip, ist, jst, ik, ib
    character(len=1024) :: temp_string
    real(real64) :: df_dv, dummy(3)
    real(real64), allocatable :: x(:, :)
    type(states_elec_t), pointer :: chi_out
    push_sub(target_chi_velocity)
 
    chi_out => opt_control_point_qs(qcchi_out)
 
    !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
 
    safe_allocate(x(1:ions%natoms, 1:ions%space%dim))
    do ip = 1, ions%natoms
      x(ip, :) = ions%vel(:, ip)
    end do
 
    !calculate dF/dn, which is the time-independent part of the inhomogenous term for the propagation of Chi
    df_dv = m_zero
    dummy(:) = m_zero
    tg%rho(:) = m_zero
    do ist = 1, ions%natoms
      do jst=1, gr%box%dim
        temp_string = tg%vel_der_array(ist, jst)
        call parse_array(temp_string, x, 'v')
        call conv_to_c_string(temp_string)
        call parse_expression(df_dv, dummy(1), 1, dummy(1:3), dummy(1), dummy(1), temp_string)
        tg%rho(:) = tg%rho(:) + df_dv*tg%grad_local_pot(ist,:,jst)/ions%mass(ist)
      end do
    end do
 
    safe_deallocate_a(x)
    nullify(chi_out)
    pop_sub(target_chi_velocity)
  end subroutine target_chi_velocity
 
 
  ! ---------------------------------------------------------------
  subroutine target_inh_velocity(tg, psi, gr, kpoints, time, inh, iter)
    class(target_velocity_t), intent(inout) :: tg
    type(states_elec_t),      intent(inout) :: psi
    type(grid_t),             intent(in)    :: gr
    type(kpoints_t),          intent(in)    :: kpoints
    real(real64),             intent(in)    :: time
    type(states_elec_t),      intent(inout) :: inh
    integer,                  intent(in)    :: iter
 
    integer :: ik, ist, ip, idim
    complex(real64), allocatable :: zpsi(:)
 
    push_sub(target_inh_velocity)
 
    safe_allocate(zpsi(1:gr%np))
 
    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
 
    safe_deallocate_a(zpsi)
 
    pop_sub(target_inh_velocity)
  end subroutine target_inh_velocity
  !----------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine target_tdcalc_velocity(tg, namespace, space, hm, gr, ions, ext_partners, psi, time, max_time)
    class(target_velocity_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
 
    complex(real64), allocatable :: opsi(:, :), zpsi(:, :)
    integer :: iatom, ik, ist, idim
    real(real64) :: dt
    push_sub(target_tdcalc_velocity)
 
    tg%td_fitness(time) = m_zero
 
    safe_allocate(zpsi(1:gr%np_part, 1))
    safe_allocate(opsi(1:gr%np_part, 1))
    opsi = m_z0
    ! WARNING This does not work for spinors.
    do iatom = 1, ions%natoms
      ions%tot_force(:, iatom) = hm%ep%fii(1:gr%box%dim, iatom)
      do ik = 1, psi%nik
        do ist = 1, psi%nst
          do idim = 1, gr%box%dim
            call states_elec_get_state(psi, gr, ist, ik, zpsi)
            opsi(1:gr%np, 1) = tg%grad_local_pot(iatom, 1:gr%np, idim)*zpsi(1:gr%np, 1)
            ions%tot_force(idim, iatom) = ions%tot_force(idim, iatom) &
              + real(psi%occ(ist, ik)*zmf_dotp(gr, psi%d%dim, opsi, zpsi), real64)
          end do
        end do
      end do
    end do
    safe_deallocate_a(opsi)
    safe_deallocate_a(zpsi)
 
    dt = tg%dt
    if ((time == 0) .or. (time == max_time)) dt = tg%dt * m_half
    do iatom = 1, ions%natoms
      ions%vel(:, iatom) = ions%vel(:, iatom) + ions%tot_force(:, iatom) * dt / ions%mass(iatom)
    end do
 
    pop_sub(target_tdcalc_velocity)
  end subroutine target_tdcalc_velocity
  ! ----------------------------------------------------------------------
 
end module target_velocity_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: