mix.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2014 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!! Copyright (C) 2024 A. Buccheri
!!
!! 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 mix_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use nl_operator_oct_m
  use namespace_oct_m
  use parser_oct_m
  use mixing_preconditioner_oct_m
  use profiling_oct_m
  use restart_oct_m
  use solvers_oct_m
  use space_oct_m
  use stencil_cube_oct_m
  use types_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                     &
    mix_t,                      &
    mix_init,                   &
    mix_clear,                  &
    mix_end,                    &
    mix_dump,                   &
    mix_load,                   &
    mixing,                     &
    dmixing,                    &
    zmixing,                    &
    mix_coefficient,            &
    mix_scheme,                 &
    mix_d3,                     &
    mix_get_field,              &
    mixfield_t,                 &
    mixfield_init,              &
    mixfield_clear,             &
    mixfield_end,               &
    mixfield_set_vin,           &
    mixfield_set_vout,          &
    mixfield_get_vnew,          &
    mix_add_auxmixfield
 
  ! A better design would be to have mixfield_t and mixfield_ptr_t separate from mix_t:
  ! mix_t would contain constants, and should be immutable with only getters.
  ! mixfield_type encapsulates the data used in mixing, and has bound procedures for
  ! operating on this data.
 
  type mixfield_t
    private
    ! Real versions of the data
    real(real64), allocatable :: ddf(:, :, :)
    real(real64), allocatable :: ddv(:, :, :)
    real(real64), allocatable :: df_old(:, :)
    real(real64), allocatable :: dvin_old(:, :)
    real(real64), allocatable :: dvin(:, :)
    real(real64), allocatable :: dvout(:, :)
    real(real64), allocatable :: dvnew(:, :)
    real(real64), allocatable :: dresidual(:, :)
 
    ! Complex versions of the data
    complex(real64), allocatable :: zdf(:, :, :)
    complex(real64), allocatable :: zdv(:, :, :)
    complex(real64), allocatable :: zf_old(:, :)
    complex(real64), allocatable :: zvin_old(:, :)
    complex(real64), allocatable :: zvin(:, :)
    complex(real64), allocatable :: zvout(:, :)
    complex(real64), allocatable :: zvnew(:, :)
    complex(real64), allocatable :: zresidual(:, :)
 
    type(type_t) :: func_type
    integer :: d1, d2, d3
    !                          d1 is mesh%np, d2 vary depending on where this called
    !                          d3 is for SCF history of a quantity.
    logical :: mix_spin_density_matrix
    real(real64) :: last_residue
  end type mixfield_t
 
  type mixfield_ptr_t
    private
    type(mixfield_t), pointer :: p
  end type mixfield_ptr_t
 
  integer, parameter :: MAX_AUXMIXFIELD = 5
 
  type mix_t
    private
    type(derivatives_t), pointer :: der
 
    integer :: scheme
    real(real64), allocatable  :: coeff(:)
    real(real64)   :: adaptive_damping
    real(real64)   :: residual_coeff
    integer :: iter
    !                             !< include the iterations done in previous calculations.
    integer, public :: ns
    integer, public :: ns_restart
    integer :: interval
 
    integer :: last_ipos
    !                                This should have been an attribute of type(mixfield_t)
 
    ! Fields
    type(mixfield_t) :: mixfield
    integer :: nauxmixfield
    type(mixfield_ptr_t) :: auxmixfield(MAX_AUXMIXFIELD)
    integer :: ipos
    !                             This should have been an attribute of type(mixfield_ptr_t)
 
    ! Preconditioning
    logical               :: kerker
    real(real64)          :: kerker_factor
    logical               :: precondition
    type(nl_operator_t)   :: preconditioner
    !                                         [GPAW](https://wiki.fysik.dtu.dk/gpaw/documentation/densitymix/densitymix.html)
 
    logical               :: mix_spin_density_matrix
  contains
    procedure :: compute_residuals_aux_field
 
  end type mix_t
 
  interface mixfield_set_vin
    module procedure dmixfield_set_vin, zmixfield_set_vin
  end interface mixfield_set_vin
 
  interface mixfield_set_vout
    module procedure dmixfield_set_vout, zmixfield_set_vout
  end interface mixfield_set_vout
 
  interface mixfield_get_vnew
    module procedure dmixfield_get_vnew, zmixfield_get_vnew
  end interface mixfield_get_vnew
 
contains
 
  subroutine mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
    type(mix_t),                   intent(out) :: smix
    type(namespace_t),             intent(in)  :: namespace
    class(space_t),                intent(in)  :: space
    type(derivatives_t), target,   intent(in)  :: der
    integer,                       intent(in)  :: d1
    integer,                       intent(in)  :: d2
    integer,             optional, intent(in)  :: def_
    type(type_t),        optional, intent(in)  :: func_type_
    character(len=*),    optional, intent(in)  :: prefix_
 
    integer :: def, ii
    character(len=32) :: prefix
    type(type_t) :: func_type
    real(real64) :: coeff
 
    push_sub(mix_init)
 
    smix%der => der
 
    def = option__mixingscheme__broyden
    if (present(def_)) def = def_
    if (present(func_type_)) then
      func_type = func_type_
    else
      func_type = type_float
    end if
    prefix = ''
    if (present(prefix_)) prefix = prefix_
 
    call messages_obsolete_variable(namespace, 'TypeOfMixing', 'MixingScheme')
 
    !%Variable MixingScheme
    !%Type integer
    !%Default broyden
    !%Section SCF::Mixing
    !%Description
    !% The scheme used to produce, at each iteration in the self-consistent cycle
    !% that attempts to solve the Kohn-Sham equations, the input density from the value
    !% of the input and output densities of previous iterations.
    !%Option linear 0
    !% Simple linear mixing.
    !%Option broyden 2
    !% Broyden scheme [C. G Broyden, <i>Math. Comp.</i> <b>19</b>, 577 (1965);
    !% D. D. Johnson, <i>Phys. Rev. B</i> <b>38</b>, 12807 (1988)].
    !% The scheme is slightly adapted, see the comments in the code.
    !% For complex functions (e.g. Sternheimer with <tt>EMEta</tt> > 0), we use the generalization
    !% with a complex dot product.
    !%Option diis 9
    !% Direct inversion in the iterative subspace (diis)
    !% scheme [P. Pulay, <i>Chem. Phys. Lett.</i>, <b>73</b>, 393
    !% (1980)] as described in [G. Kresse, and J. Furthmueller,
    !% <i>Phys. Rev. B</i> <b>54</b>, 11169 (1996)].
    !%Option broyden_adaptive 11
    !% Broyden mixing with adaptive damping. The Broyden update follows
    !% [C. G. Broyden, <i>Math. Comp.</i> <b>19</b>, 577 (1965);
    !% D. D. Johnson, <i>Phys. Rev. B</i> <b>38</b>, 12807 (1988)].
    !% The adaptive damping logic is inspired by
    !% [M. F. Herbst and A. Levitt, <i>J. Comput. Phys.</i> <b>459</b>, 111127 (2022)]
    !% and by the adaptive-coefficient analysis of
    !% [M. Novak et al., <i>Comput. Phys. Commun.</i> <b>292</b>, 108865 (2023)].
    !%End
    call parse_variable(namespace, trim(prefix)//'MixingScheme', def, smix%scheme)
    if (.not. varinfo_valid_option('MixingScheme', smix%scheme)) then
      call messages_input_error(namespace, 'MixingScheme', 'invalid option')
    end if
    call messages_print_var_option("MixingScheme", smix%scheme, namespace=namespace)
 
    !%Variable MixingPreconditioner
    !%Type logical
    !%Default false
    !%Section SCF::Mixing
    !%Description
    !% (Experimental) If set to yes, Octopus will use a preconditioner
    !% for the mixing operator.
    !% This preconditioner is disabled for systems with dimension other than 3.
    !%End
    call parse_variable(namespace, trim(prefix)+'MixingPreconditioner', .false., smix%precondition)
    if (der%dim /= 3) smix%precondition = .false.
    if (smix%precondition) call messages_experimental('MixingPreconditioner', namespace=namespace)
 
    !%Variable Mixing
    !%Type float
    !%Default 0.3
    !%Section SCF::Mixing
    !%Description
    !% The linear, Broyden, adaptive Broyden, and DIIS scheme depend on a "mixing parameter", set by this variable.
    !% Must be 0 < <tt>Mixing</tt> <= 1.
    !% For adaptive Broyden this is the maximum damping, which may be reduced automatically.
    !%End
    call parse_variable(namespace, trim(prefix)+'Mixing', 0.3_real64, coeff)
    if (coeff <= m_zero .or. coeff > m_one) then
      call messages_input_error(namespace, 'Mixing', 'Value should be positive and smaller than one')
    end if
    safe_allocate(smix%coeff(1:d2))
    smix%coeff = coeff
    smix%adaptive_damping = m_one
 
    if (d2 > 1) then
      !%Variable MixingSpinDensityMatrix
      !%Type logical
      !%Default false
      !%Section SCF::Mixing
      !%Description
      !% If set to yes, Octopus mixes the components up, down, and real and imaginary parts of
      !% the up down (for spinor calculations) of the density or potential.
      !% By default, Octopus mixes the magnetization and the charge density with a different weight,
      !% where the magnetization mixing coefficient is given by <tt>MixingMagnetization</tt>
      !%End
      call parse_variable(namespace, trim(prefix)+'MixingSpinDensityMatrix', .false., smix%mix_spin_density_matrix)
 
      !%Variable MixingMagnetization
      !%Type float
      !%Default 1.5
      !%Section SCF::Mixing
      !%Description
      !% Same as <tt>Mixing</tt>, but for the "magnetization".
      !% The default value is taken from D. D. Johnson, <i>Phys. Rev. B</i> <b>38</b>, 12807 (1988).
      !%End
      call parse_variable(namespace, trim(prefix)+'MixingMagnetization', 1.5_real64, coeff)
 
      if (.not. smix%mix_spin_density_matrix) smix%coeff(2:d2) = coeff
    else
      smix%mix_spin_density_matrix = .true.
    end if
 
    !%Variable MixingResidual
    !%Type float
    !%Default 0.05
    !%Section SCF::Mixing
    !%Description
    !% In the DIIS mixing it is benefitial to include a bit of
    !% residual into the mixing. This parameter controls this amount.
    !%End
    call parse_variable(namespace, trim(prefix)+'MixingResidual', 0.05_real64, smix%residual_coeff)
    if (smix%residual_coeff <= m_zero .or. smix%residual_coeff > m_one) then
      call messages_input_error(namespace, 'MixingResidual', 'Value should be positive and smaller than one.')
    end if
 
    !%Variable MixNumberSteps
    !%Type integer
    !%Default 4
    !%Section SCF::Mixing
    !%Description
    !% In the Broyden scheme, the new input density or potential is constructed
    !% from the values of the densities/potentials of a given number of previous iterations.
    !% This number is set by this variable. Must be greater than 1.
    !%End
    if (smix%scheme /= option__mixingscheme__linear) then
      call parse_variable(namespace, trim(prefix)//'MixNumberSteps', 4, smix%ns)
      if (smix%ns <= 1) call messages_input_error(namespace, 'MixNumberSteps')
    else
      smix%ns = 0
    end if
 
    !%Variable MixingRestart
    !%Type integer
    !%Default 20
    !%Section SCF::Mixing
    !%Description
    !% In the Broyden scheme, the mixing is restarted after
    !% the number of iterations given by this variable.
    !% Set this to zero to disable restarting the mixing.
    !%End
    if (smix%scheme /= option__mixingscheme__linear) then
      call parse_variable(namespace, trim(prefix)//'MixingRestart', 20, smix%ns_restart)
      if (smix%ns_restart < 0) call messages_input_error(namespace, 'MixingRestart')
    else
      smix%ns_restart = 0
    end if
 
    write(message(1), '(A,I4,A,I4,A)') "Info: Mixing uses ", smix%ns, " steps and restarts after ", &
      smix%ns_restart, " steps."
    call messages_info(1, namespace=namespace)
 
    !%Variable MixInterval
    !%Type integer
    !%Default 1
    !%Section SCF::Mixing
    !%Description
    !% When this variable is set to a value different than 1 (the
    !% default) a combined mixing scheme will be used, with MixInterval
    !% - 1 steps of linear mixing followed by 1 step of the selected
    !% mixing. For the moment this variable only works with DIIS mixing.
    !%End
    call parse_variable(namespace, trim(prefix)//'MixInterval', 1, smix%interval)
    if (smix%interval < 1) call messages_input_error(namespace, 'MixInterval', 'MixInterval must be larger or equal than 1')
 
    smix%iter = 0
 
    !%Variable MixingKerker
    !%Type logical
    !%Default false
    !%Section SCF::Mixing
    !%Description
    !% (Experimental) If set to yes, Octopus will use the Kerker preconditioner
    !% for the mixing operator applied in the linear and Broyden mixing schemes.
    !% When using preconditioning with the Broyden scheme, the user should increase the <tt>Mixing</tt> from its default
    !% to between 0.5 and 0.8, for optimal results.
    !% The implementation follows Yoshinori Shiihara et al., Modelling Simul. Mater. Sci. Eng. 16 (2008), 035004
    !%End
    call parse_variable(namespace, trim(prefix)+'MixingKerker', .false., smix%kerker)
    if(smix%kerker) call messages_experimental('MixingKerker')
 
    !%Variable MixingKerkerFactor
    !%Type float
    !%Default 1.0
    !%Section SCF::Mixing
    !%Description
    !% The screening factor, $q_0$ in the Kerker preconditioner in units of inverse length.
    !% For small wave vectors, q, screening behaves like $\frac{Aq}{q_0}$, where A is the mixing coefficient.
    !% An optimal value is determined empirically for a given system and mixing coeffient, however the default
    !% is a reasonable choice in many cases.
    !%End
    call parse_variable(namespace, trim(prefix)+'MixingKerkerFactor', 1._real64, smix%kerker_factor)
 
    smix%nauxmixfield = 0
    do ii = 1,max_auxmixfield
      nullify(smix%auxmixfield(ii)%p)
    end do
 
    call mixfield_init(smix, smix%mixfield, d1, d2, smix%ns, func_type)
    call mix_clear(smix)
    if (smix%precondition) call init_preconditioner()
 
    pop_sub(mix_init)
 
  contains
 
    subroutine init_preconditioner()
 
      integer :: ns, maxp, ip, is
      real(real64), parameter :: weight = 50.0_real64
 
      ! This the mixing preconditioner from GPAW:
      !
      !   https://wiki.fysik.dtu.dk/gpaw/documentation/densitymix/densitymix.html
      !
 
      assert(.not. der%mesh%use_curvilinear)
      assert(der%dim == 3)
 
      call nl_operator_init(smix%preconditioner, "Mixing preconditioner")
      call stencil_cube_get_lapl(smix%preconditioner%stencil, der%dim, 1)
      call nl_operator_build(space, der%mesh, smix%preconditioner, der%mesh%np, const_w = .not. der%mesh%use_curvilinear)
 
      ns = smix%preconditioner%stencil%size
 
      if (smix%preconditioner%const_w) then
        maxp = 1
      else
        maxp = der%mesh%np
      end if
 
      do ip = 1, maxp
 
        do is = 1, ns
          select case (sum(abs(smix%preconditioner%stencil%points(1:der%dim, is))))
          case (0)
            smix%preconditioner%w(is, ip) = m_one + weight/8.0_real64
          case (1)
            smix%preconditioner%w(is, ip) = weight/16.0_real64
          case (2)
            smix%preconditioner%w(is, ip) = weight/32.0_real64
          case (3)
            smix%preconditioner%w(is, ip) = weight/64.0_real64
          case default
            assert(.false.)
          end select
 
        end do
      end do
 
      call nl_operator_allocate_gpu_buffers(smix%preconditioner)
      call nl_operator_update_gpu_buffers(smix%preconditioner)
      call nl_operator_output_weights(smix%preconditioner)
 
    end subroutine init_preconditioner
 
  end subroutine mix_init
 
 
  ! ---------------------------------------------------------
  subroutine mix_clear(smix)
    type(mix_t),             intent(inout) :: smix
 
    push_sub(mix_clear)
 
    call mixfield_clear(smix%scheme, smix%mixfield)
 
    smix%iter = 0
    smix%last_ipos = 0
    smix%adaptive_damping = m_one
 
    pop_sub(mix_clear)
  end subroutine mix_clear
 
 
  ! ---------------------------------------------------------
  subroutine mix_end(smix)
    type(mix_t), intent(inout) :: smix
 
    integer :: ii
 
    push_sub(mix_end)
 
    if (smix%precondition) call nl_operator_end(smix%preconditioner)
 
    call mixfield_end(smix, smix%mixfield)
 
    smix%nauxmixfield = 0
    do ii = 1,max_auxmixfield
      nullify(smix%auxmixfield(ii)%p)
    end do
 
    safe_deallocate_a(smix%coeff)
 
    pop_sub(mix_end)
  end subroutine mix_end
 
 
  ! ---------------------------------------------------------
  subroutine mix_dump(namespace, restart, smix, mesh, ierr)
    type(namespace_t), intent(in)  :: namespace
    class(restart_t),  intent(in)  :: restart
    type(mix_t),       intent(in)  :: smix
    class(mesh_t),     intent(in)  :: mesh
    integer,           intent(out) :: ierr
 
    integer :: iunit, id2, id3, err, err2(4)
    character(len=40) :: lines(7)
    character(len=80) :: filename
 
    push_sub(mix_dump)
 
    ierr = 0
 
    if (restart%skip()) then
      pop_sub(mix_dump)
      return
    end if
 
    message(1) = "Debug: Writing mixing restart."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    ! functions to be written need to be compatible with the mesh
    assert(mesh%np == smix%mixfield%d1)
 
    ! First we write some information about the mixing
    iunit = restart%open('mixing')
    if(restart%do_i_write()) then
      write(lines(1), '(a11,i1)')  'scheme=    ', smix%scheme
      ! Number of global mesh points have to be written, not only smix%d1
      write(lines(2), '(a11,i10)') 'd1=        ', mesh%np_global
      write(lines(3), '(a11,i10)') 'd2=        ', smix%mixfield%d2
      write(lines(4), '(a11,i10)') 'd3=        ', smix%mixfield%d3
      write(lines(5), '(a11,i10)') 'iter=      ', smix%iter
      write(lines(6), '(a11,i10)') 'ns=        ', smix%ns
      write(lines(7), '(a11,i10)') 'last_ipos= ', smix%last_ipos
      call restart%write(iunit, lines, 7, err)
      ierr = ierr + err
    end if
    call restart%close(iunit)
 
    ! Now we write the different functions.
    ! These are not needed when using linear mixing, so we will make sure we skip this step in that case.
    err2 = 0
    if (smix%scheme /= option__mixingscheme__linear) then
      do id2 = 1, smix%mixfield%d2
        do id3 = 1, smix%mixfield%d3
 
          write(filename,'(a3,i2.2,i2.2)') 'df_', id2, id3
          if (smix%mixfield%func_type == type_float) then
            call restart%write_mesh_function(filename, mesh, smix%mixfield%ddf(1:mesh%np, id2, id3), err)
          else
            call restart%write_mesh_function(filename, mesh, smix%mixfield%zdf(1:mesh%np, id2, id3), err)
          end if
          err2(1) = err2(1) + err
 
          write(filename,'(a3,i2.2,i2.2)') 'dv_', id2, id3
          if (smix%mixfield%func_type == type_float) then
            call restart%write_mesh_function(filename, mesh, smix%mixfield%ddv(1:mesh%np, id2, id3), err)
          else
            call restart%write_mesh_function(filename, mesh, smix%mixfield%zdv(1:mesh%np, id2, id3), err)
          end if
          err2(2) = err2(2) + err
 
        end do
 
        write(filename,'(a6,i2.2)') 'f_old_', id2
        if (smix%mixfield%func_type == type_float) then
          call restart%write_mesh_function(filename, mesh, smix%mixfield%df_old(1:mesh%np, id2), err)
        else
          call restart%write_mesh_function(filename, mesh, smix%mixfield%zf_old(1:mesh%np, id2), err)
        end if
        err2(3) = err2(3) + err
 
        write(filename,'(a8,i2.2)') 'vin_old_', id2
        if (smix%mixfield%func_type == type_float) then
          call restart%write_mesh_function(filename, mesh, smix%mixfield%dvin_old(1:mesh%np, id2), err)
        else
          call restart%write_mesh_function(filename, mesh, smix%mixfield%zvin_old(1:mesh%np, id2), err)
        end if
        err2(4) = err2(4) + err
 
      end do
 
      if (err2(1) /= 0) ierr = ierr + 2
      if (err2(2) /= 0) ierr = ierr + 4
      if (err2(3) /= 0) ierr = ierr + 8
      if (err2(4) /= 0) ierr = ierr + 16
    end if
 
    message(1) = "Debug: Writing mixing restart done."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(mix_dump)
  end subroutine mix_dump
 
 
  !---------------------------------------------------------
  subroutine mix_load(namespace, restart, smix, mesh, ierr)
    type(namespace_t), intent(in)    :: namespace
    type(restart_t),   intent(in)    :: restart
    type(mix_t),       intent(inout) :: smix
    class(mesh_t),     intent(in)    :: mesh
    integer,           intent(out)   :: ierr
 
    integer :: iunit, err, err2(4)
    integer :: scheme, d1, d2, d3, ns
    integer :: id2, id3
    character(len=11)  :: str
    character(len=80)  :: filename
    character(len=256) :: lines(7)
 
    push_sub(mix_load)
 
    ierr = 0
 
    if (restart%skip()) then
      ierr = -1
      pop_sub(mix_load)
      return
    end if
 
    message(1) = "Debug: Reading mixing restart."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    ! First we read some information about the mixing
    iunit = restart%open('mixing')
    call restart%read(iunit, lines, 7, err)
    if (err /= 0) then
      ierr = ierr + 1
    else
      read(lines(1), *) str, scheme
      read(lines(2), *) str, d1
      read(lines(3), *) str, d2
      read(lines(4), *) str, d3
      read(lines(5), *) str, smix%iter
      read(lines(6), *) str, ns
      read(lines(7), *) str, smix%last_ipos
    end if
    call restart%close(iunit)
 
 
    if (ierr == 0) then
      ! We can only use the restart information if the mixing scheme and the number of steps used remained the same
      if (scheme /= smix%scheme .or. ns /= smix%ns) then
        message(1) = "The mixing scheme from the restart data is not the same as the one used in the current calculation."
        call messages_warning(1, namespace=namespace)
        ierr = ierr + 2
      end if
 
      ! Check the dimensions of the arrays to be read
      if (mesh%np_global /= d1 .or. mesh%np /= smix%mixfield%d1 .or. d2 /= smix%mixfield%d2) then
        message(1) = "The dimensions of the arrays from the mixing restart data"
        message(2) = "are not the same as the ones used in this calculation."
        call messages_warning(2, namespace=namespace)
        ierr = ierr + 4
      end if
    end if
 
 
    ! Now we read the different functions.
    ! Note that we may have more or less functions than the ones needed (d3 /= smix%d3)
    if (ierr == 0) then
      if (smix%scheme /= option__mixingscheme__linear) then
        err2 = 0
        do id2 = 1, smix%mixfield%d2
          do id3 = 1, smix%mixfield%d3
 
            write(filename,'(a3,i2.2,i2.2)') 'df_', id2, id3
            if (smix%mixfield%func_type == type_float) then
              call restart%read_mesh_function(filename, mesh, smix%mixfield%ddf(1:mesh%np, id2, id3), err)
            else
              call restart%read_mesh_function(filename, mesh, smix%mixfield%zdf(1:mesh%np, id2, id3), err)
            end if
            if (err /= 0) err2(1) = err2(1) + 1
 
            write(filename,'(a3,i2.2,i2.2)') 'dv_', id2, id3
            if (smix%mixfield%func_type == type_float) then
              call restart%read_mesh_function(filename, mesh, smix%mixfield%ddv(1:mesh%np, id2, id3), err)
            else
              call restart%read_mesh_function(filename, mesh, smix%mixfield%zdv(1:mesh%np, id2, id3), err)
            end if
            if (err /= 0) err2(2) = err2(2) + 1
 
          end do
 
          write(filename,'(a6,i2.2)') 'f_old_', id2
          if (smix%mixfield%func_type == type_float) then
            call restart%read_mesh_function(filename, mesh, smix%mixfield%df_old(1:mesh%np, id2), err)
          else
            call restart%read_mesh_function(filename, mesh, smix%mixfield%zf_old(1:mesh%np, id2), err)
          end if
          if (err /= 0) err2(3) = err2(3) + 1
 
          write(filename,'(a8,i2.2)') 'vin_old_', id2
          if (smix%mixfield%func_type == type_float) then
            call restart%read_mesh_function(filename, mesh, smix%mixfield%dvin_old(1:mesh%np, id2), err)
          else
            call restart%read_mesh_function(filename, mesh, smix%mixfield%zvin_old(1:mesh%np, id2), err)
          end if
          if (err /= 0) err2(4) = err2(4) + 1
 
        end do
 
        if (err2(1) /= 0) ierr = ierr + 8
        if (err2(2) /= 0) ierr = ierr + 16
        if (err2(3) /= 0) ierr = ierr + 32
        if (err2(4) /= 0) ierr = ierr + 64
      end if
    end if
 
    if (ierr /= 0) then
      ! Something went wront, so make sure we start from scratch
      call mix_clear(smix)
    end if
 
    if (ierr == 0 .and. smix%scheme == option__mixingscheme__broyden_adaptive) then
      smix%adaptive_damping = m_one
      smix%mixfield%last_residue = m_zero
    end if
 
 
    message(1) = "Debug: Reading mixing restart done."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(mix_load)
  end subroutine mix_load
 
  real(real64) pure function mix_coefficient(this) result(coefficient)
    type(mix_t), intent(in) :: this
 
    coefficient = this%coeff(1)
  end function mix_coefficient
 
  integer pure function mix_scheme(this) result(scheme)
    type(mix_t), intent(in) :: this
 
    scheme = this%scheme
  end function mix_scheme
 
  integer pure function mix_d3(this)
    type(mix_t), intent(in) :: this
 
    mix_d3 = this%mixfield%d3
  end function mix_d3
 
  subroutine mix_get_field(this, mixfield)
    type(mix_t), target,  intent(in) :: this
    type(mixfield_t), pointer, intent(out) :: mixfield
 
    mixfield => this%mixfield
  end subroutine mix_get_field
 
  subroutine mixing(namespace, smix)
    type(namespace_t), intent(in)    :: namespace
    type(mix_t),       intent(inout) :: smix
 
    push_sub(mixing)
 
    if (smix%mixfield%func_type == type_float) then
      call dmixing(namespace, smix, smix%mixfield%dvin, smix%mixfield%dvout, smix%mixfield%dvnew)
    else
      call zmixing(namespace, smix, smix%mixfield%zvin, smix%mixfield%zvout, smix%mixfield%zvnew)
    end if
 
    pop_sub(mixing)
  end subroutine mixing
 
  subroutine mix_add_auxmixfield(namespace, smix, mixfield)
    type(namespace_t), intent(in)         :: namespace
    type(mix_t),       intent(inout)      :: smix
    type(mixfield_t),  target, intent(in) :: mixfield
 
    push_sub(mix_add_auxmixfield)
 
    smix%nauxmixfield = smix%nauxmixfield + 1
    smix%auxmixfield(smix%nauxmixfield)%p => mixfield
 
    if (smix%scheme == option__mixingscheme__diis) then
      message(1) = 'Mixing scheme DIIS is not implemented for auxiliary mixing fields'
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(mix_add_auxmixfield)
  end subroutine mix_add_auxmixfield
 
  subroutine mixfield_init(smix, mixfield, d1, d2, d3, func_type)
    type(mix_t),      intent(in)    :: smix
    type(mixfield_t), intent(inout) :: mixfield
    integer,          intent(in)    :: d1, d2, d3
    type(type_t),     intent(in)    :: func_type
 
    push_sub(mixfield_init)
 
    mixfield%d1 = d1
    mixfield%d2 = d2
    mixfield%d3 = d3
 
    mixfield%func_type = func_type
 
    mixfield%mix_spin_density_matrix = smix%mix_spin_density_matrix
    mixfield%last_residue = m_zero
 
    if (smix%scheme /= option__mixingscheme__linear) then
      if (mixfield%func_type == type_float) then
        safe_allocate(     mixfield%ddf(1:d1, 1:d2, 1:d3))
        safe_allocate(     mixfield%ddv(1:d1, 1:d2, 1:d3))
        safe_allocate(mixfield%dvin_old(1:d1, 1:d2))
        safe_allocate(  mixfield%df_old(1:d1, 1:d2))
      else
        safe_allocate(     mixfield%zdf(1:d1, 1:d2, 1:d3))
        safe_allocate(     mixfield%zdv(1:d1, 1:d2, 1:d3))
        safe_allocate(mixfield%zvin_old(1:d1, 1:d2))
        safe_allocate(  mixfield%zf_old(1:d1, 1:d2))
      end if
    end if
 
    if (mixfield%func_type == type_float) then
      safe_allocate(mixfield%dvin(1:d1, 1:d2))
      safe_allocate(mixfield%dvout(1:d1, 1:d2))
      safe_allocate(mixfield%dvnew(1:d1, 1:d2))
      safe_allocate(mixfield%dresidual(1:d1, 1:d2))
    else
      safe_allocate(mixfield%zvin(1:d1, 1:d2))
      safe_allocate(mixfield%zvout(1:d1, 1:d2))
      safe_allocate(mixfield%zvnew(1:d1, 1:d2))
      safe_allocate(mixfield%zresidual(1:d1, 1:d2))
    end if
 
    pop_sub(mixfield_init)
  end subroutine mixfield_init
 
  subroutine mixfield_end(smix, mixfield)
    type(mix_t),      intent(inout) :: smix
    type(mixfield_t), intent(inout) :: mixfield
 
    push_sub(mixfield_end)
 
    ! Arrays got allocated for all mixing schemes, except linear mixing
    if (smix%scheme /= option__mixingscheme__linear) then
      if (mixfield%func_type == type_float) then
        safe_deallocate_a(mixfield%ddf)
        safe_deallocate_a(mixfield%ddv)
        safe_deallocate_a(mixfield%dvin_old)
        safe_deallocate_a(mixfield%df_old)
      else
        safe_deallocate_a(mixfield%zdf)
        safe_deallocate_a(mixfield%zdv)
        safe_deallocate_a(mixfield%zvin_old)
        safe_deallocate_a(mixfield%zf_old)
      end if
    end if
 
    if (mixfield%func_type == type_float) then
      safe_deallocate_a(mixfield%dvin)
      safe_deallocate_a(mixfield%dvout)
      safe_deallocate_a(mixfield%dvnew)
      safe_deallocate_a(mixfield%dresidual)
    else
      safe_deallocate_a(mixfield%zvin)
      safe_deallocate_a(mixfield%zvout)
      safe_deallocate_a(mixfield%zvnew)
      safe_deallocate_a(mixfield%zresidual)
    end if
 
    pop_sub(mixfield_end)
  end subroutine mixfield_end
 
  subroutine mixfield_clear(scheme, mixfield)
    integer,             intent(in) :: scheme
    type(mixfield_t), intent(inout) :: mixfield
    integer :: d1, d2, d3
 
    push_sub(mixfield_clear)
 
    d1 = mixfield%d1
    d2 = mixfield%d2
    d3 = mixfield%d3
 
    if (scheme /= option__mixingscheme__linear) then
      if (mixfield%func_type == type_float) then
        assert(allocated(mixfield%ddf))
        mixfield%ddf(1:d1, 1:d2, 1:d3) = m_zero
        mixfield%ddv(1:d1, 1:d2, 1:d3) = m_zero
        mixfield%dvin_old(1:d1, 1:d2) = m_zero
        mixfield%df_old(1:d1, 1:d2) = m_zero
      else
        assert(allocated(mixfield%zdf))
        mixfield%zdf(1:d1, 1:d2, 1:d3) = m_z0
        mixfield%zdv(1:d1, 1:d2, 1:d3) = m_z0
        mixfield%zvin_old(1:d1, 1:d2) = m_z0
        mixfield%zf_old(1:d1, 1:d2) = m_z0
      end if
    end if
 
    if (mixfield%func_type == type_float) then
      mixfield%dvin(1:d1, 1:d2)  = m_zero
      mixfield%dvout(1:d1, 1:d2) = m_zero
      mixfield%dvnew(1:d1, 1:d2) = m_zero
      mixfield%dresidual(1:d1, 1:d2) = m_zero
    else
      mixfield%zvin(1:d1, 1:d2)  = m_z0
      mixfield%zvout(1:d1, 1:d2) = m_z0
      mixfield%zvnew(1:d1, 1:d2) = m_z0
      mixfield%zresidual(1:d1, 1:d2) = m_z0
    end if
 
    pop_sub(mixfield_clear)
  end subroutine mixfield_clear
 
  subroutine compute_residuals_aux_field(this)
    class(mix_t), intent(inout) :: this
 
    type(mixfield_t),   pointer :: aux_field
    integer                     :: dims(2)
    integer                     :: i
 
    push_sub(compute_residuals_aux_field)
 
    do i = 1, this%nauxmixfield
      aux_field => this%auxmixfield(i)%p
      dims = [aux_field%d1, aux_field%d2]
 
      if (aux_field%func_type == type_float) then
        call dcompute_residual(this, dims, aux_field%dvin, aux_field%dvout, aux_field%dresidual)
      else
        call zcompute_residual(this, dims, aux_field%zvin, aux_field%zvout, aux_field%zresidual)
      end if
 
    end do
    nullify(aux_field)
 
    pop_sub(compute_residuals_aux_field)
 
  end subroutine compute_residuals_aux_field
 
  subroutine linear_mixing_aux_field(smix)
    type(mix_t), intent(inout) :: smix
 
    type(mixfield_t), pointer :: field
    integer :: i
 
    push_sub(linear_mixing_aux_field)
 
    call smix%compute_residuals_aux_field()
 
    do i = 1, smix%nauxmixfield
      field => smix%auxmixfield(i)%p
 
      if (field%func_type == type_float) then
        call dmixing_linear(field%d1, field%d2, smix%coeff, field%dvin, field%dresidual, field%dvnew)
      else
        call zmixing_linear(field%d1, field%d2, smix%coeff, field%zvin, field%zresidual, field%zvnew)
      end if
 
    end do
 
    nullify(field)
 
    pop_sub(linear_mixing_aux_field)
  end subroutine linear_mixing_aux_field
 
#include "undef.F90"
#include "real.F90"
 
#include "mix_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
 
#include "mix_inc.F90"
 
end module mix_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: