photon_mode_mf.F90

Go to the documentation of this file.

!! Copyright (C) 2017 Johannes Flick
!! Copyright (C) 2021 Davis Welakuh (merged)
!!
!! 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.
 
!!   v_mxc = - omega*q_alpha*(lambda*r) + (lambda*dipole)*(lambda*r)
!!   see PRL 115, 093001 (2015) and PRL 121, 113002 (2018)
#include "global.h"
 
module photon_mode_mf_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use io_function_oct_m
  use io_oct_m
  use ions_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use photon_mode_oct_m
  use profiling_oct_m
  use ps_oct_m
  use restart_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
 
 
  implicit none
 
  private
 
  public ::                &
    mf_t,                  &
    mf_init,               &
    mf_end,                &
    mf_calc,               &
    mf_photons_dump,       &
    mf_photons_load
  type mf_t
    real(real64),   pointer     :: vmf(:)
    real(real64),   allocatable :: dipole(:)
    real(real64),   allocatable :: dipole_former(:)
    complex(real64),   allocatable :: integral(:)
    real(real64),   allocatable :: pt_q(:), pt_p(:)
    real(real64),   allocatable :: pt_q_former(:)
    real(real64)         :: time_former
    real(real64),   pointer     :: fmf(:)  !! meanfield force
    logical              :: has_restart
 
  end type mf_t
 
contains
 
  subroutine mf_init(this, gr, st, ions, pt_mode)
    type(mf_t),          intent(out)  :: this
    type(grid_t),        intent(in)   :: gr
    type(states_elec_t), intent(in)   :: st
    type(ions_t),        intent(in)   :: ions
    type(photon_mode_t), intent(in)   :: pt_mode
 
    integer :: ions_dim
 
    push_sub(mf_init)
 
    ions_dim = gr%box%dim
 
    safe_allocate(this%vmf(1:gr%np))
    safe_allocate(this%dipole(1:ions_dim))
    safe_allocate(this%dipole_former(1:ions_dim))
 
    safe_allocate(this%integral(1:pt_mode%nmodes))
    safe_allocate(this%pt_q(1:pt_mode%nmodes))
    safe_allocate(this%pt_p(1:pt_mode%nmodes))
    safe_allocate(this%pt_q_former(1:pt_mode%nmodes))
    safe_allocate(this%fmf(1:ions_dim))
 
    this%vmf = m_zero
    this%has_restart = .false.
 
 
    this%dipole = - ions%dipole() - st%dipole(gr)
 
    this%dipole_former = m_zero
    this%integral = m_zero
 
    if (pt_mode%has_q0_p0) then
      this%pt_q(1:pt_mode%nmodes) = pt_mode%pt_coord_q0(1:pt_mode%nmodes)
      this%pt_p(1:pt_mode%nmodes) = pt_mode%pt_momen_p0(1:pt_mode%nmodes)
    else
      this%pt_q = m_zero
      this%pt_p = m_zero
    end if
 
    this%pt_q_former = m_zero
    this%time_former = m_zero
 
    this%fmf(1:ions_dim) = m_zero
 
    pop_sub(mf_init)
  end subroutine mf_init
 
!------------------------------------------
 
  subroutine mf_end(this)
    type(mf_t), intent(inout) :: this
 
    push_sub(mf_end)
 
    safe_deallocate_p(this%vmf)
    safe_deallocate_a(this%dipole_former)
    safe_deallocate_a(this%dipole)
 
    safe_deallocate_a(this%integral)
    safe_deallocate_a(this%pt_p)
    safe_deallocate_a(this%pt_q)
    safe_deallocate_a(this%pt_q_former)
 
    safe_deallocate_p(this%fmf)
 
    pop_sub(mf_end)
  end subroutine mf_end
 
!------------------------------------------
 
  subroutine mf_calc(this, gr, st, ions, pt_mode, time)
    type(mf_t),          intent(inout)    :: this
    type(grid_t),        intent(inout)    :: gr
    type(states_elec_t), intent(inout)    :: st
    type(ions_t),        intent(in)       :: ions
    type(photon_mode_t), intent(inout)    :: pt_mode
    real(real64),        intent(in)       :: time
 
    real(real64) :: lambda_pol_dipole, lambda_pol_dipole_former
    complex(real64) :: integrand
    real(real64) :: q0, p0
    integer :: ii, jj, ions_dim
    logical, save :: first = .true.
 
    push_sub(mf_calc)
 
    ions_dim = gr%box%dim
 
    if (.not. (first .and. this%has_restart)) then
      do ii = 1, pt_mode%nmodes
        this%pt_q_former(ii) = this%pt_q(ii)
      end do
 
      this%dipole_former = this%dipole
      this%dipole = - ions%dipole() - st%dipole(gr)
 
    end if
 
    this%vmf(1:gr%np) = m_zero
    this%fmf(1:ions_dim) = m_zero
 
    do ii = 1, pt_mode%nmodes
      lambda_pol_dipole = m_zero
      lambda_pol_dipole_former = m_zero
      do jj = 1, ions_dim
        lambda_pol_dipole = lambda_pol_dipole +  &
          pt_mode%lambda(ii)*pt_mode%pol(ii, jj)*this%dipole(jj)
        lambda_pol_dipole_former = lambda_pol_dipole_former + &
          pt_mode%lambda(ii)*pt_mode%pol(ii, jj)*this%dipole_former(jj)
      end do
 
      if (.not.(first .and. this%has_restart)) then
        if (abs(time) <= m_epsilon) then
          this%integral(ii) = m_zero
        else if (abs(this%integral(ii)) <= m_epsilon) then
          this%integral(ii) = m_half*lambda_pol_dipole*exp(-m_zi*pt_mode%omega(ii)*(this%time_former))
        else
          this%integral(ii) = this%integral(ii) + lambda_pol_dipole*exp(-m_zi*pt_mode%omega(ii)*(this%time_former))
        end if
      end if
 
      integrand = -(this%integral(ii)*exp(m_zi*pt_mode%omega(ii)*(time)*pt_mode%mu))* &
        (time*pt_mode%mu - this%time_former)
 
      this%pt_q(ii) = aimag(integrand)
      this%pt_p(ii) = pt_mode%omega(ii)*real(integrand)
 
      if (pt_mode%has_q0_p0) then
        q0 = pt_mode%pt_coord_q0(ii)*cos(pt_mode%omega(ii)*(time)*pt_mode%mu)
        q0 = q0 + pt_mode%pt_momen_p0(ii)/pt_mode%omega(ii)*sin(pt_mode%omega(ii)*(time)*pt_mode%mu)
        this%pt_q(ii) = this%pt_q(ii) + q0
 
        p0 = -pt_mode%pt_coord_q0(ii)*pt_mode%omega(ii)*sin(pt_mode%omega(ii)*(time)*pt_mode%mu)
        p0 = p0 + pt_mode%pt_momen_p0(ii)*cos(pt_mode%omega(ii)*(time)*pt_mode%mu)
        this%pt_p(ii) = this%pt_p(ii) + p0
      end if
 
      if (.not.allocated(pt_mode%pol_dipole)) then
        call photon_mode_compute_dipoles(pt_mode, gr)
      end if
 
      ! we need the negative sign due to the electric pol_dipole
      this%vmf(1:gr%np) = this%vmf(1:gr%np) - &
        m_half*((lambda_pol_dipole + lambda_pol_dipole_former) + pt_mode%omega(ii)* &
        (this%pt_q(ii) + this%pt_q_former(ii)))*(pt_mode%lambda(ii)*pt_mode%pol_dipole(1:gr%np,ii))
      do jj = 1, ions_dim
        this%fmf(jj) = this%fmf(jj) - pt_mode%omega(ii)*pt_mode%lambda(ii)* &
          pt_mode%pol(ii, jj)*(this%pt_q(ii) + lambda_pol_dipole/pt_mode%omega(ii)) !minus?
      end do
    end do
 
    if (first .and. this%has_restart) then
      first = .false.
    end if
 
    this%time_former = time*pt_mode%mu
 
    pop_sub(mf_calc)
  end subroutine mf_calc
 
! ---------------------------------------------------------
 
  subroutine mf_photons_dump(restart, this, gr, dt, pt_mode, ierr)
    type(restart_t), intent(in)  :: restart
    type(mf_t),      intent(in)  :: this
    type(grid_t),        intent(in)    :: gr
    real(real64),        intent(in)    :: dt
    type(photon_mode_t), intent(in)    :: pt_mode
    integer,         intent(out) :: ierr
 
    character(len=80), allocatable :: lines(:)
    integer :: iunit, err, jj, ions_dim
 
    push_sub(mf_photons_dump)
    ions_dim = gr%box%dim
 
    safe_allocate(lines(1:2*ions_dim + 4))
 
    ierr = 0
 
    iunit = restart%open('photon_mf')
    write(lines(1), '(a10,2x,es19.12)') 'pt_integral_real', real(this%integral(1))
    write(lines(2), '(a10,2x,es19.12)') 'pt_integral_aimag', aimag(this%integral(1))
    write(lines(3), '(a10,2x,es19.12)') 'pt_q_former', this%pt_q_former(1)
    write(lines(4), '(a10,2x,es19.12)') 'pt_time_former', this%time_former - dt*pt_mode%mu
    do jj = 1, ions_dim
      write(lines(4 + jj), '(a10,2x,es19.12)') 'dipole', this%dipole(jj)
    end do
    do jj = 1, ions_dim
      write(lines(4 + ions_dim + jj), '(a10,2x,es19.12)') 'dipole_former', this%dipole_former(jj)
    end do
    call restart%write(iunit, lines, 2*ions_dim + 4, err)
    if (err /= 0) ierr = ierr + 1
    call restart%close(iunit)
 
    safe_deallocate_a(lines)
 
    pop_sub(mf_photons_dump)
  end subroutine mf_photons_dump
 
! ---------------------------------------------------------
 
  subroutine mf_photons_load(restart, this, gr, ierr)
    type(restart_t),   intent(in)    :: restart
    type(mf_t),        intent(inout) :: this
    type(grid_t),      intent(in)    :: gr
    integer,           intent(out)   :: ierr
 
    integer :: err, iunit, jj, ions_dim
    character(len=128), allocatable :: lines(:)
    character(len=7) :: dummy
    real(real64), allocatable :: rr(:)
 
    push_sub(mf_photons_load)
 
    ierr = 0
    ions_dim = gr%box%dim
 
    if (restart%skip()) then
      ierr = -1
      pop_sub(mf_photons_load)
      return
    end if
 
    message(1) = "Debug: Reading Photons restart."
    call messages_info(1, namespace=restart%namespace, debug_only=.true.)
 
    safe_allocate(rr(1:2*ions_dim + 4))
    safe_allocate(lines(1:2*ions_dim + 4))
    iunit = restart%open('photon_mf')
    call restart%read(iunit, lines, 2*ions_dim + 4, err)
    if (err /= 0) then
      ierr = ierr + 1
    else
      do jj = 1, 2*ions_dim + 4
        read(lines(jj),'(a10,2x,es19.12)') dummy, rr(jj)
      end do
 
      this%integral(1) = cmplx(rr(1), rr(2), real64)
      this%pt_q_former(1) = rr(3)
      this%time_former = rr(4)
      do jj = 1, ions_dim
        this%dipole(jj) = rr(jj + 4)
      end do
      do jj = 1, ions_dim
        this%dipole_former(jj) = rr(jj + ions_dim + 4)
      end do
      this%has_restart = .true.
 
    end if
    call restart%close(iunit)
 
    message(1) = "Debug: Reading Photons restart done."
    call messages_info(1, namespace=restart%namespace, debug_only=.true.)
 
    safe_deallocate_a(rr)
    safe_deallocate_a(lines)
 
    pop_sub(mf_photons_load)
  end subroutine mf_photons_load
 
! ---------------------------------------------------------
 
end module photon_mode_mf_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: