pcm_eom.F90

Go to the documentation of this file.

!! Copyright (C) 2017 Gabriel Gil, Stefano Corni, Silvio Pipolo, Carlo Andrea Rozzi
!!
!! 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_eom_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  implicit none
 
  private
  public ::                  &
    pcm_charges_propagation, &
    pcm_eom_enough_initial,  &
    pcm_eom_end,             &
    pcm_tessera_t,           &
    debye_param_t,           &
    drude_param_t,           &
    pcm_debye_model,         &
    pcm_drude_model,         &
    pcm_nuclei,              &
    pcm_electrons,           &
    pcm_external_potential,  &
    pcm_external_plus_kick,  &
    pcm_kick
  save
 
  type :: pcm_tessera_t
    real(real64) :: point(1:3)
    real(real64) :: area
    real(real64) :: normal(1:3)
    real(real64) :: r_sphere
  end type pcm_tessera_t
 
  type(pcm_tessera_t), allocatable :: cts_act(:)
  integer :: nts_act
 
  type :: debye_param_t
    real(real64) :: eps_0
    real(real64) :: eps_d
    real(real64) :: tau
  end type debye_param_t
 
  type :: drude_param_t
    real(real64) :: aa
    real(real64) :: gm
    real(real64) :: w0
  end type drude_param_t
 
  integer, parameter ::  &
    PCM_DEBYE_MODEL = 1, &
    pcm_drude_model = 2
 
  integer, parameter ::         &
    PCM_ELECTRONS          = 0, &
    pcm_nuclei             = 1, &
    pcm_external_potential = 2, &
    pcm_kick               = 3, &
    pcm_external_plus_kick = 4
 
 
  type(debye_param_t) :: deb
  type(drude_param_t) :: drl
 
 
  integer :: which_eom
 
  integer :: which_eps
 
  real(real64) :: dt
 
  !  polarization charges and variation of it in the previous iteration:
  real(real64), allocatable :: q_tp(:), dq_tp(:)
  real(real64), allocatable :: qext_tp(:), dqext_tp(:)
  real(real64), allocatable :: qkick_tp(:), dqkick_tp(:)
 
  real(real64), allocatable :: pot_tp(:)
  real(real64), allocatable :: potext_tp(:)
 
  !  See Chem.Phys.Lett. 429 (2006) 310-316 for Velocity-Verlet (VV) algorithm...
  real(real64) :: f1, f2, f3, f4, f5
  !  analogous to force in the equation of motion for the pol.charges, prev. iter.
  real(real64), allocatable :: force_tp(:)
  real(real64), allocatable :: force_qext_tp(:)
  real(real64), allocatable :: force_qkick_tp(:)
 
  ! In Ref.1 - J.Phys.Chem.A 2015, 119, 5405-5416... - Debye dielectric func. (eps)
  ! In Ref.2 - In J. Phys. Chem. C 2016, 120, 28774−28781... - Drude-Lorentz eps
  real(real64), allocatable :: cals(:,:), cald(:,:)
  real(real64), allocatable :: eigv(:), eigt(:,:)
  real(real64), allocatable :: sm12(:,:), sp12(:,:)
  real(real64), allocatable :: matq0(:,:), matqd(:,:)
  real(real64), allocatable :: matq0_lf(:,:), matqd_lf(:,:)
  !                                               !! \f$ \tilde{Q} \f$ and R matrices from Eq.(38)-(39) (Ref.1), respectively
  !                                               !! Q_f and Q_{\omega} from Eq.(17)-(16) (Ref.2), respectively
  real(real64), allocatable :: matqv(:,:), matqq(:,:)
  real(real64), allocatable :: matqv_lf(:,:)
  !                                               !! Q^{IEF(d)}_d, \f$ \tilde{Q} \f$ and R matrices enter the EOM in eq.(37), Ref.1
  !                                               !! Q_f and Q_{\omega} matrices enter the EOM in eq.(15), Ref.2
  !                                               !! N.B.: matrices R (in case of Debye) and Q_{\omega} (in case of Drude-Lorentz),
  !                                               !! i.e., matqq in this implementation,
  !                                               !! are the same in the EOMs for polarization charges due to the solute or ext.pot.
 
  logical :: enough_initial = .false.
 
contains
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine pcm_charges_propagation(q_t, pot_t, this_dt, this_cts_act, input_asc, this_eom, &
    this_eps, namespace, this_deb, this_drl)
    save
    real(real64),                  intent(out) :: q_t(:)
    real(real64),                  intent(in)  :: pot_t(:)
    real(real64),                  intent(in)  :: this_dt
    type(pcm_tessera_t),           intent(in)  :: this_cts_act(:)
    logical,                       intent(in)  :: input_asc
    integer,                       intent(in)  :: this_eom
    !                                                      !! or due to external potential ('external')
    integer,                       intent(in)  :: this_eps
    !                                                      !!(PCM_DEBYE_MODEL) or Drude-Lorentz (PCM_DRUDE_MODEL)
    type(namespace_t),             intent(in)  :: namespace
    type(debye_param_t), optional, intent(in)  :: this_deb
    type(drude_param_t), optional, intent(in)  :: this_drl
 
    logical :: firsttime = .true.
    logical :: initial_electron = .true.
    logical :: initial_external = .true.
    logical :: initial_kick = .true.
 
    push_sub(pcm_charges_propagation)
 
    which_eom = this_eom
    if (which_eom /= pcm_electrons .and. which_eom /= pcm_external_potential .and. &
      which_eom /= pcm_external_plus_kick .and. which_eom /= pcm_kick) then
      message(1) = "pcm_charges_propagation: EOM evolution of PCM charges can only be due to solute electrons"
      message(2) = "or external potential (including the kick)."
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (firsttime) then
      dt = this_dt
      nts_act = size(this_cts_act)
      if (size(q_t) /= nts_act) then
        message(1) = "pcm_charges_propagation: Number of tesserae do not coincide with size of PCM charges array."
        call messages_fatal(1, namespace=namespace)
      end if
 
      safe_allocate(cts_act(1:nts_act))
      cts_act = this_cts_act
 
      which_eps = this_eps
      select case (which_eps)
      case (pcm_debye_model)
        if (present(this_deb)) then
          deb = this_deb
        else
          message(1) = "pcm_charges_propagation: EOM-PCM error. Debye dielectric function requires three parameters."
          call messages_fatal(1, namespace=namespace)
        end if
      case (pcm_drude_model)
        if (present(this_drl)) then
          drl = this_drl
        else
          message(1) = "pcm_charges_propagation: EOM-PCM error. Drude-Lorentz dielectric function requires three parameters."
          call messages_fatal(1, namespace=namespace)
        end if
      case default
        message(1) = "pcm_charges_propagation: EOM-PCM error. Only Debye or Drude-Lorent dielectric models are allowed."
        call messages_fatal(1, namespace=namespace)
      end select
 
      if (abs(deb%tau) <= m_epsilon) then
        message(1) = "pcm_charges_propagation: EOM-PCM error. Debye EOM-PCM require a non-null Debye relaxation time."
        call messages_fatal(1, namespace=namespace)
      end if
      firsttime = .false.
    end if
 
 
    if (input_asc) then
      select case (which_eps)
      case (pcm_debye_model)
        ! initialize pcm charges due to electrons, external potential or kick
        call pcm_charges_from_input_file(q_t, pot_t, namespace)
 
        pop_sub(pcm_charges_propagation)
        return
      case default
        message(1) = "pcm_charges_propagation: EOM-PCM error. Only Debye EOM-PCM can startup from input charges."
        call messages_fatal(1, namespace=namespace)
      end select
    end if
 
    if ((initial_electron .and. which_eom == pcm_electrons) .or. &
      (initial_external .and. (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick)) .or. &
      (initial_kick     .and. which_eom == pcm_kick)) then
 
      ! initialize pcm matrices
      call pcm_bem_init(namespace)
 
      ! initialize pcm charges due to electrons, external potential or kick
      call init_charges(q_t, pot_t, namespace)
 
      if (initial_electron .and. which_eom == pcm_electrons) initial_electron = .false.
      if (initial_external .and. (which_eom == pcm_external_potential .or. &
        which_eom == pcm_external_plus_kick)) initial_external = .false.
      if (initial_kick     .and. which_eom == pcm_kick) initial_kick = .false.
 
    else
 
      ! propagate pcm charges due to electrons or external potential (including possible kick)
      if (which_eps == pcm_debye_model) then
        call pcm_ief_prop_deb(q_t, pot_t)
      else if (which_eps == pcm_drude_model) then
        call pcm_ief_prop_vv_ief_drl(q_t, pot_t)
      end if
 
    end if
 
    pop_sub(pcm_charges_propagation)
  end subroutine pcm_charges_propagation
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine pcm_charges_from_input_file(q_t, pot_t, namespace)
    real(real64), intent(out) :: q_t(:)
    real(real64), intent(in)  :: pot_t(:)
    type(namespace_t), intent(in) :: namespace
 
    real(real64) :: aux1(3)
    integer :: aux2
    integer :: asc_unit
    integer :: ia
 
    push_sub(pcm_charges_from_input_file)
 
    if (which_eom == pcm_electrons) then
      safe_allocate(q_tp(1:nts_act))
      asc_unit = io_open(pcm_dir//'ASC_e.dat', namespace, action='read')
      do ia = 1, nts_act
        read(asc_unit,*) aux1, q_t(ia), aux2
      end do
      q_tp = q_t
 
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick) then
      safe_allocate(qext_tp(1:nts_act))
      asc_unit = io_open(pcm_dir//'ASC_ext.dat', namespace, action='read')
      do ia = 1, nts_act
        read(asc_unit,*) aux1, q_t(ia), aux2
      end do
      qext_tp = q_t
 
    else if (which_eom == pcm_kick) then
      safe_allocate(qkick_tp(1:nts_act))
      asc_unit = io_open(pcm_dir//'ASC_kick.dat', namespace, action='read')
      do ia = 1, nts_act
        read(asc_unit,*) aux1, q_t(ia), aux2
      end do
      qkick_tp = q_t
    end if
    call io_close(asc_unit)
 
    safe_allocate(pot_tp(1:nts_act))
    pot_tp = pot_t
 
    pop_sub(pcm_charges_from_input_file)
  end subroutine pcm_charges_from_input_file
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine init_charges(q_t,pot_t, namespace)
    real(real64),      intent(out) :: q_t(:)
    real(real64),      intent(in)  :: pot_t(:)
    type(namespace_t), intent(in)  :: namespace
 
    push_sub(init_charges)
 
    if (which_eom == pcm_electrons) then
      ! Here we consider the potential at any earlier time equal to the initial potential.
      ! Therefore, we can suppose that the solvent is already in equilibrium with the initial solute potential.
      ! The latter is valid when starting the propagation from the ground state but does not hold in general.
      message(1) = 'EOM-PCM for solvent polarization due to solute electrons considers that you start from a ground state run.'
      call messages_warning(1, namespace=namespace)
 
      safe_allocate(pot_tp(1:nts_act))
      pot_tp = pot_t
 
      ! applying the static IEF-PCM response matrix (corresponging to epsilon_0) to the initial potential
      q_t = matmul(matq0, pot_t)
 
      safe_allocate(q_tp(1:nts_act))
      q_tp = q_t
 
      if (which_eps == pcm_drude_model) then
        safe_allocate(dq_tp(1:nts_act))
        safe_allocate(force_tp(1:nts_act))
        dq_tp = m_zero
        force_tp = m_zero
      end if
 
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick) then
      ! Here (instead) we consider zero the potential at any earlier time.
      ! Therefore, the solvent is not initially in equilibrium with the external potential unless its initial value is zero.
 
      safe_allocate(potext_tp(1:nts_act))
      potext_tp = pot_t
 
      ! applying the dynamic IEF-PCM response matrix (for epsilon_d) to the initial potential
      q_t = matmul(matqd_lf, pot_t)
 
      safe_allocate(qext_tp(1:nts_act))
      qext_tp = q_t
 
      if (which_eps == pcm_drude_model) then
        safe_allocate(dqext_tp(1:nts_act))
        safe_allocate(force_qext_tp(1:nts_act))
        dqext_tp = m_zero
        force_qext_tp = m_zero
      end if
 
    else if (which_eom == pcm_kick) then
 
      if (which_eps == pcm_drude_model) then
        safe_allocate(dqkick_tp(1:nts_act))
        safe_allocate(force_qkick_tp(1:nts_act))
        dqkick_tp = m_zero
        force_qkick_tp = m_zero
      end if
 
      if (which_eps == pcm_drude_model) then
        dqkick_tp = matmul(matqv_lf, pot_t)*dt
      else
        q_t = matmul(matqv_lf - matmul(matqq, matqd_lf), pot_t)
      end if
 
      safe_allocate(qkick_tp(1:nts_act))
      qkick_tp = q_t
 
    end if
 
    ! initializing Velocity-Verlet algorithm for the integration of EOM-PCM for Drude-Lorentz
    if (which_eps == pcm_drude_model) call init_vv_propagator
 
    pop_sub(init_charges)
  end subroutine init_charges
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine pcm_ief_prop_deb(q_t, pot_t)
    real(real64), intent(out) :: q_t(:)
    real(real64), intent(in)  :: pot_t(:)
 
    push_sub(pcm_ief_prop_deb)
 
    if (which_eom == pcm_electrons) then
      ! Eq.(47) in S. Corni, S. Pipolo and R. Cammi, J.Phys.Chem.A 2015, 119, 5405-5416.
      q_t(:) = q_tp(:) - dt*matmul(matqq, q_tp)   + &
        dt*matmul(matqv, pot_tp) + &
        matmul(matqd, pot_t - pot_tp)
 
      q_tp = q_t
 
      pot_tp = pot_t
 
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick) then
      ! analogous to Eq.(47) ibid. with different matrices
      q_t(:) = qext_tp(:) - dt*matmul(matqq, qext_tp)   + &
        dt*matmul(matqv_lf, potext_tp) + &
        matmul(matqd_lf, pot_t - potext_tp)             ! N.B. matqq
 
      qext_tp = q_t
 
      potext_tp = pot_t
 
    else if (which_eom == pcm_kick) then
      ! simplifying for kick
      q_t(:) = qkick_tp(:) - dt*matmul(matqq, qkick_tp)    ! N.B. matqq
 
      qkick_tp = q_t
 
    end if
 
    pop_sub(pcm_ief_prop_deb)
  end subroutine pcm_ief_prop_deb
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine init_vv_propagator
    push_sub(init_vv_propagator)
 
    f1 = dt*(m_one - dt*m_half*drl%gm)
    f2 = dt*dt*m_half
    f3 = m_one - dt*drl%gm*(m_one - dt*m_half*drl%gm)
    f4 = m_half*dt
    f5 = drl%gm*f2
 
    pop_sub(init_vv_propagator)
  end subroutine init_vv_propagator
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine pcm_ief_prop_vv_ief_drl(q_t, pot_t)
    real(real64), intent(out) :: q_t(:)
    real(real64), intent(in)  :: pot_t(:)
 
    real(real64) :: force(nts_act)
 
    push_sub(pcm_ief_prop_vv_ief_drl)
 
    if (which_eom == pcm_electrons) then
      ! From Eq.(15) S. Pipolo and S. Corni, J.Phys.Chem.C 2016, 120, 28774-28781.
      ! Using the scheme in Eq.(21) and (17) of E. Vanden-Eijnden, G. Ciccotti, Chem.Phys.Lett. 429 (2006) 310-316.
      q_t = q_tp + f1*dq_tp + f2*force_tp
      force = -matmul(matqq, q_t) + matmul(matqv, pot_t)
      dq_tp = f3*dq_tp + f4*(force+force_tp) -f5*force_tp
      force_tp = force
      q_tp = q_t
 
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick) then
      ! analagous to Eq.(15) ibid. with different matrices
      q_t = qext_tp + f1*dqext_tp + f2*force_qext_tp
      force = -matmul(matqq, q_t) + matmul(matqv_lf, pot_t) ! N.B. matqq
      dqext_tp = f3*dqext_tp + f4*(force + force_qext_tp) -f5*force_qext_tp
      force_qext_tp = force
      q_tp = q_t
 
    else if (which_eom == pcm_kick) then
      ! simplifying for kick
      q_t = qkick_tp + f1*dqkick_tp + f2*force_qkick_tp
      force = -matmul(matqq, q_t)                         ! N.B. matqq
      dqkick_tp = f3*dqkick_tp + f4*(force + force_qkick_tp) -f5*force_qkick_tp
      force_qkick_tp = force
      q_tp = q_t
 
    end if
 
    pot_tp = pot_t
 
    pop_sub(pcm_ief_prop_vv_ief_drl)
  end subroutine pcm_ief_prop_vv_ief_drl
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine pcm_bem_init(namespace)
    type(namespace_t), intent(in) :: namespace
    integer :: itess, jtess
    integer :: pcmmat0_unit, pcmmatd_unit
 
    push_sub(pcm_bem_init)
 
    if (which_eom == pcm_electrons) then
      safe_allocate(matq0(1:nts_act, 1:nts_act))
      safe_allocate(matqd(1:nts_act, 1:nts_act))
      safe_allocate(matqv(1:nts_act, 1:nts_act))
      if (.not. allocated(matqq)) then
        safe_allocate(matqq(1:nts_act, 1:nts_act))
      end if
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick .or. which_eom == pcm_kick) then
      safe_allocate(matq0_lf(1:nts_act, 1:nts_act)) ! not used yet
      safe_allocate(matqd_lf(1:nts_act, 1:nts_act))
      safe_allocate(matqv_lf(1:nts_act, 1:nts_act))
      if (.not. allocated(matqq)) then
        safe_allocate(matqq(1:nts_act, 1:nts_act))
      end if
    end if
    call do_pcm_propmat(namespace)
 
    if (which_eom == pcm_electrons) then
      pcmmat0_unit = io_open(pcm_dir//'pcm_matrix_static_from_eom.out', namespace, action='write')
      pcmmatd_unit = io_open(pcm_dir//'pcm_matrix_dynamic_from_eom.out', namespace, action='write')
      do jtess = 1, nts_act
        do itess = 1, nts_act
          write(pcmmat0_unit,*) matq0(itess, jtess)
          write(pcmmatd_unit,*) matqd(itess, jtess)
        end do
      end do
      call io_close(pcmmat0_unit)
      call io_close(pcmmatd_unit)
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick) then
      pcmmat0_unit = io_open(pcm_dir//'pcm_matrix_static_lf_from_eom.out', namespace, action='write')
      pcmmatd_unit = io_open(pcm_dir//'pcm_matrix_dynamic_lf_from_eom.out', namespace, action='write')
      do jtess = 1, nts_act
        do itess = 1, nts_act
          write(pcmmat0_unit,*) matq0_lf(itess, jtess)
          write(pcmmatd_unit,*) matqd_lf(itess, jtess)
        end do
      end do
      call io_close(pcmmat0_unit)
      call io_close(pcmmatd_unit)
    end if
 
    pop_sub(pcm_bem_init)
  end subroutine pcm_bem_init
 
  !------------------------------------------------------------------------------------------------------------------------------
 
  subroutine pcm_eom_end()
    push_sub(pcm_eom_end)
 
    ! pcm charges
    safe_deallocate_a(q_tp)
    safe_deallocate_a(qext_tp)
    safe_deallocate_a(qkick_tp)
 
    ! increment pcm charges
    safe_deallocate_a(dq_tp)
    safe_deallocate_a(dqext_tp)
    safe_deallocate_a(dqkick_tp)
 
    ! force-like term in pcm-eom equation
    safe_deallocate_a(force_tp)
    safe_deallocate_a(force_qext_tp)
    safe_deallocate_a(force_qkick_tp)
 
    ! pcm-eom bem matrices
    safe_deallocate_a(matq0)
    safe_deallocate_a(matqd)
    safe_deallocate_a(matq0_lf)
    safe_deallocate_a(matqd_lf)
    safe_deallocate_a(matqv)
    safe_deallocate_a(matqv_lf)
    safe_deallocate_a(matqq)
 
    safe_deallocate_a(pot_tp)
    safe_deallocate_a(cts_act)
 
    call deallocate_ts_matrix()
 
    pop_sub(pcm_eom_end)
  end subroutine pcm_eom_end
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine do_pcm_propmat(namespace)
    type(namespace_t), intent(in) :: namespace
    save
    integer :: i, j
    real(real64), allocatable :: scr4(:,:), scr1(:,:)
    real(real64), allocatable :: scr2(:,:), scr3(:,:)
    real(real64), allocatable :: fact1(:), fact2(:)
    real(real64), allocatable :: Kdiag0(:), Kdiagd(:)
    real(real64) :: sgn,sgn_lf, fac_eps0, fac_epsd
    real(real64) :: temp
 
    logical :: firsttime = .true.
 
    push_sub(do_pcm_propmat)
 
    sgn = m_one ! In the case of NP this is -1 because the normal to the cavity is always pointing outward
    !             and the field acreated by a positive unit charge outside the cavity is directed inward.
 
    sgn_lf = m_one
    ! ''local field'' differ from ''standard'' PCM response matrix in some sign changes
    if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick .or. which_eom == pcm_kick) sgn_lf = -m_one
 
    safe_allocate(cals(1:nts_act, 1:nts_act))
    safe_allocate(cald(1:nts_act, 1:nts_act))
    safe_allocate(kdiag0(1:nts_act))
    safe_allocate(kdiagd(1:nts_act))
    safe_allocate(fact1(1:nts_act))
    safe_allocate(fact2(1:nts_act))
 
    ! generate Calderon S and D matrices
    do i = 1, nts_act
      do j = 1, nts_act
        call green_s(i, j, temp)
        cals(i, j) = temp
        call green_d(i, j, temp)
        cald(i, j) = temp
      end do
    end do
    if (firsttime) then
      call allocate_ts_matrix()
      call do_ts_matrix(namespace)
      firsttime = .false.
    end if
    safe_deallocate_a(cals)
    safe_deallocate_a(cald)
 
    safe_allocate(scr4(1:nts_act, 1:nts_act))
    safe_allocate(scr1(1:nts_act, 1:nts_act))
    safe_allocate(scr2(1:nts_act, 1:nts_act))
    safe_allocate(scr3(1:nts_act, 1:nts_act))
 
    if (which_eps == pcm_debye_model) then
      if (.not. is_close(deb%eps_0, m_one)) then
        fac_eps0 = (deb%eps_0 + m_one)/(deb%eps_0 - m_one)
        kdiag0(:) = sgn_lf*(m_two*m_pi - sgn*sgn_lf*eigv(:))/(m_two*m_pi*fac_eps0 - sgn*eigv(:)) ! Eq.(14) with eps_0 in Ref.1
      else
        kdiag0(:) = m_zero
      end if
      if (.not. is_close(deb%eps_d, m_one)) then
        fac_epsd = (deb%eps_d + m_one)/(deb%eps_d - m_one)
        kdiagd(:) = sgn_lf*(m_two*m_pi - sgn*sgn_lf*eigv(:))/(m_two*m_pi*fac_epsd - sgn*eigv(:)) ! Eq.(14) with eps_d, ibid.
      else
        kdiagd(:) = m_zero
      end if
      fact1(:) = ((m_two*m_pi - sgn*eigv(:))*deb%eps_0 + m_two*m_pi + eigv(:))/ &       ! inverse of Eq.(32), ibid.
        ((m_two*m_pi - sgn*eigv(:))*deb%eps_d + m_two*m_pi + eigv(:))/deb%tau
      fact2(:) = kdiag0(:)*fact1(:)                                 ! tau^{-1}K_0 in Eq.(38), ibid.
 
    else if (which_eps == pcm_drude_model) then
      kdiagd(:) = m_zero                                         ! from Eq.(10) up in Ref.2
      fact2(:) = (m_two*m_pi - sgn*eigv(:))*drl%aa/(m_four*m_pi) ! Eq.(10) down
      do i = 1, nts_act
        if (fact2(i) < m_zero)fact2(i) = m_zero ! check out
      end do
      if (abs(drl%w0) <= m_epsilon) drl%w0 = 1.0e-8_real64                       ! check out
      fact1(:) = fact2(:) + drl%w0*drl%w0                                       ! Eq.(19), ibid.
      fact2(:) = sgn_lf*(m_two*m_pi - sgn*sgn_lf*eigv(:))*drl%aa/(m_four*m_pi)  ! Eq.(10) down, local field analogous
      kdiag0(:) = fact2(:)/fact1(:)                                             ! from Eq.(10) up, ibid.
    end if
 
    scr3 = matmul(sm12, eigt)
    scr2 = matmul(transpose(eigt), sp12)
    scr4 = transpose(scr3)
    do i = 1, nts_act
      scr1(:,i) = scr3(:, i)*fact1(i)
    end do
    matqq = matmul(scr1, scr2) ! Eq.(39) in Ref.1 and Eq.(16) in Ref.2
    do i = 1,nts_act
      scr1(:,i) = scr3(:, i)*fact2(i)
    end do
    if (which_eom == pcm_electrons) then
      matqv = -matmul(scr1, scr4) ! Eq.(38) in Ref.1 and Eq.(17) in Ref.2
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick .or. which_eom == pcm_kick) then
      matqv_lf = -matmul(scr1, scr4) ! local field analogous
    end if
    do i = 1, nts_act
      scr1(:,i) = scr3(:,i)*kdiag0(i)
    end do
    if (which_eom == pcm_electrons) then
      matq0 = -matmul(scr1, scr4) ! from Eq.(14) and (18) for eps_0 in Ref.1
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick .or. which_eom == pcm_kick) then
      matq0_lf = -matmul(scr1, scr4) ! local field analogous (not used yet)
    end if
    do i = 1, nts_act
      scr1(:,i) = scr3(:, i)*kdiagd(i)
    end do
    if (which_eom == pcm_electrons) then
      matqd = -matmul(scr1, scr4) ! from Eq.(14) and (18) for eps_d, ibid.
    else if (which_eom == pcm_external_potential .or. which_eom == pcm_external_plus_kick .or. which_eom == pcm_kick) then
      matqd_lf = -matmul(scr1, scr4) ! local field analogous
    end if
 
    safe_deallocate_a(scr4)
    safe_deallocate_a(scr1)
    safe_deallocate_a(scr2)
    safe_deallocate_a(scr3)
    safe_deallocate_a(fact1)
    safe_deallocate_a(fact2)
    safe_deallocate_a(kdiag0)
    safe_deallocate_a(kdiagd)
 
    if (enough_initial) then
      call deallocate_ts_matrix()
    end if
 
    pop_sub(do_pcm_propmat)
  end subroutine do_pcm_propmat
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine green_d(i, j, value)
    integer, intent(in)  :: i, j
    real(real64),   intent(out) :: value
 
    real(real64) :: dist,diff(3)
 
    push_sub(green_d)
 
    if (i /= j) then
      diff = cts_act(i)%point - cts_act(j)%point
      dist = norm2(diff)
      value = dot_product(cts_act(j)%normal, diff)/dist**3 ! Eq.(5) in Refs.1-2
    else
      value = -1.0694_real64*sqrt(m_four*m_pi*cts_act(i)%area)
      value = value/(m_two*cts_act(i)%r_sphere)/cts_act(i)%area ! diagonal part is a bit cumbersome
    end if
 
    pop_sub(green_d)
  end subroutine green_d
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine green_s(i, j, value)
    integer, intent(in):: i,j
    real(real64),   intent(out) :: value
 
    real(real64) :: dist,diff(3)
 
    push_sub(green_s)
 
    if (i /= j) then
      diff = cts_act(i)%point - cts_act(j)%point
      dist = norm2(diff)
      value = m_one/dist ! Eq.(5) in Refs.1-2
    else
      value = 1.0694_real64*sqrt(m_four*m_pi/cts_act(i)%area) ! diagonal part is a bit cumbersome
    end if
 
    pop_sub(green_s)
  end subroutine green_s
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine allocate_ts_matrix()
    push_sub(allocate_ts_matrix)
 
    safe_allocate(eigv(1:nts_act))
    safe_allocate(eigt(1:nts_act, 1:nts_act))
    safe_allocate(sm12(1:nts_act, 1:nts_act))
    safe_allocate(sp12(1:nts_act, 1:nts_act))
 
    pop_sub(allocate_ts_matrix)
  end subroutine allocate_ts_matrix
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine deallocate_ts_matrix()
    push_sub(deallocate_ts_matrix)
 
    safe_deallocate_a(eigv)
    safe_deallocate_a(eigt)
    safe_deallocate_a(sm12)
    safe_deallocate_a(sp12)
 
    pop_sub(deallocate_ts_matrix)
  end subroutine deallocate_ts_matrix
 
  !------------------------------------------------------------------------------------------------------------------------------
  subroutine do_ts_matrix(namespace)
    type(namespace_t), intent(in) :: namespace
    integer :: i, j
    integer :: info, lwork, liwork
    real(real64), allocatable :: scr1(:,:), scr2(:,:), eigt_t(:,:)
    real(real64) :: sgn
    character jobz, uplo
    integer, allocatable :: iwork(:)
    real(real64), allocatable :: work(:)
 
    push_sub(do_ts_matrix)
 
    safe_allocate(scr1(1:nts_act, 1:nts_act))
    safe_allocate(scr2(1:nts_act, 1:nts_act))
    safe_allocate(eigt_t(1:nts_act, 1:nts_act))
    safe_allocate(work(1:1 + 6*nts_act + 2*nts_act*nts_act))
    safe_allocate(iwork(1:3 + 5*nts_act))
 
    sgn = m_one
 
    jobz = 'V'
    uplo = 'U'
    lwork = 1 + 6*nts_act + 2*nts_act*nts_act
    liwork = 3 + 5*nts_act
    eigt = cals
    call dsyevd(jobz, uplo, nts_act, eigt, nts_act, eigv, work, lwork, iwork, liwork, info)
    do i = 1, nts_act
      if (eigv(i) < m_zero) then
        write(message(1),*) "Eigenvalue ", i, " of S when constructing the TS matrix is negative!"
        write(message(2),*) "I put it to 1e-8"
        call messages_warning(2, namespace=namespace)
        eigv(i) = 1.0e-8_real64
      end if
      scr1(:,i) = eigt(:,i)*sqrt(eigv(i))
    end do
    eigt_t = transpose(eigt)
    sp12 = matmul(scr1, eigt_t)          ! building S^{1/2} to be used in R (Ref.1) and Q_{\omega} (Ref.2)
    do i = 1, nts_act
      scr1(:, i) = eigt(:, i)/sqrt(eigv(i))
    end do
    sm12 = matmul(scr1, eigt_t)          ! building S^{-1/2} to be used in R and \tilde{Q} (Ref.1) and Q_{\omega} and Q_f (Ref.2)
    do i = 1, nts_act
      scr1(:,i) = cald(:, i)*cts_act(i)%area
    end do
    scr2 = matmul(matmul(sm12, scr1), sp12)   ! Eq.(10) in Ref.1 (paragraph after Eq.(9), Ref.2)
    do j = 1, nts_act
      do i = 1, nts_act
        eigt(i, j) = m_half*(scr2(i, j) + scr2(j, i)) ! re-symmetrizing for numerical reasons
      end do
    end do
    call dsyevd(jobz,uplo,nts_act,eigt,nts_act,eigv,work,lwork, &
      iwork,liwork,info)        ! obtaining \Lambda (eigv) and T (eigt), Eq.(10), ibid.
 
    safe_deallocate_a(work)
    safe_deallocate_a(iwork)
    safe_deallocate_a(scr1)
    safe_deallocate_a(scr2)
    safe_deallocate_a(eigt_t)
 
    pop_sub(do_ts_matrix)
  end subroutine do_ts_matrix
 
  ! -----------------------------------------------------------------------------
  subroutine pcm_eom_enough_initial(not_yet_called)
    logical, intent(out) :: not_yet_called
 
    push_sub(pcm_eom_enough_initial)
 
    enough_initial = .true.
    not_yet_called = .false.
 
    pop_sub(pcm_eom_enough_initial)
  end subroutine pcm_eom_enough_initial
 
end module pcm_eom_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: