propagator_expmid.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 propagator_expmid_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use ext_partner_list_oct_m
  use grid_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use ks_potential_oct_m
  use lda_u_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use potential_interpolation_oct_m
  use propagator_base_oct_m
  use propagation_ops_elec_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_group_oct_m
  use v_ks_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                        &
    exponential_midpoint,          &
    exponential_midpoint_predictor
 
contains
 
  ! ---------------------------------------------------------
  subroutine exponential_midpoint(hm, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    type(hamiltonian_elec_t), target, intent(inout) :: hm
    type(namespace_t),                intent(in)    :: namespace
    type(electron_space_t),           intent(in)    :: space
    type(grid_t),             target, intent(inout) :: gr
    type(states_elec_t),      target, intent(inout) :: st
    type(propagator_base_t),  target, intent(inout) :: tr
    real(real64),                     intent(in)    :: time
    real(real64),                     intent(in)    :: dt
    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
 
    type(gauge_field_t), pointer  :: gfield
 
    push_sub(propagator_dt.exponential_midpoint)
 
    call hm%ks_pot%interpolate_potentials(tr%vks_old, 3, time, dt, time - dt/m_two)
 
    !move the ions to time 'time - dt/2'
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time - m_half*dt, m_half*dt, save_pos = .true.)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, &
        m_half*dt, time, save_gf = .true.)
    end if
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt*m_half)
 
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt)
 
    !restore to time 'time - dt'
    call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions)
 
    if(associated(gfield)) then
      call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners)
    end if
 
    pop_sub(propagator_dt.exponential_midpoint)
  end subroutine exponential_midpoint
 
  ! ---------------------------------------------------------
  subroutine exponential_midpoint_predictor(hm, ks, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    type(hamiltonian_elec_t), target, intent(inout) :: hm
    type(v_ks_t),                     intent(inout) :: ks
    type(namespace_t),                intent(in)    :: namespace
    type(electron_space_t),           intent(in)    :: space
    type(grid_t),             target, intent(inout) :: gr
    type(states_elec_t),      target, intent(inout) :: st
    type(propagator_base_t),  target, intent(inout) :: tr
    real(real64),                     intent(in)    :: time
    real(real64),                     intent(in)    :: dt
    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
 
    type(gauge_field_t), pointer  :: gfield
    type(states_elec_t) :: st0
 
    push_sub(propagator_dt.exponential_midpoint_predictor)
 
    if (hm%theory_level /= independent_particles) then
      ! store copy of initial states
      call states_elec_copy(st0, st)
 
      ! get states, density and potential at time - dt/2
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
      ! compute potential, also with states being predicted at time - dt/2 to get orbital-dependent
      ! features correct (e.g. MGGA)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
 
      ! restore initial states to apply H(time - dt/2)
      call states_elec_group_copy(st0%d, st0%group, st%group)
      call states_elec_end(st0)
    end if
 
    ! move the ions to time time - dt/2
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time - m_half*dt, m_half*dt, save_pos = .true.)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, &
        m_half*dt, time, save_gf = .true.)
    end if
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt*m_half)
 
    ! apply Hamiltonian with dt
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt)
 
    !restore to time 'time - dt'
    call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions)
 
    if(associated(gfield)) then
      call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners)
    end if
 
    pop_sub(propagator_dt.exponential_midpoint_predictor)
  end subroutine exponential_midpoint_predictor
end module propagator_expmid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: