qs_mixing_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
MODULE qs_mixing_utils
   USE bibliography,                    ONLY: sundararaman2017,&
                                              cite_reference
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_dbcsr_api,                    ONLY: dbcsr_create,&
                                              dbcsr_p_type,&
                                              dbcsr_set,&
                                              dbcsr_type,&
                                              dbcsr_type_symmetric
   USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
   USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set
   USE distribution_1d_types,           ONLY: distribution_1d_type
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: z_zero
   USE message_passing,                 ONLY: mp_para_env_type
   USE pw_types,                        ONLY: pw_c1d_gs_type
   USE qs_density_mixing_types,         ONLY: broyden_mixing_nr,&
                                              cp_1d_z_p_type,&
                                              gspace_mixing_nr,&
                                              mixing_storage_type,&
                                              modified_broyden_mixing_nr,&
                                              multisecant_mixing_nr,&
                                              new_pulay_mixing_nr,&
                                              pulay_mixing_nr
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
   USE qs_rho_atom_types,               ONLY: rho_atom_type
   USE qs_rho_types,                    ONLY: qs_rho_get,&
                                              qs_rho_type
   USE qs_scf_methods,                  ONLY: cp_sm_mix
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils'
   INTEGER, PARAMETER, PRIVATE          :: max_dftb_scc_vars = 1, &
                                           max_xtb_scc_vars = 28
 
   PUBLIC :: mixing_allocate, mixing_init, mixing_density, charge_mixing_init, &
             self_consistency_check
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param rho_ao ...
!> \param p_delta ...
!> \param para_env ...
!> \param p_out ...
!> \param delta ...
! **************************************************************************************************
   SUBROUTINE self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao, p_delta
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: p_out
      REAL(kind=dp), INTENT(INOUT)                       :: delta
 
      CHARACTER(len=*), PARAMETER :: routinen = 'self_consistency_check'
 
      INTEGER                                            :: handle, ic, ispin, nimg, nspins
      REAL(kind=dp)                                      :: tmp
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_q, p_in
 
      CALL timeset(routinen, handle)
 
      NULLIFY (matrix_q, p_in)
 
      cpassert(ASSOCIATED(p_out))
      NULLIFY (matrix_q, p_in)
      p_in => rho_ao
      matrix_q => p_delta
      nspins = SIZE(p_in, 1)
      nimg = SIZE(p_in, 2)
 
      ! Compute the difference (p_out - p_in)and check convergence
      delta = 0.0_dp
      DO ispin = 1, nspins
         DO ic = 1, nimg
            CALL dbcsr_set(matrix_q(ispin, ic)%matrix, 0.0_dp)
            CALL cp_sm_mix(m1=p_out(ispin, ic)%matrix, m2=p_in(ispin, ic)%matrix, &
                           p_mix=1.0_dp, delta=tmp, para_env=para_env, &
                           m3=matrix_q(ispin, ic)%matrix)
            delta = max(tmp, delta)
         END DO
      END DO
      CALL timestop(handle)
 
   END SUBROUTINE self_consistency_check
 
! **************************************************************************************************
!> \brief  allocation needed when density mixing is used
!> \param qs_env ...
!> \param mixing_method ...
!> \param p_mix_new ...
!> \param p_delta ...
!> \param nspins ...
!> \param mixing_store ...
!> \par History
!>      05.2009 created [MI]
!>      08.2014 kpoints [JGH]
!>      02.2015 changed input to make g-space mixing available in linear scaling SCF [Patrick Seewald]
!>      02.2019 charge mixing [JGH]
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      INTEGER                                            :: mixing_method
      TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
         POINTER                                         :: p_mix_new, p_delta
      INTEGER, INTENT(IN)                                :: nspins
      TYPE(mixing_storage_type), POINTER                 :: mixing_store
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'mixing_allocate'
 
      INTEGER                                            :: handle, i, ia, iat, ic, ikind, ispin, &
                                                            max_shell, na, natom, nbuffer, nel, &
                                                            nimg, nkind
      LOGICAL                                            :: charge_mixing
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
      TYPE(dbcsr_type), POINTER                          :: refmatrix
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(distribution_1d_type), POINTER                :: distribution_1d
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_orb
      TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom
 
      CALL timeset(routinen, handle)
 
      NULLIFY (matrix_s, dft_control, sab_orb, refmatrix, rho_atom)
      CALL get_qs_env(qs_env, &
                      sab_orb=sab_orb, &
                      matrix_s_kp=matrix_s, &
                      dft_control=dft_control)
 
      charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
                      .OR. dft_control%qs_control%semi_empirical
      refmatrix => matrix_s(1, 1)%matrix
      nimg = dft_control%nimages
 
      !   *** allocate p_mix_new ***
      IF (PRESENT(p_mix_new)) THEN
         IF (.NOT. ASSOCIATED(p_mix_new)) THEN
            CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
            DO i = 1, nspins
               DO ic = 1, nimg
                  ALLOCATE (p_mix_new(i, ic)%matrix)
                  CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
                                    name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
                  CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
                  CALL dbcsr_set(p_mix_new(i, ic)%matrix, 0.0_dp)
               END DO
            END DO
         END IF
      END IF
 
      !   *** allocate p_delta ***
      IF (PRESENT(p_delta)) THEN
         IF (mixing_method >= gspace_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(p_delta)) THEN
               CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
               DO i = 1, nspins
                  DO ic = 1, nimg
                     ALLOCATE (p_delta(i, ic)%matrix)
                     CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
                                       name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
                     CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
                     CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
                  END DO
               END DO
            END IF
            cpassert(ASSOCIATED(mixing_store))
         END IF
      END IF
 
      IF (charge_mixing) THEN
         !   *** allocate buffer for TB SCC variable mixing ***
         IF (mixing_method >= gspace_mixing_nr) THEN
            cpassert(.NOT. mixing_store%gmix_p)
            IF (dft_control%qs_control%dftb) THEN
               max_shell = max_dftb_scc_vars
            ELSE IF (dft_control%qs_control%xtb) THEN
               max_shell = max_xtb_scc_vars
            ELSE
               cpabort('UNKNOWN METHOD')
            END IF
            nbuffer = mixing_store%nbuffer
            mixing_store%ncall = 0
            CALL get_qs_env(qs_env, local_particles=distribution_1d)
            nkind = SIZE(distribution_1d%n_el)
            na = sum(distribution_1d%n_el(:))
            IF (ASSOCIATED(mixing_store%atlist)) DEALLOCATE (mixing_store%atlist)
            ALLOCATE (mixing_store%atlist(na))
            mixing_store%nat_local = na
            mixing_store%max_shell = max_shell
            ia = 0
            DO ikind = 1, nkind
               nel = distribution_1d%n_el(ikind)
               DO iat = 1, nel
                  ia = ia + 1
                  mixing_store%atlist(ia) = distribution_1d%list(ikind)%array(iat)
               END DO
            END DO
            IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
            ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
            IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
            ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
         END IF
         IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
            IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
            ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
            mixing_store%pulay_matrix = 0.0_dp
         ELSE IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
            IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
            ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
            mixing_store%abroy = 0.0_dp
            IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
            ALLOCATE (mixing_store%wbroy(nbuffer))
            mixing_store%wbroy = 0.0_dp
            IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
            ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
            mixing_store%dfbroy = 0.0_dp
            IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
            ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
            mixing_store%ubroy = 0.0_dp
         ELSE IF (mixing_method == multisecant_mixing_nr) THEN
            cpabort("multisecant_mixing not available")
         END IF
      ELSE
         !   *** allocate buffer for gspace mixing ***
         IF (mixing_method >= gspace_mixing_nr) THEN
            nbuffer = mixing_store%nbuffer
            mixing_store%ncall = 0
            mixing_store%ncall_p = 0
            IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
               ALLOCATE (mixing_store%rhoin(nspins))
               DO ispin = 1, nspins
                  NULLIFY (mixing_store%rhoin(ispin)%cc)
               END DO
            END IF
 
            IF (mixing_store%gmix_p .AND. dft_control%qs_control%gapw) THEN
               CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
               natom = SIZE(rho_atom)
               IF (.NOT. ASSOCIATED(mixing_store%paw)) THEN
                  ALLOCATE (mixing_store%paw(natom))
                  mixing_store%paw = .false.
                  ALLOCATE (mixing_store%cpc_h_in(natom, nspins))
                  ALLOCATE (mixing_store%cpc_s_in(natom, nspins))
                  DO ispin = 1, nspins
                     DO iat = 1, natom
                        NULLIFY (mixing_store%cpc_h_in(iat, ispin)%r_coef)
                        NULLIFY (mixing_store%cpc_s_in(iat, ispin)%r_coef)
                     END DO
                  END DO
               END IF
            END IF
         END IF
 
         !   *** allocare res_buffer if needed
         IF (mixing_method >= pulay_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
               ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
               DO ispin = 1, nspins
                  DO i = 1, nbuffer
                     NULLIFY (mixing_store%res_buffer(i, ispin)%cc)
                  END DO
               END DO
            END IF
         END IF
 
         !   *** allocate pulay
         IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
               ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
            END IF
 
            IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
               ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
               DO ispin = 1, nspins
                  DO i = 1, nbuffer
                     NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
                  END DO
               END DO
            END IF
            IF (mixing_store%gmix_p) THEN
               IF (dft_control%qs_control%gapw) THEN
                  IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer)) THEN
                     ALLOCATE (mixing_store%cpc_h_in_buffer(nbuffer, natom, nspins))
                     ALLOCATE (mixing_store%cpc_s_in_buffer(nbuffer, natom, nspins))
                     ALLOCATE (mixing_store%cpc_h_res_buffer(nbuffer, natom, nspins))
                     ALLOCATE (mixing_store%cpc_s_res_buffer(nbuffer, natom, nspins))
                     DO ispin = 1, nspins
                        DO iat = 1, natom
                           DO i = 1, nbuffer
                              NULLIFY (mixing_store%cpc_h_in_buffer(i, iat, ispin)%r_coef)
                              NULLIFY (mixing_store%cpc_s_in_buffer(i, iat, ispin)%r_coef)
                              NULLIFY (mixing_store%cpc_h_res_buffer(i, iat, ispin)%r_coef)
                              NULLIFY (mixing_store%cpc_s_res_buffer(i, iat, ispin)%r_coef)
                           END DO
                        END DO
                     END DO
                  END IF
               END IF
            END IF
 
         END IF
         !   *** allocate broyden buffer ***
         IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
               ALLOCATE (mixing_store%rhoin_old(nspins))
               DO ispin = 1, nspins
                  NULLIFY (mixing_store%rhoin_old(ispin)%cc)
               END DO
            END IF
            IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
               ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
               ALLOCATE (mixing_store%last_res(nspins))
               DO ispin = 1, nspins
                  DO i = 1, nbuffer
                     NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
                  END DO
                  NULLIFY (mixing_store%last_res(ispin)%cc)
               END DO
            END IF
            IF (mixing_method == modified_broyden_mixing_nr) THEN
               IF (.NOT. ASSOCIATED(mixing_store%weight)) THEN
                  ALLOCATE (mixing_store%weight(nbuffer, nspins))
               END IF
               mixing_store%weight = 0.0_dp
            END IF
            IF (mixing_store%gmix_p) THEN
 
               IF (dft_control%qs_control%gapw) THEN
                  IF (.NOT. ASSOCIATED(mixing_store%cpc_h_old)) THEN
                     ALLOCATE (mixing_store%cpc_h_old(natom, nspins))
                     ALLOCATE (mixing_store%cpc_s_old(natom, nspins))
                     DO ispin = 1, nspins
                        DO iat = 1, natom
                           NULLIFY (mixing_store%cpc_h_old(iat, ispin)%r_coef)
                           NULLIFY (mixing_store%cpc_s_old(iat, ispin)%r_coef)
                        END DO
                     END DO
                  END IF
                  IF (.NOT. ASSOCIATED(mixing_store%dcpc_h_in)) THEN
                     ALLOCATE (mixing_store%dcpc_h_in(nbuffer, natom, nspins))
                     ALLOCATE (mixing_store%dcpc_s_in(nbuffer, natom, nspins))
                     ALLOCATE (mixing_store%cpc_h_lastres(natom, nspins))
                     ALLOCATE (mixing_store%cpc_s_lastres(natom, nspins))
                     DO ispin = 1, nspins
                        DO iat = 1, natom
                           DO i = 1, nbuffer
                              NULLIFY (mixing_store%dcpc_h_in(i, iat, ispin)%r_coef)
                              NULLIFY (mixing_store%dcpc_s_in(i, iat, ispin)%r_coef)
                           END DO
                           NULLIFY (mixing_store%cpc_h_lastres(iat, ispin)%r_coef)
                           NULLIFY (mixing_store%cpc_s_lastres(iat, ispin)%r_coef)
                        END DO
                     END DO
                  END IF
               END IF
            END IF
         END IF
 
         !   *** allocate multisecant buffer ***
         IF (mixing_method == multisecant_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(mixing_store%norm_res_buffer)) THEN
               ALLOCATE (mixing_store%norm_res_buffer(nbuffer, nspins))
            END IF
         END IF
 
         IF (mixing_method == multisecant_mixing_nr) THEN
            IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
               ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
               DO ispin = 1, nspins
                  DO i = 1, nbuffer
                     NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
                  END DO
               END DO
            END IF
         END IF
 
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE mixing_allocate
 
! **************************************************************************************************
!> \brief Allocate or resize a mixing vector, optionally clearing its contents.
!> \param buffer Vector or array of vectors.
!> \param n Local vector length, including any auxiliary coordinates.
!> \param clear Whether to clear the vector; Pulay retains history across SCF runs.
! **************************************************************************************************
   ELEMENTAL SUBROUTINE mixing_buffer_allocate(buffer, n, clear)
      TYPE(cp_1d_z_p_type), INTENT(INOUT)                :: buffer
      INTEGER, INTENT(IN)                                :: n
      LOGICAL, INTENT(IN)                                :: clear
 
      IF (ASSOCIATED(buffer%cc)) THEN
         IF (SIZE(buffer%cc) /= n) DEALLOCATE (buffer%cc)
      END IF
      IF (.NOT. ASSOCIATED(buffer%cc)) ALLOCATE (buffer%cc(n))
      IF (clear) buffer%cc = z_zero
   END SUBROUTINE mixing_buffer_allocate
 
! **************************************************************************************************
!> \brief View the physical density, appending auxiliary coordinates on their owning rank.
!> \param mixing_store Owns the optional enlarged work buffer.
!> \param rho_g Physical plane-wave density; its storage and grid are unchanged.
!> \param density Vectors passed to the existing mixing algorithms.
!> \param auxiliary Replicated density coordinates in electrons, scaled to Fourier-density units.
! **************************************************************************************************
   SUBROUTINE mixing_density(mixing_store, rho_g, density, auxiliary)
      TYPE(mixing_storage_type), POINTER                 :: mixing_store
      TYPE(pw_c1d_gs_type), INTENT(IN)                   :: rho_g(:)
      TYPE(cp_1d_z_p_type), ALLOCATABLE, INTENT(OUT)     :: density(:)
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: auxiliary(:, :)
 
      INTEGER                                            :: ispin, ng
 
      IF (mixing_store%naux > 0) THEN
         cpassert(PRESENT(auxiliary))
         cpassert(SIZE(auxiliary, 1) == mixing_store%naux)
         cpassert(SIZE(auxiliary, 2) == SIZE(rho_g))
      END IF
      ng = SIZE(rho_g(1)%array)
      ALLOCATE (density(SIZE(rho_g)))
      DO ispin = 1, SIZE(rho_g)
         density(ispin)%cc => rho_g(ispin)%array
         IF (ALLOCATED(mixing_store%aux_buffer)) THEN
            density(ispin)%cc => mixing_store%aux_buffer(:, ispin)
            density(ispin)%cc(:ng) = rho_g(ispin)%array
            density(ispin)%cc(ng + 1:) = auxiliary(:, ispin)/rho_g(ispin)%pw_grid%vol
         END IF
      END DO
   END SUBROUTINE mixing_density
 
! **************************************************************************************************
!> \brief  initialiation needed when gspace mixing is used
!> \param mixing_method ...
!> \param rho ...
!> \param mixing_store ...
!> \param para_env ...
!> \param rho_atom ...
!> \param auxiliary Replicated real density coordinates, included once in the mixing metric.
!> \par History
!>      05.2009 created [MI]
!> \author MI
! **************************************************************************************************
   SUBROUTINE mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom, auxiliary)
      INTEGER, INTENT(IN)                                :: mixing_method
      TYPE(qs_rho_type), POINTER                         :: rho
      TYPE(mixing_storage_type), POINTER                 :: mixing_store
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
         POINTER                                         :: rho_atom
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: auxiliary(:, :)
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'mixing_init'
 
      INTEGER :: handle, iat, ib, icoef1, icoef2, ig, ig1, ig_count, iproc, ispin, n1, n2, natom, &
         nbuffer, ncoef_h1, ncoef_h2, ncoef_s1, ncoef_s2, ng, nimg, nspin, nvec
      REAL(dp)                                           :: bconst, beta, cpc_h_tmp, cpc_s_tmp, &
                                                            fdamp, g2max, g2min, kmin, qk, qkappa, &
                                                            qm
      REAL(dp), DIMENSION(:), POINTER                    :: g2
      REAL(dp), DIMENSION(:, :), POINTER                 :: g_vec
      TYPE(cp_1d_z_p_type), ALLOCATABLE                  :: density(:)
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
      TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
 
      CALL timeset(routinen, handle)
 
      NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
      CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
 
      nspin = SIZE(rho_g)
      ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
      nimg = SIZE(rho_ao_kp, 2)
      IF (PRESENT(auxiliary) .OR. mixing_store%naux > 0) THEN
         mixing_store%ncall = 0
         mixing_store%ncall_p = 0
      END IF
      mixing_store%naux = 0
      IF (PRESENT(auxiliary)) mixing_store%naux = SIZE(auxiliary, 1)
      nvec = ng
      IF (ALLOCATED(mixing_store%aux_buffer)) DEALLOCATE (mixing_store%aux_buffer)
      IF (mixing_store%naux > 0) THEN
         IF (para_env%mepos == 0) THEN
            nvec = ng + mixing_store%naux
            ALLOCATE (mixing_store%aux_buffer(nvec, nspin))
         END IF
      END IF
      CALL mixing_density(mixing_store, rho_g, density, auxiliary)
      mixing_store%ig_max = ng
      g2 => rho_g(1)%pw_grid%gsq
      g_vec => rho_g(1)%pw_grid%g
 
      IF (mixing_store%max_gvec_exp > 0._dp) THEN
         DO ig = 1, ng
            IF (g2(ig) > mixing_store%max_g2) THEN
               mixing_store%ig_max = ig
               EXIT
            END IF
         END DO
      END IF
 
      IF (ASSOCIATED(mixing_store%kerker_factor)) DEALLOCATE (mixing_store%kerker_factor)
      ALLOCATE (mixing_store%kerker_factor(nvec))
      IF (ASSOCIATED(mixing_store%kerker_factor_mag)) DEALLOCATE (mixing_store%kerker_factor_mag)
      ALLOCATE (mixing_store%kerker_factor_mag(nvec))
      IF (ASSOCIATED(mixing_store%special_metric)) DEALLOCATE (mixing_store%special_metric)
      ALLOCATE (mixing_store%special_metric(nvec))
      beta = mixing_store%beta
      kmin = 0.1_dp
      mixing_store%kerker_factor = 1.0_dp
      mixing_store%kerker_factor_mag = 1.0_dp
      mixing_store%special_metric = 1.0_dp
      ig1 = 1
 
      IF (.NOT. (mixing_method == new_pulay_mixing_nr)) THEN
         IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
         DO ig = ig1, mixing_store%ig_max
            mixing_store%kerker_factor(ig) = max(g2(ig)/(g2(ig) + beta*beta), kmin)
            IF (mixing_store%beta_mag > 1.0e-10_dp) THEN
               mixing_store%kerker_factor_mag(ig) = max(g2(ig)/(g2(ig) + mixing_store%beta_mag*mixing_store%beta_mag), kmin)
            ELSE
               ! beta_mag ~ 0 means no Kerker screening on the magnetization channel
               mixing_store%kerker_factor_mag(ig) = 1.0_dp
            END IF
            mixing_store%special_metric(ig) = &
               1.0_dp + 50.0_dp/8.0_dp*( &
               1.0_dp + cos(g_vec(1, ig)) + cos(g_vec(2, ig)) + cos(g_vec(3, ig)) + &
               cos(g_vec(1, ig))*cos(g_vec(2, ig)) + &
               cos(g_vec(2, ig))*cos(g_vec(3, ig)) + &
               cos(g_vec(1, ig))*cos(g_vec(3, ig)) + &
               cos(g_vec(1, ig))*cos(g_vec(2, ig))*cos(g_vec(3, ig)))
         END DO
      ELSE
         CALL cite_reference(sundararaman2017)
         qkappa = mixing_store%qkappa
         qk = mixing_store%qk
         qm = mixing_store%qm
         IF (rho_g(1)%pw_grid%have_g0) ig1 = 1
         DO ig = ig1, mixing_store%ig_max
            mixing_store%kerker_factor(ig) = (g2(ig) + qkappa*qkappa)/(g2(ig) + qkappa*qkappa + qk*qk)
            mixing_store%special_metric(ig) = (g2(ig) + qkappa*qkappa + qm*qm)/(g2(ig) + qkappa*qkappa)
         END DO
      END IF
      nbuffer = mixing_store%nbuffer
      CALL mixing_buffer_allocate(mixing_store%rhoin, nvec, .false.)
      IF (ASSOCIATED(mixing_store%rhoin_buffer)) CALL mixing_buffer_allocate(mixing_store%rhoin_buffer, nvec, .false.)
      IF (ASSOCIATED(mixing_store%res_buffer)) CALL mixing_buffer_allocate(mixing_store%res_buffer, nvec, .false.)
      DO ispin = 1, nspin
         mixing_store%rhoin(ispin)%cc = density(ispin)%cc
      END DO
      IF (nspin == 2) THEN
         mixing_store%rhoin(1)%cc = density(1)%cc + density(2)%cc
         mixing_store%rhoin(2)%cc = density(1)%cc - density(2)%cc
      END IF
      IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
         DO ispin = 1, nspin
            mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
         END DO
      END IF
 
      IF (mixing_store%gmix_p) THEN
         IF (PRESENT(rho_atom)) THEN
            natom = SIZE(rho_atom)
            DO ispin = 1, nspin
               DO iat = 1, natom
                  IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
                     mixing_store%paw(iat) = .true.
                     n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
                     n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
                     IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
                        IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
                           ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
                        END IF
                        mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
                        mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
                     END IF
                  END IF
               END DO
            END DO
            IF (nspin == 2 .AND. ASSOCIATED(mixing_store%cpc_s_in)) THEN
               DO iat = 1, natom
                  IF (mixing_store%paw(iat)) THEN
                     cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 1)%r_coef))
                     cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 2)%r_coef))
                     cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 1)%r_coef))
                     cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 2)%r_coef))
 
                     ncoef_h1 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 1)
                     ncoef_h2 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 2)
                     ncoef_s1 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 1)
                     ncoef_s2 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 2)
                     cpassert(ncoef_h1 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 1))
                     cpassert(ncoef_h2 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 2))
                     cpassert(ncoef_s1 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 1))
                     cpassert(ncoef_s2 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 2))
 
                     DO icoef2 = 1, ncoef_h2
                        DO icoef1 = 1, ncoef_h1
                           cpc_h_tmp = mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2)
                           mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2) = &
                              cpc_h_tmp + mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
                           mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2) = &
                              cpc_h_tmp - mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
                        END DO
                     END DO
 
                     DO icoef2 = 1, ncoef_s2
                        DO icoef1 = 1, ncoef_s1
                           cpc_s_tmp = mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2)
                           mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2) = &
                              cpc_s_tmp + mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
                           mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2) = &
                              cpc_s_tmp - mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
                        END DO
                     END DO
                  END IF
               END DO
            END IF
         END IF
      END IF
 
      IF (mixing_method == gspace_mixing_nr) THEN
      ELSE IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
         IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
            DO ispin = 1, nspin
               DO iat = 1, natom
                  IF (mixing_store%paw(iat)) THEN
                     n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
                     n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
                     IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
                        DO ib = 1, nbuffer
                           ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
                        END DO
                     END IF
                     DO ib = 1, nbuffer
                        mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
                        mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
                        mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
                        mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
                     END DO
                  END IF
               END DO
            END DO
         END IF
      ELSE IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
         CALL mixing_buffer_allocate(mixing_store%rhoin_old, nvec, .true.)
         CALL mixing_buffer_allocate(mixing_store%drho_buffer, nvec, .true.)
         CALL mixing_buffer_allocate(mixing_store%last_res, nvec, .true.)
         IF (mixing_store%gmix_p) THEN
            IF (PRESENT(rho_atom)) THEN
               DO ispin = 1, nspin
                  DO iat = 1, natom
                     IF (mixing_store%paw(iat)) THEN
                        n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
                        n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
                        IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
                           ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
                        END IF
                        mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
                        mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
                        IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
                           DO ib = 1, nbuffer
                              ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
                              ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
                           END DO
                           ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
                           ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
                        END IF
                        DO ib = 1, nbuffer
                           mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
                           mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
                        END DO
                        mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
                        mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
                     END IF
                  END DO
               END DO
            END IF
         END IF
 
         IF (ASSOCIATED(mixing_store%p_metric)) DEALLOCATE (mixing_store%p_metric)
         ALLOCATE (mixing_store%p_metric(nvec))
         mixing_store%p_metric = 1.0_dp
         bconst = mixing_store%bconst
         g2min = 1.e30_dp
         DO ig = 1, ng
            IF (g2(ig) > 1.e-10_dp) g2min = min(g2min, g2(ig))
         END DO
         g2max = -1.e30_dp
         DO ig = 1, ng
            g2max = max(g2max, g2(ig))
         END DO
         CALL para_env%min(g2min)
         CALL para_env%max(g2max)
         ! fdamp/g2 varies between (bconst-1) and 0
         ! i.e. p_metric varies between bconst and 1
         ! fdamp = (bconst-1.0_dp)*g2min
         fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
         DO ig = 1, ng
            mixing_store%p_metric(ig) = (g2(ig) + fdamp)/max(g2(ig), 1.e-10_dp)
         END DO
         IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
         ! Keep the original G cutoff; auxiliary coordinates have unit metric.
         IF (nvec > ng) mixing_store%p_metric(mixing_store%ig_max + 1:ng) = 0.0_dp
      ELSE IF (mixing_method == multisecant_mixing_nr) THEN
         IF (ASSOCIATED(mixing_store%ig_global_index)) DEALLOCATE (mixing_store%ig_global_index)
         ALLOCATE (mixing_store%ig_global_index(nvec))
         mixing_store%ig_global_index = 0
         ig_count = 0
         DO iproc = 0, para_env%num_pe - 1
            IF (para_env%mepos == iproc) THEN
               DO ig = 1, nvec
                  ig_count = ig_count + 1
                  mixing_store%ig_global_index(ig) = ig_count
               END DO
            END IF
            CALL para_env%bcast(ig_count, iproc)
         END DO
      END IF
 
      IF (nvec > ng) mixing_store%ig_max = nvec
      CALL timestop(handle)
 
   END SUBROUTINE mixing_init
 
! **************************************************************************************************
!> \brief initialiation needed when charge mixing is used
!> \param mixing_store ...
!> \par History
!>      02.2019 created [JGH]
!> \author JGH
! **************************************************************************************************
   ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
      TYPE(mixing_storage_type), INTENT(INOUT)           :: mixing_store
 
      mixing_store%ncall = 0
      mixing_store%tb_scc_mixer_error = 0.0_dp
      mixing_store%tb_scc_mixer_step = 0
      IF (ASSOCIATED(mixing_store%acharge)) mixing_store%acharge = 0.0_dp
      IF (ASSOCIATED(mixing_store%dacharge)) mixing_store%dacharge = 0.0_dp
 
   END SUBROUTINE charge_mixing_init
 
END MODULE qs_mixing_utils