td_calc.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 td_calc_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use elec_matrix_elements_oct_m
  use ext_partner_list_oct_m
  use forces_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_base_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lasers_oct_m
  use loct_math_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_calc_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::       &
    td_calc_tacc, &
    td_calc_tvel, &
    td_calc_ionch
 
contains
 
! ---------------------------------------------------------
! ---------------------------------------------------------
  subroutine td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, time)
    type(namespace_t),        intent(in)  :: namespace
    class(space_t),           intent(in)  :: space
    type(grid_t),             intent(in)  :: gr
    type(ions_t),             intent(in)  :: ions
    type(partner_list_t),     intent(in)  :: ext_partners
    type(states_elec_t),      intent(in)  :: st
    type(hamiltonian_elec_t), intent(in)  :: hm
    real(real64),             intent(in)  :: time
    real(real64),             intent(out) :: acc(:)
 
    real(real64) :: field(space%dim), x(space%dim)
    complex(real64), allocatable :: zpsi(:, :), hzpsi(:,:), hhzpsi(:,:), xzpsi(:,:,:), vnl_xzpsi(:,:)
    integer  :: j, k, ik, ist, idim
    type(lasers_t), pointer :: lasers
 
    push_sub(td_calc_tacc)
 
    ! See #683
    if (family_is_mgga(hm%xc%family) .or. family_is_hybrid(hm%xc)) then
      call messages_not_implemented("Dipole acceleration with mGGAs and hybrid functionals")
    end if
 
    ! The term i<[V_l,p]> + i<[V_nl,p]> may be considered as equal but opposite to the
    ! force exerted by the electrons on the ions. COMMENT: This has to be thought about.
    ! Maybe we are forgetting something....
    call total_force_calculate(space, gr, ions, hm%ep, st, hm%kpoints, acc, hm%lda_u_level)
 
    ! Adds the laser contribution : i<[V_laser, p]>
    ! WARNING: this ignores the possibility of non-electric td external fields.
    field = m_zero
    lasers => list_get_lasers(ext_partners)
    if(associated(lasers)) then
      do j = 1, lasers%no_lasers
        call laser_electric_field(lasers%lasers(j), field(1:space%dim), time, 0.001_real64)
        acc(1:space%dim) = acc(1:space%dim) - st%qtot*field(1:space%dim)
      end do
    end if
 
    if (.not. hm%ep%non_local) then
      pop_sub(td_calc_tacc)
      return
    end if
 
    ! And now, i<[H,[V_nl,x]]>
    x = m_zero
    safe_allocate(zpsi(1:gr%np_part, 1:st%d%dim))
    safe_allocate(hzpsi(1:gr%np_part, 1:st%d%dim))
    safe_allocate(hhzpsi(1:3, 1:gr%np_part))
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ist = st%st_start, st%st_end
 
        call states_elec_get_state(st, gr, ist, ik, zpsi)
 
        call zhamiltonian_elec_apply_single(hm, namespace, gr, zpsi, hzpsi, ist, ik)
 
        safe_allocate(xzpsi(1:gr%np_part, 1:st%d%dim, 1:3))
        safe_allocate(vnl_xzpsi(1:gr%np_part, 1:st%d%dim))
        xzpsi = m_z0
        do k = 1, gr%np
          do j = 1, space%dim
            xzpsi(k, 1:st%d%dim, j) = gr%x(k, j)*zpsi(k, 1:st%d%dim)
          end do
        end do
 
        do j = 1, space%dim
          call zhamiltonian_elec_apply_single(hm, namespace, gr, xzpsi(:, :, j), vnl_xzpsi, ist, ik, &
            terms = term_non_local_potential)
 
          do idim = 1, st%d%dim
            x(j) = x(j) - 2*st%occ(ist, ik)*real(zmf_dotp(gr, hzpsi(1:gr%np, idim), vnl_xzpsi(:, idim)), real64)
          end do
        end do
 
        xzpsi = m_z0
        do k = 1, gr%np
          do j = 1, space%dim
            xzpsi(k, 1:st%d%dim, j) = gr%x(k, j)*hzpsi(k, 1:st%d%dim)
          end do
        end do
 
        do j = 1, space%dim
          call zhamiltonian_elec_apply_single(hm, namespace, gr, xzpsi(:, :, j), vnl_xzpsi, ist, ik, &
            terms = term_non_local_potential)
 
          do idim = 1, st%d%dim
            x(j) = x(j) + 2*st%occ(ist, ik)*real(zmf_dotp(gr, zpsi(:, idim), vnl_xzpsi(:, idim)), real64)
          end do
        end do
        safe_deallocate_a(xzpsi)
        safe_deallocate_a(vnl_xzpsi)
 
      end do
    end do
    safe_deallocate_a(hzpsi)
    safe_deallocate_a(hhzpsi)
 
    if (st%parallel_in_states) then
      call comm_allreduce(st%mpi_grp, x)
    end if
    acc = acc + x
 
    pop_sub(td_calc_tacc)
  end subroutine td_calc_tacc
 
! ---------------------------------------------------------
! ---------------------------------------------------------
  subroutine td_calc_tvel(elec_me, kpoints, vel)
    type(elec_matrix_elements_t), intent(in)  :: elec_me
    type(kpoints_t),              intent(in)  :: kpoints
    real(real64),                 intent(out) :: vel(:)
 
    real(real64), allocatable :: momentum(:,:,:)
    integer :: st_start, st_end
 
    push_sub(td_calc_tvel)
 
    st_start = elec_me%states%st_start
    st_end = elec_me%states%st_end
 
    safe_allocate(momentum(1:elec_me%space%dim, st_start:st_end, elec_me%states%d%kpt%start:elec_me%states%d%kpt%end))
    call elec_me%momentum_me(kpoints, momentum)
 
    momentum(1:elec_me%space%dim, st_start:st_end, 1) = &
      sum(momentum(1:elec_me%space%dim, st_start:st_end, elec_me%states%d%kpt%start:elec_me%states%d%kpt%end), 3)
    momentum(1:elec_me%space%dim, 1, 1) = sum(momentum(1:elec_me%space%dim, st_start:st_end, 1), 2)
    vel = momentum(:, 1, 1)
 
    safe_deallocate_a(momentum)
 
    pop_sub(td_calc_tvel)
  end subroutine td_calc_tvel
 
! ---------------------------------------------------------
! ---------------------------------------------------------
  subroutine td_calc_ionch(mesh, st, ch, Nch)
    type(mesh_t),        intent(in)    :: mesh
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: nch
    real(real64),        intent(out)   :: ch(0:nch)
 
    integer :: ik, ist, ii, jj, idim, nid
    real(real64)   :: prod, prod0
    real(real64), allocatable :: n(:), nnot(:)
    complex(real64),  allocatable :: zpsi(:)
    !combinations
    integer :: next
    type(c_ptr) :: c
    integer, allocatable :: idx0(:), idx(:), idxref(:)
 
    push_sub(td_calc_ionch)
 
    safe_allocate(   n(1: nch))
    safe_allocate(nnot(1: nch))
    safe_allocate(zpsi(1:mesh%np))
 
    n(:)   = m_zero
    nnot(:)= m_zero
 
 
    ii = 1
    do ik = 1, st%nik
      do ist = 1, st%nst
        do idim = 1, st%d%dim
 
          if (st%st_start <= ist .and. ist <= st%st_end .and. &
            st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then
            call states_elec_get_state(st, mesh, idim, ist, ik, zpsi)
            n(ii) = zmf_nrm2(mesh, zpsi)
            nnot(ii) = m_one - n(ii)**2
          end if
 
          ii = ii + 1
 
        end do
      end do
    end do
 
    if (st%parallel_in_states) then
      call comm_allreduce(st%mpi_grp, n)
      call comm_allreduce(st%mpi_grp, nnot)
    end if
 
! print * ,mpi_world%rank, "N    =", N(:)
! print * ,mpi_world%rank, "Nnot =", Nnot(:)
 
    ch(:) = m_zero
 
    nid = nch - 1
    safe_allocate(idx(0:nid))
    safe_allocate(idxref(0:nid))
    idxref = (/ (ii, ii = 0, nid) /)
    do ii = 0, nch
!     print *, "Nch= ", Nch, "ii", ii
      call loct_combination_init(c, nch, ii)
      if (ii == 0) then
        safe_allocate(idx0(0:0))
      else
        safe_allocate(idx0(0:ii-1))
      end if
!     print *,"size(idx0,1)=",size(idx0,1)
      if (debug%info) then
        call messages_write("P(")
        call messages_write(ii)
        call messages_write(") = 0")
        call messages_new_line()
      end if
      ! Loop over combinations
      do
        prod  = m_one
        prod0 = m_one
        if (debug%info) then
          call messages_write("P(")
          call messages_write(ii)
          call messages_write(") += ")
        end if
        call loct_get_combination(c, idx0)
        idx(:) = idxref(:)
        do jj = 0, ii-1
          idx(idx0(jj))= -1
          if (debug%info) then
            call messages_write(" No(")
            call messages_write(idx0(jj)+1)
            call messages_write(") ")
          end if
          prod0 = prod0 * nnot(idx0(jj)+1)
        end do
 
        do jj = 0, nid
          if (idx(jj) >= 0) then
            if (debug%info) then
              call messages_write(" N(")
              call messages_write(idx(jj)+1)
              call messages_write(") ")
            end if
            prod = prod * n(idx(jj)+1)
          end if
        end do
 
        if (debug%info) call messages_new_line()
 
        ch(ii) = ch(ii) + prod*prod0
 
        call loct_combination_next(c, next)
        if (next /= 0) exit
      end do
      safe_deallocate_a(idx0)
      call loct_combination_end(c)
 
      if (debug%info) call messages_info()
    end do
 
    safe_deallocate_a(idx)
    safe_deallocate_a(idxref)
 
 
    safe_deallocate_a(n)
    safe_deallocate_a(nnot)
    safe_deallocate_a(zpsi)
 
    pop_sub(td_calc_ionch)
  end subroutine td_calc_ionch
 
 
end module td_calc_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: