ks_potential.F90

Go to the documentation of this file.

!! Copyright (C) 2024 N. Tancogne-Dejean
!!
!! 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 ks_potential_oct_m
  use accel_oct_m
  use atom_oct_m
  use debug_oct_m
  use global_oct_m
  use hamiltonian_elec_base_oct_m
  use io_function_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mix_oct_m
  use namespace_oct_m
  use output_low_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_dim_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                   &
    ks_potential_t,           &
    xc_copied_potentials_t,   &
    xc_copied_potentials_end, &
    vtau_set_vout,            &
    vtau_set_vin,             &
    vtau_get_vnew
 
  integer, public, parameter ::        &
    INDEPENDENT_PARTICLES = 2,         &
    hartree               = 1,         &
    hartree_fock          = 3,         &
    kohn_sham_dft         = 4,         &
    generalized_kohn_sham_dft = 5,     &
    rdmft                 = 7
 
 
  type ks_potential_t
    private
    integer :: np, nspin
    integer :: theory_level = independent_particles
    logical :: needs_vtau = .false.
 
    real(real64), public, allocatable :: vhartree(:)
    real(real64), public, allocatable :: vxc(:,:)
    real(real64), public, allocatable :: vhxc(:,:)
    !                                                !! other contributions: photon xc potential, constrained DFT, ...
    !                                                !! Check v_ks_calc for more details
    real(real64), allocatable :: vtau(:,:)
 
    type(accel_mem_t) :: vtau_accel
 
  contains
    procedure :: init => ks_potential_init
    procedure :: end => ks_potential_end
    procedure :: set_vtau => ks_potential_set_vtau
    procedure :: add_vhxc => ks_potential_add_vhxc
    procedure :: dmult_vhxc => dks_potential_mult_vhxc
    procedure :: zmult_vhxc => zks_potential_mult_vhxc
    procedure :: load => ks_potential_load_vhxc
    procedure :: dump => ks_potential_dump_vhxc
    procedure :: init_interpolation => ks_potential_init_interpolation
    procedure :: run_zero_iter => ks_potential_run_zero_iter
    procedure :: interpolation_new => ks_potential_interpolation_new
    procedure :: get_interpolated_potentials => ks_potential_get_interpolated_potentials
    procedure :: set_interpolated_potentials => ks_potential_set_interpolated_potentials
    procedure :: interpolate_potentials => ks_potential_interpolate_potentials
    procedure :: dapply_vtau_psi => dks_potential_apply_vtau_psi
    procedure :: zapply_vtau_psi => zks_potential_apply_vtau_psi
    procedure :: dcurrent_mass_renormalization => dks_potential_current_mass_renormalization
    procedure :: zcurrent_mass_renormalization => zks_potential_current_mass_renormalization
    procedure :: output_potentials => ks_potential_output_potentials
    procedure :: store_potentials => ks_potential_store_copy
    procedure :: restore_potentials => ks_potential_restore_copy
    procedure :: check_convergence => ks_potential_check_convergence
    procedure :: perform_interpolation => ks_potential_perform_interpolation
    procedure :: mix_potentials => ks_potential_mix_potentials
  end type ks_potential_t
 
  ! We sometimes need to store a copy of  the potentials, say for interpolation
  ! This is used to store the corresponding data, still keeping them private
  type xc_copied_potentials_t
    private
    real(real64), public, allocatable :: vhxc(:,:)
    real(real64), allocatable :: vtau(:,:)
  contains
    procedure :: copy_vhxc_to_buffer => xc_copied_potentials_copy_vhxc_to_buffer
    final :: xc_copied_potentials_end
  end type xc_copied_potentials_t
 
contains
 
  subroutine ks_potential_init(this, np, np_part, nspin, theory_level, needs_vtau)
    class(ks_potential_t), intent(inout) :: this
    integer,               intent(in)    :: np, np_part
    integer,               intent(in)    :: nspin
    integer,               intent(in)    :: theory_level
    logical,               intent(in)    :: needs_vtau
 
    push_sub(ks_potential_init)
 
    this%theory_level = theory_level
    this%needs_vtau = needs_vtau
    this%np = np
    this%nspin = nspin
 
    ! In the case of spinors, vxc_11 = hm%vxc(:, 1), vxc_22 = hm%vxc(:, 2), Re(vxc_12) = hm%vxc(:. 3);
    ! Im(vxc_12) = hm%vxc(:, 4)
    safe_allocate(this%vhxc(1:np, 1:nspin))
    this%vhxc(1:np, 1:nspin) = m_zero
 
    if (theory_level /= independent_particles) then
 
      safe_allocate(this%vhartree(1:np_part))
      this%vhartree=m_zero
 
      safe_allocate(this%vxc(1:np, 1:nspin))
      this%vxc=m_zero
 
      if (needs_vtau) then
        safe_allocate(this%vtau(1:np, 1:nspin))
        this%vtau=m_zero
        if (accel_is_enabled()) then
          call accel_create_buffer(this%vtau_accel, accel_mem_read_only, type_float, accel_padded_size(np)*nspin)
        end if
 
      end if
    end if
 
    pop_sub(ks_potential_init)
  end subroutine ks_potential_init
 
  subroutine ks_potential_end(this)
    class(ks_potential_t), intent(inout) :: this
 
    push_sub(ks_potential_end)
 
    safe_deallocate_a(this%vhxc)
    safe_deallocate_a(this%vhartree)
    safe_deallocate_a(this%vxc)
    safe_deallocate_a(this%vtau)
 
    if (this%needs_vtau .and. accel_is_enabled()) then
      call accel_release_buffer(this%vtau_accel)
    end if
 
    pop_sub(ks_potential_end)
  end subroutine ks_potential_end
 
  subroutine ks_potential_set_vtau(this, vtau)
    class(ks_potential_t), intent(inout) :: this
    real(real64),          intent(in)    :: vtau(:,:)
 
    push_sub(ks_potential_set_vtau)
 
    assert(this%needs_vtau)
 
    call lalg_copy(this%np, this%nspin, vtau, this%vtau)
 
    call ks_potential_update_vtau_buffer(this)
 
    pop_sub(ks_potential_set_vtau)
  end subroutine ks_potential_set_vtau
 
  subroutine ks_potential_update_vtau_buffer(this)
    class(ks_potential_t), intent(inout) :: this
 
    integer :: offset, ispin
 
    push_sub(ks_potential_update_vtau_buffer)
 
    assert(this%needs_vtau)
 
    if (accel_is_enabled()) then
      offset = 0
      do ispin = 1, this%nspin
        call accel_write_buffer(this%vtau_accel, this%np, this%vtau(:, ispin), offset = offset)
        offset = offset + accel_padded_size(this%np)
      end do
    end if
 
    pop_sub(ks_potential_update_vtau_buffer)
  end subroutine ks_potential_update_vtau_buffer
 
  subroutine ks_potential_add_vhxc(this, pot, nspin)
    class(ks_potential_t), intent(in)    :: this
    real(real64),          intent(inout) :: pot(:,:)
    integer, optional,     intent(in)    :: nspin
 
    integer :: ispin, ip, nspin_
 
    push_sub(ks_potential_add_vhxc)
 
    assert(not_in_openmp())
    nspin_ = optional_default(nspin, this%nspin)
    assert(size(pot, dim=1) >= this%np)
    assert(size(pot, dim=2) >= nspin_)
 
    !TODO: We could skip this for IPA, and also use axpy here
 
    !$omp parallel private(ip, ispin)
    do ispin = 1, nspin_
      !$omp do simd schedule(static)
      do ip = 1, this%np
        pot(ip, ispin) = pot(ip, ispin) + this%vhxc(ip, ispin)
      end do
      !$omp end do simd nowait
    end do
    !$omp end parallel
 
    pop_sub(ks_potential_add_vhxc)
  end subroutine ks_potential_add_vhxc
 
  ! -----------------------------------------------------------------
  subroutine ks_potential_dump_vhxc(this, restart, space, mesh, ierr)
    class(ks_potential_t), intent(in)  :: this
    type(restart_t),       intent(in)  :: restart
    class(space_t),        intent(in)  :: space
    class(mesh_t),         intent(in)  :: mesh
    integer,               intent(out) :: ierr
 
    integer :: iunit, err, err2(2), isp
    character(len=MAX_PATH_LEN) :: filename
    character(len=100) :: lines(2)
 
    push_sub(ks_potential_dump_vhxc)
 
    ierr = 0
 
    if (restart_skip(restart) .or. this%theory_level == independent_particles) then
      pop_sub(ks_potential_dump_vhxc)
      return
    end if
 
    if (debug%info) then
      message(1) = "Debug: Writing Vhxc restart."
      call messages_info(1)
    end if
 
    !write the different components of the Hartree+XC potential
    iunit = restart_open(restart, 'vhxc')
    lines(1) = '#     #spin    #nspin    filename'
    lines(2) = '%vhxc'
    call restart_write(restart, iunit, lines, 2, err)
    if (err /= 0) ierr = ierr + 1
 
    err2 = 0
    do isp = 1, this%nspin
      filename = get_filename_with_spin('vhxc', this%nspin, isp)
      write(lines(1), '(i8,a,i8,a)') isp, ' | ', this%nspin, ' | "'//trim(adjustl(filename))//'"'
      call restart_write(restart, iunit, lines, 1, err)
      if (err /= 0) err2(1) = err2(1) + 1
 
      call restart_write_mesh_function(restart, space, filename, mesh, this%vhxc(:,isp), err)
      if (err /= 0) err2(2) = err2(2) + 1
 
    end do
    if (err2(1) /= 0) ierr = ierr + 2
    if (err2(2) /= 0) ierr = ierr + 4
 
    lines(1) = '%'
    call restart_write(restart, iunit, lines, 1, err)
    if (err /= 0) ierr = ierr + 4
 
    ! MGGAs and hybrid MGGAs have an extra term that also needs to be dumped
    if (this%needs_vtau) then
      lines(1) = '#     #spin    #nspin    filename'
      lines(2) = '%vtau'
      call restart_write(restart, iunit, lines, 2, err)
      if (err /= 0) ierr = ierr + 8
 
      err2 = 0
      do isp = 1, this%nspin
        filename = get_filename_with_spin('vtau', this%nspin, isp)
        write(lines(1), '(i8,a,i8,a)') isp, ' | ', this%nspin, ' | "'//trim(adjustl(filename))//'"'
        call restart_write(restart, iunit, lines, 1, err)
        if (err /= 0) err2(1) = err2(1) + 16
 
        call restart_write_mesh_function(restart, space, filename, mesh, this%vtau(:,isp), err)
        if (err /= 0) err2(1) = err2(1) + 1
 
      end do
      if (err2(1) /= 0) ierr = ierr + 32
      if (err2(2) /= 0) ierr = ierr + 64
 
      lines(1) = '%'
      call restart_write(restart, iunit, lines, 1, err)
      if (err /= 0) ierr = ierr + 128
    end if
 
    call restart_close(restart, iunit)
 
    if (debug%info) then
      message(1) = "Debug: Writing Vhxc restart done."
      call messages_info(1)
    end if
 
    pop_sub(ks_potential_dump_vhxc)
  end subroutine ks_potential_dump_vhxc
 
 
  ! ---------------------------------------------------------
  subroutine ks_potential_load_vhxc(this, restart, space, mesh, ierr)
    class(ks_potential_t), intent(inout) :: this
    type(restart_t),       intent(in)    :: restart
    class(space_t),        intent(in)    :: space
    class(mesh_t),         intent(in)    :: mesh
    integer,               intent(out)   :: ierr
 
    integer :: err, err2, isp
    character(len=MAX_PATH_LEN) :: filename
 
    push_sub(ks_potential_load_vhxc)
 
    ierr = 0
 
    if (restart_skip(restart) .or. this%theory_level == independent_particles) then
      ierr = -1
      pop_sub(ks_potential_load_vhxc)
      return
    end if
 
    if (debug%info) then
      message(1) = "Debug: Reading Vhxc restart."
      call messages_info(1)
    end if
 
    err2 = 0
    do isp = 1, this%nspin
      filename = get_filename_with_spin('vhxc', this%nspin, isp)
 
      call restart_read_mesh_function(restart, space, filename, mesh, this%vhxc(:,isp), err)
      if (err /= 0) err2 = err2 + 1
 
    end do
    if (err2 /= 0) ierr = ierr + 1
 
    ! MGGAs and hybrid MGGAs have an extra term that also needs to be read
    err2 = 0
    if (this%needs_vtau) then
      do isp = 1, this%nspin
        filename = get_filename_with_spin('vtau', this%nspin, isp)
 
        call restart_read_mesh_function(restart, space, filename, mesh, this%vtau(:,isp), err)
        if (err /= 0) err2 = err2 + 1
 
      end do
 
      if (err2 /= 0) ierr = ierr + 2
    end if
 
    if (debug%info) then
      message(1) = "Debug: Reading Vhxc restart done."
      call messages_info(1)
    end if
 
    pop_sub(ks_potential_load_vhxc)
  end subroutine ks_potential_load_vhxc
 
  subroutine ks_potential_init_interpolation(this, vksold, order)
    class(ks_potential_t),           intent(in)    :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    integer, optional,               intent(in)    :: order
 
    push_sub(ks_potential_init_interpolation)
 
    call potential_interpolation_init(vksold, this%np, this%nspin, this%needs_vtau, order=order)
 
    pop_sub(ks_potential_init_interpolation)
  end subroutine ks_potential_init_interpolation
 
  subroutine ks_potential_run_zero_iter(this, vksold)
    class(ks_potential_t),           intent(in)    :: this
    type(potential_interpolation_t), intent(inout) :: vksold
 
    push_sub(ks_potential_run_zero_iter)
 
    call potential_interpolation_run_zero_iter(vksold, this%np, this%nspin, this%vhxc, vtau=this%vtau)
 
    pop_sub(ks_potential_run_zero_iter)
  end subroutine ks_potential_run_zero_iter
 
  subroutine ks_potential_interpolation_new(this, vksold, current_time, dt)
    class(ks_potential_t),           intent(inout) :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    real(real64),                    intent(in)    :: current_time
    real(real64),                    intent(in)    :: dt
 
    push_sub(ks_potential_interpolation_new)
 
    call potential_interpolation_new(vksold, this%np, this%nspin, current_time, dt, this%vhxc, vtau=this%vtau)
 
    pop_sub(ks_potential_interpolation_new)
  end subroutine ks_potential_interpolation_new
 
  subroutine ks_potential_get_interpolated_potentials(this, vksold, history, storage)
    class(ks_potential_t),           intent(inout) :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    integer,                         intent(in)    :: history
    type(xc_copied_potentials_t),  optional, intent(inout) :: storage
 
    if (this%theory_level == independent_particles) return
 
    push_sub(ks_potential_get_interpolated_potentials)
 
    if (.not. present(storage)) then
      if (this%needs_vtau) then
        call potential_interpolation_get(vksold, this%np, this%nspin, history, this%vhxc, vtau = this%vtau)
        call ks_potential_update_vtau_buffer(this)
      else
        call potential_interpolation_get(vksold, this%np, this%nspin, history, this%vhxc)
      end if
    else
      call ks_potential_storage_allocate(this, storage)
      if (this%needs_vtau) then
        call potential_interpolation_get(vksold, this%np, this%nspin, history, storage%vhxc, vtau = storage%vtau)
      else
        call potential_interpolation_get(vksold, this%np, this%nspin, history, storage%vhxc)
      end if
    end if
 
    pop_sub(ks_potential_get_interpolated_potentials)
  end subroutine ks_potential_get_interpolated_potentials
 
  subroutine ks_potential_set_interpolated_potentials(this, vksold, history)
    class(ks_potential_t),           intent(inout) :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    integer,                         intent(in)    :: history
 
    if (this%theory_level == independent_particles) return
 
    push_sub(ks_potential_set_interpolated_potentials)
 
    if (this%needs_vtau) then
      call potential_interpolation_set(vksold, this%np, this%nspin, history, this%vhxc, vtau = this%vtau)
    else
      call potential_interpolation_set(vksold, this%np, this%nspin, history, this%vhxc)
    end if
 
    pop_sub(ks_potential_set_interpolated_potentials)
  end subroutine ks_potential_set_interpolated_potentials
 
 
  subroutine ks_potential_interpolate_potentials(this, vksold, order, current_time, dt, interpolation_time)
    class(ks_potential_t),           intent(inout) :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    integer,                         intent(in)    :: order
    real(real64),                    intent(in)    :: current_time
    real(real64),                    intent(in)    :: dt
    real(real64),                    intent(in)    :: interpolation_time
 
    if (this%theory_level == independent_particles) return
 
    push_sub(ks_potential_interpolate_potentials)
 
    if (this%needs_vtau) then
      call potential_interpolation_interpolate(vksold, order, current_time, dt, interpolation_time, &
        this%vhxc, vtau = this%vtau)
      call ks_potential_update_vtau_buffer(this)
    else
      call potential_interpolation_interpolate(vksold, order, current_time, dt, interpolation_time, &
        this%vhxc)
    end if
 
    pop_sub(ks_potential_interpolate_potentials)
  end subroutine ks_potential_interpolate_potentials
 
  ! ---------------------------------------------------------
  subroutine vtau_set_vout(field, this)
    type(mixfield_t),     intent(inout) :: field
    type(ks_potential_t), intent(in)    :: this
 
    push_sub(vtau_set_vout)
 
    call mixfield_set_vout(field, this%vtau)
 
    pop_sub(vtau_set_vout)
  end subroutine vtau_set_vout
 
  ! ---------------------------------------------------------
  subroutine vtau_set_vin(field, this)
    type(mixfield_t),     intent(inout) :: field
    type(ks_potential_t), intent(in)    :: this
 
    push_sub(vtau_set_vin)
 
    call mixfield_set_vin(field, this%vtau)
 
    pop_sub(vtau_set_vin)
  end subroutine vtau_set_vin
 
  ! ---------------------------------------------------------
  subroutine vtau_get_vnew(field, this)
    type(mixfield_t),     intent(in)    :: field
    type(ks_potential_t), intent(inout) :: this
 
    push_sub(vtau_get_vnew)
 
    call mixfield_get_vnew(field, this%vtau)
 
    call ks_potential_update_vtau_buffer(this)
 
    pop_sub(vtau_get_vnew)
  end subroutine vtau_get_vnew
 
  subroutine ks_potential_output_potentials(this, namespace, how, dir, space, mesh, pos, atoms, grp)
    class(ks_potential_t), intent(in)    :: this
    type(namespace_t),     intent(in)    :: namespace
    integer(int64),        intent(in)    :: how
    character(len=*),      intent(in)    :: dir
    class(space_t),        intent(in)    :: space
    class(mesh_t),         intent(in)    :: mesh
    real(real64),     optional, intent(in)  :: pos(:,:)
    type(atom_t),     optional, intent(in)  :: atoms(:)
    type(mpi_grp_t),  optional, intent(in)  :: grp
 
    character(len=MAX_PATH_LEN) :: fname
    integer :: is, err
 
    push_sub(ks_potential_output_potentials)
 
    call dio_function_output(how, dir, 'vh', namespace, &
      space, mesh, this%vhartree, units_out%energy, err, pos=pos, atoms=atoms, grp = grp)
 
    do is = 1, this%nspin
      fname = get_filename_with_spin('vxc', this%nspin, is)
      call dio_function_output(how, dir, fname, namespace, space, &
        mesh, this%vxc(:, is), units_out%energy, err, pos=pos, atoms=atoms, grp = grp)
 
      if (this%needs_vtau) then
        fname = get_filename_with_spin('vtau', this%nspin, is)
        call dio_function_output(how, dir, fname, namespace, space, &
          mesh, this%vtau(:, is), units_out%energy, err, pos=pos, atoms=atoms, grp = grp)
      end if
    end do
 
    pop_sub(ks_potential_output_potentials)
  end subroutine ks_potential_output_potentials
 
  subroutine ks_potential_storage_allocate(this, copy)
    class(ks_potential_t),       intent(in) :: this
    type(xc_copied_potentials_t),  intent(inout) :: copy
 
    push_sub(ks_potential_storage_allocate)
 
    if (.not. allocated(copy%vhxc)) then
      safe_allocate(copy%vhxc(1:this%np, 1:this%nspin))
    end if
    if (this%needs_vtau) then
      if (.not. allocated(copy%vtau)) then
        safe_allocate(copy%vtau(1:this%np, 1:this%nspin))
      end if
    end if
 
    pop_sub(ks_potential_storage_allocate)
  end subroutine ks_potential_storage_allocate
 
  subroutine ks_potential_store_copy(this, copy)
    class(ks_potential_t),       intent(in) :: this
    type(xc_copied_potentials_t),  intent(inout) :: copy
 
    if (this%theory_level == independent_particles) return
 
    push_sub(ks_potential_store_copy)
 
    call ks_potential_storage_allocate(this, copy)
    call lalg_copy(this%np, this%nspin, this%vhxc, copy%vhxc)
    if (this%needs_vtau) then
      call lalg_copy(this%np,  this%nspin, this%vtau, copy%vtau)
    end if
 
    pop_sub(ks_potential_store_copy)
  end subroutine ks_potential_store_copy
 
  subroutine ks_potential_restore_copy(this, copy)
    class(ks_potential_t),       intent(inout) :: this
    type(xc_copied_potentials_t),  intent(in) :: copy
 
    if (this%theory_level == independent_particles) return
 
    push_sub(ks_potential_restore_copy)
 
    assert(allocated(copy%vhxc))
    call lalg_copy(this%np, this%nspin, copy%vhxc, this%vhxc)
    if (this%needs_vtau) then
      assert(allocated(copy%vtau))
      call this%set_vtau(copy%vtau)
    end if
 
    pop_sub(ks_potential_restore_copy)
  end subroutine ks_potential_restore_copy
 
  subroutine xc_copied_potentials_copy_vhxc_to_buffer(this, np, nspin, pnp, buffer)
    class(xc_copied_potentials_t), intent(in)    :: this
    integer(int64),              intent(in)    :: np
    integer,                     intent(in)    :: nspin
    integer(int64),              intent(in)    :: pnp
    type(accel_mem_t),           intent(inout) :: buffer
 
    integer :: ispin
 
    push_sub(xc_copied_potentials_copy_vhxc_to_buffer)
 
    do ispin = 1, nspin
      call accel_write_buffer(buffer, np, this%vhxc(:, ispin), offset = (ispin - 1)*pnp)
    end do
 
    pop_sub(xc_copied_potentials_copy_vhxc_to_buffer)
  end subroutine xc_copied_potentials_copy_vhxc_to_buffer
 
  subroutine xc_copied_potentials_end(this)
    type(xc_copied_potentials_t), intent(inout) :: this
 
    safe_deallocate_a(this%vhxc)
    safe_deallocate_a(this%vtau)
  end subroutine xc_copied_potentials_end
 
  real(real64) function ks_potential_check_convergence(this, copy, mesh, rho, qtot) result(diff)
    class(ks_potential_t),       intent(in) :: this
    type(xc_copied_potentials_t),  intent(in) :: copy
    class(mesh_t),               intent(in) :: mesh
    real(real64),                intent(in) :: rho(:,:)
    real(real64),                intent(in) :: qtot
 
    real(real64), allocatable :: diff_pot(:)
    integer :: ip, is
 
    push_sub(ks_potential_check_convergence)
 
    safe_allocate(diff_pot(1:this%np))
 
    !$omp parallel
    !$omp do
    do ip = 1, this%np
      diff_pot(ip) = abs(copy%vhxc(ip, 1) - this%vhxc(ip, 1))*rho(ip, 1)
    end do
    do is = 2, this%nspin
      !$omp do
      do ip = 1, this%np
        diff_pot(ip) = diff_pot(ip) + abs(copy%vhxc(ip, is) - this%vhxc(ip, is))*rho(ip, is)
      end do
    end do
    !$omp end parallel
    diff = dmf_integrate(mesh, diff_pot) / qtot
 
    safe_deallocate_a(diff_pot)
 
    pop_sub(ks_potential_check_convergence)
  end function ks_potential_check_convergence
 
  subroutine ks_potential_perform_interpolation(this, vksold, times, current_time)
    class(ks_potential_t),           intent(inout) :: this
    type(potential_interpolation_t), intent(inout) :: vksold
    real(real64),                    intent(in)    :: times(:)
    real(real64),                    intent(in)    :: current_time
 
    push_sub(ks_potential_perform_interpolation)
 
    assert(size(times) == 3)
 
    call interpolate(times, vksold%v_old(:, :, 1:3), current_time, vksold%v_old(:, :, 0))
    if (this%needs_vtau) then
      call interpolate(times, vksold%vtau_old(:, :, 1:3), current_time, vksold%vtau_old(:, :, 0))
    end if
 
    pop_sub(ks_potential_perform_interpolation)
  end subroutine ks_potential_perform_interpolation
 
  subroutine ks_potential_mix_potentials(this, vold, dt)
    class(ks_potential_t),       intent(in)    :: this
    type(xc_copied_potentials_t),  intent(inout) :: vold
    real(real64),                intent(in)    :: dt
 
    integer :: ispin, ip
 
    push_sub(ks_potential_mix_potentials)
 
    assert(not_in_openmp())
 
    !$omp parallel private(ip)
    do ispin = 1, this%nspin
      !$omp do simd
      do ip = 1, this%np
        vold%vhxc(ip, ispin) =  m_half*dt*(this%vhxc(ip, ispin) - vold%vhxc(ip, ispin))
      end do
 
      if (this%needs_vtau) then
        !$omp do simd
        do ip = 1, this%np
          vold%vtau(ip, ispin) =  m_half*dt*(this%vtau(ip, ispin) - vold%vtau(ip, ispin))
        end do
      end if
    end do
    !$omp end parallel
 
    pop_sub(ks_potential_mix_potentials)
  end subroutine ks_potential_mix_potentials
 
 
#include "undef.F90"
#include "complex.F90"
#include "ks_potential_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "ks_potential_inc.F90"
 
end module ks_potential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: