phonons_lr.F90

Go to the documentation of this file.

!! Copyright (C) 2007-2012 Xavier Andrade, David Strubbe
!!
!! 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 phonons_lr_oct_m
  use born_charges_oct_m
  use debug_oct_m
  use epot_oct_m
  use forces_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kdotp_oct_m
  use kdotp_calc_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use linear_response_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_basic_oct_m
  use namespace_oct_m
  use parser_oct_m
  use perturbation_ionic_oct_m
  use profiling_oct_m
  use restart_oct_m
  use smear_oct_m
  use space_oct_m
  use species_pot_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_restart_oct_m
  use sternheimer_oct_m
  use electrons_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use vibrations_oct_m
  use v_ks_oct_m
 
  implicit none
 
  private
  public :: &
    phonons_lr_run,    &
    phn_nm_wfs_tag,    &
    phn_wfs_tag,       &
    phn_rho_tag,       &
    axsf_mode_output
 
contains
 
  ! ---------------------------------------------------------
  subroutine phonons_lr_run(system, from_scratch)
    class(*),        intent(inout) :: system
    logical,         intent(in)    :: from_scratch
 
    push_sub(phonons_lr_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = vib_modes not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call phonons_lr_run_legacy(system, from_scratch)
    end select
 
    pop_sub(phonons_lr_run)
  end subroutine phonons_lr_run
 
  ! ---------------------------------------------------------
  subroutine phonons_lr_run_legacy(sys, fromscratch)
    type(electrons_t), target, intent(inout) :: sys
    logical,                   intent(in)    :: fromscratch
 
    type(sternheimer_t) :: sh
    type(lr_t)          :: lr(1:1), kdotp_lr(sys%space%dim)
    type(vibrations_t)  :: vib
    class(perturbation_ionic_t), pointer  :: pert
 
    type(ions_t),     pointer :: ions
    type(states_elec_t),   pointer :: st
    type(grid_t),     pointer :: gr
 
    integer :: natoms, ndim, iatom, idir, jatom, jdir, imat, jmat, iunit_restart, ierr, start_mode, ispin
    complex(real64), allocatable :: force_deriv(:,:)
    real(real64),    allocatable :: grad_rho_core(:,:), drho_core(:,:), dvxc_bare(:,:,:)
    real(real64),    allocatable :: ddrho_tot(:,:), ddvxc_resp(:,:,:)
    complex(real64), allocatable :: zdrho_tot(:,:), zdvxc_resp(:,:,:)
    character(len=80) :: str_tmp
    character(len=300) :: line(1)
    type(born_charges_t) :: born
    logical :: normal_mode_wfs, do_infrared, symmetrize
    type(restart_t) :: restart_load, restart_dump, kdotp_restart, gs_restart
 
    push_sub(phonons_lr_run_legacy)
 
    !some shortcuts
 
    ions => sys%ions
    st  => sys%st
    gr  => sys%gr
 
    if (sys%hm%pcm%run_pcm) then
      call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
    end if
 
    if (sys%space%is_periodic()) then
      call messages_not_implemented('linear-response vib_modes for periodic systems', namespace=sys%namespace)
    end if
 
    if (sys%kpoints%use_symmetries) then
      call messages_experimental("CalculationMode = vib_modes with k-points symmetries", namespace=sys%namespace)
    end if
 
    if (sys%kpoints%reduced%npoints /= sys%kpoints%full%npoints) then
      call messages_experimental("CalculationMode = vib_modes with reduced k-grid", namespace=sys%namespace)
    end if
 
    !%Variable CalcNormalModeWfs
    !%Type logical
    !%Default false
    !%Section Linear Response::Vibrational Modes
    !%Description
    !% If set to true, the response wavefunctions for each normal mode will be calculated
    !% and written in directory <tt>restart/vib_modes/phn_nm_wfs_XXXXX</tt>.
    !% This part is time-consuming and not parallel, but not needed for most purposes.
    !%End
    call parse_variable(sys%namespace, 'CalcNormalModeWfs', .false., normal_mode_wfs)
 
    !%Variable CalcInfrared
    !%Type logical
    !%Default true
    !%Section Linear Response::Vibrational Modes
    !%Description
    !% If set to true, infrared intensities (and born charges) will be calculated
    !% and written in <tt>vib_modes/infrared</tt>.
    !%End
    call parse_variable(sys%namespace, 'CalcInfrared', .true., do_infrared)
 
    !%Variable SymmetrizeDynamicalMatrix
    !%Type logical
    !%Default true
    !%Section Linear Response::Vibrational Modes
    !%Description
    !% If set to true, all entries of the dynamical matrix will be calculated and then
    !% the matrix will be symmetrized to enforce <math>D_{ij} = D_{ji}</math>. If set to false,
    !% only the upper half of the matrix will be calculated.
    !%End
    call parse_variable(sys%namespace, 'SymmetrizeDynamicalMatrix', .true., symmetrize)
 
    ! replaced by properly saving and reading the dynamical matrix
    call messages_obsolete_variable(sys%namespace, 'UseRestartDontSolve')
 
    natoms = ions%natoms
    ndim = sys%space%dim
 
    call gs_restart%init(sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=gr, exact=.true.)
    if (ierr == 0) then
      call states_elec_look_and_load(gs_restart, sys%namespace, sys%space, st, sys%gr, sys%kpoints, st%restart_fixed_occ)
      call gs_restart%end()
    else
      message(1) = "Previous gs calculation is required."
      call messages_fatal(1, namespace=sys%namespace)
    end if
 
    ! read kdotp wavefunctions if necessary (for IR intensities)
    if (sys%space%is_periodic() .and. do_infrared) then
      message(1) = "Reading kdotp wavefunctions for periodic directions."
      call messages_info(1, namespace=sys%namespace)
 
      call kdotp_restart%init(sys%namespace, restart_kdotp, restart_type_load, sys%mc, ierr, mesh=gr)
      if (ierr /= 0) then
        message(1) = "Unable to read kdotp wavefunctions."
        message(2) = "Previous kdotp calculation required."
        call messages_fatal(2, namespace=sys%namespace)
      end if
 
      do idir = 1, sys%space%periodic_dim
        call lr_init(kdotp_lr(idir))
        call lr_allocate(kdotp_lr(idir), sys%st, sys%gr)
 
        ! load wavefunctions
        str_tmp = trim(kdotp_wfs_tag(idir))
        call kdotp_restart%open_dir(wfs_tag_sigma(sys%namespace, str_tmp, 1), ierr)
        if (ierr == 0) then
          call states_elec_load(kdotp_restart, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
            sys%st%restart_fixed_occ, ierr=ierr, lr=kdotp_lr(idir))
        end if
        call kdotp_restart%close_dir()
 
        if (ierr /= 0) then
          message(1) = "Unable to read kdotp wavefunctions from '"//trim(wfs_tag_sigma(sys%namespace, str_tmp, 1))//"'."
          message(2) = "Previous kdotp calculation required."
          call messages_fatal(2, namespace=sys%namespace)
        end if
      end do
      call kdotp_restart%end()
    end if
 
    message(1) = 'Info: Setting up Hamiltonian for linear response.'
    call messages_info(1, namespace=sys%namespace)
 
    call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
    call sternheimer_init(sh, sys%namespace, sys%space, sys%gr, sys%st, sys%hm, sys%ks, sys%mc, &
      wfs_are_cplx = states_are_complex(st))
 
    if (sys%hm%ep%nlcc .and. .not. sh%add_fxc()) then
      message(1) = "Linear-response vib_modes with non-linear core corrections requires the xc kernel."
      message(2) = "HamiltonianVariation must include V_xc."
      call messages_fatal(2, namespace=sys%namespace)
    end if
 
    call vibrations_init(vib, ions%space, ions%natoms, ions%mass, "lr", sys%namespace)
 
    call epot_precalc_local_potential(sys%hm%ep, sys%namespace, sys%gr, sys%ions)
 
    if (do_infrared) then
      call born_charges_init(born, sys%namespace, ions%natoms, st%val_charge, st%qtot, ndim)
    end if
 
    call lr_init(lr(1))
    call lr_allocate(lr(1), st, gr)
 
    call restart_dump%init(sys%namespace, restart_vib_modes, restart_type_dump, sys%mc, ierr, mesh=gr)
    call restart_load%init(sys%namespace, restart_vib_modes, restart_type_load, sys%mc, ierr, mesh=gr)
 
    !CALCULATE
 
    ! the ionic contribution, see second term in Eq. 85 in Baroni et al. RMP, 73, 515 (2001)
    call build_ionic_dyn_matrix()
 
    ! Compute the  -<phi0 | v2 | phi0> term
    if (states_are_real(st)) then
      call dionic_pert_matrix_elements_2(sys%gr, sys%namespace, sys%space, sys%ions, sys%hm, 1, st, vib, vib%dyn_matrix)
    else
      call zionic_pert_matrix_elements_2(sys%gr, sys%namespace, sys%space, sys%ions, sys%hm, 1, st, vib, vib%dyn_matrix)
    end if
 
    if (fromscratch) then
      start_mode = 1
    else
      call phonons_load(restart_load, vib, start_mode)
    end if
 
    ! Delete, if fromScratch, or trying to open it failed and there is something wrong with it.
    if (start_mode == 1) call restart_dump%rm('restart')
 
    ! Output the first start_mode-1 modes
    do imat = 1, start_mode - 1
      call vibrations_out_dyn_matrix_row(vib, imat)
    end do
 
 
    ! The remaining term is <\psi_1 | v^{(1)} | \psi_0> + <\psi_0 | v^{(1)} | \psi_1>
    ! We first get \psi_1 from Sternheimer, and then compute the term using X(forces_derivative)
    ! Note that in order to avoid getting v^{(1)}, we perform an integration by part  in this routine
    pert => perturbation_ionic_t(sys%namespace, ions)
 
    if (sys%hm%ep%nlcc) then
      safe_allocate(grad_rho_core(1:gr%np, 1:ndim))
      safe_allocate(drho_core(1:gr%np, 1:st%d%nspin))
      safe_allocate(dvxc_bare(1:gr%np, 1:st%d%nspin, 1:1))
      if (states_are_real(st)) then
        safe_allocate(ddrho_tot(1:gr%np, 1:st%d%nspin))
        safe_allocate(ddvxc_resp(1:gr%np, 1:st%d%nspin, 1:1))
      else
        safe_allocate(zdrho_tot(1:gr%np, 1:st%d%nspin))
        safe_allocate(zdvxc_resp(1:gr%np, 1:st%d%nspin, 1:1))
      end if
    end if
 
    do imat = start_mode, vib%num_modes
      iatom = vibrations_get_atom(vib, imat)
      idir  = vibrations_get_dir(vib, imat)
 
      write(message(1),'(a,i5,a,a1,a)') &
        "Calculating response to displacement of atom ", iatom, " in ", index2axis(idir), "-direction."
      call messages_info(1, namespace=sys%namespace)
 
      ! the converged wfns for the previous mode are probably not a good starting point
      call lr_zero(lr(1), st)
 
      if (.not. fromscratch) then
        message(1) = "Loading restart wavefunctions for linear response."
        call messages_info(1, namespace=sys%namespace)
        call restart_load%open_dir(wfs_tag_sigma(sys%namespace, phn_wfs_tag(iatom, idir), 1), ierr)
        if (ierr == 0) then
          call states_elec_load(restart_load, sys%namespace, sys%space, st, sys%gr, sys%kpoints, &
            sys%st%restart_fixed_occ, ierr=ierr, lr = lr(1))
        end if
        if (ierr /= 0) then
          message(1) = "Unable to read response wavefunctions from '"//&
            trim(wfs_tag_sigma(sys%namespace, phn_wfs_tag(iatom, idir), 1))//"'."
          call messages_warning(1, namespace=sys%namespace)
        end if
        call restart_load%close_dir()
      end if
 
      call pert%setup_atom(iatom)
      call pert%setup_dir(idir)
 
      if (sys%hm%ep%nlcc) then
        ! The bare perturbation includes the xc response to the rigid displacement of the core density:
        ! d rho_core / d R_{iatom, idir} = -d_idir rho_core, spread over the spin channels
        call species_get_nlcc_grad_bare(ions%atom(iatom)%species, ions%space, ions%latt, ions%pos(:, iatom), &
          gr, grad_rho_core)
        drho_core(1:gr%np, 1:st%d%nspin) = m_zero
        do ispin = 1, st%d%spin_channels
          drho_core(1:gr%np, ispin) = -grad_rho_core(1:gr%np, idir) / st%d%spin_channels
        end do
        call dsternheimer_calc_hvar(sh, sys%namespace, gr, sys%hm, lr, 1, dvxc_bare, &
          lr_rho = drho_core, exclude_hartree = .true.)
        call pert%set_nlcc_response(dvxc_bare(:, :, 1))
      end if
 
      ! We now solve the Sternheimber equation in order to get the first-order change in the wavefunction
      ! This is then used to get the force derivatives
      safe_allocate(force_deriv(1:ndim, 1:natoms))
      if (states_are_real(st)) then
 
        call dsternheimer_solve(sh, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, sys%mc, &
          lr, 1, m_zero, pert, restart_dump, phn_rho_tag(iatom, idir), phn_wfs_tag(iatom, idir))
 
        if (sys%hm%ep%nlcc) then
          ! the NLCC term of the force derivatives needs the xc response to the
          ! composite (valence + core) response density
          ddrho_tot(1:gr%np, 1:st%d%nspin) = lr(1)%ddl_rho(1:gr%np, 1:st%d%nspin) + drho_core(1:gr%np, 1:st%d%nspin)
          call dsternheimer_calc_hvar(sh, sys%namespace, gr, sys%hm, lr, 1, ddvxc_resp, &
            lr_rho = ddrho_tot, exclude_hartree = .true.)
          call dforces_derivative(gr, sys%namespace, sys%space, ions, sys%hm%ep, st, sys%kpoints, lr(1), lr(1), force_deriv, &
            sys%hm%lda_u_level, vxc_response = ddvxc_resp(:, :, 1))
        else
          call dforces_derivative(gr, sys%namespace, sys%space, ions, sys%hm%ep, st, sys%kpoints, lr(1), lr(1), force_deriv, &
            sys%hm%lda_u_level)
        end if
 
      else
 
        call zsternheimer_solve(sh, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, sys%mc, &
          lr, 1, m_z0, pert, restart_dump, phn_rho_tag(iatom, idir), phn_wfs_tag(iatom, idir))
 
        if (sys%hm%ep%nlcc) then
          zdrho_tot(1:gr%np, 1:st%d%nspin) = lr(1)%zdl_rho(1:gr%np, 1:st%d%nspin) + drho_core(1:gr%np, 1:st%d%nspin)
          call zsternheimer_calc_hvar(sh, sys%namespace, gr, sys%hm, lr, 1, zdvxc_resp, &
            lr_rho = zdrho_tot, exclude_hartree = .true.)
          call zforces_derivative(gr, sys%namespace, sys%space, ions, sys%hm%ep, st, sys%kpoints, lr(1), lr(1), force_deriv, &
            sys%hm%lda_u_level, vxc_response = zdvxc_resp(:, :, 1))
        else
          call zforces_derivative(gr, sys%namespace, sys%space, ions, sys%hm%ep, st, sys%kpoints, lr(1), lr(1), force_deriv, &
            sys%hm%lda_u_level)
        end if
 
      end if
 
      do jmat = 1, vib%num_modes
        if (.not. symmetrize .and. jmat < imat) then
          vib%dyn_matrix(jmat, imat) = vib%dyn_matrix(imat, jmat)
          cycle
        end if
 
        jatom = vibrations_get_atom(vib, jmat)
        jdir  = vibrations_get_dir(vib, jmat)
 
        vib%dyn_matrix(jmat, imat) = vib%dyn_matrix(jmat, imat) + real(force_deriv(jdir, jatom), real64)
        vib%dyn_matrix(jmat, imat) = vib%dyn_matrix(jmat, imat) * vibrations_norm_factor(vib, iatom, jatom)
      end do
      safe_deallocate_a(force_deriv)
 
      call vibrations_out_dyn_matrix_row(vib, imat)
 
      if (do_infrared) then
        if (states_are_real(st)) then
          call dphonons_lr_infrared(gr, ions, st, lr(1), kdotp_lr, imat, iatom, idir, vib%infrared)
        else
          call zphonons_lr_infrared(gr, ions, st, lr(1), kdotp_lr, imat, iatom, idir, vib%infrared)
        end if
      end if
 
      iunit_restart = restart_dump%open('restart', position='append')
      ! open and close makes sure output is not buffered
      do jmat = 1, vib%num_modes
        write(line(1), *) jmat, imat, vib%dyn_matrix(jmat, imat)
        call restart_dump%write(iunit_restart, line, 1, ierr)
        if (ierr /= 0) then
          message(1) = "Could not write restart information."
          call messages_warning(1, namespace=sys%namespace)
        end if
      end do
      write(line(1), *) imat, (vib%infrared(imat, idir), idir = 1, ndim)
      call restart_dump%write(iunit_restart, line, 1, ierr)
      if (ierr /= 0) then
        message(1) = "Could not write restart information."
        call messages_warning(1, namespace=sys%namespace)
      end if
      call restart_dump%close(iunit_restart)
 
      message(1) = ""
      call messages_info(1, namespace=sys%namespace)
    end do
 
    safe_deallocate_p(pert)
 
    safe_deallocate_a(grad_rho_core)
    safe_deallocate_a(drho_core)
    safe_deallocate_a(dvxc_bare)
    safe_deallocate_a(ddrho_tot)
    safe_deallocate_a(ddvxc_resp)
    safe_deallocate_a(zdrho_tot)
    safe_deallocate_a(zdvxc_resp)
 
    if (symmetrize) call vibrations_symmetrize_dyn_matrix(vib)
    call vibrations_diag_dyn_matrix(vib)
    call vibrations_output(vib)
    call axsf_mode_output(vib, ions, gr, sys%namespace)
 
    if (do_infrared) then
      if (sys%space%is_periodic() .and. .not. smear_is_semiconducting(st%smear)) then
        message(1) = "Cannot calculate infrared intensities for periodic system with smearing (i.e. without a gap)."
        call messages_info(1, namespace=sys%namespace)
      else
        call born_from_infrared(vib, born)
        call born_output_charges(born, ions%atom, ions%charge, ions%natoms, sys%namespace, &
          ndim, vib_modes_dir, write_real = .true.)
        call calc_infrared()
      end if
 
      call born_charges_end(born)
    end if
 
    if (normal_mode_wfs) then
      message(1) = "Calculating response wavefunctions for normal modes."
      call messages_info(1, namespace=sys%namespace)
      if (states_are_real(st)) then
        call dphonons_lr_wavefunctions(lr(1), sys%namespace, sys%space, st, sys%gr, sys%kpoints, vib, restart_load, &
          restart_dump)
      else
        call zphonons_lr_wavefunctions(lr(1), sys%namespace, sys%space, st, sys%gr, sys%kpoints, vib, restart_load, &
          restart_dump)
      end if
    end if
 
    !DESTRUCT
 
    call lr_dealloc(lr(1))
    call vibrations_end(vib)
    call sternheimer_end(sh)
    call states_elec_deallocate_wfns(st)
    if (sys%space%is_periodic() .and. do_infrared) then
      do idir = 1, sys%space%periodic_dim
        call lr_dealloc(kdotp_lr(idir))
      end do
    end if
    call restart_load%end()
    call restart_dump%end()
 
    pop_sub(phonons_lr_run_legacy)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine build_ionic_dyn_matrix()
      real(real64) :: term, weight, xi(1:ndim), dx(1:ndim), r2
 
      push_sub(phonons_lr_run_legacy.build_ionic_dyn_matrix)
 
      assert(.not. ions%space%is_periodic())
 
      vib%dyn_matrix(:,:) = m_zero
 
      do iatom = 1, natoms
        xi = ions%pos(:, iatom)
 
        do jatom = 1, natoms
          if(iatom == jatom) cycle
 
          dx = xi - ions%pos(:, jatom)
          r2 = dot_product(dx, dx)
 
          weight = ions%charge(iatom) * ions%charge(jatom) /(sqrt(r2)**3)
 
          do idir = 1, ndim
            do jdir = 1, ndim
 
              term = weight * (ddelta(idir, jdir) - m_three*dx(idir)*dx(jdir)/r2)
 
              ! The force is given by F_I = Z_I \sum_K Z_K (R_I-R_K)/|R_I-R_K|^3
              ! There is therefore a diagonal term from the case J=I, with the sum over K
              vib%dyn_matrix(vibrations_get_index(vib, iatom, jdir), vibrations_get_index(vib, iatom, idir)) = &
                vib%dyn_matrix(vibrations_get_index(vib, iatom, jdir), vibrations_get_index(vib, iatom, idir)) + term
 
              vib%dyn_matrix(vibrations_get_index(vib, jatom, jdir), vibrations_get_index(vib, iatom, idir)) = &
                vib%dyn_matrix(vibrations_get_index(vib, jatom, jdir), vibrations_get_index(vib, iatom, idir)) - term
            end do
          end do
        end do
      end do
 
      pop_sub(phonons_lr_run_legacy.build_ionic_dyn_matrix)
    end subroutine build_ionic_dyn_matrix
 
    ! ---------------------------------------------------------
    subroutine calc_infrared()
 
      integer :: iunit_ir
      real(real64) :: lir(1:sys%space%dim+1)
 
      push_sub(phonons_lr_run_legacy.calc_infrared)
 
      iunit_ir = io_open(vib_modes_dir//'infrared', sys%namespace, action='write')
 
      ! dP/dQ over a dimensionless normal-mode eigenvector is a charge, so the columns
      ! below are in units of the elementary charge and do not depend on UnitsOutput.
      write(iunit_ir, '(a)', advance = 'no') '#   freq ['//trim(units_abbrev(unit_invcm))//']'
      do idir = 1, ndim
        write(iunit_ir, '(a14)', advance = 'no') '<' // index2axis(idir) // '> [e]'
      end do
      write(iunit_ir, '(a14)') 'average [e]'
 
      do iatom = 1, natoms
        do idir = 1, ndim
 
          imat = vibrations_get_index(vib, iatom, idir)
 
          write(iunit_ir, '(f17.8)', advance = 'no') units_from_atomic(unit_invcm, vib%freq(imat))
          do jdir = 1, ndim
            lir(jdir) = dot_product(vib%infrared(:, jdir), vib%normal_mode(:, imat))
            write(iunit_ir, '(f14.5)', advance = 'no') units_from_atomic(unit_one, lir(jdir))
          end do
 
          lir(ndim+1) = norm2(lir(1:ndim))/sqrt(real(ndim, real64) )
          write(iunit_ir, '(f17.8)') units_from_atomic(unit_one, lir(ndim + 1))
        end do
      end do
 
      call io_close(iunit_ir)
      pop_sub(phonons_lr_run_legacy.calc_infrared)
    end subroutine calc_infrared
 
  end subroutine phonons_lr_run_legacy
 
 
  ! ---------------------------------------------------------
  subroutine born_from_infrared(vib, born)
    type(vibrations_t),   intent(in)    :: vib
    type(born_charges_t), intent(inout) :: born
 
    integer :: imat, idir, iatom
 
    push_sub(born_from_infrared)
 
    do imat = 1, vib%num_modes
      idir = vibrations_get_dir(vib, imat)
      iatom = vibrations_get_atom(vib, imat)
      born%charge(1:vib%ndim, idir, iatom) = -vib%infrared(imat, 1:vib%ndim)
    end do
 
    pop_sub(born_from_infrared)
  end subroutine born_from_infrared
 
 
  ! ---------------------------------------------------------
  character(len=100) function phn_rho_tag(iatom, dir) result(str)
    integer, intent(in) :: iatom, dir
 
    push_sub(phn_rho_tag)
 
    write(str, '(a,i4.4,a,i1)') 'phn_rho_', iatom, '_',  dir
 
    pop_sub(phn_rho_tag)
 
  end function phn_rho_tag
 
 
  ! ---------------------------------------------------------
  character(len=100) function phn_wfs_tag(iatom, dir) result(str)
    integer, intent(in) :: iatom, dir
 
    push_sub(phn_wfs_tag)
 
    write(str, '(a,i4.4,a,a)') "phn_wfs_", iatom, "_", index2axis(dir)
 
    pop_sub(phn_wfs_tag)
 
  end function phn_wfs_tag
 
 
  ! ---------------------------------------------------------
  character(len=100) function phn_nm_wfs_tag(inm) result(str)
    integer, intent(in) :: inm
 
    push_sub(phn_nm_wfs_tag)
 
    write(str, '(a,i5.5)') "phn_nm_wfs_", inm
 
    pop_sub(phn_nm_wfs_tag)
 
  end function phn_nm_wfs_tag
 
 
  ! ---------------------------------------------------------
  subroutine axsf_mode_output(this, ions, mesh, namespace)
    type(vibrations_t), intent(in) :: this
    type(ions_t),       intent(in) :: ions
    class(mesh_t),      intent(in) :: mesh
    type(namespace_t),  intent(in) :: namespace
 
    integer :: iunit, iatom, idir, imat, jmat
    real(real64), allocatable :: forces(:,:)
    character(len=2) :: suffix
 
    if (.not. mpi_world%is_root()) return
 
    push_sub(axsf_mode_output)
 
    ! for some reason, direct usage of this%suffix gives an odd result
    suffix = vibrations_get_suffix(this)
    iunit = io_open(vib_modes_dir//'normal_modes_'//suffix//'.axsf', namespace, action='write')
 
    write(iunit, '(a,i6)') 'ANIMSTEPS ', this%num_modes
    safe_allocate(forces(1:ions%space%dim, 1:ions%natoms))
    do imat = 1, this%num_modes
      do jmat = 1, this%num_modes
        iatom = vibrations_get_atom(this, jmat)
        idir  = vibrations_get_dir(this, jmat)
        forces(idir, iatom) = this%normal_mode(jmat, imat)
      end do
      call write_xsf_geometry(iunit, ions%space, ions%latt, ions%pos, ions%atom, mesh, forces = forces, index = imat)
    end do
    safe_deallocate_a(forces)
    call io_close(iunit)
 
    pop_sub(axsf_mode_output)
  end subroutine axsf_mode_output
 
  ! ---------------------------------------------------------
  subroutine phonons_load(restart, vib, start_mode)
    type(restart_t),    intent(in)    :: restart
    type(vibrations_t), intent(inout) :: vib
    integer,            intent(out)   :: start_mode
 
    integer :: iunit, ierr, imode, jmode, imode_read, jmode_read
    character(len=120) :: line(1)
 
    push_sub(phonons_load)
 
    iunit = restart%open('restart')
    if (iunit /= -1) then
      imode_loop: do imode = 1, vib%num_modes
        do jmode = 1, vib%num_modes
          call restart%read(iunit, line, 1, ierr)
          if (ierr /= 0) exit imode_loop
          read(line(1), fmt=*, iostat=ierr) jmode_read, imode_read, vib%dyn_matrix(jmode, imode)
          if (imode_read /= imode) then
            write(message(1),'(a,i9,a,i9)') "Corruption of restart data: row ", imode, " is labeled as ", imode_read
            call messages_fatal(1)
          end if
          if (jmode_read /= jmode) then
            write(message(1),'(a,i9,a,i9)') "Corruption of restart data: column ", jmode, " is labeled as ", jmode_read
            call messages_fatal(1)
          end if
        end do
 
        call restart%read(iunit, line, 1, ierr)
        if (ierr /= 0) exit
 
        start_mode = imode + 1
 
        read(line(1), fmt=*, iostat=ierr) imode_read, vib%infrared(imode, 1:vib%ndim)
        if (imode_read /= imode) then
          write(message(1),'(a,i9,a,i9)') "Corruption of restart data: infrared row ", imode, " is labeled as ", imode_read
          call messages_fatal(1)
        end if
      end do imode_loop
 
      write(message(1),'(a,i9,a,i9)') 'Info: Read saved dynamical-matrix rows for ', &
        start_mode - 1, ' modes out of ', vib%num_modes
      call messages_info(1)
 
      call restart%close(iunit)
    else
      start_mode = 1
 
      message(1) = "Could not open restart file 'restart'. Starting from scratch."
      call messages_warning(1)
    end if
 
    pop_sub(phonons_load)
  end subroutine phonons_load
 
#include "complex.F90"
#include "phonons_lr_inc.F90"
 
#include "undef.F90"
 
#include "real.F90"
#include "phonons_lr_inc.F90"
 
end module phonons_lr_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: