pcm_potential.F90

Go to the documentation of this file.

!! Copyright (C) 2014 Alain Delgado Gran, Carlo Andrea Rozzi, Stefano Corni, Gabriel Gil
!!
!! 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 pcm_potential_oct_m
  use box_minimum_oct_m
  use debug_oct_m
  use ext_partner_list_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use index_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use kick_oct_m
  use lasers_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_interpolation_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use par_vec_oct_m
  use pcm_oct_m
  use pcm_eom_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use species_oct_m
  use space_oct_m
  use varinfo_oct_m
 
  ! to output debug info
  use unit_oct_m
  use unit_system_oct_m
  use io_function_oct_m
  use mesh_function_oct_m
 
  implicit none
 
  private
 
  public ::                 &
    pcm_hartree_potential
 
contains
 
  ! -----------------------------------------------------------------------------
  subroutine pcm_hartree_potential(pcm, space, mesh, psolver, ext_partners, vhartree, density, pcm_corr, &
    kick, time)
    type(pcm_t),         intent(inout) :: pcm
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(poisson_t),     intent(inout) :: psolver
    type(partner_list_t),intent(in)    :: ext_partners
    real(real64),        intent(in)    :: vhartree(:)
    real(real64),        intent(in)    :: density(:)
    real(real64),        intent(out)   :: pcm_corr
    type(kick_t), optional, intent(in) :: kick
    real(real64), optional,     intent(in)    :: time
 
    real(real64), allocatable :: potx(:)
    complex(real64), allocatable :: kick_eval(:)
    real(real64), allocatable :: kick_real(:)
    integer :: ii
 
    logical :: kick_time
    type(lasers_t), pointer :: lasers
 
    push_sub(pcm_hartree_potential)
 
    if (.not. pcm%run_pcm .or. .not. pcm_update(pcm)) then
      pcm_corr = m_zero
      pop_sub(pcm_hartree_potential)
      return
    end if
 
    if (pcm%solute) then
      call pcm_calc_pot_rs(pcm, mesh, psolver, v_h = vhartree, time_present = present(time))
    end if
 
    if (pcm%localf .and. present(time)) then
      lasers => list_get_lasers(ext_partners)
      if (associated(lasers) .and. present(kick)) then 
        safe_allocate(potx(1:mesh%np_part))
        safe_allocate(kick_eval(1:mesh%np_part))
        safe_allocate(kick_real(1:mesh%np_part))
        potx = m_zero
        kick_eval = m_zero
        do ii = 1, lasers%no_lasers
          call laser_potential(lasers%lasers(ii), mesh, potx, time)
        end do
        kick_real = m_zero
        kick_time = ((pcm%iter-1)*pcm%dt <= kick%time) .and. (pcm%iter*pcm%dt > kick%time)
        if (kick_time) then
          call kick_function_get(space, mesh, kick, kick_eval, 1, to_interpolate = .true.)
          kick_eval = kick%delta_strength * kick_eval
          kick_real = real(kick_eval, real64)
        end if
        call pcm_calc_pot_rs(pcm, mesh, psolver, v_ext = potx, kick = -kick_real, &
          time_present = present(time), kick_time = kick_time)
        safe_deallocate_a(potx)
        safe_deallocate_a(kick_eval)
        safe_deallocate_a(kick_real)
      else if (associated(lasers) .and. .not.present(kick)) then 
        safe_allocate(potx(1:mesh%np_part))
        potx = m_zero
        do ii = 1, lasers%no_lasers
          call laser_potential(lasers%lasers(ii), mesh, potx, time)
        end do
        call pcm_calc_pot_rs(pcm, mesh, psolver, v_ext = potx, time_present = present(time))
        safe_deallocate_a(potx)
      else if (.not.associated(lasers) .and. present(kick)) then 
        safe_allocate(kick_eval(1:mesh%np_part))
        safe_allocate(kick_real(1:mesh%np_part))
        kick_eval = m_zero
        kick_real = m_zero
        kick_time =((pcm%iter-1)*pcm%dt <= kick%time) .and. (pcm%iter*pcm%dt > kick%time)
        if (kick_time) then
          call kick_function_get(space, mesh, kick, kick_eval, 1, to_interpolate = .true.)
          kick_eval = kick%delta_strength * kick_eval
          kick_real = real(kick_eval, real64)
        end if
        call pcm_calc_pot_rs(pcm, mesh, psolver, kick = -kick_real, &
          time_present = present(time), kick_time = kick_time)
        safe_deallocate_a(kick_eval)
        safe_deallocate_a(kick_real)
      end if
 
      ! Calculating the PCM term renormalizing the sum of the single-particle energies
      ! to keep the idea of pcm_corr... but it will be added later on
      pcm_corr = dmf_dotp( mesh, density, pcm%v_e_rs + pcm%v_n_rs + pcm%v_ext_rs)
    else
      ! Calculating the PCM term renormalizing the sum of the single-particle energies
      pcm_corr = dmf_dotp( mesh, density, pcm%v_e_rs + pcm%v_n_rs)
    end if
 
    call messages_write(' PCM potential updated')
    call messages_new_line()
    call messages_write(' PCM update iteration counter: ')
    call messages_write(pcm%iter)
    call messages_info(debug_only=.true.)
 
    pop_sub(pcm_hartree_potential)
 
  end subroutine pcm_hartree_potential
 
end module pcm_potential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: