static_pol.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"
#define RESTART_FILE 'dipoles'
 
module static_pol_oct_m
  use born_charges_oct_m
  use debug_oct_m
  use elf_oct_m
  use em_resp_oct_m
  use global_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 lcao_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 profiling_oct_m
  use restart_oct_m
  use scf_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use electrons_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
 
  implicit none
 
  private
  public :: &
    static_pol_run
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine static_pol_run(system, from_scratch)
    class(*),        intent(inout) :: system
    logical,         intent(in)    :: from_scratch
 
    push_sub(static_pol_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = static_pol not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call static_pol_run_legacy(system, from_scratch)
    class default
      assert(.false.)
    end select
 
    pop_sub(static_pol_run)
  end subroutine static_pol_run
 
  ! ---------------------------------------------------------
  subroutine static_pol_run_legacy(sys, fromScratch)
    type(electrons_t),    intent(inout) :: sys
    logical,              intent(in)    :: fromScratch
 
    type(scf_t) :: scfv
    integer :: iunit, ios, i_start, ii, jj, is, isign, ierr, read_count, verbosity
    real(real64) :: e_field, e_field_saved
    real(real64), allocatable :: Vpsl_save(:), trrho(:), dipole(:, :, :)
    real(real64), allocatable :: elf(:,:), lr_elf(:,:), elfd(:,:), lr_elfd(:,:)
    real(real64), allocatable :: lr_rho(:,:), lr_rho2(:,:), gs_rho(:,:), tmp_rho(:,:)
    real(real64) :: center_dipole(1:sys%space%dim), diag_dipole(1:sys%space%dim), ionic_dipole(1:sys%space%dim), &
      print_dipole(1:sys%space%dim)
    type(born_charges_t) :: born_charges
    logical :: calc_born, start_density_is_zero_field, write_restart_densities, calc_diagonal, verbose
    logical :: diagonal_done, center_written, fromScratch_local, field_written
    character(len=MAX_PATH_LEN) :: fname, dir_name
    character(len=120) :: line(1)
    character :: sign_char
    type(restart_t) :: gs_restart, restart_load, restart_dump
 
    push_sub(static_pol_run_legacy)
 
    if (sys%hm%pcm%run_pcm) then
      call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
    end if
 
    if (sys%kpoints%use_symmetries) then
      call messages_experimental("KPoints symmetries with CalculationMode = em_resp", namespace=sys%namespace)
    end if
 
    call init_()
 
    ! load wavefunctions
    call restart_init(gs_restart, sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr, exact=.true.)
    if (ierr == 0) then
      call states_elec_load(gs_restart, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
        fixed_occ=.false., ierr=ierr)
    end if
    if (ierr /= 0) then
      message(1) = "Unable to read wavefunctions."
      call messages_fatal(1, namespace=sys%namespace)
    end if
    call restart_end(gs_restart)
 
    if (sys%space%is_periodic()) then
      message(1) = "Electric field cannot be applied to a periodic system (currently)."
      call messages_fatal(1, namespace=sys%namespace)
    end if
 
    ! set up Hamiltonian
    message(1) = 'Info: Setting up Hamiltonian.'
    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, calc_eigenval = .false.) ! we read them from restart
 
    ! Allocate the dipole
    safe_allocate(dipole(1:sys%space%dim, 1:sys%space%dim, 1:2))
    dipole = m_zero
 
    i_start = 1
    center_written = .false.
    diagonal_done = .false.
    field_written = .false.
 
    call restart_init(restart_dump, sys%namespace, restart_em_resp_fd, restart_type_dump, sys%mc, ierr, mesh=sys%gr)
 
    if (.not. fromscratch) then
      call restart_init(restart_load, sys%namespace, restart_em_resp_fd, restart_type_load, sys%mc, ierr, mesh=sys%gr)
      if (ierr == 0) then
        iunit = restart_open(restart_load, restart_file)
      else
        iunit = -1
      end if
 
      if (iunit /= -1) then
        ! Finds out how many dipoles have already been written.
        rewind(iunit)
 
        read(iunit, fmt=*, iostat = ios) e_field_saved
        field_written = (ios == 0)
 
        read(iunit, fmt=*, iostat = ios) (center_dipole(jj), jj = 1, sys%space%dim)
        center_written = (ios == 0)
 
        do ii = 1, sys%space%dim
          read(iunit, fmt=*, iostat = ios) ((dipole(ii, jj, isign), jj = 1, sys%space%dim), isign = 1, 2)
          if (ios /= 0) exit
          i_start = i_start + 1
        end do
 
        read(iunit, fmt=*, iostat = ios) (diag_dipole(jj), jj = 1, sys%space%dim)
        diagonal_done = (ios == 0)
 
        call restart_close(restart_load, iunit)
      end if
 
      ! if saved dipoles used a different e_field, we cannot use them
      if (field_written .and. abs(e_field_saved - e_field) > m_epsilon) then
        message(1) = "Saved dipoles are from a different electric field, cannot use them."
        call messages_warning(1, namespace=sys%namespace)
        center_written = .false.
        diagonal_done = .false.
        i_start = 1
      end if
 
      read_count = (i_start - 1) * 2
      if (center_written) read_count = read_count + 1
      if (diagonal_done)  read_count = read_count + 1
      write(message(1),'(a,i5,a)') "Using ", read_count, " dipole(s) from file."
      call messages_info(1, namespace=sys%namespace)
    end if
 
    if (sys%outp%what(option__output__density) .or. &
      sys%outp%what(option__output__pol_density)) then
      if (i_start > 2 .and. calc_diagonal) then
        i_start = 2
        diagonal_done = .false.
        !FIXME: take derivatives between yz and z (not y) so can restart from only last (z) calc
      end if
    end if
 
    if (i_start == 1) then
      ! open new file, erase old data, write e_field
      iunit = restart_open(restart_dump, restart_file, status='replace')
      write(line(1), fmt='(e20.12)') e_field
      call restart_write(restart_dump, iunit, line, 1, ierr)
      if (ierr /= 0) then
        message(1) = "Unsuccessful write of electric field."
        call messages_warning(1, namespace=sys%namespace)
      end if
      call restart_close(restart_dump, iunit)
      center_written = .false.
    end if
 
    ! Save local potential
    safe_allocate(vpsl_save(1:sys%gr%np))
    vpsl_save = sys%hm%ep%Vpsl
 
    ! Allocate the trrho to contain the trace of the density.
    safe_allocate(trrho(1:sys%gr%np))
    safe_allocate(gs_rho(1:sys%gr%np, 1:sys%st%d%nspin))
    safe_allocate(tmp_rho(1:sys%gr%np, 1:sys%st%d%nspin))
    trrho = m_zero
    gs_rho = m_zero
 
    call output_init_()
    call scf_init(scfv, sys%namespace, sys%gr, sys%ions, sys%st, sys%mc, sys%hm, sys%space)
    call born_charges_init(born_charges, sys%namespace, sys%ions%natoms, sys%st%val_charge, &
      sys%st%qtot, sys%space%dim)
 
    ! now calculate the dipole without field
 
    sys%hm%ep%vpsl(1:sys%gr%np) = vpsl_save(1:sys%gr%np)
    call sys%hm%update(sys%gr, sys%namespace, sys%space, sys%ext_partners)
 
    write(message(1), '(a)')
    write(message(2), '(a)') 'Info: Calculating dipole moment for zero field.'
    call messages_info(2, namespace=sys%namespace)
    call scf_run(scfv, sys%namespace, sys%space, sys%mc, sys%gr, sys%ions, sys%ext_partners, sys%st, &
      sys%ks, sys%hm, verbosity = verbosity)
 
    gs_rho(1:sys%gr%np, 1:sys%st%d%nspin) = sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin)
    trrho = m_zero
    do is = 1, sys%st%d%spin_channels
      trrho(1:sys%gr%np) = trrho(1:sys%gr%np) + gs_rho(1:sys%gr%np, is)
    end do
 
    ! calculate dipole
    do jj = 1, sys%space%dim
      center_dipole(jj) = dmf_moment(sys%gr, trrho, jj, 1)
    end do
 
    ! Writes the dipole to file
    if (.not. center_written) then
      iunit = restart_open(restart_dump, restart_file, position='append')
      write(line(1), fmt='(6e20.12)') (center_dipole(jj), jj = 1, sys%space%dim)
      call restart_write(restart_dump, iunit, line, 1, ierr)
      if (ierr /= 0) then
        message(1) = "Unsuccessful write of center dipole."
        call messages_warning(1, namespace=sys%namespace)
      end if
      call restart_close(restart_dump, iunit)
    end if
 
    if (mpi_grp_is_root(mpi_world)) then
      ionic_dipole(1:sys%space%dim) = sys%ions%dipole()
      print_dipole(1:sys%space%dim) = center_dipole(1:sys%space%dim) + ionic_dipole(1:sys%space%dim)
      call output_dipole(print_dipole, sys%space%dim, namespace=sys%namespace)
    end if
 
    do ii = i_start, sys%space%dim
      do isign = 1, 2
        write(message(1), '(a)')
        write(message(2), '(a,f6.4,5a)') 'Info: Calculating dipole moment for field ', &
          units_from_atomic(units_out%force, (-1)**isign * e_field), ' ', &
          trim(units_abbrev(units_out%force)), ' in the ', index2axis(ii), '-direction.'
        call messages_info(2, namespace=sys%namespace)
        ! there would be an extra factor of -1 in here that is for the electronic charge
        ! except that we treat electrons as positive
 
        sys%hm%ep%vpsl(1:sys%gr%np) = vpsl_save(1:sys%gr%np) + (-1)**isign * sys%gr%x(1:sys%gr%np, ii) * e_field
        call sys%hm%update(sys%gr, sys%namespace, sys%space, sys%ext_partners)
 
        if (isign == 1) then
          sign_char = '+'
        else
          sign_char = '-'
        end if
 
        write(dir_name,'(a)') "field_"//index2axis(ii)//sign_char
        fromscratch_local = fromscratch
 
        if (.not. fromscratch) then
          call restart_open_dir(restart_load, trim(dir_name), ierr)
          if (ierr == 0) then
            call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
              fixed_occ=.false., ierr=ierr)
          end if
          call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
          if (ierr /= 0) fromscratch_local = .true.
          call restart_close_dir(restart_load)
        end if
 
        if (fromscratch_local) then
          if (start_density_is_zero_field) then
            sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin) = gs_rho(1:sys%gr%np, 1:sys%st%d%nspin)
            call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
          else
            call lcao_run(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, &
              sys%hm, lmm_r = scfv%lmm_r)
          end if
        end if
 
        call scf_mix_clear(scfv)
        call scf_run(scfv, sys%namespace, sys%space, sys%mc, sys%gr, sys%ions, sys%ext_partners, sys%st, &
          sys%ks, sys%hm, verbosity = verbosity)
 
        trrho = m_zero
        do is = 1, sys%st%d%spin_channels
          trrho(1:sys%gr%np) = trrho(1:sys%gr%np) + sys%st%rho(1:sys%gr%np, is)
        end do
 
        ! calculate dipole
        do jj = 1, sys%space%dim
          dipole(ii, jj, isign) = dmf_moment(sys%gr, trrho, jj, 1)
        end do
 
        if (mpi_grp_is_root(mpi_world)) then
          print_dipole(1:sys%space%dim) = dipole(ii, 1:sys%space%dim, isign) + ionic_dipole(1:sys%space%dim)
          call output_dipole(print_dipole, sys%space%dim, namespace=sys%namespace)
        end if
 
        call output_cycle_()
 
        if (write_restart_densities) then
          call restart_open_dir(restart_dump, trim(dir_name), ierr)
          if (ierr == 0) then
            call states_elec_dump(restart_dump, sys%space, sys%st, sys%gr, sys%kpoints, ierr)
          end if
          call restart_close_dir(restart_dump)
          if (ierr /= 0) then
            message(1) = 'Unable to write states wavefunctions.'
            call messages_warning(1, namespace=sys%namespace)
          end if
        end if
      end do
 
      ! Writes the dipole to file
      iunit = restart_open(restart_dump, restart_file, position='append')
      write(line(1), '(6e20.12)') ((dipole(ii, jj, isign), jj = 1, sys%space%dim), isign = 1, 2)
      call restart_write(restart_dump, iunit, line, 1, ierr)
      if (ierr /= 0) then
        message(1) = "Unsuccessful write of dipole."
        call messages_warning(1, namespace=sys%namespace)
      end if
      call restart_close(restart_dump, iunit)
    end do
 
    if (.not. diagonal_done .and. calc_diagonal) then
      write(message(1), '(a)')
      write(message(2), '(a,f6.4,3a, f6.4, 3a)') 'Info: Calculating dipole moment for field ', &
        units_from_atomic(units_out%force, e_field), ' ', &
        trim(units_abbrev(units_out%force)), ' in the '//index2axis(2)//'-direction plus ', &
        units_from_atomic(units_out%force, e_field), ' ', &
        trim(units_abbrev(units_out%force)), ' in the '//index2axis(3)//'-direction.'
      call messages_info(2, namespace=sys%namespace)
 
      sys%hm%ep%vpsl(1:sys%gr%np) = vpsl_save(1:sys%gr%np) &
        - (sys%gr%x(1:sys%gr%np, 2) + sys%gr%x(1:sys%gr%np, 3)) * e_field
      call sys%hm%update(sys%gr, sys%namespace, sys%space, sys%ext_partners)
 
      if (isign == 1) then
        sign_char = '+'
      else
        sign_char = '-'
      end if
 
      fromscratch_local = fromscratch
 
      if (.not. fromscratch) then
        call restart_open_dir(restart_load, "field_yz+", ierr)
        if (ierr == 0) then
          call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
            fixed_occ=.false., ierr=ierr)
        end if
        call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
        if (ierr /= 0) fromscratch_local = .true.
        call restart_close_dir(restart_load)
      end if
 
      if (fromscratch_local) then
        if (start_density_is_zero_field) then
          sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin) = gs_rho(1:sys%gr%np, 1:sys%st%d%nspin)
          call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
        else
          call lcao_run(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, &
            sys%hm, lmm_r = scfv%lmm_r)
        end if
      end if
 
      call scf_mix_clear(scfv)
      call scf_run(scfv, sys%namespace, sys%space, sys%mc, sys%gr, sys%ions, sys%ext_partners, sys%st, &
        sys%ks, sys%hm, verbosity = verbosity)
 
      trrho = m_zero
      do is = 1, sys%st%d%spin_channels
        trrho(1:sys%gr%np) = trrho(1:sys%gr%np) + sys%st%rho(1:sys%gr%np, is)
      end do
 
      ! calculate dipole
      do jj = 1, sys%space%dim
        diag_dipole(jj) = dmf_moment(sys%gr, trrho, jj, 1)
      end do
 
      if (mpi_grp_is_root(mpi_world)) then
        print_dipole(1:sys%space%dim) = diag_dipole(1:sys%space%dim) + ionic_dipole(1:sys%space%dim)
        call output_dipole(print_dipole, sys%space%dim, namespace=sys%namespace)
      end if
 
      ! Writes the dipole to file
      iunit = restart_open(restart_dump, restart_file, position='append')
      write(line(1), fmt='(3e20.12)') (diag_dipole(jj), jj = 1, sys%space%dim)
      call restart_write(restart_dump, iunit, line, 1, ierr)
      if (ierr /= 0) then
        message(1) = "Unsuccessful write of dipole."
        call messages_warning(1, namespace=sys%namespace)
      end if
      call restart_close(restart_dump, iunit)
 
      if (write_restart_densities) then
        call restart_open_dir(restart_dump, "field_yz+", ierr)
        if (ierr == 0) then
          call states_elec_dump(restart_dump, sys%space, sys%st, sys%gr, sys%kpoints, ierr)
        end if
        call restart_close_dir(restart_dump)
        if (ierr /= 0) then
          message(1) = 'Unable to write states wavefunctions.'
          call messages_warning(1, namespace=sys%namespace)
        end if
      end if
 
    end if
 
    if (.not. fromscratch) call restart_end(restart_load)
    call scf_end(scfv)
    call output_end_()
    call born_charges_end(born_charges)
 
    safe_deallocate_a(vpsl_save)
    safe_deallocate_a(trrho)
    safe_deallocate_a(gs_rho)
    safe_deallocate_a(tmp_rho)
    safe_deallocate_a(dipole)
    call end_()
    pop_sub(static_pol_run_legacy)
 
  contains
 
    ! ---------------------------------------------------------
    subroutine init_()
      push_sub(static_pol_run_legacy.init_)
 
      call states_elec_allocate_wfns(sys%st, sys%gr)
 
      call messages_obsolete_variable(sys%namespace, "EMStaticField", "EMStaticElectricField")
      !%Variable EMStaticElectricField
      !%Type float
      !%Default 0.01 a.u.
      !%Section Linear Response::Static Polarization
      !%Description
      !% Magnitude of the static electric field used to calculate the static polarizability,
      !% if <tt>ResponseMethod = finite_differences</tt>.
      !%End
      call parse_variable(sys%namespace, 'EMStaticElectricField', 0.01_real64, e_field, units_inp%force)
      if (e_field <= m_zero) then
        write(message(1), '(a,e14.6,a)') "Input: '", e_field, "' is not a valid EMStaticElectricField."
        message(2) = '(Must have EMStaticElectricField > 0)'
        call messages_fatal(2, namespace=sys%namespace)
      end if
 
      ! variable defined in em_resp
      call parse_variable(sys%namespace, 'EMCalcbornCharges', .false., calc_born)
      if (calc_born) call messages_experimental("Calculation of born effective charges", namespace=sys%namespace)
 
      !%Variable EMStartDensityIsZeroField
      !%Type logical
      !%Default true
      !%Section Linear Response::Static Polarization
      !%Description
      !% Use the charge density from the zero-field calculation as the starting density for
      !% SCF calculations with applied fields. For small fields, this will be fastest.
      !% If there is trouble converging with larger fields, set to false,
      !% to initialize the calculation for each field from scratch, as specified by the LCAO variables.
      !% Only applies if <tt>ResponseMethod = finite_differences</tt>.
      !%End
      call parse_variable(sys%namespace, 'EMStartDensityIsZeroField', .true., start_density_is_zero_field)
 
      !%Variable EMCalcDiagonalField
      !%Type logical
      !%Default true
      !%Section Linear Response::Static Polarization
      !%Description
      !% Calculate <i>yz</i>-field for <math>\beta_{xyz}</math> hyperpolarizability, which is sometimes harder to converge.
      !% Only applies if <tt>ResponseMethod = finite_differences</tt>.
      !%End
      call parse_variable(sys%namespace, 'EMCalcDiagonalField', .true., calc_diagonal)
 
      !%Variable EMWriteRestartDensities
      !%Type logical
      !%Default true
      !%Section Linear Response::Static Polarization
      !%Description
      !% Write density after each calculation for restart, rather than just the resulting electronic dipole moment.
      !% Only applies if <tt>ResponseMethod = finite_differences</tt>. Restarting from calculations at smaller
      !% fields can be helpful if there are convergence problems.
      !%End
      call parse_variable(sys%namespace, 'EMWriteRestartDensities', .true., write_restart_densities)
 
      !%Variable EMVerbose
      !%Type logical
      !%Default false
      !%Section Linear Response::Static Polarization
      !%Description
      !% Write full SCF output.
      !% Only applies if <tt>ResponseMethod = finite_differences</tt>.
      !%End
      call parse_variable(sys%namespace, 'EMVerbose', .false., verbose)
 
      if (verbose) then
        verbosity = verb_full
      else
        verbosity = verb_compact
      end if
 
      pop_sub(static_pol_run_legacy.init_)
    end subroutine init_
 
    ! ---------------------------------------------------------
    subroutine end_()
 
      push_sub(end_)
 
      call states_elec_deallocate_wfns(sys%st)
 
      pop_sub(end_)
    end subroutine end_
 
 
    !-------------------------------------------------------------
    subroutine output_init_()
      push_sub(output_init_)
 
      !allocate memory for what we want to output
      if (sys%outp%what(option__output__density) .or. &
        sys%outp%what(option__output__pol_density)) then
        safe_allocate(lr_rho(1:sys%gr%np, 1:sys%st%d%nspin))
        safe_allocate(lr_rho2(1:sys%gr%np, 1:sys%st%d%nspin))
      end if
 
      if (sys%outp%what(option__output__elf)) then
        safe_allocate(    elf(1:sys%gr%np, 1:sys%st%d%nspin))
        safe_allocate( lr_elf(1:sys%gr%np, 1:sys%st%d%nspin))
        safe_allocate(   elfd(1:sys%gr%np, 1:sys%st%d%nspin))
        safe_allocate(lr_elfd(1:sys%gr%np, 1:sys%st%d%nspin))
      end if
 
      pop_sub(output_init_)
    end subroutine output_init_
 
 
    !-------------------------------------------------------------
    subroutine output_cycle_()
      integer :: iatom
      type(unit_t) :: fn_unit
 
      push_sub(output_cycle_)
 
      ! BORN CHARGES
      if (calc_born) then
        do iatom = 1, sys%ions%natoms
          if (isign == 1) then
            ! temporary assignment for use in next cycle when isign == 2
            born_charges%charge(ii, 1:sys%space%dim, iatom) = sys%ions%tot_force(:, iatom)
          else
            born_charges%charge(ii, 1:sys%space%dim, iatom) = &
              (sys%ions%tot_force(:, iatom) - born_charges%charge(ii, 1:sys%space%dim, iatom)) &
              / (m_two*e_field)
            born_charges%charge(ii, ii, iatom) = born_charges%charge(ii, ii, iatom) + sys%ions%charge(iatom)
            ! since the efield is applied in the SCF calculation by just altering the external potential felt by the electrons,
            ! the ionic force due to the efield is not included in the forces returned by the SCF run, and so the ionic
            ! contribution to the born charge must be added by hand here
          end if
        end do
      end if
 
      !DENSITY AND POLARIZABILITY DENSITY
      if (sys%outp%what(option__output__density) .or. &
        sys%outp%what(option__output__pol_density)) then
 
        if (isign == 1 .and. ii == 2) then
          tmp_rho(1:sys%gr%np, 1:sys%st%d%nspin) = sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin)
          ! for use in off-diagonal non-linear densities
        end if
 
        if (isign == 1) then
          ! temporary assignment for use in next cycle when isign == 2
          lr_rho(1:sys%gr%np, 1:sys%st%d%nspin) = sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin)
        else
          lr_rho2(1:sys%gr%np, 1:sys%st%d%nspin) = &
            -(sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin) + lr_rho(1:sys%gr%np, 1:sys%st%d%nspin) - &
            2 * gs_rho(1:sys%gr%np, 1:sys%st%d%nspin)) / e_field**2
 
          lr_rho(1:sys%gr%np, 1:sys%st%d%nspin) = &
            (sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin) - lr_rho(1:sys%gr%np, 1:sys%st%d%nspin)) / (m_two*e_field)
 
          !write
          do is = 1, sys%st%d%nspin
            if (sys%outp%what(option__output__density)) then
              fn_unit = units_out%length**(1-sys%space%dim) / units_out%energy
              write(fname, '(a,i1,2a)') 'fd_density-sp', is, '-', index2axis(ii)
              call dio_function_output(sys%outp%how(option__output__density), em_resp_fd_dir, trim(fname),&
                sys%namespace, sys%space, sys%gr, lr_rho(:, is), fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
 
              ! save the trouble of writing many copies of each density, since ii,jj = jj,ii
              fn_unit = units_out%length**(2-sys%space%dim) / units_out%energy**2
              do jj = ii, sys%space%dim
                write(fname, '(a,i1,4a)') 'fd2_density-sp', is, '-', index2axis(ii), '-', index2axis(jj)
                call dio_function_output(sys%outp%how(option__output__density), em_resp_fd_dir, trim(fname),&
                  sys%namespace, sys%space, sys%gr, lr_rho2(:, is), fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
              end do
            end if
 
            if (sys%outp%what(option__output__pol_density)) then
              do jj = ii, sys%space%dim
                fn_unit = units_out%length**(2-sys%space%dim) / units_out%energy
                write(fname, '(a,i1,4a)') 'alpha_density-sp', is, '-', index2axis(ii), '-', index2axis(jj)
                call dio_function_output(sys%outp%how(option__output__pol_density), em_resp_fd_dir, trim(fname), &
                  sys%namespace, sys%space, sys%gr, -sys%gr%x(:, jj) * lr_rho(:, is), &
                  fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
 
                fn_unit = units_out%length**(3-sys%space%dim) / units_out%energy**2
                write(fname, '(a,i1,6a)') 'beta_density-sp', is, '-', index2axis(ii), &
                  '-', index2axis(ii), '-', index2axis(jj)
                call dio_function_output(sys%outp%how(option__output__pol_density), em_resp_fd_dir, trim(fname), &
                  sys%namespace, sys%space, sys%gr, -sys%gr%x(:, jj) * lr_rho2(:, is), &
                  fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
              end do
            end if
          end do
 
        end if
      end if
 
      !ELF
      if (sys%outp%what(option__output__elf)) then
 
        if (isign == 1) then
          call elf_calc(sys%space, sys%st, sys%gr, sys%kpoints, elf, elfd)
        else
          call elf_calc(sys%space, sys%st, sys%gr, sys%kpoints, lr_elf, lr_elfd)
 
          !numerical derivative
          lr_elf(1:sys%gr%np, 1:sys%st%d%nspin) = &
            ( lr_elf(1:sys%gr%np, 1:sys%st%d%nspin) -  elf(1:sys%gr%np, 1:sys%st%d%nspin)) / (m_two * e_field)
          lr_elfd(1:sys%gr%np, 1:sys%st%d%nspin) = &
            (lr_elfd(1:sys%gr%np, 1:sys%st%d%nspin) - elfd(1:sys%gr%np, 1:sys%st%d%nspin)) / (m_two * e_field)
 
          !write
          do is = 1, sys%st%d%nspin
            write(fname, '(a,i1,2a)') 'lr_elf-sp', is, '-', index2axis(ii)
            call dio_function_output(sys%outp%how(option__output__elf), em_resp_fd_dir, trim(fname),&
              sys%namespace, sys%space, sys%gr, lr_elf(:, is), unit_one, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
            write(fname, '(a,i1,2a)') 'lr_elf_D-sp', is, '-', index2axis(ii)
            call dio_function_output(sys%outp%how(option__output__elf), em_resp_fd_dir, trim(fname),&
              sys%namespace, sys%space, sys%gr, lr_elfd(:, is), unit_one, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
          end do
        end if
 
      end if
 
      pop_sub(output_cycle_)
    end subroutine output_cycle_
 
 
    !-------------------------------------------------------------
    subroutine output_end_()
      real(real64) :: alpha(sys%space%dim, sys%space%dim)
      complex(real64) :: beta(sys%space%dim, sys%space%dim, sys%space%dim)
      integer :: iunit, idir
      type(unit_t) :: fn_unit
      real(real64) :: freq_factor(3)
 
      push_sub(output_end_)
 
      call io_mkdir(em_resp_fd_dir, sys%namespace)
 
      if ((sys%outp%what(option__output__density) .or. &
        sys%outp%what(option__output__pol_density)) .and. calc_diagonal) then
        lr_rho2(1:sys%gr%np, 1:sys%st%d%nspin) = &
          -(sys%st%rho(1:sys%gr%np, 1:sys%st%d%nspin) - lr_rho(1:sys%gr%np, 1:sys%st%d%nspin) &
          - tmp_rho(1:sys%gr%np, 1:sys%st%d%nspin) + gs_rho(1:sys%gr%np, 1:sys%st%d%nspin)) / e_field**2
 
        do is = 1, sys%st%d%nspin
          if (sys%outp%what(option__output__density)) then
            fn_unit = units_out%length**(2-sys%space%dim) / units_out%energy**2
            write(fname, '(a,i1,a)') 'fd2_density-sp', is, '-y-z'
            call dio_function_output(sys%outp%how(option__output__density), em_resp_fd_dir, trim(fname),&
              sys%namespace, sys%space, sys%gr, lr_rho2(:, is), fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
          end if
 
          if (sys%outp%what(option__output__pol_density)) then
            fn_unit = units_out%length**(3-sys%space%dim) / units_out%energy**2
            write(fname, '(a,i1,a)') 'beta_density-sp', is, '-x-y-z'
            call dio_function_output(sys%outp%how(option__output__pol_density), em_resp_fd_dir, trim(fname),&
              sys%namespace, sys%space, sys%gr, -sys%gr%x(:, 1) * lr_rho2(:, is), &
              fn_unit, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
          end if
        end do
      end if
 
      if (mpi_grp_is_root(mpi_world)) then ! output pol file
        iunit = io_open(em_resp_fd_dir//'alpha', sys%namespace, action='write')
        write(iunit, '(3a)') '# Polarizability tensor [', trim(units_abbrev(units_out%polarizability)), ']'
 
        alpha(1:sys%space%dim, 1:sys%space%dim) = (dipole(1:sys%space%dim, 1:sys%space%dim, 1) - &
          dipole(1:sys%space%dim, 1:sys%space%dim, 2)) / (m_two * e_field)
 
        beta = m_zero
 
        do idir = 1, sys%space%dim
          beta(1:sys%space%dim, idir, idir) = &
            -(dipole(idir, 1:sys%space%dim, 1) + dipole(idir, 1:sys%space%dim, 2) - &
            m_two * center_dipole(1:sys%space%dim)) / e_field**2
          beta(idir, 1:sys%space%dim, idir) = beta(1:sys%space%dim, idir, idir)
          beta(idir, idir, 1:sys%space%dim) = beta(1:sys%space%dim, idir, idir)
        end do
 
        assert(sys%space%dim == 3)
 
        if (calc_diagonal) then
          beta(1, 2, 3) = -(diag_dipole(1) - dipole(2, 1, 1) - dipole(3, 1, 1) + center_dipole(1)) / e_field**2
        else
          beta(1, 2, 3) = m_zero
        end if
 
        beta(2, 3, 1) = beta(1, 2, 3)
        beta(3, 1, 2) = beta(1, 2, 3)
        beta(3, 2, 1) = beta(1, 2, 3)
        beta(1, 3, 2) = beta(1, 2, 3)
        beta(2, 1, 3) = beta(1, 2, 3)
 
        call output_tensor(alpha, sys%space%dim, units_out%polarizability, iunit=iunit)
        call io_close(iunit)
 
        freq_factor(1:3) = m_zero ! for compatibility with em_resp version
        call out_hyperpolarizability(sys%gr%box, beta, freq_factor(1:3), .true., em_resp_fd_dir, sys%namespace)
 
        if (calc_born) then
          call born_output_charges(born_charges, sys%ions%atom, sys%ions%charge, sys%ions%natoms, &
            sys%namespace, sys%space%dim, em_resp_fd_dir, states_are_real(sys%st))
        end if
      end if
 
      if (sys%outp%what(option__output__density) .or. &
        sys%outp%what(option__output__pol_density)) then
        safe_deallocate_a(lr_rho)
        safe_deallocate_a(lr_rho2)
      end if
 
      if (sys%outp%what(option__output__elf)) then
        safe_deallocate_a(lr_elf)
        safe_deallocate_a(elf)
        safe_deallocate_a(lr_elfd)
        safe_deallocate_a(elfd)
      end if
 
      pop_sub(output_end_)
    end subroutine output_end_
 
  end subroutine static_pol_run_legacy
 
end module static_pol_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: