lcao.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module lcao_oct_m
  use atom_oct_m
  use atomic_orbital_oct_m
  use batch_oct_m
  use blacs_oct_m
  use blacs_proc_grid_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use ks_potential_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use lapack_oct_m
  use loct_oct_m
  use magnetic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_debug_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use quickrnd_oct_m
  use scalapack_oct_m
  use smear_oct_m
  use space_oct_m
  use species_oct_m
  use species_pot_oct_m
  use species_factory_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
  use states_elec_io_oct_m
  use submesh_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use v_ks_oct_m
  use varinfo_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
  use xc_sic_oct_m
 
  implicit none
 
  private
  public ::             &
    lcao_t,             &
    lcao_init,          &
    lcao_init_orbitals, &
    lcao_wf,            &
    lcao_run,           &
    lcao_end,           &
    lcao_is_available,  &
    lcao_num_orbitals
 
  integer, parameter ::        &
    INITRHO_PARAMAGNETIC  = 1, &
    initrho_ferromagnetic = 2, &
    initrho_random        = 3, &
    initrho_userdef       = 77
 
  type lcao_t
    private
    integer                 :: mode
    logical                 :: complex_ylms
    logical                 :: initialized
    integer                 :: norbs
    integer                 :: maxorbs
    integer,    allocatable :: atom(:)
    integer,    allocatable :: level(:)
    integer,    allocatable :: ddim(:)
    logical                 :: alternative
    logical                 :: derivative
    integer,    allocatable :: cst(:, :)
    integer,    allocatable :: ck(:, :)
    real(4),    allocatable :: dbuff_single(:, :, :, :)
    complex(4), allocatable :: zbuff_single(:, :, :, :)
    real(real64),      allocatable :: dbuff(:, :, :, :)
    complex(real64),      allocatable :: zbuff(:, :, :, :)
    logical                 :: save_memory
    logical                 :: initialized_orbitals
    real(real64)            :: orbital_scale_factor
 
    logical,    allocatable :: is_empty(:)
    !                                      !< This occurs in domain parallelization if the atom is not
    !                                      !< in the local domain
 
    logical               :: keep_orb
    real(real64), allocatable :: radius(:)
    real(real64)          :: lapdist
    integer               :: mult
    integer               :: maxorb
    integer, allocatable  :: basis_atom(:)
    integer, allocatable  :: basis_orb(:)
    integer, allocatable  :: atom_orb_basis(:, :)
    integer, allocatable  :: norb_atom(:)
    logical               :: parallel
    integer               :: lsize(1:2)
    integer               :: nproc(1:2)
    integer               :: myroc(1:2)
    integer               :: desc(1:BLACS_DLEN)
    logical, allocatable  :: calc_atom(:)
    real(real64)          :: diag_tol
    type(submesh_t), allocatable :: sphere(:)
    type(batch_t),   allocatable :: orbitals(:)
    logical,         allocatable :: is_orbital_initialized(:)
 
    integer                 :: gmd_opt
  end type lcao_t
 
 
 
  real(real64), parameter, public :: default_eigenval = 1.0e10_real64
 
contains
 
  ! ---------------------------------------------------------
  subroutine lcao_init(this, namespace, space, gr, ions, st, st_start)
    type(lcao_t),           intent(out) :: this
    type(namespace_t),      intent(in)  :: namespace
    type(electron_space_t), intent(in)  :: space
    type(grid_t),           intent(in)  :: gr
    type(ions_t),           intent(in)  :: ions
    type(states_elec_t),    intent(in)  :: st
    integer,                intent(in)  :: st_start
 
    integer :: ia, n, iorb, jj, maxj, idim
    integer :: ii, ll, mm, norbs, ii_tmp
    integer :: mode_default
    real(real64)   :: max_orb_radius, maxradius
    integer :: iunit_o
 
    push_sub(lcao_init)
 
    this%initialized = .true.
    this%gmd_opt = initrho_userdef
 
    ! The initial LCAO calculation is done by default if we have species representing atoms.
    ! Otherwise, it is not the default value and has to be enforced in the input file.
    mode_default = option__lcaostart__lcao_states
    if (ions%only_user_def) mode_default = option__lcaostart__lcao_none
 
    !%Variable LCAOStart
    !%Type integer
    !%Section SCF::LCAO
    !%Description
    !% Before starting a SCF calculation, <tt>Octopus</tt> can perform
    !% a linear combination of atomic orbitals (LCAO) calculation.
    !% These can provide <tt>Octopus</tt> with a good set
    !% of initial wavefunctions and with a new guess for the density.
    !% (Up to the current version, only a minimal basis set is used.)
    !% The default is <tt>lcao_states</tt> if at least one species representing an atom is present.
    !% The default is <tt>lcao_none</tt> if all species are <tt>species_user_defined</tt>,
    !% <tt>species_charge_density</tt>, <tt>species_from_file</tt>, or <tt>species_jellium_slab</tt>.
    !%
    !% The initial guess densities for LCAO are taken from the atomic orbitals for pseudopotential species;
    !% from the natural charge density for <tt>species_charge_density</tt>, <tt>species_point</tt>,
    !% <tt>species_jellium</tt>, and <tt>species_jellium_slab</tt>;
    !% or uniform for <tt>species_full_delta</tt>, <tt>species_full_gaussian</tt>,
    !% <tt>species_user_defined</tt>, or <tt>species_from_file</tt>.
    !% Pseudopotential species use the pseudo-wavefunctions as orbitals, full-potential atomic species
    !% (<tt>species_full_delta</tt> and <tt>species_full_gaussian</tt>) use hydrogenic wavefunctions, and
    !% others use harmonic-oscillator wavefunctions.
    !%
    !% The LCAO implementation currently does not simultaneously support k-point parallelism and states parallelism.
    !% You must select one or the other, or alternatively use lcao_none (random wavefunctions) as an initial guess.
    !%
    !% Note: Some pseudopotential files (CPI, FHI for example) do not
    !% contain full information about the orbitals. In this case,
    !% Octopus generates the starting density from the normalized
    !% square root of the local potential. If no orbitals are
    !% available at all from the pseudopotential files, Octopus will
    !% not be able to perform an LCAO and the initial states will be
    !% randomized.
    !%
    !%Option lcao_none 0
    !% Do not perform a LCAO calculation before the SCF cycle. Instead use random wavefunctions.
    !%Option lcao_states 2
    !% Do a LCAO calculation before the SCF cycle and use the resulting wavefunctions as
    !% initial wavefunctions without changing the guess density.
    !% This will speed up the convergence of the eigensolver during the first SCF iterations.
    !%Option lcao_full 3
    !% Do a LCAO calculation before the SCF cycle and use the LCAO wavefunctions to build a new
    !% guess density and a new KS potential.
    !% Using the LCAO density as a new guess density may improve the convergence, but can
    !% also slow it down or yield wrong results (especially for spin-polarized calculations).
    !%Option lcao_states_batch 4
    !% (Experimental) Same as lcao_states, but faster for GPUs, large systems and parallel in states.
    !% It is not working for spinors, however.
    !% This is not used for unocc calculations at the moment.
    !%End
    call parse_variable(namespace, 'LCAOStart', mode_default, this%mode)
    if (.not. varinfo_valid_option('LCAOStart', this%mode)) call messages_input_error(namespace, 'LCAOStart')
 
    call messages_print_var_option('LCAOStart', this%mode, namespace=namespace)
 
    ! LCAO alternative not implemented for st_start > 1
    if (this%mode == option__lcaostart__lcao_states_batch .and. st_start > 1) then
      message(1) = "LCAOStart = lcao_states_batch not compatible with this run."
      message(2) = "Please use LCAOStart = lcao_states instead"
      call messages_fatal(2, namespace=namespace)
    end if
 
    !TODO: Remove the alternative flag
    this%alternative = this%mode == option__lcaostart__lcao_states_batch
 
    if (this%mode == option__lcaostart__lcao_none) then
      pop_sub(lcao_init)
      return
    end if
 
    call messages_obsolete_variable(namespace, 'LCAOAlternative')
 
 
    ! DAS: For spinors, you will always get magnetization in (1, 0, 0) direction, and the
    ! eigenvalues will be incorrect. This is due to confusion between spins and spinors in the code.
    if (st%d%ispin == spinors .and. this%alternative) then
      message(1) = "LCAOStart = lcao_states_batch is not working for spinors."
      call messages_fatal(1, namespace=namespace)
    end if
    if (space%is_periodic() .and. this%alternative) then
      call messages_experimental("LCAOStart = lcao_states_batch in periodic systems", namespace=namespace)
      ! specifically, if you get the message about submesh radius > box size, results will probably be totally wrong.
    end if
 
    !%Variable LCAOComplexYlms
    !%Type logical
    !%Default false
    !%Section SCF::LCAO
    !%Description
    !% If set to true, and using complex states, complex spherical harmonics will be used, <i>i.e.</i>
    !% with <math>e^{\pm i m \phi}</math>.
    !% If false, real spherical harmonics with <math>\sin(m \phi)</math> or <math>\cos(m \phi)</math> are used.
    !% This variable will make it more likely to get states that are eigenvectors of the <math>L_z</math>
    !% operator, with a definite angular momentum.
    !%End
 
    if (states_are_complex(st)) then
      call parse_variable(namespace, 'LCAOComplexYlms', .false., this%complex_ylms)
    else
      this%complex_ylms = .false.
    end if
 
    !%Variable LCAOSaveMemory
    !%Type logical
    !%Default false
    !%Section SCF::LCAO
    !%Description
    !% If set to true, the LCAO will allocate extra memory needed in single precision instead of
    !% double precision.
    !%End
    call parse_variable(namespace, 'LCAOSaveMemory', .false., this%save_memory)
 
 
    if (debug%info .and. mpi_world%is_root()) then
      call io_mkdir('debug/lcao', namespace)
      iunit_o = io_open('debug/lcao/orbitals', namespace, action='write')
      write(iunit_o,'(7a6)') 'iorb', 'atom', 'level', 'i', 'l', 'm', 'spin'
    end if
 
    if (.not. this%alternative) then
 
      !%Variable LCAOScaleFactor
      !%Type float
      !%Default 1.0
      !%Section SCF::LCAO
      !%Description
      !% The coordinates of the atomic orbitals used by the LCAO
      !% procedure will be rescaled by the value of this variable. 1.0 means no rescaling.
      !%End
      call parse_variable(namespace, 'LCAOScaleFactor', m_one, this%orbital_scale_factor)
      call messages_print_var_value('LCAOScaleFactor', this%orbital_scale_factor, namespace=namespace)
 
      !%Variable LCAOMaximumOrbitalRadius
      !%Type float
      !%Default 20.0 a.u.
      !%Section SCF::LCAO
      !%Description
      !% The LCAO procedure will ignore orbitals that have an
      !% extent greater that this value.
      !%End
      call parse_variable(namespace, 'LCAOMaximumOrbitalRadius', 20.0_real64, max_orb_radius, unit = units_inp%length)
      call messages_print_var_value('LCAOMaximumOrbitalRadius', max_orb_radius, units_out%length, namespace=namespace)
 
      ! count the number of orbitals available
      maxj = 0
      this%maxorbs = 0
      do ia = 1, ions%natoms
        maxj = max(maxj, ions%atom(ia)%species%get_niwfs())
        this%maxorbs = this%maxorbs + ions%atom(ia)%species%get_niwfs()
      end do
 
      this%maxorbs = this%maxorbs*st%d%dim
 
      if (this%maxorbs == 0) then
        call messages_write('The are no atomic orbitals available, cannot do LCAO.')
        call messages_warning(namespace=namespace)
        this%mode = option__lcaostart__lcao_none
        pop_sub(lcao_init)
        return
      end if
 
      ! generate tables to know which indices each atomic orbital has
 
      safe_allocate( this%atom(1:this%maxorbs))
      safe_allocate(this%level(1:this%maxorbs))
      safe_allocate( this%ddim(1:this%maxorbs))
 
      safe_allocate(this%is_empty(1:this%maxorbs))
      this%is_empty = .false.
 
      ! this is determined by the stencil we are using and the spacing
      this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
 
      ! calculate the radius of each orbital
      safe_allocate(this%radius(1:ions%natoms))
 
      do ia = 1, ions%natoms
        norbs = ions%atom(ia)%species%get_niwfs()
 
        maxradius = m_zero
        do iorb = 1, norbs
          call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
          ! For all-electron species, we need to use the principal quantum number
          if(ions%atom(ia)%species%is_full()) call ions%atom(ia)%species%get_iwf_n( iorb, 1, ii)
          maxradius = max(maxradius, ions%atom(ia)%species%get_iwf_radius( ii, is = 1))
        end do
 
        maxradius = min(maxradius, m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
 
        this%radius(ia) = maxradius
      end do
 
 
      ! Each atom provides niwfs pseudo-orbitals (this number is given in
      ! ions%atom(ia)%species%niwfs for atom number ia). This number is
      ! actually multiplied by two in case of spin-unrestricted or spinors
      ! calculations.
      !
      ! The pseudo-orbitals are placed in order in the following way (Natoms
      ! is the total number of atoms).
      !
      ! n = 1 => first orbital of atom 1,
      ! n = 2 => first orbital of atom 2.
      ! n = 3 => first orbital of atom 3.
      ! ....
      ! n = Natoms => first orbital of atom Natoms
      ! n = Natoms + 1 = > second orbital of atom 1
      ! ....
      !
      ! If at some point in this loop an atom pseudo cannot provide the corresponding
      ! orbital (because the niws orbitals have been exhausted), it moves on to the following
      ! atom.
      !
      ! In the spinors case, it changes a bit:
      !
      ! n = 1 => first spin-up orbital of atom 1, assigned to the spin-up component of the spinor.
      ! n = 2 => first spin-down orbital of atom 1, assigned to the spin-down component of the spinor.
      ! n = 3 => first spin-up orbital of atom 2, assigned to the spin-up component of the spinor.
 
      iorb = 1
      do jj = 1, maxj
        do ia = 1, ions%natoms
          do idim = 1,st%d%dim
            if (jj > ions%atom(ia)%species%get_niwfs()) cycle
            call ions%atom(ia)%species%get_iwf_ilm(jj, idim, ii, ll, mm)
            ! For all-electron species, we need to use the principal quantum number
            if(ions%atom(ia)%species%is_full()) then
              ii_tmp = ii
              call ions%atom(ia)%species%get_iwf_n( ii_tmp, 1, ii)
            end if
            if (this%orbital_scale_factor*ions%atom(ia)%species%get_iwf_radius( ii, is = 1) >= max_orb_radius) cycle
 
            this%atom(iorb) = ia
            this%level(iorb) = jj
            this%ddim(iorb) = idim
 
            if (debug%info .and. mpi_world%is_root()) then
              write(iunit_o,'(7i6)') iorb, this%atom(iorb), this%level(iorb), ii, ll, mm, this%ddim(iorb)
            end if
 
            iorb = iorb + 1
          end do
        end do
      end do
 
      if (debug%info .and. mpi_world%is_root()) then
        call io_close(iunit_o)
      end if
 
      ! some orbitals might have been removed because of their radii
      if (this%maxorbs /= iorb - 1) then
        call messages_write('Info: ')
        call messages_write(this%maxorbs - iorb + 1)
        call messages_write(' of ')
        call messages_write(this%maxorbs)
        call messages_write(' orbitals cannot be used for the LCAO calculation,')
        call messages_new_line()
        call messages_write('      their radii exceeds LCAOMaximumOrbitalRadius (')
        call messages_write(max_orb_radius, units = units_out%length, fmt = '(f6.1)')
        call messages_write(').')
        call messages_warning(namespace=namespace)
 
        this%maxorbs = iorb - 1
      end if
 
      if (this%maxorbs < st%nst) then
        call messages_write('Cannot do LCAO for all states because there are not enough atomic orbitals.')
        call messages_new_line()
 
        call messages_write('Required: ')
        call messages_write(st%nst)
        call messages_write('. Available: ')
        call messages_write(this%maxorbs)
        call messages_write('. ')
        call messages_write(st%nst - this%maxorbs)
        call messages_write(' orbitals will be randomized.')
        call messages_warning(namespace=namespace)
      end if
 
      !%Variable LCAODimension
      !%Type integer
      !%Section SCF::LCAO
      !%Description
      !% (Only applies if <tt>LCAOStart /= lcao_states_batch</tt>.)
      !% Before starting the SCF cycle, an initial LCAO calculation can be performed
      !% in order to obtain reasonable initial guesses for spin-orbitals and densities.
      !% For this purpose, the code calculates a number of atomic orbitals.
      !% The number available for a species described by a pseudopotential is all the
      !% orbitals up the maximum angular momentum in the pseudopotential, minus any orbitals that
      !% are found to be unbound. For non-pseudopotential species, the number is equal to
      !% twice the valence charge.
      !% The default dimension for the LCAO basis
      !% set will be the sum of all these numbers, or twice the number of required orbitals
      !% for the full calculation, whichever is less.
      !%
      !% This dimension however can be changed by making use of this
      !% variable. Note that <tt>LCAODimension</tt> cannot be smaller than the
      !% number of orbitals needed in the full calculation -- if
      !% <tt>LCAODimension</tt> is smaller, it will be silently increased to meet
      !% this requirement. In the same way, if <tt>LCAODimension</tt> is larger
      !% than the available number of atomic orbitals, it will be
      !% reduced. If you want to use the largest possible number, set
      !% <tt>LCAODimension</tt> to a negative number.
      !%End
      call parse_variable(namespace, 'LCAODimension', 0, n)
 
      if (n > 0 .and. n <= st%nst .and. st%nst <= this%maxorbs) then
        ! n was made too small
        this%norbs = st%nst
      else if (n > st%nst .and. n <= this%maxorbs) then
        ! n is a reasonable value
        this%norbs = n
      else if (n == 0) then
        ! using the default
        this%norbs = min(this%maxorbs, 2*st%nst)
      else
        ! n was negative, or greater than maxorbs
        this%norbs = this%maxorbs
      end if
 
      assert(this%norbs <= this%maxorbs)
 
      safe_allocate(this%cst(1:this%norbs, 1:st%d%spin_channels))
      safe_allocate(this%ck(1:this%norbs, 1:st%d%spin_channels))
      this%initialized_orbitals = .false.
    else
      call lcao2_init()
    end if
 
    pop_sub(lcao_init)
 
  contains
 
    subroutine lcao2_init()
      integer :: iatom, iorb, norbs
      real(real64)   :: maxradius
      integer :: ibasis
#ifdef HAVE_SCALAPACK
      integer :: jatom, jorb, jbasis, ilbasis, jlbasis, proc(1:2), info, nbl
#endif
      push_sub(lcao_init.lcao2_init)
 
      call messages_write('Info: Using LCAO batched implementation.')
      call messages_info(namespace=namespace)
 
      call messages_experimental('LCAO alternative implementation', namespace=namespace)
 
      !%Variable LCAOKeepOrbitals
      !%Type logical
      !%Default yes
      !%Section SCF::LCAO
      !%Description
      !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>.
      !% If set to yes (the default) Octopus keeps atomic orbitals in
      !% memory during the LCAO procedure. If set to no, the orbitals
      !% are generated each time that they are needed, increasing
      !% computational time but saving memory.
      !%
      !% When set to yes, Octopus prints the amount of memory per node
      !% that is required to store the orbitals.
      !%
      !%End
      call parse_variable(namespace, 'LCAOKeepOrbitals', .true., this%keep_orb)
 
      !%Variable LCAOExtraOrbitals
      !%Type logical
      !%Default false
      !%Section SCF::LCAO
      !%Description
      !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>, and all species are pseudopotentials.
      !% (experimental) If this variable is set to yes, the LCAO
      !% procedure will add an extra set of numerical orbitals (by
      !% using the derivative of the radial part of the original
      !% orbitals). Note that this corresponds roughly to adding orbitals
      !% with higher principal quantum numbers, but the same angular momentum.
      !% This option may cause problems for unoccupied states since you may miss
      !% some lower-lying states which correspond to higher angular momenta instead
      !% of higher principal quantum number.
      !%End
      call parse_variable(namespace, 'LCAOExtraOrbitals', .false., this%derivative)
 
      ! DAS: if you calculate the Na atom this way, spin-polarized, with just one unoccupied state,
      ! you will obtain states (up and down) which are actually the 10th states if you start with
      ! random wavefunctions! We really need to implement taking the derivative of the angular part
      ! instead to be sure of getting decent results!
 
      if (this%derivative) then
        call messages_experimental('LCAO extra orbitals', namespace=namespace)
 
        if (st%nst * st%smear%el_per_state > st%qtot) then
          message(1) = "Lower-lying empty states may be missed with LCAOExtraOrbitals."
          call messages_warning(1, namespace=namespace)
        end if
      end if
 
      !%Variable LCAODiagTol
      !%Type float
      !%Default 1e-10
      !%Section SCF::LCAO
      !%Description
      !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>.
      !% The tolerance for the diagonalization of the LCAO Hamiltonian.
      !%End
      call parse_variable(namespace, 'LCAODiagTol', 1e-10_real64, this%diag_tol)
 
      if (this%derivative) then
        this%mult = 2
      else
        this%mult = 1
      end if
 
      safe_allocate(this%sphere(1:ions%natoms))
      safe_allocate(this%orbitals(1:ions%natoms))
      safe_allocate(this%is_orbital_initialized(1:ions%natoms))
      this%is_orbital_initialized = .false.
 
      safe_allocate(this%norb_atom(1:ions%natoms))
 
      this%maxorb = 0
      this%norbs = 0
      do iatom = 1, ions%natoms
        this%norb_atom(iatom) = this%mult*ions%atom(iatom)%species%get_niwfs()
        this%maxorb = max(this%maxorb, ions%atom(iatom)%species%get_niwfs())
        this%norbs = this%norbs + ions%atom(iatom)%species%get_niwfs()
      end do
 
      this%maxorb = this%maxorb*this%mult
      this%norbs = this%norbs*this%mult
 
      safe_allocate(this%basis_atom(1:this%norbs))
      safe_allocate(this%basis_orb(1:this%norbs))
      safe_allocate(this%atom_orb_basis(1:ions%natoms, 1:this%maxorb))
 
      ! Initialize the mapping between indices
 
      ibasis = 0
      do iatom = 1, ions%natoms
        norbs = ions%atom(iatom)%species%get_niwfs()
 
        do iorb = 1, this%mult*norbs
          ibasis = ibasis + 1
          this%atom_orb_basis(iatom, iorb) = ibasis
          this%basis_atom(ibasis) = iatom
          this%basis_orb(ibasis) = iorb
 
          ! no stored spin index in alternative mode
          if (debug%info .and. mpi_world%is_root()) then
            call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
            write(iunit_o,'(7i6)') ibasis, iatom, iorb, ii, ll, mm, 1
          end if
        end do
      end do
 
      if (debug%info .and. mpi_world%is_root()) then
        call io_close(iunit_o)
      end if
 
      ! this is determined by the stencil we are using and the spacing
      this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
 
      ! calculate the radius of each orbital
      safe_allocate(this%radius(1:ions%natoms))
 
      do iatom = 1, ions%natoms
        norbs = ions%atom(iatom)%species%get_niwfs()
 
        maxradius = m_zero
        do iorb = 1, norbs
          call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
          ! For all-electron species, we need to use the principal quantum number
          if(ions%atom(iatom)%species%is_full()) call ions%atom(iatom)%species%get_iwf_n( iorb, 1, ii)
          maxradius = max(maxradius, ions%atom(iatom)%species%get_iwf_radius( ii, is = 1))
        end do
 
        if (this%derivative) maxradius = maxradius + this%lapdist
 
        maxradius = min(maxradius, m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
 
        this%radius(iatom) = maxradius
      end do
 
      safe_allocate(this%calc_atom(1:ions%natoms))
      this%calc_atom = .true.
 
      ! initialize parallel data
#ifndef HAVE_SCALAPACK
      this%parallel = .false.
#else
      this%parallel = (st%parallel_in_states .or. gr%parallel_in_domains) &
        .and. .not. blacs_proc_grid_null(st%dom_st_proc_grid)
 
      if (this%parallel) then
        nbl = min(16, this%norbs)
 
        ! The size of the distributed matrix in each node
        this%lsize(1) = max(1, numroc(this%norbs, nbl, st%dom_st_proc_grid%myrow, 0, st%dom_st_proc_grid%nprow))
        this%lsize(2) = max(1, numroc(this%norbs, nbl, st%dom_st_proc_grid%mycol, 0, st%dom_st_proc_grid%npcol))
 
        this%nproc(1) = st%dom_st_proc_grid%nprow
        this%nproc(2) = st%dom_st_proc_grid%npcol
        this%myroc(1) = st%dom_st_proc_grid%myrow
        this%myroc(2) = st%dom_st_proc_grid%mycol
 
        call descinit(this%desc(1), this%norbs, this%norbs, nbl, nbl, 0, 0, &
          st%dom_st_proc_grid%context, this%lsize(1), info)
 
        if (info /= 0) then
          write(message(1), '(a,i6)') 'descinit for BLACS failed with error code ', info
          call messages_fatal(1, namespace=namespace)
        end if
 
        this%calc_atom = .false.
        do iatom = 1, ions%natoms
          ibasis = this%atom_orb_basis(iatom, 1)
 
          do jatom = 1, ions%natoms
            jbasis = this%atom_orb_basis(jatom, 1)
 
            do iorb = 1, this%norb_atom(iatom)
              do jorb = 1, this%norb_atom(jatom)
                call lcao_local_index(this,  ibasis - 1 + iorb,  jbasis - 1 + jorb, &
                  ilbasis, jlbasis, proc(1), proc(2))
 
                this%calc_atom(this%basis_atom(jbasis)) = &
                  this%calc_atom(this%basis_atom(jbasis)) .or. proc(2) == this%myroc(2)
 
              end do
            end do
 
          end do
        end do
 
      end if
#endif
 
      pop_sub(lcao_init.lcao2_init)
    end subroutine lcao2_init
 
  end subroutine lcao_init
 
 
  ! ---------------------------------------------------------
  subroutine lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r, known_lower_bound)
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,        optional, intent(in)    :: st_start
    real(real64),   optional, intent(in)    :: lmm_r
    logical,        optional, intent(out)   :: known_lower_bound
 
    type(lcao_t) :: lcao
    integer :: st_start_random, required_min_nst
    logical :: lcao_done
    logical :: is_orbital_dependent
 
    push_sub(lcao_run)
 
    if (present(known_lower_bound)) known_lower_bound = .false.
 
    if (present(st_start)) then
      ! If we are doing unocc calculation, do not mess with the correct eigenvalues
      ! of the occupied states.
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        calc_eigenval=.not. present(st_start), calc_current=.false.)
 
      assert(st_start >= 1)
      if (st_start > st%nst) then ! nothing to be done in LCAO
        if (present(known_lower_bound)) known_lower_bound = .true.
        pop_sub(lcao_run)
        return
      end if
    end if
 
    call profiling_in('LCAO_RUN')
 
    call lcao_init(lcao, namespace, space, gr, ions, st, optional_default(st_start, 1))
 
    call lcao_init_orbitals(lcao, namespace, st, gr, ions, start = st_start)
 
    ! By default, we want to use vxc for the LCAO.
    ! However, we can only do this if vxc depends on the density only
    ! For cases like MGGA or hybrids, OEP, we cannot do this
    is_orbital_dependent = (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
      .or. (ks%theory_level == kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
      .or. (ks%theory_level == generalized_kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
      .or. ks%sic%level == sic_pz_oep)
 
    if (.not. present(st_start)) then
      call lcao_guess_density(lcao, namespace, st, gr, hm, ions, st%qtot, st%d%ispin, st%rho)
 
      if (st%d%ispin > unpolarized) then
        assert(present(lmm_r))
        call write_magnetic_moments(gr, st, ions, gr%der%boundaries, lmm_r, namespace=namespace)
      end if
 
      ! set up Hamiltonian (we do not call v_ks_h_setup here because we do not want to
      ! overwrite the guess density)
      message(1) = 'Info: Setting up Hamiltonian.'
      call messages_info(1, namespace=namespace)
 
      ! get the effective potential (we don`t need the eigenvalues yet)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.false., &
        calc_current=.false., calc_energy=.false., force_semilocal=is_orbital_dependent)
      ! eigenvalues have nevertheless to be initialized to something
      ! This value must be larger that the highest eigenvalue from LCAO to get the correct occupations
      st%eigenval = default_eigenval
    end if
 
    lcao_done = .false.
 
    ! after initialized, can check that LCAO is possible
    if (lcao_is_available(lcao)) then
      lcao_done = .true.
 
      if (present(st_start)) then
        write(message(1),'(a,i8,a)') 'Performing LCAO for states ', st_start, ' and above'
        call messages_info(1, namespace=namespace)
      end if
 
      call lcao_wf(lcao, st, gr, ions, hm, namespace, start = st_start)
 
      ! In some rare cases, like bad pseudopotentials with unbound orbitals,
      ! we might not have enough orbitals to go up to the Fermi energy
      ! In this case, we cannot set the eigenvales to a huge value, as
      ! the smearing will not converge. We set them to zero then
      select case (st%d%ispin)
      case (unpolarized)
        required_min_nst = int(st%qtot/2)
      case (spin_polarized)
        required_min_nst = int(st%qtot/2)
      case (spinors)
        required_min_nst = int(st%qtot)
      end select
      if (st%smear%method /= smear_fixed_occ .and. st%smear%method /= smear_semiconductor) then
        if (lcao%norbs <= required_min_nst .and. lcao%norbs < st%nst) then
          st%eigenval(lcao%norbs+1:,:) = m_zero
        end if
      end if
 
      if (.not. present(st_start)) then
        call states_elec_fermi(st, namespace, gr)
        call states_elec_write_eigenvalues(min(st%nst, lcao%norbs), st, space, hm%kpoints, namespace=namespace)
 
        ! Update the density and the Hamiltonian
        if (lcao%mode == option__lcaostart__lcao_full) then
          call v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, &
            calc_eigenval = .false., calc_current=.false.)
          if (st%d%ispin > unpolarized) then
            assert(present(lmm_r))
            call write_magnetic_moments(gr, st, ions, gr%der%boundaries, lmm_r, namespace=namespace)
          end if
        end if
      end if
 
      if (present(known_lower_bound) .and. lcao%norbs >= st%nst) known_lower_bound = .true.
    end if
 
    if (.not. lcao_done .or. lcao%norbs < st%nst) then
 
      if (lcao_done) then
        st_start_random = lcao%norbs + 1
      else
        st_start_random = 1
      end if
      if (present(st_start)) st_start_random = max(st_start, st_start_random)
 
      if (st_start_random > 1) then
        write(message(1),'(a,i8,a)') 'Generating random wavefunctions for states ', st_start_random, ' and above'
        call messages_info(1, namespace=namespace)
      end if
 
      ! Randomly generate the initial wavefunctions.
      call states_elec_generate_random(st, gr, hm%kpoints, ist_start_ = st_start_random, normalized = .false.)
 
      ! For spinors, we rotate the random states in the local frame of the magnetization
      ! Unless the user specified the spin directions using InitialSpins
      if (.not. st%fixed_spins) then
        call rotate_random_states_to_local_frame(st, gr, st_start_random, lcao%gmd_opt)
      end if
 
      call messages_write('Orthogonalizing wavefunctions.')
      call messages_info(namespace=namespace)
      call states_elec_orthogonalize(st, namespace, gr)
 
      if (.not. lcao_done) then
        ! If we are doing unocc calculation, do not mess with the correct eigenvalues and occupations
        ! of the occupied states.
        call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
          calc_eigenval=.not. present(st_start), calc_current=.false., force_semilocal=is_orbital_dependent) ! get potentials
        if (.not. present(st_start)) then
          call states_elec_fermi(st, namespace, gr) ! occupations
        end if
 
      end if
 
    else if (present(st_start)) then
 
      if (st_start > 1) then
        call messages_write('Orthogonalizing wavefunctions.')
        call messages_info(namespace=namespace)
        call states_elec_orthogonalize(st, namespace, gr)
      end if
 
    end if
 
    call lcao_end(lcao)
 
 
    call profiling_out('LCAO_RUN')
    pop_sub(lcao_run)
  end subroutine lcao_run
 
  ! ---------------------------------------------------------
  subroutine lcao_end(this)
    type(lcao_t), intent(inout) :: this
 
    push_sub(lcao_end)
 
    safe_deallocate_a(this%calc_atom)
    safe_deallocate_a(this%norb_atom)
    safe_deallocate_a(this%basis_atom)
    safe_deallocate_a(this%basis_orb)
    safe_deallocate_a(this%atom_orb_basis)
    safe_deallocate_a(this%radius)
    safe_deallocate_a(this%sphere)
    safe_deallocate_a(this%orbitals)
 
    safe_deallocate_a(this%atom)
    safe_deallocate_a(this%level)
    safe_deallocate_a(this%ddim)
    safe_deallocate_a(this%cst)
    safe_deallocate_a(this%ck)
    safe_deallocate_a(this%dbuff_single)
    safe_deallocate_a(this%zbuff_single)
    safe_deallocate_a(this%dbuff)
    safe_deallocate_a(this%zbuff)
 
    this%initialized = .false.
    pop_sub(lcao_end)
  end subroutine lcao_end
 
 
  ! ---------------------------------------------------------
  subroutine lcao_wf(this, st, gr, ions, hm, namespace, start)
    type(lcao_t),             intent(inout) :: this
    type(states_elec_t),      intent(inout) :: st
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(namespace_t),        intent(in)    :: namespace
    integer, optional,        intent(in)    :: start
 
    integer :: start_
 
    assert(this%initialized)
 
    call profiling_in("LCAO")
    push_sub(lcao_wf)
 
    start_ = 1
    if (present(start)) start_ = start
 
    if (this%alternative) then
      if (states_are_real(st)) then
        call dlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
      else
        call zlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
      end if
    else
      if (states_are_real(st)) then
        call dlcao_wf(this, st, gr, ions, hm, namespace, start_)
      else
        call zlcao_wf(this, st, gr, ions, hm, namespace, start_)
      end if
    end if
    pop_sub(lcao_wf)
    call profiling_out("LCAO")
  end subroutine lcao_wf
 
 
  ! ---------------------------------------------------------
  logical function lcao_is_available(this) result(available)
    type(lcao_t), intent(in) :: this
 
    push_sub(lcao_is_available)
 
    available = this%initialized .and. this%mode /= option__lcaostart__lcao_none &
      .and. this%norbs > 0
 
    pop_sub(lcao_is_available)
  end function lcao_is_available
 
 
  ! ---------------------------------------------------------
  integer function lcao_num_orbitals(this) result(norbs)
    type(lcao_t), intent(in) :: this
 
    push_sub(lcao_num_orbitals)
    norbs = this%norbs
 
    pop_sub(lcao_num_orbitals)
  end function lcao_num_orbitals
 
  ! ---------------------------------------------------------
 
  subroutine lcao_local_index(this, ig, jg, il, jl, prow, pcol)
    type(lcao_t), intent(in)  :: this
    integer,      intent(in)  :: ig
    integer,      intent(in)  :: jg
    integer,      intent(out) :: il
    integer,      intent(out) :: jl
    integer,      intent(out) :: prow
    integer,      intent(out) :: pcol
 
    ! no PUSH_SUB, called too often
#ifdef HAVE_SCALAPACK
    call infog2l(ig, jg, this%desc(1), this%nproc(1), this%nproc(2), this%myroc(1), this%myroc(2), &
      il, jl, prow, pcol)
#else
    il = ig
    jl = jg
    prow = 0
    pcol = 0
#endif
 
  end subroutine lcao_local_index
 
  ! ---------------------------------------------------------
 
  subroutine lcao_alt_end_orbital(this, iatom)
    type(lcao_t),   intent(inout) :: this
    integer,        intent(in)    :: iatom
 
    push_sub(lcao_alt_end_orbital)
 
    if (this%is_orbital_initialized(iatom)) then
      call this%orbitals(iatom)%end()
      this%is_orbital_initialized(iatom) = .false.
    end if
 
    pop_sub(lcao_alt_end_orbital)
 
  end subroutine lcao_alt_end_orbital
 
  ! ---------------------------------------------------------
 
  subroutine lcao_atom_density(this, st, mesh, ions, iatom, spin_channels, rho)
    type(lcao_t),             intent(inout) :: this
    type(states_elec_t),      intent(in)    :: st
    class(mesh_t),            intent(in)    :: mesh
    type(ions_t),     target, intent(in)    :: ions
    integer,                  intent(in)    :: iatom
    integer,                  intent(in)    :: spin_channels
    real(real64),             intent(inout) :: rho(:, :)
 
    real(real64), allocatable :: dorbital(:, :)
    complex(real64), allocatable :: zorbital(:, :)
    real(real64), allocatable :: factors(:)
    real(real64) :: factor, aa
    integer :: iorb, ip, ii, ll, mm, ispin
    type(ps_t), pointer :: ps
    logical :: use_stored_orbitals
 
    push_sub(lcao_atom_density)
 
    rho = m_zero
 
    use_stored_orbitals = ions%atom(iatom)%species%is_ps() &
      .and. states_are_real(st) .and. spin_channels == 1 .and. lcao_is_available(this) &
      .and. st%d%dim == 1 .and. .not. ions%space%is_periodic()
 
 
    ! we can use the orbitals we already have calculated
    if (use_stored_orbitals) then
      !There is no periodic copies here, so this will not work for periodic systems
      assert(.not. ions%space%is_periodic())
 
      select type(spec=>ions%atom(iatom)%species)
      class is(pseudopotential_t)
        ps => spec%ps
      class default
        assert(.false.)
      end select
 
      if (.not. this%alternative) then
 
        if (states_are_real(st)) then
          safe_allocate(dorbital(1:mesh%np, 1:st%d%dim))
        else
          safe_allocate(zorbital(1:mesh%np, 1:st%d%dim))
        end if
 
        do iorb = 1, this%norbs
          if (iatom /= this%atom(iorb)) cycle
 
          call ions%atom(iatom)%species%get_iwf_ilm(this%level(iorb), 1, ii, ll, mm)
          factor = ps%conf%occ(ii, 1)/(m_two*ll + m_one)
 
          if (states_are_real(st)) then
            call dget_ao(this, st, mesh, ions, iorb, 1, dorbital, use_psi = .true.)
            !$omp parallel do
            do ip = 1, mesh%np
              rho(ip, 1) = rho(ip, 1) + factor*dorbital(ip, 1)**2
            end do
          else
            call zget_ao(this, st, mesh, ions, iorb, 1, zorbital, use_psi = .true.)
            !$omp parallel do
            do ip = 1, mesh%np
              rho(ip, 1) = rho(ip, 1) + factor*abs(zorbital(ip, 1))**2
            end do
          end if
 
        end do
 
        safe_deallocate_a(dorbital)
        safe_deallocate_a(zorbital)
 
      else
 
        ! for simplicity, always use real ones here.
        call dlcao_alt_get_orbital(this, this%sphere(iatom), ions, 1, iatom, this%norb_atom(iatom))
 
        ! the extra orbitals with the derivative are not relevant here, hence we divide by this%mult
        safe_allocate(factors(1:this%norb_atom(iatom)/this%mult))
 
        do iorb = 1, this%norb_atom(iatom)/this%mult
          call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
          factors(iorb) = ps%conf%occ(ii, 1)/(m_two*ll + m_one)
        end do
 
        !$omp parallel do private(ip, aa, iorb) if(.not. this%sphere(iatom)%overlap)
        do ip = 1, this%sphere(iatom)%np
          aa = m_zero
          do iorb = 1, this%norb_atom(iatom)/this%mult
            aa = aa + factors(iorb)*this%orbitals(iatom)%dff_linear(ip, iorb)**2
          end do
          rho(this%sphere(iatom)%map(ip), 1) = rho(this%sphere(iatom)%map(ip), 1) + aa
        end do
 
        safe_deallocate_a(factors)
 
      end if
 
    else
      call species_atom_density(ions%atom(iatom)%species, ions%namespace, ions%space, ions%latt, &
        ions%pos(:, iatom), mesh, spin_channels, rho)
    end if
 
    ! The above code can sometimes return negative values of the density. Here we avoid introducing
    ! them in the calculation of v_s, mixing, ...
    do ispin = 1, spin_channels
      !$omp parallel do simd
      do ip = 1, mesh%np
        rho(ip, ispin) = max(rho(ip, ispin), m_zero)
      end do
    end do
 
    pop_sub(lcao_atom_density)
  end subroutine lcao_atom_density
 
  ! ---------------------------------------------------------
  subroutine lcao_guess_density(this, namespace, st, gr, hm, ions, qtot, ispin, rho)
    type(lcao_t),             intent(inout) :: this
    type(namespace_t),        intent(in)    :: namespace
    type(states_elec_t),      intent(in)    :: st
    type(grid_t),             intent(in)    :: gr
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(ions_t),             intent(in)    :: ions
    real(real64),             intent(in)    :: qtot
    integer,                  intent(in)    :: ispin
    real(real64), contiguous, intent(out)   :: rho(:, :)
 
    integer :: ia, is, idir,  ip, m_dim
    integer(int64), save :: iseed = splitmix64_321
    type(block_t) :: blk
    real(real64) :: rr, rnd, phi, theta, lmag, n1, n2
    real(real64), allocatable :: atom_rho(:,:), mag(:,:)
    real(real64), parameter :: tol_min_mag = 1.0e-20_real64
 
    push_sub(lcao_guess_density)
 
    if (st%d%spin_channels == 1) then
      this%gmd_opt = initrho_paramagnetic
    else
      !%Variable GuessMagnetDensity
      !%Type integer
      !%Default ferromagnetic
      !%Section SCF::LCAO
      !%Description
      !% The guess density for the SCF cycle is just the sum of all the atomic densities.
      !% When performing spin-polarized or non-collinear-spin calculations this option sets
      !% the guess magnetization density.
      !%
      !% For anti-ferromagnetic configurations, the <tt>user_defined</tt> option should be used.
      !%
      !% Note that if the <tt>paramagnetic</tt> option is used, the final ground state will also be
      !% paramagnetic, but the same is not true for the other options.
      !%Option paramagnetic 1
      !% Magnetization density is zero.
      !%Option ferromagnetic 2
      !% Magnetization density is the sum of the atomic magnetization densities.
      !%Option random 3
      !% Each atomic magnetization density is randomly rotated.
      !%Option user_defined 77
      !% The atomic magnetization densities are rotated so that the magnetization
      !% vector has the same direction as a vector provided by the user. In this case,
      !% the <tt>AtomsMagnetDirection</tt> block has to be set.
      !%End
      call parse_variable(namespace, 'GuessMagnetDensity', initrho_ferromagnetic, this%gmd_opt)
      if (.not. varinfo_valid_option('GuessMagnetDensity', this%gmd_opt)) call messages_input_error(namespace, 'GuessMagnetDensity')
      call messages_print_var_option('GuessMagnetDensity', this%gmd_opt, namespace=namespace)
    end if
 
    if (parse_is_defined(namespace, 'GuessMagnetDensity') .and. (hm%theory_level == hartree_fock &
      .or. hm%theory_level == generalized_kohn_sham_dft) .and. st%d%spin_channels > 1) then
      message(1) = "GuessMagnetDensity cannot be used for Hartree-Fock and generalized Kohn-Sham calculation."
      message(2) = "Please perform a LDA or GGA calculation first and restart from this calculation."
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (this%gmd_opt == initrho_userdef) then
      !%Variable AtomsMagnetDirection
      !%Type block
      !%Section SCF::LCAO
      !%Description
      !% This option is only used when <tt>GuessMagnetDensity</tt> is
      !% set to <tt>user_defined</tt>. It provides a direction for the
      !% magnetization vector of each atom when building the guess
      !% density. In order to do that, the user should specify the
      !% coordinates of a vector that has the desired direction and
      !% norm.  Note that it is necessary to maintain the ordering in
      !% which the species were defined in the coordinates
      !% specifications.
      !%
      !% For spin-polarized calculations, the vectors should have only
      !% one component; for non-collinear-spin calculations, they
      !% should have three components. If the norm of the vector is greater
      !% than the number of valence electrons in the atom, it will be rescaled
      !% to this number, which is the maximum possible magnetization.
      !%End
      if (parse_block(namespace, 'AtomsMagnetDirection', blk) < 0) then
        message(1) = "AtomsMagnetDirection block is not defined."
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (parse_block_n(blk) /= ions%natoms) then
        message(1) = "The number of rows in the AtomsMagnetDirection block does not equal the number of atoms."
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (ispin == spin_polarized) then
        m_dim = 1
      elseif(ispin == spinors) then
        m_dim = 3
      endif
 
      safe_allocate(mag(1:m_dim, 1:ions%natoms))
      do ia = 1, ions%natoms
        !Read from AtomsMagnetDirection block
        do idir = 1, m_dim
          call parse_block_float(blk, ia-1, idir-1, mag(idir, ia))
          if (abs(mag(idir, ia)) < tol_min_mag) mag(idir, ia) = m_zero
        end do
      end do
      call parse_block_end(blk)
    end if
 
    rho = m_zero
 
    safe_allocate(atom_rho(1:gr%np, 1:st%d%spin_channels))
    select case (this%gmd_opt)
    case (initrho_paramagnetic)
      do ia = ions%atoms_dist%start, ions%atoms_dist%end
        call lcao_atom_density(this, st, gr, ions, ia, st%d%spin_channels, atom_rho)
        call lalg_axpy(gr%np, st%d%spin_channels, m_one, atom_rho, rho)
      end do
 
      if (st%d%spin_channels == 2) then
        !$omp parallel do
        do ip = 1, gr%np
          rho(ip, 1) = m_half*(rho(ip, 1) + rho(ip, 2))
          rho(ip, 2) = rho(ip, 1)
        end do
      end if
 
    case (initrho_ferromagnetic)
      do ia = ions%atoms_dist%start, ions%atoms_dist%end
        call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho(1:gr%np, 1:2))
        rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
      end do
 
    case (initrho_random) ! Randomly oriented spins
      do ia = ions%atoms_dist%start, ions%atoms_dist%end
        call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho)
 
        ! For charge neutral atoms, the magnetization density will always be zero
        ! In order to still make a random spin structure, we then reassign the charge in
        ! one spin channel to have a 1 \nu_B for this atom, with a random direction
        !
        ! An example where this is needed is a Xe3 cluster where we put an excess charge
        ! Each atom is neutral but the full system has a net magnetiation
        n1 = dmf_integrate(gr, atom_rho(:, 1))
        n2 = dmf_integrate(gr, atom_rho(:, 2))
 
        if (is_close(n1, n2)) then
          lmag = m_one
 
          call lalg_axpy(gr%np, (lmag - n1 + n2)/m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
          call lalg_scal(gr%np, (n1 + n2 - lmag)/m_two/n2, atom_rho(:, 2))
        end if
 
 
        if (ispin == spin_polarized) then
          call quickrnd(iseed, rnd)
          rnd = rnd - m_half
          if (rnd > m_zero) then
            rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
          else
            rho(1:gr%np, 1) = rho(1:gr%np, 1) + atom_rho(1:gr%np, 2)
            rho(1:gr%np, 2) = rho(1:gr%np, 2) + atom_rho(1:gr%np, 1)
          end if
        elseif (ispin == spinors) then
          call quickrnd(iseed, phi)
          call quickrnd(iseed, theta)
          phi = phi*m_two*m_pi
          theta = theta*m_pi*m_half
 
          call accumulate_rotated_density(gr, rho, atom_rho, theta, phi)
        end if
      end do
 
    case (initrho_userdef) ! User-defined
      do ia = ions%atoms_dist%start, ions%atoms_dist%end
        !Get atomic density
        call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho)
 
        !Scale magnetization density
        n1 = dmf_integrate(gr, atom_rho(:, 1))
        n2 = dmf_integrate(gr, atom_rho(:, 2))
 
        lmag = norm2(mag(:, ia))
        if (lmag > n1 + n2) then
          mag = mag*(n1 + n2)/lmag
          lmag = n1 + n2
        elseif (abs(lmag) <= m_epsilon) then
          if (abs(n1 - n2) <= m_epsilon) then
            call lalg_axpy(gr%np, 2, m_one, atom_rho, rho)
          else
            !$omp parallel do simd
            do ip = 1, gr%np
              atom_rho(ip, 1) = m_half*(atom_rho(ip, 1) + atom_rho(ip, 2))
              rho(ip, 1) = rho(ip, 1) + atom_rho(ip, 1)
              rho(ip, 2) = rho(ip, 2) + atom_rho(ip, 1)
            end do
          end if
          cycle
        end if
 
        if (.not. is_close(n1 - n2, lmag) .and. abs(n2) > m_epsilon) then
          if (n1 - n2 < lmag) then
            call lalg_axpy(gr%np, (lmag - n1 + n2)/m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
            call lalg_scal(gr%np, (n1 + n2 - lmag)/m_two/n2, atom_rho(:, 2))
          elseif (n1 - n2 > lmag) then
            call lalg_axpy(gr%np, (n1 - n2 - lmag)/m_two/n1, atom_rho(:, 1), atom_rho(:, 2))
            call lalg_scal(gr%np, (n1 + n2 + lmag)/m_two/n1, atom_rho(:, 1))
          end if
        end if
 
        !Rotate magnetization density
        if (ispin == spin_polarized) then
 
          if (mag(1, ia) > m_zero) then
            call lalg_axpy(gr%np, 2, m_one, atom_rho, rho)
          else
            call lalg_axpy(gr%np, m_one, atom_rho(:,2), rho(:,1))
            call lalg_axpy(gr%np, m_one, atom_rho(:,1), rho(:,2))
          end if
 
        elseif (ispin == spinors) then
 
          call get_angles_from_magnetization(mag(:, ia), lmag, theta, phi)
          call accumulate_rotated_density(gr, rho, atom_rho, theta, phi)
        end if
      end do
 
    end select
 
 
    if (ions%atoms_dist%parallel) then
      do is = 1, st%d%nspin
        call lalg_copy(gr%np, rho(:,is), atom_rho(:,1))
        call ions%atoms_dist%mpi_grp%allreduce(atom_rho(1, 1), rho(1, is), gr%np, mpi_double_precision, mpi_sum)
      end do
    end if
 
    ! we now renormalize the density (necessary if we have a charged system)
    rr = integrated_charge_density(gr, st, rho)
    write(message(1),'(a,f13.6)')'Info: Unnormalized total charge = ', rr
    call messages_info(1, namespace=namespace)
 
    ! We only renormalize if the density is not zero
    if (abs(rr) > m_epsilon) then
      call lalg_scal(gr%np, st%d%nspin, qtot/rr, rho)
      rr = integrated_charge_density(gr, st, rho)
      write(message(1),'(a,f13.6)')'Info: Renormalized total charge = ', rr
      call messages_info(1, namespace=namespace)
    end if
 
    ! Symmetrize the density if needed
    if (st%symmetrize_density) then
      do is = 1, st%d%nspin
        call dgrid_symmetrize_scalar_field(gr, rho(:, is))
      end do
    end if
 
    safe_deallocate_a(atom_rho)
    safe_deallocate_a(mag)
 
    pop_sub(lcao_guess_density)
  end subroutine lcao_guess_density
 
  real(real64) function integrated_charge_density(gr, st, rho) result(rr)
    type(grid_t),             intent(in) :: gr
    type(states_elec_t),      intent(in) :: st
    real(real64),             intent(in) :: rho(:,:)
 
    integer :: is
 
    rr = m_zero
    do is = 1, st%d%spin_channels
      rr = rr + dmf_integrate(gr, rho(:, is), reduce = .false.)
    end do
 
    call gr%allreduce(rr)
  end function integrated_charge_density
 
  ! ---------------------------------------------------------
  subroutine accumulate_rotated_density(mesh, rho, atom_rho, theta, phi)
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(inout) :: rho(:,:)
    real(real64),        intent(in)    :: atom_rho(:,:)
    real(real64),        intent(in)    :: theta, phi
 
    integer :: ip
    real(real64) :: half_theta
 
    push_sub(accumulate_rotated_density)
 
    half_theta = m_half * theta
 
    !$omp parallel do simd
    do ip = 1, mesh%np
      rho(ip, 1) = rho(ip, 1) + cos(half_theta)**2*atom_rho(ip, 1) + sin(half_theta)**2*atom_rho(ip, 2)
      rho(ip, 2) = rho(ip, 2) + sin(half_theta)**2*atom_rho(ip, 1) + cos(half_theta)**2*atom_rho(ip, 2)
      rho(ip, 3) = rho(ip, 3) + cos(half_theta)*sin(half_theta)*cos(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
      rho(ip, 4) = rho(ip, 4) - cos(half_theta)*sin(half_theta)*sin(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
    end do
 
    pop_sub(accumulate_rotated_density)
  end subroutine accumulate_rotated_density
 
  ! ---------------------------------------------------------
 
  subroutine lcao_init_orbitals(this, namespace, st, gr, ions, start)
    type(lcao_t),        intent(inout) :: this
    type(namespace_t),   intent(in)    :: namespace
    type(states_elec_t), intent(inout) :: st
    type(grid_t),        intent(in)    :: gr
    type(ions_t),        intent(in)    :: ions
    integer, optional,   intent(in)    :: start
 
    if (.not. lcao_is_available(this)) return
 
    push_sub(lcao_init_orbitals)
 
    if (.not. this%alternative) then
      if (states_are_real(st)) then
        call dinit_orbitals(this, namespace, st, gr, ions, start)
      else
        call zinit_orbitals(this, namespace, st, gr, ions, start)
      end if
    else
      if (states_are_real(st)) then
        call dlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
      else
        call zlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
      end if
 
    end if
 
    pop_sub(lcao_init_orbitals)
  end subroutine lcao_init_orbitals
 
  subroutine get_angles_from_magnetization(mag, lmag, theta, phi)
    real(real64), intent(in)  :: mag(3)
    real(real64), intent(in)  :: lmag
    real(real64), intent(out) :: theta
    real(real64), intent(out) :: phi
 
    assert(lmag > m_zero)
 
    theta = acos(max(-1.0_real64, min(1.0_real64, mag(3)/lmag)))
 
    if (abs(mag(1)) <= m_epsilon) then
      if (abs(mag(2)) <= m_epsilon) then
        phi = m_zero
      elseif (mag(2) < m_zero) then
        phi = m_pi*m_three*m_half
      else
        phi = m_pi*m_half
      end if
    else
 
      phi = atan2(mag(2), mag(1))
      !note: atan2 returns the correct phase. The correction from wikipedia on spherical coordinates is not necessary.
 
    end if
 
  end subroutine get_angles_from_magnetization
 
 
  subroutine rotate_random_states_to_local_frame(st, gr, ist_start, gmd_opt)
    type(states_elec_t), intent(inout) :: st
    type(grid_t),        intent(in)    :: gr
    integer,             intent(in)    :: ist_start
    integer,             intent(in)    :: gmd_opt
 
    integer :: ik, ist, ip
    real(real64), allocatable :: md(:,:), up(:)
    complex(real64), allocatable :: down(:), psi(:,:)
    real(real64) :: lmag, theta, phi, norm, mm(3)
 
    if (st%d%ispin /= spinors) return
    if (gmd_opt == initrho_random) return
 
    push_sub(rotate_random_states_to_local_frame)
 
    safe_allocate(psi(1:gr%np, 1:st%d%dim))
 
    select case (gmd_opt)
    case(initrho_paramagnetic)
      ! Random spinors should have no net magnetization
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = ist_start, st%st_end
          call states_elec_get_state(st, gr, ist, ik, psi)
          !$omp parallel do private(norm)
          do ip = 1, gr%np
            psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
            psi(ip, 2) = psi(ip, 1)
          end do
          call zmf_normalize(gr, st%d%dim, psi)
          call states_elec_set_state(st, gr, ist, ik, psi)
        end do
      end do
 
    case(initrho_ferromagnetic)
      ! Random spinors magnetization is aligned along the z axis
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = ist_start, st%st_end
          call states_elec_get_state(st, gr, ist, ik, psi)
          !$omp parallel do
          do ip = 1, gr%np
            psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
            psi(ip, 2) = m_zero
          end do
          call zmf_normalize(gr, st%d%dim, psi)
          call states_elec_set_state(st, gr, ist, ik, psi)
        end do
      end do
 
    case(initrho_userdef)
      ! Random spinors magnetization is aligned along the local frame
      safe_allocate(md(1:gr%np, 1:3))
      safe_allocate(up(1:gr%np))
      safe_allocate(down(1:gr%np))
      call magnetic_density(gr, st%d, st%rho, md)
 
      mm(1) = dmf_integrate(gr, md(:, 1), reduce = .false.)
      mm(2) = dmf_integrate(gr, md(:, 2), reduce = .false.)
      mm(3) = dmf_integrate(gr, md(:, 3), reduce = .false.)
 
      if (gr%parallel_in_domains) then
        call gr%allreduce(mm)
      end if
 
      ! If the density has a net direction, we align globally the spinors to this direction
      ! This allows to generate random pure states
      lmag = norm2(mm)
      if (lmag > 1.0e-2_real64) then
        call get_angles_from_magnetization(mm, lmag, theta, phi)
        up(:) = cos(theta/m_two)
        down(:) = exp(-m_zi * phi) * sin(theta/m_two)
      else
        ! We do a rotation in the local frame
        ! Note that this produces mixed states, so the norm of the magnetization is smaller than 1
        !
        !$omp parallel do private(lmag, theta, phi)
        do ip = 1, gr%np
          lmag = norm2(md(ip, 1:3))
          call get_angles_from_magnetization(md(ip,:), lmag, theta, phi)
          up(ip) = cos(theta/m_two)
          down(ip) = exp(-m_zi * phi) * sin(theta/m_two)
        end do
        safe_deallocate_a(md)
      end if
 
      ! Rotate the spinors
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = ist_start, st%st_end
          call states_elec_get_state(st, gr, ist, ik, psi)
          call zmf_normalize(gr, st%d%dim, psi)
          !$omp parallel do private(norm)
          do ip = 1, gr%np
            norm = sqrt(real(psi(ip, 1)*conjg(psi(ip, 1)) + psi(ip, 2)*conjg(psi(ip, 2))))
            psi(ip, 1) = norm * up(ip)
            psi(ip, 2) = norm * down(ip)
          end do
          call states_elec_set_state(st, gr, ist, ik, psi)
        end do
      end do
      safe_deallocate_a(up)
      safe_deallocate_a(down)
 
    end select
    safe_deallocate_a(psi)
 
    pop_sub(rotate_random_states_to_local_frame)
  end subroutine rotate_random_states_to_local_frame
 
#include "undef.F90"
#include "real.F90"
#include "lcao_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "lcao_inc.F90"
 
 
end module lcao_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: