target_excited.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_excited_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use excited_states_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, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use output_oct_m
  use output_low_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
  use states_elec_restart_oct_m
  use target_low_oct_m
  use td_oct_m
  use types_oct_m
 
  implicit none
 
  private
  public ::                &
    target_init_excited,   &
    target_j1_excited,     &
    target_output_excited, &
    target_chi_excited
 
 
contains
 
 
  ! ----------------------------------------------------------------------
  subroutine target_init_excited(mesh, namespace, space, tg, td, restart, kpoints)
    class(mesh_t),     intent(in)    :: mesh
    type(namespace_t), intent(in)    :: namespace
    class(space_t),    intent(in)    :: space
    type(target_t),    intent(inout) :: tg
    type(td_t),        intent(in)    :: td
    type(restart_t),   intent(in)    :: restart
    type(kpoints_t),   intent(in)    :: kpoints
 
    integer :: ierr, nik, dim
 
    push_sub(target_init_excited)
 
    message(1) =  'Info: TargetOperator is a linear combination of Slater determinants.'
    call messages_info(1, namespace=namespace)
 
    tg%move_ions = td%ions_dyn%ions_move()
    tg%dt = td%dt
 
    call states_elec_look(restart, nik, dim, tg%st%nst, ierr)
    if (ierr /= 0) then
      message(1) = "Unable to read states information."
      call messages_fatal(1, namespace=namespace)
    end if
    tg%st%st_start = 1
    tg%st%st_end   = tg%st%nst
 
    safe_deallocate_a(tg%st%occ)
    safe_deallocate_a(tg%st%eigenval)
    safe_deallocate_a(tg%st%node)
 
    safe_allocate(     tg%st%occ(1:tg%st%nst, 1:tg%st%nik))
    safe_allocate(tg%st%eigenval(1:tg%st%nst, 1:tg%st%nik))
    safe_allocate(    tg%st%node(1:tg%st%nst))
    if (tg%st%d%ispin == spinors) then
      safe_deallocate_a(tg%st%spin)
      safe_allocate(tg%st%spin(1:3, 1:tg%st%nst, 1:tg%st%nik))
    end if
    call states_elec_allocate_wfns(tg%st, mesh, type_cmplx)
    tg%st%node(:)  = 0
 
    call states_elec_load(restart, namespace, space, tg%st, mesh, kpoints, ierr)
    if (ierr /= 0) then
      message(1) = "Unable to read wavefunctions."
      call messages_fatal(1, namespace=namespace)
    end if
 
    call excited_states_init(tg%est, tg%st, "oct-excited-state-target", namespace)
 
    pop_sub(target_init_excited)
  end subroutine target_init_excited
 
 
  ! ----------------------------------------------------------------------
  subroutine target_output_excited(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_excited)
 
    call io_mkdir(trim(dir), namespace)
    call output_states(outp, namespace, space, trim(dir)//'/st', tg%est%st, gr, ions, hm, -1)
    call excited_states_output(tg%est, trim(dir), namespace)
 
    pop_sub(target_output_excited)
  end subroutine target_output_excited
  ! ----------------------------------------------------------------------
 
 
  ! ----------------------------------------------------------------------
  real(real64) function target_j1_excited(tg, namespace, gr, psi) result(j1)
    type(target_t),      intent(in) :: tg
    type(namespace_t),   intent(in)    :: namespace
    type(grid_t),        intent(in) :: gr
    type(states_elec_t), intent(in) :: psi
 
    push_sub(target_j1_excited)
 
    j1 = abs(zstates_elec_mpdotp(namespace, gr, tg%est, psi))**2
 
    pop_sub(target_j1_excited)
  end function target_j1_excited
 
 
  ! ----------------------------------------------------------------------
  subroutine target_chi_excited(tg, namespace, gr, psi_in, chi_out)
    type(target_t),         intent(in)    :: tg
    type(namespace_t),      intent(in)    :: namespace
    type(grid_t),           intent(in)    :: gr
    type(states_elec_t),    intent(in)    :: psi_in
    type(states_elec_t),    intent(inout) :: chi_out
 
    complex(real64), allocatable :: ci(:), di(:), mat(:, :, :), mm(:, :, :, :), mk(:, :), lambda(:, :)
    complex(real64), allocatable :: zpsi(:, :), zchi(:, :)
    integer :: ik, ist, jst, ia, ib, n_pairs, nst, kpoints, jj, idim, ip
    push_sub(target_chi_excited)
 
    n_pairs = tg%est%n_pairs
    kpoints = psi_in%nik
    nst = psi_in%nst
 
 
    safe_allocate(zpsi(1:gr%np, 1:psi_in%d%dim))
    safe_allocate(zchi(1:gr%np, 1:psi_in%d%dim))
    safe_allocate(ci(1:n_pairs))
    safe_allocate(di(1:n_pairs))
    safe_allocate(mat(1:tg%est%st%nst, 1:nst, 1:psi_in%nik))
    safe_allocate(mm(1:nst, 1:nst, 1:kpoints, 1:n_pairs))
    safe_allocate(mk(1:gr%np_part, 1:psi_in%d%dim))
    safe_allocate(lambda(1:n_pairs, 1:n_pairs))
 
    call zstates_elec_matrix(tg%est%st, psi_in, gr, mat)
 
    do ia = 1, n_pairs
      ci(ia) = tg%est%weight(ia)
      call zstates_elec_matrix_swap(mat, tg%est%pair(ia))
      mm(1:nst, 1:nst, 1:kpoints, ia) = mat(1:nst, 1:kpoints, 1:kpoints)
      di(ia) = zstates_elec_mpdotp(namespace, gr, tg%est%st, psi_in, mat)
      if (abs(di(ia)) > 1.0e-12_real64) then
        do ik = 1, kpoints
          call lalg_inverse(nst, mm(1:nst, 1:nst, ik, ia), 'dir')
        end do
      end if
      call zstates_elec_matrix_swap(mat, tg%est%pair(ia))
    end do
 
    do ia = 1, n_pairs
      do ib = 1, n_pairs
        lambda(ia, ib) = conjg(ci(ib)) * ci(ia) * conjg(di(ia)) * di(ib)
      end do
    end do
 
    select case (psi_in%d%ispin)
    case (unpolarized)
      write(message(1), '(a)') 'Internal error in target.target_chi: unpolarized.'
      call messages_fatal(1, namespace=namespace)
 
    case (spin_polarized)
      assert(chi_out%nik == 2)
 
      do ik = 1, kpoints
        do ist = chi_out%st_start, chi_out%st_end
 
          zchi(1:gr%np, 1:psi_in%d%dim) = m_z0
 
          do ia = 1, n_pairs
            if (ik /= tg%est%pair(ia)%kk) cycle
            if (abs(di(ia)) < 1.0e-12_real64) cycle
            do ib = 1, n_pairs
              if (abs(di(ib)) < 1.0e-12_real64) cycle
              mk = m_z0
 
              do jst = 1, nst
                if (jst == tg%est%pair(ib)%i) jj = tg%est%pair(ia)%a
                call states_elec_get_state(tg%est%st, gr, jj, ik, zpsi)
 
                do idim = 1, psi_in%d%dim
                  do ip = 1, gr%np
                    mk(ip, idim) = mk(ip, idim) + conjg(mm(ist, jst, ik, ib))*zpsi(ip, idim)
                  end do
                end do
              end do
 
              call lalg_axpy(gr%np_part, psi_in%d%dim, m_z1, lambda(ib, ia)*mk(:, :), zchi)
 
            end do
          end do
 
          call states_elec_set_state(chi_out, gr, ist, ik, zchi)
 
        end do
      end do
 
    case (spinors)
      assert(chi_out%nik == 1)
 
      do ist = chi_out%st_start, chi_out%st_end
 
        zchi(1:gr%np, 1:psi_in%d%dim) = m_z0
 
        do ia = 1, n_pairs
          if (abs(di(ia)) < 1.0e-12_real64) cycle
 
          do ib = 1, n_pairs
            if (abs(di(ib)) < 1.0e-12_real64) cycle
 
            mk = m_z0
            do jst = 1, nst
              if (jst == tg%est%pair(ib)%i) jj = tg%est%pair(ia)%a
              call states_elec_get_state(tg%est%st, gr, jj, ik, zpsi)
 
              do idim = 1, psi_in%d%dim
                do ip = 1, gr%np
                  mk(ip, idim) = mk(ip, idim) + conjg(mm(ist, jst, 1, ib))*zpsi(ip, idim)
                end do
              end do
            end do
 
            call lalg_axpy(gr%np_part, 2, m_z1, lambda(ib, ia)*mk(:, :), zchi)
          end do
        end do
 
        call states_elec_set_state(chi_out, gr, ist, ik, zchi)
 
      end do
 
    end select
 
    safe_deallocate_a(zpsi)
    safe_deallocate_a(zchi)
    safe_deallocate_a(ci)
    safe_deallocate_a(di)
    safe_deallocate_a(mat)
    safe_deallocate_a(mm)
    safe_deallocate_a(mk)
    safe_deallocate_a(lambda)
    pop_sub(target_chi_excited)
  end subroutine target_chi_excited
 
end module target_excited_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: