xc_ks_inversion.F90

Go to the documentation of this file.

!! Copyright (C) 2010 H. Appel, N. Helbig
!!
!! 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 xc_ks_inversion_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use eigensolver_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 interaction_partner_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lasers_oct_m
  use lalg_basic_oct_m
  use mesh_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use splines_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_io_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use xc_f03_lib_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                        &
    xc_ks_inversion_t,             &
    xc_ks_inversion_init,          &
    xc_ks_inversion_end,           &
    xc_ks_inversion_write_info,    &
    xc_ks_inversion_calc,          &
    invertks_2part,                &
    invertks_iter,                 &
    invertks_update_hamiltonian
 
  integer, public, parameter ::      &
    XC_INV_METHOD_TWO_PARTICLE = 1,  &
    xc_inv_method_vs_iter      = 2,  &
    xc_inv_method_iter_stella  = 3,  &
    xc_inv_method_iter_godby   = 4
 
  integer, public, parameter ::      &
    XC_KS_INVERSION_NONE      = 1,   &
    xc_ks_inversion_adiabatic = 2,   &
    xc_ks_inversion_td_exact  = 3
 
  integer, public, parameter ::      &
    XC_ASYMPTOTICS_NONE    = 1,      &
    xc_asymptotics_sc      = 2
 
  integer, parameter ::              &
    XC_FLAGS_NONE = 0
 
  type xc_ks_inversion_t
    private
    integer,             public :: method
    integer                     :: level = xc_ks_inversion_none
    integer                     :: asymptotics
    real(real64), allocatable          :: vhxc_previous_step(:,:)
    type(states_elec_t), public :: aux_st
    type(hamiltonian_elec_t)    :: aux_hm
    type(eigensolver_t), public :: eigensolver
    ! List with all the external partners
    ! This will become a list of interactions in the future
    type(partner_list_t) :: ext_partners
 
    integer,             public :: max_iter
  end type xc_ks_inversion_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine xc_ks_inversion_init(ks_inv, namespace, gr, ions, st, xc, mc, space, kpoints)
    type(xc_ks_inversion_t), intent(inout) :: ks_inv
    type(namespace_t),       intent(in)    :: namespace
    type(grid_t),            intent(inout) :: gr
    type(ions_t),            intent(inout) :: ions
    type(states_elec_t),     intent(in)    :: st
    type(xc_t),              intent(in)    :: xc
    type(multicomm_t),       intent(in)    :: mc
    class(space_t),          intent(in)    :: space
    type(kpoints_t),         intent(in)    :: kpoints
 
    class(lasers_t), pointer :: lasers
 
    push_sub(xc_ks_inversion_init)
 
    if(gr%parallel_in_domains) then
      call messages_not_implemented("Kohn-Sham inversion in parallel", namespace=namespace)
    end if
 
    call messages_experimental("Kohn-Sham inversion")
 
    !%Variable InvertKSmethod
    !%Type integer
    !%Default iter_godby
    !%Section Calculation Modes::Invert KS
    !%Description
    !% Selects whether the exact two-particle method or the iterative scheme
    !% is used to invert the density to get the KS potential.
    !%Option two_particle 1
    !% Exact two-particle scheme.
    !%Option iterative 2
    !% Iterative scheme for <math>v_s</math>.
    !%Option iter_stella 3
    !% Iterative scheme for <math>v_s</math> using Stella and Verstraete method.
    !%Option iter_godby 4
    !% Iterative scheme for <math>v_s</math> using power method from Rex Godby.
    !%End
    call parse_variable(namespace, 'InvertKSmethod', xc_inv_method_iter_godby, ks_inv%method)
 
    if (ks_inv%method < xc_inv_method_two_particle &
      .or. ks_inv%method > xc_inv_method_iter_godby) then
      call messages_input_error(namespace, 'InvertKSmethod')
      call messages_fatal(1, namespace=namespace)
    end if
 
    !%Variable KSInversionLevel
    !%Type integer
    !%Default ks_inversion_adiabatic
    !%Section Calculation Modes::Invert KS
    !%Description
    !% At what level <tt>Octopus</tt> shall handle the KS inversion.
    !%Option ks_inversion_none 1
    !% Do not compute KS inversion.
    !%Option ks_inversion_adiabatic 2
    !% Compute exact adiabatic <math>v_{xc}</math>.
    !%End
    call messages_obsolete_variable(namespace, 'KS_Inversion_Level', 'KSInversionLevel')
    call parse_variable(namespace, 'KSInversionLevel', xc_ks_inversion_adiabatic, ks_inv%level)
    if (.not. varinfo_valid_option('KSInversionLevel', ks_inv%level)) call messages_input_error(namespace, 'KSInversionLevel')
 
    !%Variable KSInversionAsymptotics
    !%Type integer
    !%Default xc_asymptotics_none
    !%Section Calculation Modes::Invert KS
    !%Description
    !% Asymptotic correction applied to <math>v_{xc}</math>.
    !%Option xc_asymptotics_none 1
    !% Do not apply any correction in the asymptotic region.
    !%Option xc_asymptotics_sc 2
    !% Applies the soft-Coulomb decay of <math>-1/\sqrt{r^2+1}</math> to <math>v_{xc}</math> in the asymptotic region.
    !%End
    call parse_variable(namespace, 'KSInversionAsymptotics', xc_asymptotics_none, ks_inv%asymptotics)
    if (ks_inv%asymptotics /= xc_asymptotics_none .and. space%dim > 1) then
      call messages_not_implemented("KSInversionAsymptotics /= xc_asymptotics_sc for dimensions higher than 1.")
    end if
 
    ! In order to get properly converged states, given a vxc, we might need to run more than
    ! one eigensolver run.
    ! This variable is documented in scf.F90
    call parse_variable(namespace, 'MaximumIter', 50, ks_inv%max_iter)
 
    if (ks_inv%level /= xc_ks_inversion_none) then
      call states_elec_copy(ks_inv%aux_st, st, exclude_wfns = .true.)
 
      ! initialize auxiliary random wavefunctions
      call states_elec_allocate_wfns(ks_inv%aux_st, gr)
      call states_elec_generate_random(ks_inv%aux_st, gr, kpoints)
 
      ! initialize densities, hamiltonian and eigensolver
      call states_elec_densities_init(ks_inv%aux_st, gr)
 
      lasers => lasers_t(namespace)
      call lasers_parse_external_fields(lasers)
      call lasers_generate_potentials(lasers, gr, space, ions%latt)
 
      if(lasers%no_lasers > 0) then
        call ks_inv%ext_partners%add(lasers)
        call lasers_check_symmetries(lasers, kpoints)
      else
        deallocate(lasers)
      end if
 
      call hamiltonian_elec_init(ks_inv%aux_hm, namespace, space, gr, ions, ks_inv%ext_partners, ks_inv%aux_st, &
        kohn_sham_dft, xc, mc, kpoints)
      call eigensolver_init(ks_inv%eigensolver, namespace, gr, ks_inv%aux_st, mc, space)
    end if
 
    pop_sub(xc_ks_inversion_init)
  end subroutine xc_ks_inversion_init
 
 
  ! ---------------------------------------------------------
  subroutine xc_ks_inversion_end(ks_inv)
    type(xc_ks_inversion_t), intent(inout) :: ks_inv
 
    type(partner_iterator_t) :: iter
    class(interaction_partner_t), pointer :: partner
 
    push_sub(xc_ks_inversion_end)
 
    if (ks_inv%level /= xc_ks_inversion_none) then
      ! cleanup
      call eigensolver_end(ks_inv%eigensolver)
      call hamiltonian_elec_end(ks_inv%aux_hm)
      call states_elec_end(ks_inv%aux_st)
 
      call iter%start(ks_inv%ext_partners)
      do while (iter%has_next())
        partner => iter%get_next()
        safe_deallocate_p(partner)
      end do
      call ks_inv%ext_partners%empty()
 
    end if
 
    pop_sub(xc_ks_inversion_end)
  end subroutine xc_ks_inversion_end
 
 
  ! ---------------------------------------------------------
  subroutine xc_ks_inversion_write_info(ks_inversion, iunit, namespace)
    type(xc_ks_inversion_t),     intent(in) :: ks_inversion
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    if (ks_inversion%level == xc_ks_inversion_none) return
 
    push_sub(xc_ks_inversion_write_info)
 
    call messages_print_var_option('KSInversionLevel', ks_inversion%level, iunit=iunit, namespace=namespace)
 
    pop_sub(xc_ks_inversion_write_info)
  end subroutine xc_ks_inversion_write_info
 
 
  ! ---------------------------------------------------------
  subroutine invertks_2part(ks_inv, target_rho, nspin, aux_hm, gr, st, eigensolver, &
    namespace, space, ext_partners)
    type(xc_ks_inversion_t),  intent(in)    :: ks_inv
    real(real64),             intent(in)    :: target_rho(:,:)
    integer,                  intent(in)    :: nspin
    type(hamiltonian_elec_t), intent(inout) :: aux_hm
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(eigensolver_t),      intent(inout) :: eigensolver
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(partner_list_t),     intent(in)    :: ext_partners
 
    integer :: asym1, asym2, ip, ispin
    integer :: np
    real(real64)   :: rr, shift
    real(real64), parameter :: smalldensity = 5e-12_real64
    real(real64), allocatable :: sqrtrho(:,:), laplace(:,:), vks(:,:), x_shifted(:)
    type(spline_t) :: spl, lapl_spl
 
    push_sub(invertks_2part)
 
    assert(.not. gr%parallel_in_domains)
 
    ! In periodic systems, the aymptotics needs to done differently.
    ! In the KLI code, I used a formula from the litterature. The same logic could  used here too - NTD
    assert(space%periodic_dim == 0)
 
    ! The treatment of the asymptotics is only valid for the 1D case
    assert(space%dim == 1)
 
    np = gr%np
 
    safe_allocate(sqrtrho(1:gr%np_part, 1:nspin))
    safe_allocate(vks(1:np, 1:nspin))
    safe_allocate(laplace(1:gr%np, 1:nspin))
 
    if (any(target_rho(:,:) < -m_epsilon)) then
      write(message(1),*) "Target density has negative points. min value = ", minval(target_rho(:,:))
      call messages_warning(1, namespace=namespace)
    end if
 
    ! Computing the square root of the density
    !$omp parallel private(ip)
    do ispin = 1, nspin
      !$omp do simd
      do ip = 1, gr%np
        sqrtrho(ip, ispin) = sqrt(target_rho(ip, ispin))
      end do
      !$omp end do simd
    end do
    !$omp end parallel
 
    ! Computing the Laplacian of \sqrt{\rho}
    if (space%dim == 1) then
      ! In the 1D case, we fit by a spline and compute the derivative with it
      do ispin = 1, nspin
        call spline_init(spl)
        call spline_init(lapl_spl)
        call spline_fit(gr%np, gr%x(:, 1), sqrtrho(:, ispin), spl)
        ! Importantly, here we use a staggered grid to force interpolation
        ! Without this, the low density region shows large oscilations, as for the case of finite differences.
        ! To avoid this, we evaluate the derivative of the cubic spline on a different grid
        ! than the one used in the calculation, fit this result and reevalute (interpolate) it on the actual grid.
        ! This is needed as a cubic spline is by construction passing by the points on which it is fitted.
        ! This leads to a much smoother result for Laplacian of the density in the low-density regions.
 
        call spline_generate_shifted_grid(spl, x_shifted)
        call spline_der2(spl, lapl_spl, m_zero, grid=x_shifted)
        safe_deallocate_a(x_shifted)
 
 
        do ip = 1, gr%np
          laplace(ip, ispin) = spline_eval(lapl_spl, gr%x(ip, 1))
        end do
        call spline_end(spl)
        call spline_end(lapl_spl)
      end do
 
    else
      do ispin = 1, nspin
        call dderivatives_lapl(gr%der, sqrtrho(:, ispin), laplace(:, ispin))
      end do
    end if
 
    ! Used to check the asymptotics
    ! Note that this does not work for open-shell systems where the up and down density are different.
    ! This could cause problems for spin-polarized H atom for instance
    asym1 = 0
    asym2 = 0
 
    ! Note that the code below does not works for grid parallelization
    ! It also assumes a special ordering of the points, which which is only guarantied in 1D
    ! A better approach would be to fix the constant with the eigenvalue (user input option or restarting)
    do ispin = 1, nspin
      !avoid division by zero and set parameters for asymptotics
      !only for 1D potentials at the moment
      !need to find a way to find all points from where asymptotics should start in 2D and 3D
      do ip = 1, int(np/2)
        if (target_rho(ip, ispin) < smalldensity) then
          vks(ip, ispin) = aux_hm%ep%vpsl(ip) + aux_hm%vhartree(ip)
          asym1 = ip
        end if
        if (target_rho(np-ip+1, ispin) < smalldensity) then
          vks(np-ip+1, ispin) = aux_hm%ep%vpsl(np-ip+1) + aux_hm%vhartree(np-ip+1)
          asym2 = np - ip + 1
        end if
      end do
 
      ! The actual expression is \epsilon + 1/2 [\nabla^2 \sqrt{\rho}]/\sqrt{\rho} - v_{ext} - v_H
      ! We first compute the first part
      !$omp parallel do private(ip)
      do ip = asym1+1, asym2-1
        vks(ip, ispin) = m_half * laplace(ip, ispin) / (sqrtrho(ip, ispin) + m_tiny)
      end do
 
      ! And then compute the remaining part
      !$omp parallel do private(ip)
      do ip = 1, gr%np
        aux_hm%vxc(ip, ispin) = vks(ip, ispin) - aux_hm%ep%vpsl(ip) - aux_hm%vhartree(ip)
      end do
    end do
 
    !ensure correct asymptotic behavior, only for 1D potentials at the moment
    !need to find a way to find all points from where asymptotics should start in 2D and 3D
    if (ks_inv%asymptotics == xc_asymptotics_sc) then
      do ispin = 1, nspin
        !$omp parallel do private(rr, ip)
        do ip = 1, asym1
          call mesh_r(gr, ip, rr)
          aux_hm%vxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
        end do
        !$omp end parallel do
 
        ! calculate constant shift for correct asymptotics and shift accordingly
        call mesh_r(gr, asym1+1, rr)
        shift  = aux_hm%vxc(asym1+1, ispin) + st%qtot/sqrt(rr**2 + m_one)
 
        !$omp parallel do simd private(ip)
        do ip = asym1+1, asym2-1
          aux_hm%vxc(ip, ispin) = aux_hm%vxc(ip, ispin) - shift
        end do
        !$omp end parallel do simd
 
 
        call mesh_r(gr, asym2-1, rr)
        shift  = aux_hm%vxc(asym2-1, ispin) + st%qtot/sqrt(rr**2 + m_one)
 
        !$omp parallel private(rr, ip)
        !$omp do simd
        do ip = 1, asym2-1
          aux_hm%vxc(ip, ispin) = aux_hm%vxc(ip, ispin) - shift
        end do
        !$omp end do simd
 
        !$omp do
        do ip = asym2, np
          call mesh_r(gr, ip, rr)
          aux_hm%vxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
        end do
        !$omp end do
        !$omp end parallel
      end do
    end if !apply asymptotic correction
 
    ! Remove a shift, imposing vxc to be continuous
    if (ks_inv%asymptotics == xc_asymptotics_none) then
      assert(asym1 > 0)
      assert(asym2 > 0)
      do ispin = 1, nspin
        ! calculate constant shift to make potential continuous
 
        ! Note that this will fail if the density does not vanish below small_density
        shift  = aux_hm%vxc(asym1+1, ispin)
        !$omp parallel do simd
        do ip = asym1+1, asym2-1
          aux_hm%vxc(ip, ispin) = aux_hm%vxc(ip, ispin) - shift
        end do
        !$omp end parallel do simd
 
        shift  = aux_hm%vxc(asym2-1, ispin)
        !$omp parallel do simd
        do ip = 1, asym2-1
          aux_hm%vxc(ip, ispin) = aux_hm%vxc(ip, ispin) - shift
        end do
        !$omp end parallel do simd
      end do
    end if
 
    ! Store the final result
    do ispin = 1, nspin
      aux_hm%vhxc(1:gr%np, ispin) = aux_hm%vxc(1:gr%np, ispin) + aux_hm%vhartree(1:gr%np)
    end do
 
 
    call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
 
    safe_deallocate_a(sqrtrho)
    safe_deallocate_a(laplace)
    safe_deallocate_a(vks)
 
    pop_sub(invertks_2part)
  end subroutine invertks_2part
 
  subroutine invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, hm, max_iter)
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(in)    :: gr
    class(space_t),           intent(in)    :: space
    type(partner_list_t),     intent(in)    :: ext_partners
    type(eigensolver_t),      intent(inout) :: eigensolver
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: max_iter
 
    logical :: converged
    integer :: iter, nst_conv
 
    push_sub(invertks_update_hamiltonian)
 
    call hm%update(gr, namespace, space, ext_partners)
 
    !Poor-man estimate of the number of occupied states. Should be improved
    nst_conv = floor(st%qtot/st%smear%el_per_state)
 
    do iter = 1, max_iter
      call eigensolver%run(namespace, gr, st, hm, 1, converged, nst_conv)
 
      !TODO: Should we allow for updating the occupations? This might be necessary for open-shell or metallic systems
      ! call states_elec_fermi(st, namespace, gr)
 
      call write_iter()
 
      if (converged) exit
    end do
 
    if (converged) then
      message(1) = 'All states converged.'
      call messages_info(1, namespace=namespace)
    else
      message(1) = 'Some of the states are not fully converged!'
      call messages_info(1, namespace=namespace)
    end if
 
    write(message(1),'(a, e17.6)') 'Criterion = ', eigensolver%tolerance
    call messages_info(1, namespace=namespace)
    call messages_new_line()
 
    ! Get the new guess density
    call density_calc(st, gr, st%rho)
 
    pop_sub(invertks_update_hamiltonian)
 
  contains
    ! ---------------------------------------------------------
    subroutine write_iter()
      character(len=50) :: str
 
      push_sub(invertks_update_hamiltonian.write_iter)
 
      write(str, '(a,i5)') 'Kohn-Sham inversion eigensolver iteration #', iter
      call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
 
      write(message(1),'(a,i6,a,i6)') 'Converged states: ', minval(eigensolver%converged(1:st%nik))
      call messages_info(1, namespace=namespace)
 
      call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, &
        eigensolver%diff, compact = .true., namespace=namespace)
 
      call messages_print_with_emphasis(namespace=namespace)
 
      pop_sub(invertks_update_hamiltonian.write_iter)
    end subroutine write_iter
  end subroutine invertks_update_hamiltonian
 
  ! ---------------------------------------------------------
  subroutine invertks_iter(ks_inv, target_rho, namespace, space, ext_partners, nspin, aux_hm, gr, &
    st, eigensolver)
    type(xc_ks_inversion_t),  intent(in)    :: ks_inv
    real(real64),             intent(in)    :: target_rho(:,:)
    type(grid_t),             intent(in)    :: gr
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: aux_hm
    type(eigensolver_t),      intent(inout) :: eigensolver
    integer,                  intent(in)    :: nspin
 
    integer :: ip, ispin, ierr, asym1, asym2
    integer :: iunit, verbosity, counter, np
    integer :: max_iter
    integer :: imax
    real(real64) :: rr, shift
    real(real64) :: alpha, beta
    real(real64) :: mu, npower, npower_in ! these constants are from Rex Godbys scheme
    real(real64) :: convdensity, diffdensity
    real(real64), allocatable :: vhxc(:,:)
    real(real64), parameter :: small_density = 5e-12_real64
 
    character(len=256) :: fname
 
    push_sub(invertks_iter)
 
    ! The treatment of the asymptotics is only valid for the 1D case
    assert(space%dim == 1)
 
    np = gr%np
 
    !%Variable InvertKSConvAbsDens
    !%Type float
    !%Default 1e-5
    !%Section Calculation Modes::Invert KS
    !%Description
    !% Absolute difference between the calculated and the target density in the KS
    !% inversion. Has to be larger than the convergence of the density in the SCF run.
    !%End
    call parse_variable(namespace, 'InvertKSConvAbsDens', 1e-5_real64, convdensity)
 
    !%Variable InvertKSStellaBeta
    !%Type float
    !%Default 1.0
    !%Section Calculation Modes::Invert KS
    !%Description
    !% residual term in Stella iterative scheme to avoid 0 denominators
    !%End
    call parse_variable(namespace, 'InvertKSStellaBeta', 1e-6_real64, beta)
 
    !%Variable InvertKSStellaAlpha
    !%Type float
    !%Default 0.05
    !%Section Calculation Modes::Invert KS
    !%Description
    !% prefactor term in iterative scheme from L Stella
    !%End
    call parse_variable(namespace, 'InvertKSStellaAlpha', 0.25_real64, alpha)
 
    !%Variable InvertKSGodbyMu
    !%Type float
    !%Default 1.0
    !%Section Calculation Modes::Invert KS
    !%Description
    !% prefactor for iterative KS inversion convergence scheme from Godby based on van Leeuwen scheme
    !%End
    call parse_variable(namespace, 'InvertKSGodbyMu', m_one, mu)
 
    !%Variable InvertKSGodbyPower
    !%Type float
    !%Default 0.05
    !%Section Calculation Modes::Invert KS
    !%Description
    !% power to which density is elevated for iterative KS inversion convergence
    !% scheme from Godby based on van Leeuwen scheme
    !%End
    call parse_variable(namespace, 'InvertKSGodbyPower', 0.05_real64, npower_in)
    npower = npower_in
 
    !%Variable InvertKSVerbosity
    !%Type integer
    !%Default 0
    !%Section Calculation Modes::Invert KS
    !%Description
    !% Selects what is output during the calculation of the KS potential.
    !%Option 0
    !% Only outputs the converged density and KS potential.
    !%Option 1
    !% Same as 0 but outputs the maximum difference to the target density in each
    !% iteration in addition.
    !%Option 2
    !% Same as 1 but outputs the density and the KS potential in each iteration in
    !% addition.
    !%End
    call parse_variable(namespace, 'InvertKSVerbosity', 0, verbosity)
    if (verbosity < 0 .or. verbosity > 2) then
      call messages_input_error(namespace, 'InvertKSVerbosity')
      call messages_fatal(1, namespace=namespace)
    end if
 
    !%Variable InvertKSMaxIter
    !%Type integer
    !%Default 200
    !%Section Calculation Modes::Invert KS
    !%Description
    !% Selects how many iterations of inversion will be done in the iterative scheme
    !%End
    call parse_variable(namespace, 'InvertKSMaxIter', 200, max_iter)
 
    safe_allocate(vhxc(1:np, 1:nspin))
    call lalg_copy(np, nspin, aux_hm%vhxc, vhxc)
 
    if (verbosity == 1 .or. verbosity == 2) then
      iunit = io_open('InvertKSconvergence', namespace, action = 'write')
    end if
 
    counter = 0
    imax = 0
 
    diffdensity = m_zero
    do ispin = 1, nspin
      do ip = 1, np
        if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity) then
          diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
          imax = ispin
        end if
      end do
    end do
 
 
    do while(diffdensity > convdensity .and. counter < max_iter)
 
      counter = counter + 1
 
      if (verbosity == 2) then
        write(fname,'(i6.6)') counter
        call dio_function_output(io_function_fill_how("AxisX"), ".", "vhxc"//fname, namespace, space, &
          gr, aux_hm%vhxc(:,1), units_out%energy, ierr)
        call dio_function_output(io_function_fill_how("AxisX"), ".", "rho"//fname, namespace, space, &
          gr, st%rho(:,1), units_out%length**(-space%dim), ierr)
      end if
 
      ! Update the Hamiltonian
      call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
 
      ! Iterative inversion with fixed parameters in Stella Verstraete method
      select case(ks_inv%method)
      case (xc_inv_method_vs_iter)
        !$omp parallel private(ip)
        do ispin = 1, nspin
          !$omp do simd
          do ip = 1, np
            vhxc(ip, ispin) = vhxc(ip, ispin) &
              + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
          end do
          !$omp end do simd
        end do
        !$omp end parallel
 
        ! adaptative iterative method, with update of alpha and beta coefficients
        ! based on residual in density
      case (xc_inv_method_iter_stella)
        beta = diffdensity*1e-3_real64 !parameter to avoid numerical problems due to small denominator
 
        ! proposition to increase convergence speed progressively
        alpha = max(0.05_real64, m_half - diffdensity*100.0_real64*0.45_real64)
        write(message(1),'(a,2E15.4,3I8, 2E15.4)') &
          ' KSinversion: diffdensity, convdensity, imax, counter, max_iter, alpha, beta ', &
          diffdensity, convdensity, imax, counter, max_iter, alpha, beta
        call messages_info(1, namespace=namespace)
 
        !$omp parallel private(ip)
        do ispin = 1, nspin
          !$omp do simd
          do ip = 1, np
            vhxc(ip, ispin) = vhxc(ip, ispin) &
              + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
          end do
          !$omp end do simd
        end do
        !$omp end parallel
 
      case (xc_inv_method_iter_godby)
!        ! below 1.e-3 start reducing power down to 0.01
!        if (diffdensity < 0.001_real64) then
!          npower = min(npower_in, diffdensity*50.0_real64)
!        end if
        write(message(1),'(a,2E15.4,3I8, 2E15.4)') &
          ' KSinversion: diffdensity, convdensity, imax, counter, max_iter, power, mu ', &
          diffdensity, convdensity, imax, counter, max_iter, npower, mu
        call messages_info(1, namespace=namespace)
 
        !$omp parallel private(ip)
        do ispin = 1, nspin
          !$omp do simd
          do ip = 1, np
            vhxc(ip, ispin) = vhxc(ip, ispin) &
              + (st%rho(ip, ispin)**npower - target_rho(ip, ispin)**npower)*mu
          end do
          !$omp end do simd
        end do
        !$omp end parallel
 
      case default
        assert(.false.)
      end select
 
      diffdensity = m_zero
      !diffdensity = maxval(abs(st%rho(1:np,1:nspin)-target_rho(1:np,1:nspin)))
 
      ! Note: The code below requires a reduction over the grid.
      ! We now have routines for this, see mesh_minmaxloc - NTD
      assert(.not. gr%parallel_in_domains)
 
      do ispin = 1, nspin
        do ip = 1, np
          if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity) then
            diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
            imax = ispin
          end if
        end do
      end do
 
      if (verbosity == 1 .or. verbosity == 2) then
        write(iunit,'(i6.6)', advance = 'no') counter
        write(iunit,'(es18.10)') diffdensity
#ifdef HAVE_FLUSH
        call flush(iunit)
#endif
      end if
 
      ! Store the final result
      call lalg_copy(np, nspin, vhxc, aux_hm%vhxc)
      do ispin = 1, nspin
        aux_hm%vxc(1:np, ispin) = vhxc(1:np, ispin) - aux_hm%vhartree(1:np)
      end do
 
    end do ! end while statement on convergence
 
    !ensure correct asymptotic behavior, only for 1D potentials at the moment
    !need to find a way to find all points from where asymptotics should start in 2 and 3D
    if (ks_inv%asymptotics == xc_asymptotics_sc) then
      !Asympotics is only for vxc
      call lalg_copy(np, nspin, aux_hm%vxc, vhxc)
 
      do ispin = 1, nspin
        ! In case the density is never smaller than 'small_density', we take the edge points
        asym1 = 1
        asym2 = np
 
        do ip = 1, int(np/2)
          if (target_rho(ip, ispin) < small_density) then
            call mesh_r(gr, ip, rr)
            vhxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
            asym1 = ip
          end if
          if (target_rho(np-ip+1, ispin) < small_density) then
            asym2 = np - ip + 1
          end if
        end do
 
        ! calculate constant shift for correct asymptotics and shift accordingly
        call mesh_r(gr, asym1+1, rr)
        shift  = vhxc(asym1+1, ispin) + st%qtot/sqrt(rr**2 + m_one)
 
        !$omp parallel do simd
        do ip = asym1+1, asym2-1
          vhxc(ip, ispin) = vhxc(ip, ispin) - shift
        end do
        !$omp end parallel do simd
 
        call mesh_r(gr, asym2-1, rr)
        shift  = vhxc(asym2-1, ispin) + st%qtot/sqrt(rr**2 + m_one)
 
        !$omp parallel private(rr)
        !$omp do simd
        do ip = 1, asym2-1
          vhxc(ip, ispin) = vhxc(ip, ispin) - shift
        end do
        !$omp end do simd
 
        !$omp do
        do ip = asym2, np
          call mesh_r(gr, ip, rr)
          vhxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
        end do
        !$omp end do
        !$omp end parallel
      end do
 
      ! See the above comment, vhxc contains vxc here
      call lalg_copy(np, nspin, vhxc, aux_hm%vxc)
      do ispin = 1, nspin
        aux_hm%vhxc(1:np, ispin) = vhxc(1:np, ispin) + aux_hm%vhartree(1:np)
      end do
    end if
 
    !TODO: check that all arrays needed by hamiltonian update are sync`d in MPI fashion
 
    !calculate final density
    call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
 
    write(message(1),'(a,I8)') "Iterative KS inversion, iterations needed:", counter
    call messages_info(1, namespace=namespace)
 
    call io_close(iunit)
 
    safe_deallocate_a(vhxc)
 
    pop_sub(invertks_iter)
  end subroutine invertks_iter
 
  ! ---------------------------------------------------------
  subroutine xc_ks_inversion_calc(ks_inversion, namespace, space, gr, hm, ext_partners, st, vxc, time)
    type(xc_ks_inversion_t),  intent(inout) :: ks_inversion
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    real(real64),             intent(inout) :: vxc(:,:)
    real(real64), optional,   intent(in)    :: time
 
    integer :: ispin
 
    if (ks_inversion%level == xc_ks_inversion_none) return
 
    push_sub(x(xc_ks_inversion_calc))
 
    if (present(time) .and. time > m_zero) then
      write(message(1),'(A,F18.12)') 'xc_ks_inversion_calc - time:', time
      call messages_info(1, namespace=namespace)
    end if
 
    ks_inversion%aux_hm%energy%intnvxc     = m_zero
    ks_inversion%aux_hm%energy%hartree     = m_zero
    ks_inversion%aux_hm%energy%exchange    = m_zero
    ks_inversion%aux_hm%energy%correlation = m_zero
 
    ks_inversion%aux_hm%vhartree = hm%vhartree
 
    if (present(time) .and. time > m_zero) then
      call lalg_copy(gr%np, st%d%nspin, ks_inversion%vhxc_previous_step, ks_inversion%aux_hm%vhxc)
 
    else
! no restart data available, start with vhxc = vh, which we know from exact input rho
      do ispin = 1, st%d%nspin
        ks_inversion%aux_hm%vxc(1:gr%np, ispin)  = m_zero !hm%ep%vpsl(:)
        ks_inversion%aux_hm%vhxc(1:gr%np, ispin) = hm%vhartree(1:gr%np) + ks_inversion%aux_hm%vxc(1:gr%np, ispin)
      end do
! TODO: restart data found. Use first KS orbital to invert equation and get starting vhxc
!      call invertks_2part(ks_inversion%aux_st%rho, st%d%nspin, ks_inversion%aux_hm, gr, &
!                         ks_inversion%aux_st, ks_inversion%eigensolver, namespace, ks_inversion%asymp)
 
    end if
    !ks_inversion%aux_hm%ep%vpsl(:)  = M_ZERO ! hm%ep%vpsl(:)
    call lalg_copy(gr%np, hm%ep%vpsl, ks_inversion%aux_hm%ep%vpsl)
 
    ! compute ks inversion, vhxc contains total KS potential
    ! Update the Hamiltonian
    call invertks_update_hamiltonian(namespace, gr, space, ext_partners, &
      ks_inversion%eigensolver, ks_inversion%aux_st, ks_inversion%aux_hm, ks_inversion%max_iter)
 
 
    ! these 2 routines need to be cleaned - they are not consistent in updating
    ! the hamiltonian, states, etc...
    select case (ks_inversion%method)
      ! adiabatic ks inversion
    case (xc_inv_method_two_particle)
      call invertks_2part(ks_inversion, ks_inversion%aux_st%rho, st%d%nspin, ks_inversion%aux_hm, gr, &
        ks_inversion%aux_st, ks_inversion%eigensolver, namespace, space, ext_partners)
 
    case (xc_inv_method_vs_iter : xc_inv_method_iter_godby)
      call invertks_iter(ks_inversion, st%rho, namespace, space, ext_partners, st%d%nspin, ks_inversion%aux_hm, gr, &
        ks_inversion%aux_st, ks_inversion%eigensolver)
    end select
 
    ! subtract Hartree potential
    ! ATTENTION: subtracts true external potential not adiabatic one
 
    do ispin = 1, st%d%nspin
      call lalg_axpy(gr%np, -m_one, hm%vhartree, ks_inversion%aux_hm%vhxc(:,ispin))
    end do
 
    vxc = ks_inversion%aux_hm%vxc
 
    ! save vhxc for next step if we are running td
    if (present(time)) then
      ks_inversion%vhxc_previous_step = ks_inversion%aux_hm%vhxc
    end if
 
    pop_sub(x(xc_ks_inversion_calc))
 
  end subroutine xc_ks_inversion_calc
 
 
end module xc_ks_inversion_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: