propagation_ops_elec.F90

Go to the documentation of this file.

!! Copyright (C) 2019 N. Tancogne-Dejean
!!
!! 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 propagation_ops_elec_oct_m
  use accel_oct_m
  use batch_oct_m
  use debug_oct_m
  use density_oct_m
  use electron_space_oct_m
  use exponential_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use interaction_partner_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lasers_oct_m
  use lda_u_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use mxll_elec_coupling_oct_m
  use namespace_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use propagator_verlet_oct_m
  use space_oct_m
  use states_elec_oct_m
  use varinfo_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                             &
    propagation_ops_elec_t,                      &
    propagation_ops_elec_update_hamiltonian,     &
    propagation_ops_elec_exp_apply,              &
    propagation_ops_elec_fuse_density_exp_apply, &
    propagation_ops_elec_move_ions,              &
    propagation_ops_elec_restore_ions,           &
    propagation_ops_elec_propagate_gauge_field,  &
    propagation_ops_elec_restore_gauge_field,    &
    propagation_ops_elec_interpolate_get,        &
    propagation_ops_do_pack,                     &
    propagation_ops_do_unpack,                   &
    propagation_ops_elec_propagate_ions_and_cell,&
    propagation_ops_finish_unpack
 
  type :: propagation_ops_elec_t
    private
    type(ion_state_t) :: ions_state
    real(real64), allocatable :: vecpot(:), vecpot_vel(:)
  end type propagation_ops_elec_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_update_hamiltonian(namespace, space, st, mesh, hm, ext_partners, time)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(states_elec_t),      intent(inout) :: st
    class(mesh_t),            intent(in)    :: mesh
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    real(real64),             intent(in)    :: time
 
 
    push_sub(propagation_ops_elec_update_hamiltonian)
 
    call profiling_in('ELEC_UPDATE_H')
 
    call calculate_mxll_dipole_field(hm, mesh, st)
    call hm%update(mesh, namespace, space, ext_partners, time = time)
    call lda_u_update_occ_matrices(hm%lda_u, namespace, mesh, st, hm%phase, hm%energy)
 
    call profiling_out('ELEC_UPDATE_H')
 
    pop_sub(propagation_ops_elec_update_hamiltonian)
  end subroutine propagation_ops_elec_update_hamiltonian
 
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_move_ions(wo, gr, hm, st, namespace, space, ions_dyn, ions, &
    ext_partners, mc, time, dt, save_pos)
    class(propagation_ops_elec_t), intent(inout) :: wo
    type(grid_t),                  intent(inout) :: gr
    type(hamiltonian_elec_t),      intent(inout) :: hm
    type(states_elec_t),           intent(inout) :: st
    type(namespace_t),             intent(in)    :: namespace
    type(electron_space_t),        intent(in)    :: space
    type(ion_dynamics_t),          intent(inout) :: ions_dyn
    type(ions_t),                  intent(inout) :: ions
    type(partner_list_t),          intent(in)    :: ext_partners
    type(multicomm_t),             intent(inout) :: mc
    real(real64),                  intent(in)    :: time
    real(real64),                  intent(in)    :: dt
    logical, optional,             intent(in)    :: save_pos
 
    real(real64) :: dt_ions
 
    push_sub(propagation_ops_elec_move_ions)
 
    call profiling_in('ELEC_MOVE_IONS')
 
 
    if (ions_dyn%is_active()) then
      dt_ions = dt * ions_dyn%ionic_scale
      if (optional_default(save_pos, .false.)) then
        call ion_dynamics_save_state(ions_dyn, ions, wo%ions_state)
      end if
 
      call propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, ions_dyn, ions, &
        mc, time, dt_ions)
 
      call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = time)
    end if
 
    call profiling_out('ELEC_MOVE_IONS')
 
    pop_sub(propagation_ops_elec_move_ions)
  end subroutine propagation_ops_elec_move_ions
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, ions_dyn, ions, &
    mc, time, dt_ions)
    type(grid_t),                  intent(inout) :: gr
    type(hamiltonian_elec_t),      intent(inout) :: hm
    type(states_elec_t),           intent(inout) :: st
    type(namespace_t),             intent(in)    :: namespace
    type(electron_space_t),        intent(in)    :: space
    type(ion_dynamics_t),          intent(inout) :: ions_dyn
    type(ions_t),                  intent(inout) :: ions
    type(multicomm_t),             intent(inout) :: mc
    real(real64),                  intent(in)    :: time
    real(real64),                  intent(in)    :: dt_ions
 
    push_sub(propagation_ops_elec_propagate_ions_and_cell)
 
    call ion_dynamics_propagate(ions_dyn, ions, time, dt_ions, namespace)
 
    ! Update lattice vectors and regenerate grid
    if (ions_dyn%cell_relax()) then
      call ion_dynamics_box_update(namespace, gr, space, ions%latt)
      call electrons_lattice_vectors_update(namespace, gr, space, hm%psolver, hm%kpoints, &
        mc, st%qtot, ions%latt)
    end if
 
    pop_sub(propagation_ops_elec_propagate_ions_and_cell)
  end subroutine propagation_ops_elec_propagate_ions_and_cell
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_restore_ions(wo, ions_dyn, ions)
    class(propagation_ops_elec_t),    intent(inout) :: wo
    type(ion_dynamics_t),    intent(inout) :: ions_dyn
    type(ions_t),            intent(inout) :: ions
 
 
    push_sub(propagation_ops_elec_restore_ions)
 
    call profiling_in('ELEC_RESTORE_IONS')
 
    if (ions_dyn%is_active()) then
      call ion_dynamics_restore_state(ions_dyn, ions, wo%ions_state)
    end if
 
    call profiling_out('ELEC_RESTORE_IONS')
 
    pop_sub(propagation_ops_elec_restore_ions)
  end subroutine propagation_ops_elec_restore_ions
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_propagate_gauge_field(wo, gfield, dt, time, save_gf)
    class(propagation_ops_elec_t), intent(inout) :: wo
    type(gauge_field_t),           intent(inout) :: gfield
    real(real64),                  intent(in)    :: dt
    real(real64),                  intent(in)    :: time
    logical,  optional,            intent(in)    :: save_gf
 
 
    push_sub(propagation_ops_elec_propagate_gauge_field)
 
    call profiling_in('ELEC_MOVE_GAUGE')
 
    if (gauge_field_is_propagated(gfield)) then
      if (optional_default(save_gf, .false.)) then
        safe_allocate(wo%vecpot(1:gfield%space%dim))
        safe_allocate(wo%vecpot_vel(1:gfield%space%dim))
        call gauge_field_get_vec_pot(gfield, wo%vecpot)
        call gauge_field_get_vec_pot_vel(gfield, wo%vecpot_vel)
      end if
      call gauge_field_do_algorithmic_operation(gfield, op_verlet_compute_acc, dt, time)
    end if
 
    call profiling_out('ELEC_MOVE_GAUGE')
 
    pop_sub(propagation_ops_elec_propagate_gauge_field)
  end subroutine propagation_ops_elec_propagate_gauge_field
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_restore_gauge_field(wo, namespace, space, hm, mesh, ext_partners)
    class(propagation_ops_elec_t), intent(inout) :: wo
    type(namespace_t),             intent(in)    :: namespace
    class(space_t),                intent(in)    :: space
    type(hamiltonian_elec_t),      intent(inout) :: hm
    class(mesh_t),                 intent(in)    :: mesh
    type(partner_list_t),          intent(in)    :: ext_partners
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(propagation_ops_elec_restore_gauge_field)
 
    call profiling_in('ELEC_RESTORE_GAUGE')
 
    gfield => list_get_gauge_field(ext_partners)
    if (associated(gfield)) then
      call gauge_field_set_vec_pot(gfield, wo%vecpot)
      call gauge_field_set_vec_pot_vel(gfield, wo%vecpot_vel)
      safe_deallocate_a(wo%vecpot)
      safe_deallocate_a(wo%vecpot_vel)
      call hm%update(mesh, namespace, space, ext_partners)
    end if
 
    call profiling_out('ELEC_RESTORE_GAUGE')
 
    pop_sub(propagation_ops_elec_restore_gauge_field)
  end subroutine propagation_ops_elec_restore_gauge_field
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_exp_apply(te, namespace, st, mesh, hm, dt, op)
    type(exponential_t),         intent(inout) :: te
    type(namespace_t),           intent(in)    :: namespace
    type(states_elec_t),         intent(inout) :: st
    class(mesh_t),               intent(in)    :: mesh
    type(hamiltonian_elec_t),    intent(inout) :: hm
    real(real64),                intent(in)    :: dt
    class(operator_t), optional, intent(in)    :: op
 
    integer :: ik, ib
 
    push_sub(propagation_ops_elec_exp_apply)
 
    call profiling_in('ELEC_EXP_APPLY')
 
    do ik = st%d%kpt%start, st%d%kpt%end
      call propagation_ops_do_pack(st, hm, st%group%block_start, ik)
      do ib = st%group%block_start, st%group%block_end
        if (ib + 1 <= st%group%block_end) call propagation_ops_do_pack(st, hm, ib+1, ik)
        call accel_set_stream(ib)
 
        call hm%phase%set_phase_corr(mesh, st%group%psib(ib, ik), async=.true.)
        if (hamiltonian_elec_inh_term(hm)) then
          call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt, &
            inh_psib = hm%inh_st%group%psib(ib, ik), op=op)
        else
          call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt, op=op)
        end if
        call hm%phase%unset_phase_corr(mesh, st%group%psib(ib, ik))
 
        call propagation_ops_do_unpack(st, hm, ib, ik)
        if (ib-1 >= st%group%block_start) call propagation_ops_finish_unpack(st, hm, ib-1, ik)
      end do
      call propagation_ops_finish_unpack(st, hm, st%group%block_end, ik)
    end do
 
    call accel_finish()
 
    call profiling_out('ELEC_EXP_APPLY')
 
    pop_sub(propagation_ops_elec_exp_apply)
 
  end subroutine propagation_ops_elec_exp_apply
 
 
  subroutine propagation_ops_elec_fuse_density_exp_apply(te, namespace, st, gr, hm, dt, dt2, op, st_td)
    type(exponential_t),           intent(inout) :: te
    type(namespace_t),             intent(in)    :: namespace
    type(states_elec_t),           intent(inout) :: st
    type(grid_t),                  intent(in)    :: gr
    type(hamiltonian_elec_t),      intent(inout) :: hm
    real(real64),                  intent(in)    :: dt
    real(real64),        optional, intent(in)    :: dt2
    class(operator_t),   optional, intent(in)    :: op
    type(states_elec_t), optional, target, intent(inout) :: st_td
    ! potential update of orbital-dependent functionals. This routine only sets the wave functions
    ! but several API expect a states instance
 
    integer                   :: ik, ib
    type(wfs_elec_t), pointer :: zpsib_dt
    type(density_calc_t)      :: dens_calc
 
    push_sub(propagation_ops_elec_fuse_density_exp_apply)
 
    call profiling_in('ELEC_FUSE_DENS_EXP_APPLY')
 
    if (present(st_td)) then
      call density_calc_init(dens_calc, st, gr, st_td%rho)
    else
      call density_calc_init(dens_calc, st, gr, st%rho)
    endif
 
    if (present(dt2) .and. .not. present(st_td)) then
      safe_allocate(zpsib_dt)
    end if
 
    do ik = st%d%kpt%start, st%d%kpt%end
      call propagation_ops_do_pack(st, hm, st%group%block_start, ik)
 
      do ib = st%group%block_start, st%group%block_end
 
        ! Pack next batch. On GPU, this overlap comms
        if (ib + 1 <= st%group%block_end) call propagation_ops_do_pack(st, hm, ib+1, ik)
        ! Switch stream back to current block for work
        call accel_set_stream(ib)
 
        call hm%phase%set_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
 
        if (present(dt2)) then
          if (present(st_td)) then
            zpsib_dt => st_td%group%psib(ib, ik)
            call propagation_ops_do_pack(st_td, hm, ib, ik, copy = .false.)
          else
            call st%group%psib(ib, ik)%copy_to(zpsib_dt)
            if (st%group%psib(ib, ik)%is_packed()) call zpsib_dt%do_pack(copy = .false.)
          end if
 
          !propagate the state to dt/2 and dt, simultaneously, with H(time - dt)
          if (hamiltonian_elec_inh_term(hm)) then
            call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, &
              deltat2 = dt2, inh_psib = hm%inh_st%group%psib(ib, ik), op = op)
          else
            call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, &
              deltat2 = dt2, op = op)
          end if
 
          call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
 
          !use the dt propagation to calculate the density
          call density_calc_accumulate(dens_calc, zpsib_dt)
 
          if (present(st_td)) then
            call hm%phase%unset_phase_corr(gr, zpsib_dt, async=.true.)
            call propagation_ops_do_unpack(st_td, hm, ib, ik)
          else
            call zpsib_dt%end()
          end if
 
        else
          !propagate the state to dt with H(time - dt)
          if (hamiltonian_elec_inh_term(hm)) then
            call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, &
              inh_psib = hm%inh_st%group%psib(ib, ik), op = op)
          else
            call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, op = op)
          end if
          call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
 
          !use the dt propagation to calculate the density
          call density_calc_accumulate(dens_calc, st%group%psib(ib, ik), async=.true.)
        end if
 
        call propagation_ops_do_unpack(st, hm, ib, ik)
        if (ib-1 >= st%group%block_start) then
          call propagation_ops_finish_unpack(st, hm, ib-1, ik)
          if (present(st_td)) call propagation_ops_finish_unpack(st_td, hm, ib-1, ik)
        end if
      end do
 
      call propagation_ops_finish_unpack(st, hm, st%group%block_end, ik)
      if (present(st_td)) call propagation_ops_finish_unpack(st_td, hm, st%group%block_end, ik)
    end do
 
    call density_calc_end(dens_calc)
 
    if (present(dt2) .and. .not. present(st_td)) then
      safe_deallocate_p(zpsib_dt)
    end if
 
    call accel_finish()
 
    call profiling_out('ELEC_FUSE_DENS_EXP_APPLY')
 
    pop_sub(propagation_ops_elec_fuse_density_exp_apply)
  end subroutine propagation_ops_elec_fuse_density_exp_apply
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_do_pack(st, hm, ib, ik, copy)
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: ib
    integer,                  intent(in)    :: ik
    logical, optional,        intent(in)    :: copy
 
    logical :: copy_
 
    push_sub(propagation_ops_do_pack)
 
    copy_ = optional_default(copy, .true.)
    if (hm%apply_packed()) then
      call accel_set_stream(ib)
      call st%group%psib(ib, ik)%do_pack(copy=copy_, async=.true.)
      if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_pack(async=.true.)
    end if
 
    pop_sub(propagation_ops_do_pack)
 
  end subroutine propagation_ops_do_pack
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_do_unpack(st, hm, ib, ik)
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: ib
    integer,                  intent(in)    :: ik
 
    push_sub(propagation_ops_do_unpack)
    if (hm%apply_packed()) then
      call accel_set_stream(ib)
      call st%group%psib(ib, ik)%do_unpack(async=.true.)
      if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_unpack(async=.true.)
    end if
    pop_sub(propagation_ops_do_unpack)
  end subroutine propagation_ops_do_unpack
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_finish_unpack(st, hm, ib, ik)
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: ib
    integer,                  intent(in)    :: ik
 
    push_sub(propagation_ops_finish_unpack)
    if (hm%apply_packed()) then
      call accel_set_stream(ib)
      call st%group%psib(ib, ik)%finish_unpack()
      if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%finish_unpack()
    end if
    pop_sub(propagation_ops_finish_unpack)
  end subroutine propagation_ops_finish_unpack
 
  ! ---------------------------------------------------------
  subroutine propagation_ops_elec_interpolate_get(hm, vks_old)
    type(hamiltonian_elec_t),        intent(inout) :: hm
    type(potential_interpolation_t), intent(inout) :: vks_old
 
    push_sub(propagation_ops_elec_interpolate_get)
 
    call hm%ks_pot%get_interpolated_potentials(vks_old, 0)
 
    pop_sub(propagation_ops_elec_interpolate_get)
 
  end subroutine propagation_ops_elec_interpolate_get
 
  ! ---------------------------------------------------------
  subroutine calculate_mxll_dipole_field(hm, mesh, st)
    type(hamiltonian_elec_t), target,intent(inout) :: hm
    class(mesh_t),                   intent(in)    :: mesh
    type(states_elec_t),             intent(in)    :: st
 
    real(real64), parameter :: density_threshold = 1.0e-8_real64
    integer :: idir
 
    if (.not. hm%mxll%calc_field_dip) return
 
    select case (hm%mxll%coupling_mode)
    case (length_gauge_dipole, multipolar_expansion)
      if (hm%mxll%add_electric_dip) then
        hm%mxll%e_field_dip = get_field_dip(hm%mxll%e_field)
      end if
      if (hm%mxll%add_magnetic_dip) then
        hm%mxll%b_field_dip = get_field_dip(hm%mxll%b_field)
      end if
    case (velocity_gauge_dipole)
      hm%mxll%vec_pot_dip = get_field_dip(hm%mxll%vec_pot)
    case default
      return
    end select
 
  contains
 
    function get_field_dip(field_full) result(field_dip)
      real(real64) :: field_dip(3)
      real(real64), contiguous, intent(in) :: field_full(:,:)
 
      real(real64), allocatable :: total_density(:), mask_density(:)
      real(real64) :: integral_mask
 
      field_dip(1:3) = m_zero
 
      select case (hm%mxll%dipole_field)
      case (dipole_average)
        safe_allocate(total_density(1:mesh%np))
        safe_allocate(mask_density(size(total_density)))
        total_density = sum(st%rho(1:mesh%np, 1:st%d%spin_channels), dim=2)
        ! This mask defines a region in which the fields are averaged
        ! The region is determined by a non-zero electronic density (density > density_threshold)
        ! The purpose of defining a region is to average the EM fields in regions overlapping with the
        ! electronic system, and not including values of the EM fields in the rest of the box.
        mask_density = merge(m_one, m_zero, total_density > density_threshold)
        integral_mask = dmf_integrate(mesh, mask_density)
        do idir = 1, mesh%box%dim
          field_dip(idir) = dmf_integrate(mesh, mask_density*field_full(:,idir))/integral_mask
        end do
        safe_deallocate_a(total_density)
        safe_deallocate_a(mask_density)
 
      case (dipole_at_com)
        if (mesh%mpi_grp%rank == hm%mxll%center_of_mass_rankmin) then
          field_dip(1:3) = field_full(hm%mxll%center_of_mass_ip,1:3)
        end if
        call mesh%allreduce(field_dip(:))
      end select
 
    end function get_field_dip
 
  end subroutine calculate_mxll_dipole_field
 
 
end module propagation_ops_elec_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: