lda_u.F90

Go to the documentation of this file.

!! Copyright (C) 2016-2020 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 lda_u_oct_m
  use accel_oct_m
  use atomic_orbital_oct_m
  use boundaries_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use distributed_oct_m
  use electron_space_oct_m
  use energy_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_base_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use loct_oct_m
  use loct_math_oct_m
  use loewdin_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use orbitalbasis_oct_m
  use orbitalset_oct_m
  use orbitalset_utils_oct_m
  use parser_oct_m
  use poisson_oct_m
  use phase_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use nonlocal_pseudopotential_oct_m
  use projector_matrix_oct_m
  use submesh_oct_m
  use symmetries_oct_m
  use types_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                          &
    lda_u_t,                         &
    lda_u_init,                      &
    dlda_u_apply,                    &
    zlda_u_apply,                    &
    lda_u_update_basis,              &
    lda_u_update_occ_matrices,       &
    lda_u_rebind_after_copy,         &
    lda_u_accel_rebuild,             &
    lda_u_end,                       &
    lda_u_build_phase_correction,    &
    lda_u_freeze_occ,                &
    lda_u_freeze_u,                  &
    dlda_u_set_occupations,          &
    zlda_u_set_occupations,          &
    dlda_u_get_occupations,          &
    zlda_u_get_occupations,          &
    dlda_u_update_potential,         &
    zlda_u_update_potential,         &
    lda_u_get_effectiveu,            &
    lda_u_set_effectiveu,            &
    lda_u_get_effectivev,            &
    lda_u_set_effectivev,            &
    dlda_u_commute_r,                &
    zlda_u_commute_r,                &
    dlda_u_commute_r_single,         &
    zlda_u_commute_r_single,         &
    dlda_u_force,                    &
    zlda_u_force,                    &
    dlda_u_rvu,                      &
    zlda_u_rvu,                      &
    lda_u_write_info,                &
    compute_acbno_u_kanamori,        &
    dcompute_dftu_energy,            &
    zcompute_dftu_energy
 
 
  integer, public, parameter ::        &
    DFT_U_NONE                    = 0, &
    dft_u_empirical               = 1, &
    dft_u_acbn0                   = 2
 
  integer, public, parameter ::        &
    DFT_U_FLL                     = 0, &
    dft_u_amf                     = 1, &
    dft_u_mix                     = 2
 
 
  type lda_u_t
    private
    integer,            public   :: level = dft_u_none 
 
    ! DFT+U basic variables
    real(real64), allocatable, public   :: dn(:,:,:,:)
    real(real64), allocatable           :: dV(:,:,:,:)
 
    ! ACBN0 variables
    complex(real64), allocatable, public   :: zn(:,:,:,:)
    complex(real64), allocatable           :: zV(:,:,:,:)
    real(real64),    allocatable, public   :: dn_alt(:,:,:,:)
    complex(real64), allocatable, public   :: zn_alt(:,:,:,:)
 
    real(real64), allocatable           :: renorm_occ(:,:,:,:,:)
 
    ! Coulomb integrales
    real(real64), allocatable, public   :: coulomb(:,:,:,:,:)
    !                                                                 !<(for the ACBN0 functional)
    complex(real64), allocatable, public   :: zcoulomb(:,:,:,:,:,:,:)
    !                                                                 !< (for the ACBN0 functional with spinors)
 
    type(orbitalbasis_t),        public :: basis
    type(orbitalset_t), pointer, public :: orbsets(:) => null() 
    integer,                     public :: norbsets = 0
 
    integer,              public :: nspins = 0
    integer,              public :: spin_channels = 0
    integer                      :: nspecies = 0
    integer,              public :: maxnorbs = 0
    integer                      :: max_np = 0
 
    logical                      :: useAllOrbitals = .false.       
    logical                      :: skipSOrbitals = .true.         
    logical                      :: freeze_occ = .false.           
    logical                      :: freeze_u = .false.             
    logical,              public :: intersite = .false.            
    real(real64)                 :: intersite_radius = m_zero      
    logical,              public :: basisfromstates = .false.      
    real(real64)                 :: acbn0_screening = m_one        
    integer, allocatable         :: basisstates(:)
    integer, allocatable         :: basisstates_os(:)
    !                                                              !! the state specified in basisstate(:)
    logical                      :: rot_inv = .false.              
    integer                      :: double_couting = dft_u_fll     
    integer                      :: sm_poisson = sm_poisson_direct 
    real(real64), allocatable    :: dc_alpha(:)
 
    type(lattice_vectors_t), pointer :: latt
 
    type(distributed_t) :: orbs_dist
    type(mpi_grp_t)     :: grp
 
    ! Intersite interaction variables
    integer, public     :: maxneighbors = 0
    real(real64), allocatable, public  :: dn_ij(:,:,:,:,:), dn_alt_ij(:,:,:,:,:), dn_alt_ii(:,:,:,:,:)
    complex(real64), allocatable, public  :: zn_ij(:,:,:,:,:), zn_alt_ij(:,:,:,:,:), zn_alt_ii(:,:,:,:,:)
 
    ! Symmetrization-related variables
    logical                   :: symmetrize_occ_matrices
    integer, allocatable      :: inv_map_symm(:,:)
    integer                   :: nsym
    real(real64), allocatable :: symm_weight(:,:,:,:)
  end type lda_u_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
    type(lda_u_t),     target, intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    integer,                   intent(in)    :: level
    type(grid_t),              intent(in)    :: gr
    type(ions_t),      target, intent(in)    :: ions
    type(states_elec_t),       intent(in)    :: st
    type(multicomm_t),         intent(in)    :: mc
    type(kpoints_t),           intent(in)    :: kpoints
 
    integer :: is, ierr
    type(block_t) :: blk
 
    push_sub(lda_u_init)
 
    assert(.not. (level == dft_u_none))
 
    call messages_print_with_emphasis(msg="DFT+U", namespace=namespace)
    if (gr%parallel_in_domains) call messages_experimental("dft+u parallel in domains", namespace=namespace)
    this%level = level
 
    this%latt => ions%latt
 
    this%grp = mc%master_grp
 
    !%Variable DFTUBasisFromStates
    !%Type logical
    !%Default no
    !%Section Hamiltonian::DFT+U
    !%Description
    !% If set to yes, Octopus will construct the localized basis from
    !% user-defined states. The states are taken at the Gamma point (or the first k-point of the
    !% states in the restart_proj folder.
    !% The states are defined via the block DFTUBasisStates
    !%End
    call parse_variable(namespace, 'DFTUBasisFromStates', .false., this%basisfromstates)
    if (this%basisfromstates) call messages_experimental("DFTUBasisFromStates", namespace=namespace)
 
    !%Variable DFTUDoubleCounting
    !%Type integer
    !%Default dft_u_fll
    !%Section Hamiltonian::DFT+U
    !%Description
    !% This variable selects which DFT+U
    !% double counting term is used.
    !%Option dft_u_fll 0
    !% (Default) The Fully Localized Limit (FLL)
    !%Option dft_u_amf 1
    !% (Experimental) Around mean field double counting, as defined in PRB 44, 943 (1991) and PRB 49, 14211 (1994).
    !%Option dft_u_mix 2
    !% (Experimental) Mixed double countind term as introduced by Petukhov et al., PRB 67, 153106 (2003).
    !% This recovers the FLL and AMF as limiting cases.
    !%End
    call parse_variable(namespace, 'DFTUDoubleCounting', dft_u_fll, this%double_couting)
    call messages_print_var_option('DFTUDoubleCounting', this%double_couting, namespace=namespace)
    if (this%double_couting /= dft_u_fll) call messages_experimental("DFTUDoubleCounting /= dft_u_ffl", namespace=namespace)
    if (st%d%ispin == spinors .and. this%double_couting /= dft_u_fll) then
      call messages_not_implemented("AMF and MIX double counting with spinors", namespace=namespace)
    end if
 
    !%Variable DFTUPoissonSolver
    !%Type integer
    !%Section Hamiltonian::DFT+U
    !%Description
    !% This variable selects which Poisson solver
    !% is used to compute the Coulomb integrals over a submesh.
    !% These are non-periodic Poisson solvers.
    !% The FFT Poisson solver with spherical cutoff is used by default.
    !%
    !%Option dft_u_poisson_direct 0
    !% Direct Poisson solver. Slow but working in all cases.
    !%Option dft_u_poisson_isf 1
    !% (Experimental) ISF Poisson solver on a submesh.
    !% This does not work for non-orthogonal cells nor domain parallelization.
    !%Option dft_u_poisson_psolver 2
    !% (Experimental) PSolver Poisson solver on a submesh.
    !% This does not work for non-orthogonal cells nor domain parallelization.
    !% Requires the PSolver external library.
    !%Option dft_u_poisson_fft 3
    !% (Default) FFT Poisson solver on a submesh.
    !% This uses the 0D periodic version of the FFT kernels.
    !%End
    call parse_variable(namespace, 'DFTUPoissonSolver', sm_poisson_fft, this%sm_poisson)
    call messages_print_var_option('DFTUPoissonSolver', this%sm_poisson, namespace=namespace)
    if (this%sm_poisson /= sm_poisson_direct .and. this%sm_poisson /= sm_poisson_fft) then
      call messages_experimental("DFTUPoissonSolver different from dft_u_poisson_direct", namespace=namespace)
      call messages_experimental("and dft_u_poisson_fft", namespace=namespace)
    end if
    if (this%sm_poisson == sm_poisson_isf) then
      if (gr%parallel_in_domains) then
        call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_isf with domain parallelization", namespace=namespace)
      end if
      if (ions%latt%nonorthogonal) then
        call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_isf with non-orthogonal cells", namespace=namespace)
      end if
    end if
    if (this%sm_poisson == sm_poisson_psolver) then
#if !(defined HAVE_PSOLVER)
      message(1) = "The PSolver Poisson solver cannot be used since the code was not compiled with the PSolver library."
      call messages_fatal(1, namespace=namespace)
#endif
      if (gr%parallel_in_domains) then
        call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_psolver with domain parallelization", namespace=namespace)
      end if
      if (ions%latt%nonorthogonal) then
        call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_psolver with non-orthogonal cells", namespace=namespace)
      end if
    end if
 
    if (this%level == dft_u_acbn0) then
      !%Variable UseAllAtomicOrbitals
      !%Type logical
      !%Default no
      !%Section Hamiltonian::DFT+U
      !%Description
      !% If set to yes, Octopus will determine the effective U for all atomic orbitals
      !% from the peusopotential. Only available with ACBN0 functional.
      !% It is strongly recommended to set AOLoewdin=yes when using the option.
      !%End
      call parse_variable(namespace, 'UseAllAtomicOrbitals', .false., this%useAllOrbitals)
      if (this%useAllOrbitals) call messages_experimental("UseAllAtomicOrbitals", namespace=namespace)
 
      !%Variable SkipSOrbitals
      !%Type logical
      !%Default no
      !%Section Hamiltonian::DFT+U
      !%Description
      !% If set to yes, Octopus will determine the effective U for all atomic orbitals
      !% from the peusopotential but s orbitals. Only available with ACBN0 functional.
      !%End
      call parse_variable(namespace, 'SkipSOrbitals', .true., this%skipSOrbitals)
      if (.not. this%SkipSOrbitals) call messages_experimental("SkipSOrbitals", namespace=namespace)
 
      !%Variable ACBN0Screening
      !%Type float
      !%Default 1.0
      !%Section Hamiltonian::DFT+U
      !%Description
      !% If set to 0, no screening will be included in the ACBN0 functional, and the U
      !% will be estimated from bare Hartree-Fock. If set to 1 (default), the full screening
      !% of the U, as defined in the ACBN0 functional, is used.
      !%End
      call parse_variable(namespace, 'ACBN0Screening', m_one, this%acbn0_screening)
      call messages_print_var_value('ACBN0Screening', this%acbn0_screening, namespace=namespace)
 
      !%Variable ACBN0RotationallyInvariant
      !%Type logical
      !%Section Hamiltonian::DFT+U
      !%Description
      !% If set to yes, Octopus will use for U and J a formula which is rotationally invariant.
      !% This is different from the original formula for U and J.
      !% This is activated by default, except in the case of spinors, as this is not yet implemented in this case.
      !%End
      call parse_variable(namespace, 'ACBN0RotationallyInvariant', st%d%ispin /= spinors, this%rot_inv)
      call messages_print_var_value('ACBN0RotationallyInvariant', this%rot_inv, namespace=namespace)
      if (this%rot_inv .and. st%d%ispin == spinors) then
        call messages_not_implemented("Rotationally invariant ACBN0 with spinors", namespace=namespace)
      end if
 
      !%Variable ACBN0IntersiteInteraction
      !%Type logical
      !%Default no
      !%Section Hamiltonian::DFT+U
      !%Description
      !% If set to yes, Octopus will determine the effective intersite interaction V
      !% Only available with ACBN0 functional.
      !% It is strongly recommended to set AOLoewdin=yes when using the option.
      !%End
      call parse_variable(namespace, 'ACBN0IntersiteInteraction', .false., this%intersite)
      call messages_print_var_value('ACBN0IntersiteInteraction', this%intersite, namespace=namespace)
      if (this%intersite) call messages_experimental("ACBN0IntersiteInteraction", namespace=namespace)
 
      if (this%intersite) then
 
        !This is a non local operator. To make this working, one probably needs to apply the
        ! symmetries to the generalized occupation matrices
        if (kpoints%use_symmetries) then
          call messages_not_implemented("Intersite interaction with kpoint symmetries", namespace=namespace)
        end if
 
        !%Variable ACBN0IntersiteCutoff
        !%Type float
        !%Section Hamiltonian::DFT+U
        !%Description
        !% The cutoff radius defining the maximal intersite distance considered.
        !% Only available with ACBN0 functional with intersite interaction.
        !%End
        call parse_variable(namespace, 'ACBN0IntersiteCutoff', m_zero, this%intersite_radius, unit = units_inp%length)
        if (abs(this%intersite_radius) < m_epsilon) then
          call messages_write("ACBN0IntersiteCutoff must be greater than 0")
          call messages_fatal(1, namespace=namespace)
        end if
 
      end if
 
    end if
 
    ! DFT+U with ultrasoft/PAW (generalized overlap) uses the S-augmented occupation
    ! <phi|S|psi~> (and the S-dressed +U operator). The following sub-paths additionally
    ! need an augmentation that is not yet implemented and are fenced here.
    if (st%requires_generalized_overlap) then
      if (this%intersite) then
        call messages_not_implemented("Intersite DFT+U with ultrasoft/PAW", namespace=namespace)
      end if
      if (this%basisfromstates) then
        call messages_not_implemented("DFTUBasisFromStates with ultrasoft/PAW", namespace=namespace)
      end if
      ! The empirical +U occupation and operator are fully S-augmented (and the orbitals are
      ! S-normalized). ACBN0 additionally needs the augmentation of the Coulomb integrals (the
      ! orbital pair density gains a Q contribution), which is not implemented, so ACBN0 is fenced.
      if (this%level == dft_u_acbn0) then
        call messages_not_implemented("ACBN0 (computed U/J) with ultrasoft/PAW; use dft_u_empirical", &
          namespace=namespace)
      end if
    end if
 
    call lda_u_write_info(this, namespace=namespace)
 
    if (.not. this%basisfromstates) then
 
      call orbitalbasis_init(this%basis, namespace, space%periodic_dim)
 
      ! Loewdin orthonormalization with ultrasoft/PAW would need the S metric
      if (st%requires_generalized_overlap .and. this%basis%orthogonalization) then
        call messages_not_implemented("AOLoewdin orthonormalization with ultrasoft/PAW", &
          namespace=namespace)
      end if
 
      if (states_are_real(st)) then
        call dorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
          this%skipSOrbitals, this%useAllOrbitals)
      else
        call zorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
          this%skipSOrbitals, this%useAllOrbitals)
      end if
      this%orbsets => this%basis%orbsets
      this%norbsets = this%basis%norbsets
      this%maxnorbs = this%basis%maxnorbs
      this%max_np = this%basis%max_np
      this%nspins = st%d%nspin
      this%spin_channels = st%d%spin_channels
      this%nspecies = ions%nspecies
 
      !We allocate the necessary ressources
      if (states_are_real(st)) then
        call dlda_u_allocate(this, st)
      else
        call zlda_u_allocate(this, st)
      end if
 
      if (gr%parallel_in_domains) then
        call distributed_init_serial(this%orbs_dist, this%norbsets)
      else
        call distributed_init(this%orbs_dist, this%norbsets, mpi_world%comm, "orbsets")
      end if
 
    else
 
      !%Variable DFTUBasisStates
      !%Type block
      !%Default none
      !%Section Hamiltonian::DFT+U
      !%Description
      !% This block starts by a line containing a single integer describing the number of
      !% orbital sets. One orbital set is a group of orbitals on which one adds a Hubbard U.
      !% Each following line of this block contains the index of a state to be used to construct the
      !% localized basis, followed by the index of the corresponding orbital set.
      !% See DFTUBasisFromStates for details.
      !%End
      if (parse_block(namespace, 'DFTUBasisStates', blk) == 0) then
        call parse_block_integer(blk, 0, 0, this%norbsets)
        this%maxnorbs = parse_block_n(blk)-1
        if (this%maxnorbs <1) then
          write(message(1),'(a,i3,a,i3)') 'DFTUBasisStates must contains at least one state.'
          call messages_fatal(1, namespace=namespace)
        end if
        safe_allocate(this%basisstates(1:this%maxnorbs))
        safe_allocate(this%basisstates_os(1:this%maxnorbs))
        do is = 1, this%maxnorbs
          call parse_block_integer(blk, is, 0, this%basisstates(is))
          call parse_block_integer(blk, is, 1, this%basisstates_os(is))
        end do
        call parse_block_end(blk)
      else
        write(message(1),'(a,i3,a,i3)') 'DFTUBasisStates must be specified if DFTUBasisFromStates=yes'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (states_are_real(st)) then
        call dorbitalbasis_build_empty(this%basis, gr, st%d%dim, this%norbsets, this%basisstates_os)
      else
        call zorbitalbasis_build_empty(this%basis, gr, st%d%dim, this%norbsets, this%basisstates_os)
      end if
 
      this%max_np = gr%np
      this%nspins = st%d%nspin
      this%spin_channels = st%d%spin_channels
      this%nspecies = 1
 
      this%orbsets => this%basis%orbsets
 
      call distributed_init_serial(this%orbs_dist, this%norbsets)
 
      !We allocate the necessary ressources
      if (states_are_real(st)) then
        call dlda_u_allocate(this, st)
      else
        call zlda_u_allocate(this, st)
      end if
 
      call lda_u_loadbasis(this, namespace, space, st, gr, mc, ierr)
      if (ierr /= 0) then
        message(1) = "Unable to load DFT+U basis from selected states."
        call messages_fatal(1)
      end if
 
    end if
 
    ! Symmetrization of the occupation matrices
    this%symmetrize_occ_matrices = st%symmetrize_density .or. kpoints%use_symmetries
    if (this%basisfromstates) this%symmetrize_occ_matrices = .false.
    if (this%symmetrize_occ_matrices) then
      call build_symmetrization_map(this, ions, gr)
    end if
 
    safe_allocate(this%dc_alpha(this%norbsets))
    this%dc_alpha = m_one
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(lda_u_init)
  end subroutine lda_u_init
 
 
  ! ---------------------------------------------------------
  subroutine lda_u_init_coulomb_integrals(this, namespace, space, gr, st, psolver)
    type(lda_u_t),     target, intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    type(grid_t),              intent(in)    :: gr
    type(states_elec_t),       intent(in)    :: st
    type(poisson_t),           intent(in)    :: psolver
 
    logical :: complex_coulomb_integrals
    integer :: ios
    integer :: norbs
 
    push_sub(lda_u_init_coulomb_integrals)
 
    norbs = this%maxnorbs
 
    if (.not. this%basisfromstates) then
 
      if (this%level == dft_u_acbn0) then
 
        complex_coulomb_integrals = .false.
        do ios = 1, this%norbsets
          if (this%orbsets(ios)%ndim  > 1) complex_coulomb_integrals = .true.
        end do
 
        if (.not. complex_coulomb_integrals) then
          write(message(1),'(a)')    'Computing the Coulomb integrals of the localized basis.'
          call messages_info(1, namespace=namespace)
          safe_allocate(this%coulomb(1:norbs,1:norbs,1:norbs,1:norbs, 1:this%norbsets))
          if (states_are_real(st)) then
            call dcompute_coulomb_integrals(this, namespace, space, gr, psolver)
          else
            call zcompute_coulomb_integrals(this, namespace, space, gr, psolver)
          end if
        else
          assert(.not. states_are_real(st))
          write(message(1),'(a)')    'Computing complex Coulomb integrals of the localized basis.'
          call messages_info(1, namespace=namespace)
          safe_allocate(this%zcoulomb(1:norbs, 1:norbs, 1:norbs, 1:norbs, 1:st%d%dim, 1:st%d%dim, 1:this%norbsets))
          call compute_complex_coulomb_integrals(this, gr, st, psolver, namespace, space)
        end if
      end if
 
    else
      write(message(1),'(a)')    'Computing the Coulomb integrals of the localized basis.'
      call messages_info(1, namespace=namespace)
      safe_allocate(this%coulomb(1:norbs, 1:norbs, 1:norbs, 1:norbs, 1:this%norbsets))
      if (states_are_real(st)) then
        call dcompute_coulomb_integrals(this, namespace, space, gr, psolver)
      else
        call zcompute_coulomb_integrals(this, namespace, space, gr, psolver)
      end if
 
    end if
 
    pop_sub(lda_u_init_coulomb_integrals)
 
  end subroutine lda_u_init_coulomb_integrals
 
 
  ! ---------------------------------------------------------
  subroutine lda_u_end(this)
    implicit none
    type(lda_u_t), intent(inout) :: this
 
    push_sub(lda_u_end)
 
    this%level = dft_u_none
 
    safe_deallocate_a(this%dn)
    safe_deallocate_a(this%zn)
    safe_deallocate_a(this%dn_alt)
    safe_deallocate_a(this%zn_alt)
    safe_deallocate_a(this%dV)
    safe_deallocate_a(this%zV)
    safe_deallocate_a(this%coulomb)
    safe_deallocate_a(this%zcoulomb)
    safe_deallocate_a(this%renorm_occ)
    safe_deallocate_a(this%dn_ij)
    safe_deallocate_a(this%zn_ij)
    safe_deallocate_a(this%dn_alt_ij)
    safe_deallocate_a(this%zn_alt_ij)
    safe_deallocate_a(this%dn_alt_ii)
    safe_deallocate_a(this%zn_alt_ii)
    safe_deallocate_a(this%basisstates)
    safe_deallocate_a(this%basisstates_os)
    safe_deallocate_a(this%dc_alpha)
    safe_deallocate_a(this%inv_map_symm)
    safe_deallocate_a(this%symm_weight)
 
    nullify(this%orbsets)
    call orbitalbasis_end(this%basis)
 
    this%max_np = 0
 
    if (.not. this%basisfromstates) then
      call distributed_end(this%orbs_dist)
    end if
 
    pop_sub(lda_u_end)
  end subroutine lda_u_end
 
  ! ---------------------------------------------------------
  subroutine lda_u_rebind_after_copy(this, ions)
    type(lda_u_t), target, intent(inout) :: this
    type(ions_t), target, intent(in)  :: ions
 
    integer :: ios, ispec
 
    push_sub(lda_u_rebind_after_copy)
 
    if (this%level == dft_u_none) then
      pop_sub(lda_u_rebind_after_copy)
      return
    end if
 
    this%latt => ions%latt
 
    if (allocated(this%basis%orbsets)) then
      this%orbsets => this%basis%orbsets
    else
      nullify(this%orbsets)
      pop_sub(lda_u_rebind_after_copy)
      return
    end if
 
    do ios = 1, this%norbsets
      if (.not. associated(this%orbsets(ios)%spec)) cycle
 
      ispec = this%orbsets(ios)%spec_index
      if (ispec < 1 .or. ispec > ions%nspecies) then
        call messages_not_implemented("lda_u_rebind_after_copy with invalid species index")
      end if
      if (.not. associated(ions%species(ispec)%s)) then
        call messages_not_implemented("lda_u_rebind_after_copy with unassociated ion species")
      end if
      this%orbsets(ios)%spec => ions%species(ispec)%s
    end do
 
    pop_sub(lda_u_rebind_after_copy)
  end subroutine lda_u_rebind_after_copy
 
  ! ---------------------------------------------------------
  subroutine lda_u_accel_rebuild(this, kpt)
    type(lda_u_t),       intent(inout) :: this
    type(distributed_t), intent(in)    :: kpt
 
    integer :: ios, ik
 
    push_sub(lda_u_accel_rebuild)
 
    if (.not. accel_is_enabled()) then
      pop_sub(lda_u_accel_rebuild)
      return
    end if
 
    if (this%level == dft_u_none .or. .not. associated(this%orbsets)) then
      pop_sub(lda_u_accel_rebuild)
      return
    end if
 
    do ios = 1, this%norbsets
      call accel_detach_buffer(this%orbsets(ios)%dbuff_orb)
      call accel_detach_buffer(this%orbsets(ios)%zbuff_orb)
      call accel_detach_buffer(this%orbsets(ios)%sphere%buff_map)
 
      if (allocated(this%orbsets(ios)%buff_eorb)) then
        do ik = lbound(this%orbsets(ios)%buff_eorb, dim=1), ubound(this%orbsets(ios)%buff_eorb, dim=1)
          call accel_detach_buffer(this%orbsets(ios)%buff_eorb(ik))
        end do
        safe_deallocate_a(this%orbsets(ios)%buff_eorb)
      end if
 
      if (allocated(this%orbsets(ios)%dorb)) then
        call dorbitalset_transfer_to_device(this%orbsets(ios), kpt, this%orbsets(ios)%allocated_on_mesh)
      else if (allocated(this%orbsets(ios)%zorb)) then
        call zorbitalset_transfer_to_device(this%orbsets(ios), kpt, this%orbsets(ios)%allocated_on_mesh)
      else
        call messages_not_implemented("lda_u_accel_rebuild with missing orbital arrays")
      end if
    end do
 
    pop_sub(lda_u_accel_rebuild)
  end subroutine lda_u_accel_rebuild
 
  ! When moving the ions, the basis must be reconstructed
  subroutine lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase, uspp_overlap)
    type(lda_u_t),     target, intent(inout) :: this
    class(space_t),            intent(in)    :: space
    type(grid_t),              intent(in)    :: gr
    type(ions_t),      target, intent(in)    :: ions
    type(states_elec_t),       intent(in)    :: st
    type(poisson_t),           intent(in)    :: psolver
    type(namespace_t),         intent(in)    :: namespace
    type(kpoints_t),           intent(in)    :: kpoints
    logical,                   intent(in)    :: has_phase
    type(nonlocal_pseudopotential_t), intent(in) :: uspp_overlap
 
    integer :: ios, maxorbs, nspin
 
    if(this%level == dft_u_none) return
 
    push_sub(lda_u_update_basis)
 
    if(.not. this%basisfromstates) then
      !We clean the orbital basis, to be able to reconstruct it
      call orbitalbasis_end(this%basis)
      nullify(this%orbsets)
 
      !We now reconstruct the basis
      if (states_are_real(st)) then
        call dorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
          this%skipSOrbitals, this%useAllOrbitals, verbose = .false.)
      else
        call zorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
          this%skipSOrbitals, this%useAllOrbitals, verbose = .false.)
      end if
      this%orbsets => this%basis%orbsets
      this%max_np = this%basis%max_np
 
      ! Ultrasoft/PAW: AONormalize=yes means <phi|S|phi> = 1, so the
      ! occupation <phi|S|psi~> is not over-counted. Done before the Bloch phases are set.
      if (st%requires_generalized_overlap .and. this%basis%normalize) then
        if (states_are_real(st)) then
          call dlda_u_normalize_orbitals_uspp(this, gr%der%mesh, uspp_overlap)
        else
          call zlda_u_normalize_orbitals_uspp(this, gr%der%mesh, uspp_overlap)
        end if
      end if
    end if
 
    !In case of intersite interaction we need to reconstruct the basis
    if (this%intersite) then
      this%maxneighbors = 0
      do ios = 1, this%norbsets
        call orbitalset_init_intersite(this%orbsets(ios), namespace, space, this%grp, ios, ions, gr%der, psolver, &
          this%orbsets, this%norbsets, this%maxnorbs, this%intersite_radius, st%d%kpt, has_phase, &
          this%sm_poisson, this%basisfromstates, this%basis%combine_j_orbitals)
        this%maxneighbors = max(this%maxneighbors, this%orbsets(ios)%nneighbors)
      end do
 
      maxorbs = this%maxnorbs
      nspin = this%nspins
 
      if (states_are_real(st)) then
        safe_deallocate_a(this%dn_ij)
        safe_allocate(this%dn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%dn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
        safe_deallocate_a(this%dn_alt_ij)
        safe_allocate(this%dn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%dn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
        safe_deallocate_a(this%dn_alt_ii)
        safe_allocate(this%dn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%dn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
      else
        safe_deallocate_a(this%zn_ij)
        safe_allocate(this%zn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%zn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
        safe_deallocate_a(this%zn_alt_ij)
        safe_allocate(this%zn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%zn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
        safe_deallocate_a(this%zn_alt_ii)
        safe_allocate(this%zn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
        this%zn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
      end if
    end if
 
    ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
    ! In case of a laser field, the phase is recomputed in hamiltonian_elec_update
    if (has_phase) then
      call lda_u_build_phase_correction(this, space, st%d, gr%der%boundaries, namespace, kpoints)
    else
      if(.not. this%basisfromstates) then
        !In case there is no phase, we perform the orthogonalization here
        if(this%basis%orthogonalization) then
          call dloewdin_orthogonalize(this%basis, st%d%kpt, namespace)
        else
          if(debug%info .and. space%is_periodic()) then
            call dloewdin_info(this%basis, st%d%kpt, namespace)
          end if
        end if
      end if
    end if
 
    ! Rebuild the Coulomb integrals
    if (allocated(this%coulomb)) then
      safe_deallocate_a(this%coulomb)
    end if
    if (allocated(this%zcoulomb)) then
      safe_deallocate_a(this%zcoulomb)
    end if
    call lda_u_init_coulomb_integrals(this, namespace, space, gr, st, psolver)
 
    pop_sub(lda_u_update_basis)
 
  end subroutine lda_u_update_basis
 
  ! Interface for the X(update_occ_matrices) routines
  subroutine lda_u_update_occ_matrices(this, namespace, mesh, st, phase, energy, uspp_overlap)
    type(lda_u_t),                 intent(inout) :: this
    type(namespace_t),             intent(in)    :: namespace
    class(mesh_t),                 intent(in)    :: mesh
    type(states_elec_t),           intent(inout) :: st
    type(phase_t),                 intent(in)    :: phase
    type(energy_t),                intent(inout) :: energy
    type(nonlocal_pseudopotential_t), intent(in) :: uspp_overlap
 
    if (this%level == dft_u_none .or. this%freeze_occ) return
    push_sub(lda_u_update_occ_matrices)
 
    if (states_are_real(st)) then
      call dupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u, uspp_overlap)
    else
      if (phase%is_allocated()) then
        call zupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u, uspp_overlap, phase)
      else
        call zupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u, uspp_overlap)
      end if
    end if
 
    pop_sub(lda_u_update_occ_matrices)
  end subroutine lda_u_update_occ_matrices
 
 
  subroutine lda_u_build_phase_correction(this, space, std, boundaries, namespace, kpoints, vec_pot, vec_pot_var)
    type(lda_u_t),                 intent(inout) :: this
    class(space_t),                intent(in)    :: space
    type(states_elec_dim_t),       intent(in)    :: std
    type(boundaries_t),            intent(in)    :: boundaries
    type(namespace_t),             intent(in)    :: namespace
    type(kpoints_t),               intent(in)    :: kpoints
    real(real64), optional,  allocatable, intent(in)    :: vec_pot(:)
    real(real64), optional,  allocatable, intent(in)    :: vec_pot_var(:, :)
 
    integer :: ios
 
    if (boundaries%spiralBC) call messages_not_implemented("DFT+U with spiral boundary conditions", &
      namespace=namespace)
 
    push_sub(lda_u_build_phase_correction)
 
    write(message(1), '(a)') 'Debug: Building the phase correction for DFT+U orbitals.'
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    do ios = 1, this%norbsets
      call orbitalset_update_phase(this%orbsets(ios), space%dim, std%kpt, kpoints, &
        (std%ispin==spin_polarized), vec_pot, vec_pot_var)
      call orbitalset_update_phase_shift(this%orbsets(ios), space%dim, std%kpt, kpoints, &
        (std%ispin==spin_polarized), vec_pot, vec_pot_var)
    end do
 
    if (.not. this%basisfromstates) then
      if (this%basis%orthogonalization) then
        call zloewdin_orthogonalize(this%basis, std%kpt, namespace)
      else
        if (debug%info .and. space%is_periodic()) call zloewdin_info(this%basis, std%kpt, namespace)
      end if
    end if
 
    pop_sub(lda_u_build_phase_correction)
 
  end subroutine lda_u_build_phase_correction
 
  ! ---------------------------------------------------------
  subroutine compute_acbno_u_kanamori(this, st, kanamori)
    type(lda_u_t), intent(in) :: this
    type(states_elec_t), intent(in) :: st
    real(real64), contiguous, intent(out) :: kanamori(:,:)
 
    if (this%nspins == 1) then
      if (states_are_real(st)) then
        call dcompute_acbno_u_kanamori_restricted(this, kanamori)
      else
        call zcompute_acbno_u_kanamori_restricted(this, kanamori)
      end if
    else
      if (states_are_real(st)) then
        call dcompute_acbno_u_kanamori(this, kanamori)
      else
        call zcompute_acbno_u_kanamori(this, kanamori)
      end if
    end if
 
 
  end subroutine compute_acbno_u_kanamori
 
  ! ---------------------------------------------------------
  subroutine lda_u_freeze_occ(this)
    type(lda_u_t),     intent(inout) :: this
 
    this%freeze_occ = .true.
  end subroutine lda_u_freeze_occ
 
  ! ---------------------------------------------------------
  subroutine lda_u_freeze_u(this)
    type(lda_u_t),     intent(inout) :: this
 
    this%freeze_u = .true.
  end subroutine lda_u_freeze_u
 
  ! ---------------------------------------------------------
  subroutine lda_u_set_effectiveu(this, Ueff)
    type(lda_u_t), intent(inout) :: this
    real(real64), contiguous, intent(in) :: ueff(:)
 
    integer :: ios
 
    push_sub(lda_u_set_effectiveu)
 
    do ios = 1,this%norbsets
      this%orbsets(ios)%Ueff = ueff(ios)
    end do
 
    pop_sub(lda_u_set_effectiveu)
  end subroutine lda_u_set_effectiveu
 
  ! ---------------------------------------------------------
  subroutine lda_u_get_effectiveu(this, Ueff)
    type(lda_u_t), intent(in) :: this
    real(real64), contiguous, intent(inout) :: ueff(:)
 
    integer :: ios
 
    push_sub(lda_u_get_effectiveu)
 
    do ios = 1,this%norbsets
      ueff(ios) = this%orbsets(ios)%Ueff
    end do
 
    pop_sub(lda_u_get_effectiveu)
  end subroutine lda_u_get_effectiveu
 
  ! ---------------------------------------------------------
  subroutine lda_u_set_effectivev(this, Veff)
    type(lda_u_t), intent(inout) :: this
    real(real64), contiguous, intent(in) :: veff(:)
 
    integer :: ios, ncount
 
    push_sub(lda_u_set_effectivev)
 
    ncount = 0
    do ios = 1, this%norbsets
      this%orbsets(ios)%V_ij(1:this%orbsets(ios)%nneighbors,0) = veff(ncount+1:ncount+this%orbsets(ios)%nneighbors)
      ncount = ncount + this%orbsets(ios)%nneighbors
    end do
 
    pop_sub(lda_u_set_effectivev)
  end subroutine lda_u_set_effectivev
 
  ! ---------------------------------------------------------
  subroutine lda_u_get_effectivev(this, Veff)
    type(lda_u_t), intent(in) :: this
    real(real64), contiguous, intent(inout) :: veff(:)
 
    integer :: ios, ncount
 
    push_sub(lda_u_get_effectivev)
 
    ncount = 0
    do ios = 1, this%norbsets
      veff(ncount+1:ncount+this%orbsets(ios)%nneighbors) = this%orbsets(ios)%V_ij(1:this%orbsets(ios)%nneighbors,0)
      ncount = ncount + this%orbsets(ios)%nneighbors
    end do
 
    pop_sub(lda_u_get_effectivev)
  end subroutine lda_u_get_effectivev
 
  ! ---------------------------------------------------------
  subroutine lda_u_write_info(this, iunit, namespace)
    type(lda_u_t),               intent(in) :: this
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    push_sub(lda_u_write_info)
 
    write(message(1), '(1x)')
    call messages_info(1, iunit=iunit, namespace=namespace)
    if (this%level == dft_u_empirical) then
      write(message(1), '(a)') "Method:"
      write(message(2), '(a)') "  [1] Dudarev et al., Phys. Rev. B 57, 1505 (1998)"
      call messages_info(2, iunit=iunit, namespace=namespace)
    else
      if (.not. this%intersite) then
        write(message(1), '(a)') "Method:"
        write(message(2), '(a)') "  [1] Agapito et al., Phys. Rev. X 5, 011006 (2015)"
      else
        write(message(1), '(a)') "Method:"
        write(message(2), '(a)') "  [1] Tancogne-Dejean, and Rubio, Phys. Rev. B 102, 155117 (2020)"
      end if
      call messages_info(2, iunit=iunit, namespace=namespace)
    end if
    write(message(1), '(a)') "Implementation:"
    write(message(2), '(a)') "  [1] Tancogne-Dejean, Oliveira, and Rubio, Phys. Rev. B 69, 245133 (2017)"
    write(message(3), '(1x)')
    call messages_info(3, iunit=iunit, namespace=namespace)
 
    pop_sub(lda_u_write_info)
 
  end subroutine lda_u_write_info
 
  ! ---------------------------------------------------------
  subroutine lda_u_loadbasis(this, namespace, space, st, mesh, mc, ierr)
    type(lda_u_t),        intent(inout) :: this
    type(namespace_t),    intent(in)    :: namespace
    class(space_t),       intent(in)    :: space
    type(states_elec_t),  intent(in)    :: st
    class(mesh_t),        intent(in)    :: mesh
    type(multicomm_t),    intent(in)    :: mc
    integer,              intent(out)   :: ierr
 
    integer :: err, wfns_file, is, ist, idim, ik, ios, iorb
    type(restart_t) :: restart_gs
    character(len=256)   :: lines(3)
    character(len=256), allocatable :: restart_file(:, :)
    logical,            allocatable :: restart_file_present(:, :)
    character(len=12)    :: filename
    character(len=1)     :: char
    character(len=50)    :: str
    type(orbitalset_t), pointer :: os
    integer, allocatable :: count(:)
    real(real64)         :: norm, center(space%dim)
    real(real64), allocatable   :: dpsi(:,:,:)
    complex(real64), allocatable   :: zpsi(:,:,:)
 
 
    push_sub(lda_u_loadbasis)
 
    ierr = 0
 
    message(1) = "Debug: Loading DFT+U basis from states."
    call messages_info(1, debug_only=.true.)
 
    call restart_gs%init(namespace, restart_proj, restart_type_load, mc, err, mesh=mesh)
 
    ! If any error occurred up to this point then it is not worth continuing,
    ! as there something fundamentally wrong with the restart files
    if (err /= 0) then
      call restart_gs%end()
      message(1) = "Error loading DFT+U basis from states, cannot proceed with the calculation"
      call messages_fatal(1)
      pop_sub(lda_u_loadbasis)
      return
    end if
    ! only supported for obf format
    if (restart_gs%file_format_states /= option__restartfileformatstates__obf) then
      message(1) = "Error: loading DFT+U basis only supported for OBF restart format"
      message(2) = "Please set RestartFileFormatStates = obf and re-run the GS calculation."
      call messages_fatal(2)
    end if
 
    ! open files to read
    wfns_file  = restart_gs%open('wfns')
    call restart_gs%read(wfns_file, lines, 2, err)
    if (err /= 0) then
      ierr = ierr - 2**5
    else if (states_are_real(st)) then
      read(lines(2), '(a)') str
      if (str(2:8) == 'Complex') then
        message(1) = "Cannot read real states from complex wavefunctions."
        call messages_fatal(1, namespace=namespace)
      else if (str(2:5) /= 'Real') then
        message(1) = "Restart file 'wfns' does not specify real/complex; cannot check compatibility."
        call messages_warning(1, namespace=namespace)
      end if
    end if
    ! complex can be restarted from real, so there is no problem.
 
    ! If any error occurred up to this point then it is not worth continuing,
    ! as there something fundamentally wrong with the restart files
    if (err /= 0) then
      call restart_gs%close(wfns_file)
      call restart_gs%end()
      message(1) = "Error loading DFT+U basis from states, cannot proceed with the calculation"
      call messages_fatal(1)
      pop_sub(lda_u_loadbasis)
      return
    end if
 
    safe_allocate(restart_file(1:st%d%dim, 1:st%nst))
    safe_allocate(restart_file_present(1:st%d%dim, 1:st%nst))
    restart_file_present = .false.
 
    ! Next we read the list of states from the files.
    ! Errors in reading the information of a specific state from the files are ignored
    ! at this point, because later we will skip reading the wavefunction of that state.
    do
      call restart_gs%read(wfns_file, lines, 1, err)
      if (err == 0) then
        read(lines(1), '(a)') char
        if (char == '%') then
          !We reached the end of the file
          exit
        else
          read(lines(1), *) ik, char, ist, char, idim, char, filename
        end if
      end if
 
      if (any(this%basisstates==ist) .and. ik == 1) then
        restart_file(idim, ist) = trim(filename)
        restart_file_present(idim, ist) = .true.
      end if
    end do
    call restart_gs%close(wfns_file)
 
    !We loop over the states we need
    safe_allocate(count(1:this%norbsets))
    count = 0
    do is = 1, this%maxnorbs
      ist = this%basisstates(is)
      ios = this%basisstates_os(is)
      count(ios) = count(ios)+1
      do idim = 1, st%d%dim
 
        if (.not. restart_file_present(idim, ist)) then
          write(message(1), '(a,i3,a)') "Cannot read states ", ist, "from the projection folder"
          call messages_fatal(1, namespace=namespace)
        end if
 
        if (states_are_real(st)) then
          call restart_gs%read_mesh_function(restart_file(idim, ist), mesh, &
            this%orbsets(ios)%dorb(:,idim,count(ios)), err)
        else
          call restart_gs%read_mesh_function(restart_file(idim, ist), mesh, &
            this%orbsets(ios)%zorb(:,idim,count(ios)), err)
        end if
        if (err /= 0) then
          message(1) = "Error loading mesh function "//trim(restart_file(idim, ist))
          call messages_fatal(1)
        end if
      end do
    end do
    safe_deallocate_a(count)
    safe_deallocate_a(restart_file)
    safe_deallocate_a(restart_file_present)
    call restart_gs%end()
 
    ! Normalize the orbitals. This is important if we use Wannier orbitals instead of KS states
    if(this%basis%normalize) then
      do ios = 1, this%norbsets
        do iorb = 1, this%orbsets(ios)%norbs
          if (states_are_real(st)) then
            norm = dmf_nrm2(mesh, st%d%dim, this%orbsets(ios)%dorb(:,:,iorb))
            call lalg_scal(mesh%np, st%d%dim, m_one/norm, this%orbsets(ios)%dorb(:,:,iorb))
          else
            norm = zmf_nrm2(mesh, st%d%dim, this%orbsets(ios)%zorb(:,:,iorb))
            call lalg_scal(mesh%np, st%d%dim, m_one/norm, this%orbsets(ios)%zorb(:,:,iorb))
          end if
        end do
      end do
    end if
 
    ! We rotate the orbitals in the complex plane to have them as close as possible to real functions
    if(states_are_complex(st) .and. st%d%dim == 1) then
      do ios = 1, this%norbsets
        do iorb = 1, this%orbsets(ios)%norbs
          call zmf_fix_phase(mesh, this%orbsets(ios)%zorb(:,1,iorb))
        end do
      end do
    end if
 
    ! We determine the center of charge by computing <w|r|w>
    ! We could also determine the spread by \Omega = <w|r^2|w> - <w|r|w>^2
    do ios = 1, this%norbsets
      if (states_are_real(st)) then
        call dorbitalset_get_center_of_mass(this%orbsets(ios), space, mesh, this%latt)
      else
        call zorbitalset_get_center_of_mass(this%orbsets(ios), space, mesh, this%latt)
      end if
    end do
 
    message(1) = "Debug: Converting the Wannier states to submeshes."
    call messages_info(1, debug_only=.true.)
 
    ! We now transfer the states to a submesh centered on the center of mass of the Wannier orbitals
    this%max_np = 0
    do ios = 1, this%norbsets
      os => this%orbsets(ios)
      center = os%sphere%center
      safe_deallocate_a(os%sphere%center)
      if (states_are_real(st)) then
        safe_allocate(dpsi(1:mesh%np, 1:os%ndim, 1:os%norbs))
        dpsi(1:mesh%np, 1:os%ndim, 1:os%norbs) = os%dorb(1:mesh%np, 1:os%ndim, 1:os%norbs)
 
        safe_deallocate_a(os%dorb)
        !We initialise the submesh corresponding to the orbital
        call submesh_init(os%sphere, space, mesh, this%latt, center, os%radius)
        safe_allocate(os%dorb(1:os%sphere%np, 1:os%ndim, 1:os%norbs))
        do iorb = 1, os%norbs
          do idim = 1, os%ndim
            call dsubmesh_copy_from_mesh(os%sphere, dpsi(:,idim,iorb), os%dorb(:,idim, iorb))
          end do
        end do
        safe_deallocate_a(dpsi)
      else
        safe_allocate(zpsi(1:mesh%np, 1:os%ndim, 1:os%norbs))
        zpsi(1:mesh%np, 1:os%ndim, 1:os%norbs) = os%zorb(1:mesh%np, 1:os%ndim, 1:os%norbs)
        safe_deallocate_a(os%zorb)
        !We initialise the submesh corresponding to the orbital
        call submesh_init(os%sphere, space, mesh, this%latt, center, os%radius)
        safe_allocate(os%zorb(1:os%sphere%np, 1:os%ndim, 1:os%norbs))
        do iorb = 1, os%norbs
          do idim = 1, os%ndim
            call zsubmesh_copy_from_mesh(os%sphere, zpsi(:,idim,iorb), os%zorb(:,idim, iorb))
          end do
        end do
        safe_deallocate_a(zpsi)
 
        safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
        safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
      end if
      os%use_submesh = .true. ! We are now on a submesh
      os%allocated_on_mesh = .false.
      this%max_np = max(this%max_np, os%sphere%np)
    end do
 
    this%basis%use_submesh = .true.
 
    ! If we use GPUs, we need to transfert the orbitals on the device
    if (accel_is_enabled() .and. st%d%dim == 1) then
      do ios = 1, this%norbsets
        os => this%orbsets(ios)
 
        os%ldorbs = max(accel_padded_size(os%sphere%np), 1)
        os%ldorbs_eorb = max(accel_padded_size(os%sphere%np), 1)
        if (states_are_real(st)) then
          call accel_create_buffer(os%dbuff_orb, accel_mem_read_only, type_float, os%ldorbs*os%norbs)
        else
          call accel_create_buffer(os%zbuff_orb, accel_mem_read_only, type_cmplx, os%ldorbs*os%norbs)
          safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
 
          do ik= st%d%kpt%start, st%d%kpt%end
            call accel_create_buffer(os%buff_eorb(ik), accel_mem_read_only, type_cmplx, os%ldorbs_eorb*os%norbs)
          end do
        end if
 
        call accel_create_buffer(os%sphere%buff_map, accel_mem_read_only, type_integer, max(os%sphere%np, 1))
        call accel_write_buffer(os%sphere%buff_map, os%sphere%np, os%sphere%map)
 
        do iorb = 1, os%norbs
          if(states_are_complex(st)) then
            call accel_write_buffer(os%zbuff_orb, os%sphere%np, os%zorb(:, 1, iorb), &
              offset = (iorb - 1)*os%ldorbs)
          else
            call accel_write_buffer(os%dbuff_orb, os%sphere%np, os%dorb(:, 1, iorb), &
              offset = (iorb - 1)*os%ldorbs)
          end if
        end do
      end do
    end if
 
 
    message(1) = "Debug: Loading DFT+U basis from states done."
    call messages_info(1, debug_only=.true.)
 
    pop_sub(lda_u_loadbasis)
  end subroutine lda_u_loadbasis
 
  subroutine build_symmetrization_map(this, ions, gr)
    type(lda_u_t),        intent(inout) :: this
    type(ions_t),         intent(in)    :: ions
    type(grid_t),         intent(in)    :: gr
 
    integer :: nsym, iop, ios, iatom, iatom_sym, ios_sym
 
    push_sub(build_symmetrization_map)
 
    nsym = ions%symm%nops
    safe_allocate(this%inv_map_symm(1:this%norbsets, 1:nsym))
    this%inv_map_symm = -1
 
    this%nsym = nsym
 
    do ios = 1, this%norbsets
      iatom = this%orbsets(ios)%iatom
      do iop = 1, nsym
        iatom_sym = ions%inv_map_symm_atoms(iatom, iop)
 
        do ios_sym = 1, this%norbsets
          if (this%orbsets(ios_sym)%iatom == iatom_sym .and. this%orbsets(ios_sym)%norbs == this%orbsets(ios)%norbs &
            .and. this%orbsets(ios_sym)%nn == this%orbsets(ios)%nn .and. this%orbsets(ios_sym)%ll == this%orbsets(ios)%ll &
            .and. is_close(this%orbsets(ios_sym)%jj, this%orbsets(ios)%jj)) then
            this%inv_map_symm(ios, iop) = ios_sym
            exit
          end if
        end do
        assert(this%inv_map_symm(ios, iop) > 0)
      end do
    end do
 
    safe_allocate(this%symm_weight(1:this%maxnorbs, 1:this%maxnorbs, 1:this%nsym, 1:this%norbsets))
    this%symm_weight = m_zero
 
    do ios = 1, this%norbsets
      ! s-orbitals
      if (this%orbsets(ios)%norbs == 1) then
        this%symm_weight(1,1, 1:this%nsym, ios) = m_one
        cycle
      end if
 
      ! Not implemented yet
      if (this%orbsets(ios)%ndim > 1) cycle
 
      call orbitals_get_symm_weight(this%orbsets(ios), ions%space, ions%latt, gr, ions%symm, this%symm_weight(:,:,:,ios))
    end do
 
    pop_sub(build_symmetrization_map)
  end subroutine build_symmetrization_map
 
  subroutine orbitals_get_symm_weight(os, space, latt, gr, symm, weight)
    type(orbitalset_t), intent(in) :: os
    type(space_t), intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    type(grid_t), intent(in) :: gr
    type(symmetries_t), intent(in) :: symm
    real(real64), contiguous, intent(inout) :: weight(:,:,:)
 
    integer :: im, imp, iop, mm
    real(real64), allocatable :: orb(:,:), orb_sym(:), ylm(:)
    type(submesh_t) :: sphere
    real(real64) :: rc, norm, origin(space%dim)
 
    push_sub(orbitals_get_symm_weight)
 
    assert(os%ndim == 1)
 
    safe_allocate(orb_sym(1:gr%np))
    safe_allocate(orb(1:gr%np, 1:os%norbs))
 
    assert(2*os%ll+1 == os%norbs)
 
    ! We generate an artificial submesh to compute the symmetries on it
    ! The radius is such that we fit in 20 points
    rc = (50.0_real64 * m_three/m_four/m_pi*product(gr%spacing))**m_third
    origin = m_zero
    call submesh_init(sphere, space, gr, latt, origin, rc)
 
    safe_allocate(ylm(1:sphere%np))
 
    ! We then compute the spherical harmonics in this submesh
    do im = 1, os%norbs
      mm = im-1-os%ll
      call loct_ylm(sphere%np, sphere%rel_x(1,1), sphere%r(1), os%ll, mm, ylm(1))
      orb(:,im) = m_zero
      call submesh_add_to_mesh(sphere, ylm, orb(:,im))
      norm = dmf_nrm2(gr, orb(:,im))
      call lalg_scal(gr%np, m_one / norm, orb(:,im))
    end do
    safe_deallocate_a(ylm)
 
    ! Then we put them on the grid, rotate them
    do im = 1, os%norbs
      do iop = 1, symm%nops
        call dgrid_symmetrize_single(gr, iop, orb(:,im), orb_sym)
 
        do imp = 1, os%norbs
          weight(im, imp, iop) = dmf_dotp(gr, orb(:,imp), orb_sym, reduce=.false.)
        end do
      end do
    end do
 
    call gr%allreduce(weight)
 
    safe_deallocate_a(orb)
    safe_deallocate_a(orb_sym)
 
    pop_sub(orbitals_get_symm_weight)
  end subroutine orbitals_get_symm_weight
 
#include "dft_u_noncollinear_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "lda_u_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "lda_u_inc.F90"
end module lda_u_oct_m