LCOV - code coverage report
Current view: top level - src - qs_diis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:744416f) Lines: 92.6 % 435 403
Test Date: 2026-09-20 02:09:09 Functions: 100.0 % 14 14

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Apply the direct inversion in the iterative subspace (DIIS) of Pulay
      10              : !>      in the framework of an SCF iteration for convergence acceleration
      11              : !> \par Literature
      12              : !>      - P. Pulay, Chem. Phys. Lett. 73, 393 (1980)
      13              : !>      - P. Pulay, J. Comput. Chem. 3, 556 (1982)
      14              : !> \par History
      15              : !>      - Changed to BLACS matrix usage (08.06.2001,MK)
      16              : !>      - rewritten to include LSD (1st attempt) (01.2003, Joost VandeVondele)
      17              : !>      - DIIS for ROKS (05.04.06,MK)
      18              : !>      - DIIS for k-points (04.2023, Augustin Bussy)
      19              : !> \author Matthias Krack (28.06.2000)
      20              : ! **************************************************************************************************
      21              : MODULE qs_diis
      22              :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_scale_and_add_fm,&
      23              :                                               cp_cfm_trace
      24              :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      25              :                                               cp_cfm_get_info,&
      26              :                                               cp_cfm_release,&
      27              :                                               cp_cfm_to_cfm,&
      28              :                                               cp_cfm_to_fm,&
      29              :                                               cp_cfm_type,&
      30              :                                               cp_fm_to_cfm
      31              :    USE cp_dbcsr_api,                    ONLY: &
      32              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
      33              :         dbcsr_set, dbcsr_transposed, dbcsr_type
      34              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot,&
      35              :                                               dbcsr_maxabs
      36              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm
      37              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      38              :                                               cp_fm_scale,&
      39              :                                               cp_fm_scale_and_add,&
      40              :                                               cp_fm_symm,&
      41              :                                               cp_fm_trace
      42              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      43              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      44              :                                               cp_fm_get_info,&
      45              :                                               cp_fm_maxabsval,&
      46              :                                               cp_fm_release,&
      47              :                                               cp_fm_set_all,&
      48              :                                               cp_fm_to_fm,&
      49              :                                               cp_fm_type
      50              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      51              :                                               cp_logger_type,&
      52              :                                               cp_to_string
      53              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      54              :                                               cp_print_key_unit_nr
      55              :    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
      56              :    USE input_section_types,             ONLY: section_vals_type
      57              :    USE kinds,                           ONLY: default_string_length,&
      58              :                                               dp
      59              :    USE local_gemm_api,                  ONLY: local_gemm_ctxt_type
      60              :    USE mathlib,                         ONLY: diag_complex,&
      61              :                                               diamat_all
      62              :    USE message_passing,                 ONLY: mp_para_env_type
      63              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      64              :    USE qs_diis_types,                   ONLY: qs_diis_buffer_type,&
      65              :                                               qs_diis_buffer_type_kp,&
      66              :                                               qs_diis_buffer_type_sparse
      67              :    USE qs_environment_types,            ONLY: get_qs_env,&
      68              :                                               qs_environment_type
      69              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      70              :                                               mo_set_type
      71              :    USE string_utilities,                ONLY: compress
      72              : #include "./base/base_uses.f90"
      73              : 
      74              :    IMPLICIT NONE
      75              : 
      76              :    PRIVATE
      77              : 
      78              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_diis'
      79              : 
      80              :    ! Public subroutines
      81              : 
      82              :    PUBLIC :: qs_diis_b_clear, &
      83              :              qs_diis_b_create, &
      84              :              qs_diis_b_step
      85              :    PUBLIC :: qs_diis_b_clear_sparse, &
      86              :              qs_diis_b_create_sparse, &
      87              :              qs_diis_b_step_4lscf
      88              :    PUBLIC :: qs_diis_b_clear_kp, &
      89              :              qs_diis_b_check_i_alloc_kp, &
      90              :              qs_diis_b_create_kp, &
      91              :              qs_diis_b_step_kp, &
      92              :              qs_diis_b_calc_err_kp, &
      93              :              qs_diis_b_info_kp
      94              : 
      95              : CONTAINS
      96              : 
      97              : ! **************************************************************************************************
      98              : !> \brief Allocates an SCF DIIS buffer
      99              : !> \param diis_buffer the buffer to create
     100              : !> \param nbuffer ...
     101              : !> \par History
     102              : !>      02.2003 created [fawzi]
     103              : !> \author fawzi
     104              : ! **************************************************************************************************
     105         4532 :    SUBROUTINE qs_diis_b_create(diis_buffer, nbuffer)
     106              : 
     107              :       TYPE(qs_diis_buffer_type), INTENT(OUT)             :: diis_buffer
     108              :       INTEGER, INTENT(in)                                :: nbuffer
     109              : 
     110              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_diis_b_create'
     111              : 
     112              :       INTEGER                                            :: handle
     113              : 
     114              : ! -------------------------------------------------------------------------
     115              : 
     116         4532 :       CALL timeset(routineN, handle)
     117              : 
     118         4532 :       NULLIFY (diis_buffer%b_matrix)
     119         4532 :       NULLIFY (diis_buffer%error)
     120         4532 :       NULLIFY (diis_buffer%param)
     121         4532 :       diis_buffer%nbuffer = nbuffer
     122         4532 :       diis_buffer%ncall = 0
     123              : 
     124         4532 :       CALL timestop(handle)
     125              : 
     126         4532 :    END SUBROUTINE qs_diis_b_create
     127              : 
     128              : ! **************************************************************************************************
     129              : !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
     130              : !>      variables and with a buffer size of nbuffer.
     131              : !> \param diis_buffer the buffer to initialize
     132              : !> \param matrix_struct the structure for the matrix of the buffer
     133              : !> \param nspin ...
     134              : !> \param scf_section ...
     135              : !> \par History
     136              : !>      - Creation (07.05.2001, Matthias Krack)
     137              : !>      - Changed to BLACS matrix usage (08.06.2001,MK)
     138              : !>      - DIIS for ROKS (05.04.06,MK)
     139              : !> \author Matthias Krack
     140              : !> \note
     141              : !>      check to allocate matrixes only when needed, using a linked list?
     142              : ! **************************************************************************************************
     143        98265 :    SUBROUTINE qs_diis_b_check_i_alloc(diis_buffer, matrix_struct, nspin, &
     144              :                                       scf_section)
     145              : 
     146              :       TYPE(qs_diis_buffer_type), INTENT(INOUT)           :: diis_buffer
     147              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     148              :       INTEGER, INTENT(IN)                                :: nspin
     149              :       TYPE(section_vals_type), POINTER                   :: scf_section
     150              : 
     151              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc'
     152              : 
     153              :       INTEGER                                            :: handle, ibuffer, ispin, nbuffer, &
     154              :                                                             output_unit
     155              :       TYPE(cp_logger_type), POINTER                      :: logger
     156              : 
     157              : ! -------------------------------------------------------------------------
     158              : 
     159        98265 :       CALL timeset(routineN, handle)
     160              : 
     161        98265 :       logger => cp_get_default_logger()
     162              : 
     163        98265 :       nbuffer = diis_buffer%nbuffer
     164              : 
     165        98265 :       IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
     166        34648 :          ALLOCATE (diis_buffer%error(nbuffer, nspin))
     167              : 
     168         7620 :          DO ispin = 1, nspin
     169        24292 :             DO ibuffer = 1, nbuffer
     170              :                CALL cp_fm_create(diis_buffer%error(ibuffer, ispin), &
     171              :                                  name="qs_diis_b%error("// &
     172              :                                  TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
     173              :                                  TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
     174        20840 :                                  matrix_struct=matrix_struct)
     175              :             END DO
     176              :          END DO
     177              :       END IF
     178              : 
     179        98265 :       IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
     180        34648 :          ALLOCATE (diis_buffer%param(nbuffer, nspin))
     181              : 
     182         7620 :          DO ispin = 1, nspin
     183        24292 :             DO ibuffer = 1, nbuffer
     184              :                CALL cp_fm_create(diis_buffer%param(ibuffer, ispin), &
     185              :                                  name="qs_diis_b%param("// &
     186              :                                  TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
     187              :                                  TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
     188        20840 :                                  matrix_struct=matrix_struct)
     189              :             END DO
     190              :          END DO
     191              :       END IF
     192              : 
     193        98265 :       IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
     194        13808 :          ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
     195       107012 :          diis_buffer%b_matrix = 0.0_dp
     196              :          output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
     197         3452 :                                             extension=".scfLog")
     198         3452 :          IF (output_unit > 0) THEN
     199              :             WRITE (UNIT=output_unit, FMT="(/,T9,A)") &
     200           20 :                "DIIS | The SCF DIIS buffer was allocated and initialized"
     201              :          END IF
     202              :          CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
     203         3452 :                                            "PRINT%DIIS_INFO")
     204              :       END IF
     205              : 
     206        98265 :       CALL timestop(handle)
     207              : 
     208        98265 :    END SUBROUTINE qs_diis_b_check_i_alloc
     209              : 
     210              : ! **************************************************************************************************
     211              : !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step.
     212              : !> \param diis_buffer ...
     213              : !> \param mo_array ...
     214              : !> \param kc ...
     215              : !> \param sc ...
     216              : !> \param delta ...
     217              : !> \param error_max ...
     218              : !> \param diis_step ...
     219              : !> \param eps_diis ...
     220              : !> \param nmixing ...
     221              : !> \param s_matrix ...
     222              : !> \param scf_section ...
     223              : !> \param roks ...
     224              : !> \par History
     225              : !>      - Creation (07.05.2001, Matthias Krack)
     226              : !>      - Changed to BLACS matrix usage (08.06.2001, MK)
     227              : !>      - 03.2003 rewamped [fawzi]
     228              : !>      - Adapted for high-spin ROKS (08.04.06,MK)
     229              : !> \author Matthias Krack
     230              : ! **************************************************************************************************
     231        98265 :    SUBROUTINE qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, &
     232              :                              diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
     233              : 
     234              :       TYPE(qs_diis_buffer_type), POINTER                 :: diis_buffer
     235              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     236              :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: kc
     237              :       TYPE(cp_fm_type), INTENT(IN)                       :: sc
     238              :       REAL(KIND=dp), INTENT(IN)                          :: delta
     239              :       REAL(KIND=dp), INTENT(OUT)                         :: error_max
     240              :       LOGICAL, INTENT(OUT)                               :: diis_step
     241              :       REAL(KIND=dp), INTENT(IN)                          :: eps_diis
     242              :       INTEGER, INTENT(IN), OPTIONAL                      :: nmixing
     243              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     244              :          POINTER                                         :: s_matrix
     245              :       TYPE(section_vals_type), POINTER                   :: scf_section
     246              :       LOGICAL, INTENT(IN), OPTIONAL                      :: roks
     247              : 
     248              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_diis_b_step'
     249              :       REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
     250              : 
     251              :       CHARACTER(LEN=2*default_string_length)             :: message
     252              :       INTEGER                                            :: handle, homo, ib, imo, ispin, jb, &
     253              :                                                             my_nmixing, nao, nb, nb1, nmo, nspin, &
     254              :                                                             output_unit
     255              :       LOGICAL                                            :: eigenvectors_discarded, my_roks
     256              :       REAL(KIND=dp)                                      :: maxocc, tmp
     257        98265 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ev, occ
     258        98265 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: occa, occb
     259        98265 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: a, b
     260              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     261              :       TYPE(cp_fm_type), POINTER                          :: c, new_errors, old_errors, parameters
     262              :       TYPE(cp_logger_type), POINTER                      :: logger
     263              : 
     264              : ! -------------------------------------------------------------------------
     265              : 
     266        98265 :       CALL timeset(routineN, handle)
     267              : 
     268        98265 :       nspin = SIZE(mo_array)
     269        98265 :       diis_step = .FALSE.
     270              : 
     271        98265 :       IF (PRESENT(roks)) THEN
     272         1018 :          my_roks = .TRUE.
     273         1018 :          nspin = 1
     274              :       ELSE
     275              :          my_roks = .FALSE.
     276              :       END IF
     277              : 
     278        98265 :       my_nmixing = 2
     279        98265 :       IF (PRESENT(nmixing)) my_nmixing = nmixing
     280              : 
     281        98265 :       NULLIFY (c, new_errors, old_errors, parameters, matrix_struct, a, b, occa, occb)
     282        98265 :       logger => cp_get_default_logger()
     283              : 
     284              :       ! Quick return, if no DIIS is requested
     285              : 
     286        98265 :       IF (diis_buffer%nbuffer < 1) THEN
     287            0 :          CALL timestop(handle)
     288              :          RETURN
     289              :       END IF
     290              : 
     291              :       CALL cp_fm_get_info(kc(1), &
     292        98265 :                           matrix_struct=matrix_struct)
     293              :       CALL qs_diis_b_check_i_alloc(diis_buffer, &
     294              :                                    matrix_struct=matrix_struct, &
     295              :                                    nspin=nspin, &
     296        98265 :                                    scf_section=scf_section)
     297              : 
     298        98265 :       error_max = 0.0_dp
     299              : 
     300        98265 :       ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
     301        98265 :       diis_buffer%ncall = diis_buffer%ncall + 1
     302        98265 :       nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
     303              : 
     304       212288 :       DO ispin = 1, nspin
     305              : 
     306              :          CALL get_mo_set(mo_set=mo_array(ispin), &
     307              :                          nao=nao, &
     308              :                          nmo=nmo, &
     309              :                          homo=homo, &
     310              :                          mo_coeff=c, &
     311              :                          occupation_numbers=occa, &
     312       114023 :                          maxocc=maxocc)
     313              : 
     314       114023 :          new_errors => diis_buffer%error(ib, ispin)
     315       114023 :          parameters => diis_buffer%param(ib, ispin)
     316              : 
     317              :          ! Copy the Kohn-Sham matrix K to the DIIS buffer
     318              : 
     319       114023 :          CALL cp_fm_to_fm(kc(ispin), parameters)
     320              : 
     321       114023 :          IF (my_roks) THEN
     322              : 
     323         3054 :             ALLOCATE (occ(nmo))
     324              : 
     325              :             CALL get_mo_set(mo_set=mo_array(2), &
     326         1018 :                             occupation_numbers=occb)
     327              : 
     328        16350 :             DO imo = 1, nmo
     329        16350 :                occ(imo) = SQRT(occa(imo) + occb(imo))
     330              :             END DO
     331              : 
     332         1018 :             CALL cp_fm_to_fm(c, sc)
     333         1018 :             CALL cp_fm_column_scale(sc, occ(1:homo))
     334              : 
     335              :             ! KC <- K*C
     336         1018 :             CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, parameters, sc, 0.0_dp, kc(ispin))
     337              : 
     338         1018 :             IF (PRESENT(s_matrix)) THEN
     339          558 :                CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
     340              :                ! SC <- S*C
     341          558 :                CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, new_errors, c, 0.0_dp, sc)
     342          558 :                CALL cp_fm_column_scale(sc, occ(1:homo))
     343              :             END IF
     344              : 
     345              :             ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
     346              :             ! or for an orthogonal basis
     347              :             ! new_errors <- KC*C^T - C*(KC)^T = K*P - P*K with S = I
     348         1018 :             CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
     349         1018 :             CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
     350              : 
     351         1018 :             DEALLOCATE (occ)
     352              : 
     353              :          ELSE
     354              : 
     355              :             ! KC <- K*C
     356       113005 :             CALL cp_fm_symm("L", "U", nao, homo, maxocc, parameters, c, 0.0_dp, kc(ispin))
     357              : 
     358       113005 :             IF (PRESENT(s_matrix)) THEN
     359              :                ! I guess that this copy can be avoided for LSD
     360        97029 :                CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, new_errors)
     361              :                ! sc <- S*C
     362        97029 :                CALL cp_fm_symm("L", "U", nao, homo, 2.0_dp, new_errors, c, 0.0_dp, sc)
     363              :                ! new_errors <- KC*(SC)^T - (SC)*(KC)^T = K*P*S - S*P*K
     364        97029 :                CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, sc, kc(ispin), 0.0_dp, new_errors)
     365        97029 :                CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), sc, -1.0_dp, new_errors)
     366              :             ELSE
     367              :                ! new_errors <- KC*(C)^T - C*(KC)^T = K*P - P*K
     368        15976 :                CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, c, kc(ispin), 0.0_dp, new_errors)
     369        15976 :                CALL parallel_gemm("N", "T", nao, nao, homo, 1.0_dp, kc(ispin), c, -1.0_dp, new_errors)
     370              :             END IF
     371              : 
     372              :          END IF
     373              : 
     374       114023 :          CALL cp_fm_maxabsval(new_errors, tmp)
     375       326311 :          error_max = MAX(error_max, tmp)
     376              : 
     377              :       END DO
     378              : 
     379              :       ! Check, if a DIIS step is appropriate
     380              : 
     381        98265 :       diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
     382              : 
     383              :       output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
     384        98265 :                                          extension=".scfLog")
     385        98265 :       IF (output_unit > 0) THEN
     386              :          WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
     387          251 :             "DIIS | Current SCF DIIS buffer size:         ", nb, &
     388          251 :             "DIIS | Maximum SCF DIIS error vector element:", error_max, &
     389          251 :             "DIIS | Current SCF convergence:              ", delta, &
     390          502 :             "DIIS | Threshold value for a DIIS step:      ", eps_diis
     391          251 :          IF (error_max < eps_diis) THEN
     392              :             WRITE (UNIT=output_unit, FMT="(T9,A)") &
     393          102 :                "DIIS | => The SCF DIIS buffer will be updated"
     394              :          ELSE
     395              :             WRITE (UNIT=output_unit, FMT="(T9,A)") &
     396          149 :                "DIIS | => No update of the SCF DIIS buffer"
     397              :          END IF
     398          251 :          IF (diis_step .AND. (error_max < eps_diis)) THEN
     399              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
     400           62 :                "DIIS | => A SCF DIIS step will be performed"
     401              :          ELSE
     402              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
     403          189 :                "DIIS | => No SCF DIIS step will be performed"
     404              :          END IF
     405              :       END IF
     406              : 
     407              :       ! Update the SCF DIIS buffer
     408              : 
     409        98265 :       IF (error_max < eps_diis) THEN
     410              : 
     411        84577 :          b => diis_buffer%b_matrix
     412              : 
     413       358895 :          DO jb = 1, nb
     414       274318 :             b(jb, ib) = 0.0_dp
     415       596230 :             DO ispin = 1, nspin
     416       321912 :                old_errors => diis_buffer%error(jb, ispin)
     417       321912 :                new_errors => diis_buffer%error(ib, ispin)
     418       321912 :                CALL cp_fm_trace(old_errors, new_errors, tmp)
     419       596230 :                b(jb, ib) = b(jb, ib) + tmp
     420              :             END DO
     421       358895 :             b(ib, jb) = b(jb, ib)
     422              :          END DO
     423              : 
     424              :       ELSE
     425              : 
     426        13688 :          diis_step = .FALSE.
     427              : 
     428              :       END IF
     429              : 
     430              :       ! Perform DIIS step
     431              : 
     432        98265 :       IF (diis_step) THEN
     433              : 
     434        62895 :          nb1 = nb + 1
     435              : 
     436       251580 :          ALLOCATE (a(nb1, nb1))
     437       188685 :          ALLOCATE (b(nb1, nb1))
     438       188685 :          ALLOCATE (ev(nb1))
     439              : 
     440              :          ! Set up the linear DIIS equation system
     441              : 
     442      2215279 :          b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
     443              : 
     444       289711 :          b(1:nb, nb1) = -1.0_dp
     445       289711 :          b(nb1, 1:nb) = -1.0_dp
     446        62895 :          b(nb1, nb1) = 0.0_dp
     447              : 
     448              :          ! Solve the linear DIIS equation system
     449              : 
     450       352606 :          ev(1:nb1) = 0.0_dp
     451        62895 :          CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
     452              : 
     453      3374123 :          a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
     454              : 
     455        62895 :          eigenvectors_discarded = .FALSE.
     456              : 
     457       352606 :          DO jb = 1, nb1
     458       352606 :             IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
     459        44256 :                IF (output_unit > 0) THEN
     460            5 :                   IF (.NOT. eigenvectors_discarded) THEN
     461              :                      WRITE (UNIT=output_unit, FMT="(T9,A)") &
     462            5 :                         "DIIS | Checking eigenvalues of the DIIS error matrix"
     463              :                   END IF
     464              :                   WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
     465            5 :                      "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
     466           10 :                      "threshold ", eigenvalue_threshold
     467            5 :                   CALL compress(message)
     468            5 :                   WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
     469            5 :                   eigenvectors_discarded = .TRUE.
     470              :                END IF
     471       259036 :                a(1:nb1, jb) = 0.0_dp
     472              :             ELSE
     473      1396578 :                a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
     474              :             END IF
     475              :          END DO
     476              : 
     477        62895 :          IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
     478              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
     479            5 :                "DIIS | The corresponding eigenvectors were discarded"
     480              :          END IF
     481              : 
     482      1655614 :          ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
     483              : 
     484              :          ! Update Kohn-Sham matrix
     485              : 
     486       136308 :          DO ispin = 1, nspin
     487        73413 :             CALL cp_fm_set_all(kc(ispin), 0.0_dp)
     488       402546 :             DO jb = 1, nb
     489       266238 :                parameters => diis_buffer%param(jb, ispin)
     490       339651 :                CALL cp_fm_scale_and_add(1.0_dp, kc(ispin), -ev(jb), parameters)
     491              :             END DO
     492              :          END DO
     493              : 
     494        62895 :          DEALLOCATE (a)
     495        62895 :          DEALLOCATE (b)
     496        62895 :          DEALLOCATE (ev)
     497              : 
     498              :       ELSE
     499              : 
     500        75980 :          DO ispin = 1, nspin
     501        40610 :             parameters => diis_buffer%param(ib, ispin)
     502        75980 :             CALL cp_fm_to_fm(parameters, kc(ispin))
     503              :          END DO
     504              : 
     505              :       END IF
     506              : 
     507              :       CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
     508        98265 :                                         "PRINT%DIIS_INFO")
     509              : 
     510        98265 :       CALL timestop(handle)
     511              : 
     512       196530 :    END SUBROUTINE qs_diis_b_step
     513              : 
     514              : ! **************************************************************************************************
     515              : !> \brief clears the buffer
     516              : !> \param diis_buffer the buffer to clear
     517              : !> \par History
     518              : !>      02.2003 created [fawzi]
     519              : !> \author fawzi
     520              : ! **************************************************************************************************
     521        16410 :    PURE SUBROUTINE qs_diis_b_clear(diis_buffer)
     522              : 
     523              :       TYPE(qs_diis_buffer_type), INTENT(INOUT)           :: diis_buffer
     524              : 
     525        16410 :       diis_buffer%ncall = 0
     526              : 
     527        16410 :    END SUBROUTINE qs_diis_b_clear
     528              : 
     529              : ! **************************************************************************************************
     530              : !> \brief Update the SCF DIIS buffer in linear scaling SCF (LS-SCF),
     531              : !>        and if appropriate does a diis step.
     532              : !> \param diis_buffer ...
     533              : !> \param qs_env ...
     534              : !> \param ls_scf_env ...
     535              : !> \param unit_nr ...
     536              : !> \param iscf ...
     537              : !> \param diis_step ...
     538              : !> \param eps_diis ...
     539              : !> \param nmixing ...
     540              : !> \param s_matrix ...
     541              : !> \param threshold ...
     542              : !> \par History
     543              : !>      - Adapted for LS-SCF (10-11-14) from qs_diis_b_step
     544              : !> \author Fredy W. Aquino
     545              : ! **************************************************************************************************
     546              : 
     547           18 :    SUBROUTINE qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, iscf, &
     548              :                                    diis_step, eps_diis, nmixing, s_matrix, threshold)
     549              : ! Note.- Input: ls_scf_env%matrix_p(ispin) , Density Matrix
     550              : !               matrix_ks (from qs_env)    , Kohn-Sham Matrix  (IN/OUT)
     551              : 
     552              :       TYPE(qs_diis_buffer_type_sparse), POINTER          :: diis_buffer
     553              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     554              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     555              :       INTEGER, INTENT(IN)                                :: unit_nr, iscf
     556              :       LOGICAL, INTENT(OUT)                               :: diis_step
     557              :       REAL(KIND=dp), INTENT(IN)                          :: eps_diis
     558              :       INTEGER, INTENT(IN), OPTIONAL                      :: nmixing
     559              :       TYPE(dbcsr_type), OPTIONAL                         :: s_matrix
     560              :       REAL(KIND=dp), INTENT(IN)                          :: threshold
     561              : 
     562              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_step_4lscf'
     563              :       REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
     564              : 
     565              :       INTEGER                                            :: handle, ib, ispin, jb, my_nmixing, nb, &
     566              :                                                             nb1, nspin
     567              :       REAL(KIND=dp)                                      :: error_max, tmp
     568           18 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ev
     569           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: a, b
     570              :       TYPE(cp_logger_type), POINTER                      :: logger
     571           18 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     572              :       TYPE(dbcsr_type)                                   :: matrix_KSerr_t, matrix_tmp
     573              :       TYPE(dbcsr_type), POINTER                          :: new_errors, old_errors, parameters
     574              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     575              : 
     576           18 :       CALL timeset(routineN, handle)
     577           18 :       IF (ls_scf_env%do_pao) THEN
     578            0 :          CPABORT("LS_SCF%LS_DIIS not compatible with PAO")
     579              :       END IF
     580           18 :       nspin = ls_scf_env%nspins
     581           18 :       diis_step = .FALSE.
     582           18 :       my_nmixing = 2
     583           18 :       IF (PRESENT(nmixing)) my_nmixing = nmixing
     584           18 :       NULLIFY (new_errors, old_errors, parameters, a, b)
     585           18 :       logger => cp_get_default_logger()
     586              :       ! Quick return, if no DIIS is requested
     587           18 :       IF (diis_buffer%nbuffer < 1) THEN
     588            0 :          CALL timestop(handle)
     589              :          RETURN
     590              :       END IF
     591              : 
     592              : ! Getting current Kohn-Sham matrix from qs_env
     593              :       CALL get_qs_env(qs_env, &
     594              :                       para_env=para_env, &
     595           18 :                       matrix_ks=matrix_ks)
     596              :       CALL qs_diis_b_check_i_alloc_sparse( &
     597              :          diis_buffer, &
     598              :          ls_scf_env, &
     599           18 :          nspin)
     600           18 :       error_max = 0.0_dp
     601              : 
     602           18 :       ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
     603           18 :       diis_buffer%ncall = diis_buffer%ncall + 1
     604           18 :       nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
     605              : ! Create scratch arrays
     606              :       CALL dbcsr_create(matrix_tmp, &
     607              :                         template=ls_scf_env%matrix_ks(1), &
     608           18 :                         matrix_type='N')
     609           18 :       CALL dbcsr_set(matrix_tmp, 0.0_dp) ! reset matrix
     610              :       CALL dbcsr_create(matrix_KSerr_t, &
     611              :                         template=ls_scf_env%matrix_ks(1), &
     612           18 :                         matrix_type='N')
     613           18 :       CALL dbcsr_set(matrix_KSerr_t, 0.0_dp) ! reset matrix
     614              : 
     615           46 :       DO ispin = 1, nspin ! ------ Loop-ispin----START
     616              : 
     617           28 :          new_errors => diis_buffer%error(ib, ispin)%matrix
     618           28 :          parameters => diis_buffer%param(ib, ispin)%matrix
     619              :          ! Copy the Kohn-Sham matrix K to the DIIS buffer
     620              :          CALL dbcsr_copy(parameters, & ! out
     621           28 :                          matrix_ks(ispin)%matrix) ! in
     622              : 
     623           28 :          IF (PRESENT(s_matrix)) THEN ! if-s_matrix ---------- START
     624              : ! Calculate Kohn-Sham error (non-orthogonal)= K*P*S-(K*P*S)^T
     625              : ! matrix_tmp = P*S
     626              :             CALL dbcsr_multiply("N", "N", &
     627              :                                 1.0_dp, ls_scf_env%matrix_p(ispin), &
     628              :                                 s_matrix, &
     629              :                                 0.0_dp, matrix_tmp, &
     630           28 :                                 filter_eps=threshold)
     631              : ! new_errors= K*P*S
     632              :             CALL dbcsr_multiply("N", "N", &
     633              :                                 1.0_dp, matrix_ks(ispin)%matrix, &
     634              :                                 matrix_tmp, &
     635              :                                 0.0_dp, new_errors, &
     636           28 :                                 filter_eps=threshold)
     637              : ! matrix_KSerr_t= transpose(K*P*S)
     638              :             CALL dbcsr_transposed(matrix_KSerr_t, &
     639           28 :                                   new_errors)
     640              : ! new_errors=K*P*S-transpose(K*P*S)
     641              :             CALL dbcsr_add(new_errors, &
     642              :                            matrix_KSerr_t, &
     643           28 :                            1.0_dp, -1.0_dp)
     644              :          ELSE ! if-s_matrix ---------- MID
     645              : ! Calculate Kohn-Sham error (orthogonal)= K*P - P*K
     646              : ! new_errors=K*P
     647              :             CALL dbcsr_multiply("N", "N", &
     648              :                                 1.0_dp, matrix_ks(ispin)%matrix, &
     649              :                                 ls_scf_env%matrix_p(ispin), &
     650              :                                 0.0_dp, new_errors, &
     651            0 :                                 filter_eps=threshold)
     652              : ! matrix_KSerr_t= transpose(K*P)
     653              :             CALL dbcsr_transposed(matrix_KSerr_t, &
     654            0 :                                   new_errors)
     655              : ! new_errors=K*P-transpose(K*P)
     656              :             CALL dbcsr_add(new_errors, &
     657              :                            matrix_KSerr_t, &
     658            0 :                            1.0_dp, -1.0_dp)
     659              :          END IF ! if-s_matrix ---------- END
     660              : 
     661           28 :          tmp = dbcsr_maxabs(new_errors)
     662           46 :          error_max = MAX(error_max, tmp)
     663              : 
     664              :       END DO ! ------ Loop-ispin----END
     665              : 
     666              :       ! Check, if a DIIS step is appropriate
     667              : 
     668           18 :       diis_step = (diis_buffer%ncall >= my_nmixing)
     669              : 
     670           18 :       IF (unit_nr > 0) THEN
     671              :          WRITE (unit_nr, '(A29,I3,A3,4(I3,A1))') &
     672            9 :             "DIIS: (ncall,nbuffer,ib,nb)=(", iscf, ")=(", &
     673           18 :             diis_buffer%ncall, ",", diis_buffer%nbuffer, ",", ib, ",", nb, ")"
     674              :          WRITE (unit_nr, '(A57,I3,A3,L1,A1,F10.8,A1,F4.2,A1,L1,A1)') &
     675            9 :             "DIIS: (diis_step,error_max,eps_diis,error_max<eps_diis)=(", &
     676            9 :             iscf, ")=(", diis_step, ",", error_max, ",", eps_diis, ",", &
     677           18 :             (error_max < eps_diis), ")"
     678              :          WRITE (unit_nr, '(A75)') &
     679            9 :             "DIIS: diis_step=T : Perform DIIS  error_max<eps_diis=T : Update DIIS buffer"
     680              :       END IF
     681              : 
     682              :       ! Update the SCF DIIS buffer
     683           18 :       IF (error_max < eps_diis) THEN
     684           18 :          b => diis_buffer%b_matrix
     685           66 :          DO jb = 1, nb
     686           48 :             b(jb, ib) = 0.0_dp
     687          124 :             DO ispin = 1, nspin
     688           76 :                old_errors => diis_buffer%error(jb, ispin)%matrix
     689           76 :                new_errors => diis_buffer%error(ib, ispin)%matrix
     690              :                CALL dbcsr_dot(old_errors, &
     691              :                               new_errors, &
     692           76 :                               tmp) ! out : < f_i | f_j >
     693          124 :                b(jb, ib) = b(jb, ib) + tmp
     694              :             END DO ! end-loop-ispin
     695           66 :             b(ib, jb) = b(jb, ib)
     696              :          END DO ! end-loop-jb
     697              :       ELSE
     698            0 :          diis_step = .FALSE.
     699              :       END IF
     700              : 
     701              :       ! Perform DIIS step
     702           18 :       IF (diis_step) THEN
     703           14 :          nb1 = nb + 1
     704           56 :          ALLOCATE (a(nb1, nb1))
     705           42 :          ALLOCATE (b(nb1, nb1))
     706           42 :          ALLOCATE (ev(nb1))
     707              :          ! Set up the linear DIIS equation system
     708          398 :          b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
     709           58 :          b(1:nb, nb1) = -1.0_dp
     710           58 :          b(nb1, 1:nb) = -1.0_dp
     711           14 :          b(nb1, nb1) = 0.0_dp
     712              :          ! Solve the linear DIIS equation system
     713           14 :          CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
     714          630 :          a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)
     715           72 :          DO jb = 1, nb1
     716           72 :             IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
     717            0 :                a(1:nb1, jb) = 0.0_dp
     718              :             ELSE
     719          308 :                a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
     720              :             END IF
     721              :          END DO ! end-loop-jb
     722              : 
     723          308 :          ev(1:nb) = MATMUL(a(1:nb, 1:nb1), b(nb1, 1:nb1))
     724              : 
     725              :          ! Update Kohn-Sham matrix
     726           14 :          IF (iscf >= ls_scf_env%iter_ini_diis) THEN ! if-iscf-to-updateKS------ START
     727              : 
     728           14 :             IF (unit_nr > 0) THEN
     729            7 :                WRITE (unit_nr, '(A40,I3)') 'DIIS: Updating Kohn-Sham matrix at iscf=', iscf
     730              :             END IF
     731              : 
     732           36 :             DO ispin = 1, nspin
     733              :                CALL dbcsr_set(matrix_ks(ispin)%matrix, & ! reset matrix
     734           22 :                               0.0_dp)
     735          106 :                DO jb = 1, nb
     736           70 :                   parameters => diis_buffer%param(jb, ispin)%matrix
     737              :                   CALL dbcsr_add(matrix_ks(ispin)%matrix, parameters, &
     738           92 :                                  1.0_dp, -ev(jb))
     739              :                END DO ! end-loop-jb
     740              :             END DO ! end-loop-ispin
     741              :          END IF ! if-iscf-to-updateKS------ END
     742              : 
     743           14 :          DEALLOCATE (a)
     744           14 :          DEALLOCATE (b)
     745           14 :          DEALLOCATE (ev)
     746              : 
     747              :       ELSE
     748           10 :          DO ispin = 1, nspin
     749            6 :             parameters => diis_buffer%param(ib, ispin)%matrix
     750              :             CALL dbcsr_copy(parameters, & ! out
     751           10 :                             matrix_ks(ispin)%matrix) ! in
     752              :          END DO ! end-loop-ispin
     753              :       END IF
     754           18 :       CALL dbcsr_release(matrix_tmp)
     755           18 :       CALL dbcsr_release(matrix_KSerr_t)
     756           18 :       CALL timestop(handle)
     757              : 
     758           18 :    END SUBROUTINE qs_diis_b_step_4lscf
     759              : 
     760              : ! **************************************************************************************************
     761              : !> \brief Allocate and initialize a DIIS buffer with a buffer size of nbuffer.
     762              : !> \param diis_buffer the buffer to initialize
     763              : !> \param ls_scf_env ...
     764              : !> \param nspin ...
     765              : !> \par History
     766              : !>      - Adapted from qs_diis_b_check_i_alloc for sparse matrices and
     767              : !>        used in LS-SCF module (ls_scf_main) (10-11-14)
     768              : !> \author Fredy W. Aquino
     769              : !> \note
     770              : !>      check to allocate matrices only when needed
     771              : ! **************************************************************************************************
     772              : 
     773           18 :    SUBROUTINE qs_diis_b_check_i_alloc_sparse(diis_buffer, ls_scf_env, &
     774              :                                              nspin)
     775              : 
     776              :       TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT)    :: diis_buffer
     777              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     778              :       INTEGER, INTENT(IN)                                :: nspin
     779              : 
     780              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_sparse'
     781              : 
     782              :       INTEGER                                            :: handle, ibuffer, ispin, nbuffer
     783              :       TYPE(cp_logger_type), POINTER                      :: logger
     784              : 
     785              : ! -------------------------------------------------------------------------
     786              : 
     787           18 :       CALL timeset(routineN, handle)
     788              : 
     789           18 :       logger => cp_get_default_logger()
     790              : 
     791           18 :       nbuffer = diis_buffer%nbuffer
     792              : 
     793           18 :       IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
     794           46 :          ALLOCATE (diis_buffer%error(nbuffer, nspin))
     795              : 
     796           10 :          DO ispin = 1, nspin
     797           34 :             DO ibuffer = 1, nbuffer
     798           24 :                ALLOCATE (diis_buffer%error(ibuffer, ispin)%matrix)
     799              : 
     800              :                CALL dbcsr_create(diis_buffer%error(ibuffer, ispin)%matrix, &
     801              :                                  template=ls_scf_env%matrix_ks(1), &
     802           30 :                                  matrix_type='N')
     803              :             END DO
     804              :          END DO
     805              :       END IF
     806              : 
     807           18 :       IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
     808           46 :          ALLOCATE (diis_buffer%param(nbuffer, nspin))
     809              : 
     810           10 :          DO ispin = 1, nspin
     811           34 :             DO ibuffer = 1, nbuffer
     812           24 :                ALLOCATE (diis_buffer%param(ibuffer, ispin)%matrix)
     813              :                CALL dbcsr_create(diis_buffer%param(ibuffer, ispin)%matrix, &
     814              :                                  template=ls_scf_env%matrix_ks(1), &
     815           30 :                                  matrix_type='N')
     816              :             END DO
     817              :          END DO
     818              :       END IF
     819              : 
     820           18 :       IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
     821           16 :          ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
     822              : 
     823          124 :          diis_buffer%b_matrix = 0.0_dp
     824              :       END IF
     825              : 
     826           18 :       CALL timestop(handle)
     827              : 
     828           18 :    END SUBROUTINE qs_diis_b_check_i_alloc_sparse
     829              : 
     830              : ! **************************************************************************************************
     831              : !> \brief clears the DIIS buffer in LS-SCF calculation
     832              : !> \param diis_buffer the buffer to clear
     833              : !> \par History
     834              : !>      10-11-14 created [FA] modified from qs_diis_b_clear
     835              : !> \author Fredy W. Aquino
     836              : ! **************************************************************************************************
     837              : 
     838            4 :    PURE SUBROUTINE qs_diis_b_clear_sparse(diis_buffer)
     839              : 
     840              :       TYPE(qs_diis_buffer_type_sparse), INTENT(INOUT)    :: diis_buffer
     841              : 
     842            4 :       diis_buffer%ncall = 0
     843              : 
     844            4 :    END SUBROUTINE qs_diis_b_clear_sparse
     845              : 
     846              : ! **************************************************************************************************
     847              : !> \brief Allocates an SCF DIIS buffer for LS-SCF calculation
     848              : !> \param diis_buffer the buffer to create
     849              : !> \param nbuffer ...
     850              : !> \par History
     851              : !>      10-11-14 created [FA] modified from qs_diis_b_create
     852              : !> \author Fredy W. Aquino
     853              : ! **************************************************************************************************
     854            4 :    PURE SUBROUTINE qs_diis_b_create_sparse(diis_buffer, nbuffer)
     855              : 
     856              :       TYPE(qs_diis_buffer_type_sparse), INTENT(OUT)      :: diis_buffer
     857              :       INTEGER, INTENT(in)                                :: nbuffer
     858              : 
     859              :       NULLIFY (diis_buffer%b_matrix)
     860              :       NULLIFY (diis_buffer%error)
     861              :       NULLIFY (diis_buffer%param)
     862            4 :       diis_buffer%nbuffer = nbuffer
     863            4 :       diis_buffer%ncall = 0
     864              : 
     865            4 :    END SUBROUTINE qs_diis_b_create_sparse
     866              : 
     867              : ! **************************************************************************************************
     868              : !> \brief Allocates an SCF DIIS buffer for k-points
     869              : !> \param diis_buffer the buffer to create
     870              : !> \param nbuffer ...
     871              : ! **************************************************************************************************
     872         2560 :    SUBROUTINE qs_diis_b_create_kp(diis_buffer, nbuffer)
     873              : 
     874              :       TYPE(qs_diis_buffer_type_kp), INTENT(OUT)          :: diis_buffer
     875              :       INTEGER, INTENT(in)                                :: nbuffer
     876              : 
     877              :       NULLIFY (diis_buffer%b_matrix)
     878              :       NULLIFY (diis_buffer%error)
     879              :       NULLIFY (diis_buffer%param)
     880              :       NULLIFY (diis_buffer%smat)
     881         2560 :       diis_buffer%nbuffer = nbuffer
     882         2560 :       diis_buffer%ncall = 0
     883              : 
     884         2560 :    END SUBROUTINE qs_diis_b_create_kp
     885              : 
     886              : ! **************************************************************************************************
     887              : !> \brief Allocate and initialize a DIIS buffer for nao*nao parameter
     888              : !>      variables and with a buffer size of nbuffer, in the k-point case
     889              : !> \param diis_buffer the buffer to initialize
     890              : !> \param matrix_struct the structure for the matrix of the buffer note: this is in the kp subgroup
     891              : !> \param nspin ...
     892              : !> \param nkp ...
     893              : !> \param scf_section ...
     894              : ! **************************************************************************************************
     895        27498 :    SUBROUTINE qs_diis_b_check_i_alloc_kp(diis_buffer, matrix_struct, nspin, nkp, scf_section)
     896              : 
     897              :       TYPE(qs_diis_buffer_type_kp), INTENT(INOUT)        :: diis_buffer
     898              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     899              :       INTEGER, INTENT(IN)                                :: nspin, nkp
     900              :       TYPE(section_vals_type), POINTER                   :: scf_section
     901              : 
     902              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_check_i_alloc_kp'
     903              : 
     904              :       INTEGER                                            :: handle, ibuffer, ikp, ispin, nbuffer, &
     905              :                                                             output_unit
     906              :       TYPE(cp_logger_type), POINTER                      :: logger
     907              : 
     908              : ! -------------------------------------------------------------------------
     909              : 
     910        27498 :       CALL timeset(routineN, handle)
     911              : 
     912        27498 :       logger => cp_get_default_logger()
     913              : 
     914        27498 :       nbuffer = diis_buffer%nbuffer
     915              : 
     916        27498 :       IF (.NOT. ASSOCIATED(diis_buffer%error)) THEN
     917        35863 :          ALLOCATE (diis_buffer%error(nbuffer, nspin, nkp))
     918              : 
     919         5949 :          DO ikp = 1, nkp
     920        10271 :             DO ispin = 1, nspin
     921        25483 :                DO ibuffer = 1, nbuffer
     922              :                   CALL cp_cfm_create(diis_buffer%error(ibuffer, ispin, ikp), &
     923              :                                      name="qs_diis_b%error("// &
     924              :                                      TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
     925              :                                      TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
     926        21610 :                                      matrix_struct=matrix_struct)
     927              :                END DO
     928              :             END DO
     929              :          END DO
     930              :       END IF
     931              : 
     932        27498 :       IF (.NOT. ASSOCIATED(diis_buffer%param)) THEN
     933        35863 :          ALLOCATE (diis_buffer%param(nbuffer, nspin, nkp))
     934              : 
     935         5949 :          DO ikp = 1, nkp
     936        10271 :             DO ispin = 1, nspin
     937        25483 :                DO ibuffer = 1, nbuffer
     938              :                   CALL cp_cfm_create(diis_buffer%param(ibuffer, ispin, ikp), &
     939              :                                      name="qs_diis_b%param("// &
     940              :                                      TRIM(ADJUSTL(cp_to_string(ibuffer)))//","// &
     941              :                                      TRIM(ADJUSTL(cp_to_string(ibuffer)))//")", &
     942        21610 :                                      matrix_struct=matrix_struct)
     943              :                END DO
     944              :             END DO
     945              :          END DO
     946              :       END IF
     947              : 
     948        27498 :       IF (.NOT. ASSOCIATED(diis_buffer%smat)) THEN
     949        10101 :          ALLOCATE (diis_buffer%smat(nkp))
     950         5949 :          DO ikp = 1, nkp
     951              :             CALL cp_cfm_create(diis_buffer%smat(ikp), &
     952              :                                name="kp_cfm_smat("//TRIM(cp_to_string(ikp))//")", &
     953         5949 :                                matrix_struct=matrix_struct)
     954              :          END DO
     955              :       END IF
     956              : 
     957        27498 :       IF (.NOT. ASSOCIATED(diis_buffer%b_matrix)) THEN
     958         8304 :          ALLOCATE (diis_buffer%b_matrix(nbuffer + 1, nbuffer + 1))
     959        64356 :          diis_buffer%b_matrix = 0.0_dp
     960              :          output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
     961         2076 :                                             extension=".scfLog")
     962         2076 :          IF (output_unit > 0) THEN
     963              :             WRITE (UNIT=output_unit, FMT="(/,T9,A)") &
     964            0 :                "DIIS | The SCF DIIS buffer was allocated and initialized"
     965              :          END IF
     966              :          CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
     967         2076 :                                            "PRINT%DIIS_INFO")
     968              :       END IF
     969              : 
     970        27498 :       CALL timestop(handle)
     971              : 
     972        27498 :    END SUBROUTINE qs_diis_b_check_i_alloc_kp
     973              : 
     974              : ! **************************************************************************************************
     975              : !> \brief clears the buffer
     976              : !> \param diis_buffer the buffer to clear
     977              : ! **************************************************************************************************
     978         3210 :    PURE SUBROUTINE qs_diis_b_clear_kp(diis_buffer)
     979              : 
     980              :       TYPE(qs_diis_buffer_type_kp), INTENT(INOUT)        :: diis_buffer
     981              : 
     982         3210 :       diis_buffer%ncall = 0
     983              : 
     984         3210 :    END SUBROUTINE qs_diis_b_clear_kp
     985              : 
     986              : ! **************************************************************************************************
     987              : !> \brief Update info about the current buffer step ib and the current number of buffers nb
     988              : !> \param diis_buffer ...
     989              : !> \param ib ...
     990              : !> \param nb ...
     991              : ! **************************************************************************************************
     992        27498 :    SUBROUTINE qs_diis_b_info_kp(diis_buffer, ib, nb)
     993              :       TYPE(qs_diis_buffer_type_kp), POINTER              :: diis_buffer
     994              :       INTEGER, INTENT(OUT)                               :: ib, nb
     995              : 
     996        27498 :       ib = MODULO(diis_buffer%ncall, diis_buffer%nbuffer) + 1
     997        27498 :       diis_buffer%ncall = diis_buffer%ncall + 1
     998        27498 :       nb = MIN(diis_buffer%ncall, diis_buffer%nbuffer)
     999              : 
    1000        27498 :    END SUBROUTINE qs_diis_b_info_kp
    1001              : 
    1002              : ! **************************************************************************************************
    1003              : !> \brief Calculate and store the error for a given k-point
    1004              : !> \param diis_buffer ...
    1005              : !> \param ib ...
    1006              : !> \param mos ...
    1007              : !> \param kc ...
    1008              : !> \param sc ...
    1009              : !> \param ispin ...
    1010              : !> \param ikp ...
    1011              : !> \param cmos caller-owned scratch; the controller preallocates it and the DIIS buffer
    1012              : !> \param real_wfn ...
    1013              : !> \param local use local GEMM on a one-rank k-point group
    1014              : !> \param gemm_ctx private local GEMM context, unused for distributed groups
    1015              : !> \note We assume that we always have an overlap matrix.
    1016              : !> TODO: do we need to pass the kp weight for the back Fourier transform?
    1017              : ! **************************************************************************************************
    1018        97664 :    SUBROUTINE qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, cmos, real_wfn, local, gemm_ctx)
    1019              :       TYPE(qs_diis_buffer_type_kp), POINTER              :: diis_buffer
    1020              :       INTEGER, INTENT(IN)                                :: ib
    1021              :       TYPE(mo_set_type), DIMENSION(:, :), POINTER        :: mos
    1022              :       TYPE(cp_cfm_type), INTENT(INOUT)                   :: kc, sc
    1023              :       INTEGER, INTENT(IN)                                :: ispin, ikp
    1024              :       TYPE(cp_cfm_type), INTENT(INOUT)                   :: cmos
    1025              :       LOGICAL, INTENT(IN), OPTIONAL                      :: real_wfn, local
    1026              :       TYPE(local_gemm_ctxt_type), INTENT(INOUT)          :: gemm_ctx
    1027              : 
    1028              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_diis_b_calc_err_kp'
    1029              : 
    1030              :       INTEGER                                            :: handle, homo, nao
    1031              :       LOGICAL                                            :: my_real_wfn
    1032              :       REAL(dp)                                           :: maxocc
    1033              :       TYPE(cp_cfm_type), POINTER                         :: new_errors, parameters, smat
    1034              :       TYPE(cp_fm_type), POINTER                          :: imos, rmos
    1035              : 
    1036        48832 :       NULLIFY (imos, rmos, parameters, new_errors, smat)
    1037              : 
    1038        48832 :       CALL timeset(routineN, handle)
    1039        48832 :       my_real_wfn = .FALSE.
    1040        48832 :       IF (PRESENT(real_wfn)) my_real_wfn = real_wfn
    1041              : 
    1042              :       !Calculate the error for this given k-point, store the KS matrix as well as the ovlp matrix
    1043              :       !All of this happens within the kp subgroups
    1044              : 
    1045              :       ! Quick return, if no DIIS is requested
    1046        48832 :       IF (diis_buffer%nbuffer < 1) THEN
    1047            0 :          CALL timestop(handle)
    1048            0 :          RETURN
    1049              :       END IF
    1050        48832 :       CALL get_mo_set(mos(1, ispin), nao=nao, homo=homo, mo_coeff=rmos, maxocc=maxocc)
    1051              : 
    1052        48832 :       IF (my_real_wfn) THEN
    1053              :          CALL cp_cfm_scale_and_add_fm(CMPLX(0.0_dp, KIND=dp), cmos, &
    1054            0 :                                       CMPLX(1.0_dp, KIND=dp), rmos)
    1055              :       ELSE
    1056        48832 :          CALL get_mo_set(mos(2, ispin), mo_coeff=imos)
    1057        48832 :          CALL cp_fm_to_cfm(rmos, imos, cmos)
    1058              :       END IF
    1059              : 
    1060        48832 :       new_errors => diis_buffer%error(ib, ispin, ikp)
    1061        48832 :       parameters => diis_buffer%param(ib, ispin, ikp)
    1062        48832 :       smat => diis_buffer%smat(ikp)
    1063              : 
    1064              :       !copy the KS and overlap matrices to the DIIS buffer
    1065        48832 :       CALL cp_cfm_to_cfm(kc, parameters)
    1066        48832 :       CALL cp_cfm_to_cfm(sc, smat)
    1067              : 
    1068              :       ! KC <- K*C; SC <- S*C. Retain the existing factor of two in SC.
    1069              :       ! new_errors <- KC*(SC)^H - SC*(KC)^H.
    1070        48832 :       IF (local) THEN
    1071       100734 :          CPASSERT(PRODUCT(parameters%matrix_struct%context%num_pe) == 1)
    1072       100734 :          CPASSERT(PRODUCT(cmos%matrix_struct%context%num_pe) == 1)
    1073       100734 :          CPASSERT(PRODUCT(kc%matrix_struct%context%num_pe) == 1)
    1074       100734 :          CPASSERT(PRODUCT(smat%matrix_struct%context%num_pe) == 1)
    1075       100734 :          CPASSERT(PRODUCT(sc%matrix_struct%context%num_pe) == 1)
    1076       100734 :          CPASSERT(PRODUCT(new_errors%matrix_struct%context%num_pe) == 1)
    1077              :          CALL gemm_ctx%gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), &
    1078              :                             parameters%local_data, SIZE(parameters%local_data, 1), &
    1079              :                             cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
    1080        33578 :                             kc%local_data, SIZE(kc%local_data, 1))
    1081              :          CALL gemm_ctx%gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), &
    1082              :                             smat%local_data, SIZE(smat%local_data, 1), &
    1083              :                             cmos%local_data, SIZE(cmos%local_data, 1), (0.0_dp, 0.0_dp), &
    1084        33578 :                             sc%local_data, SIZE(sc%local_data, 1))
    1085              :          CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
    1086              :                             sc%local_data, SIZE(sc%local_data, 1), &
    1087              :                             kc%local_data, SIZE(kc%local_data, 1), (0.0_dp, 0.0_dp), &
    1088        33578 :                             new_errors%local_data, SIZE(new_errors%local_data, 1))
    1089              :          CALL gemm_ctx%gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), &
    1090              :                             kc%local_data, SIZE(kc%local_data, 1), &
    1091              :                             sc%local_data, SIZE(sc%local_data, 1), (-1.0_dp, 0.0_dp), &
    1092        33578 :                             new_errors%local_data, SIZE(new_errors%local_data, 1))
    1093              :       ELSE
    1094        15254 :          CALL parallel_gemm("N", "N", nao, homo, nao, CMPLX(maxocc, KIND=dp), parameters, cmos, (0.0_dp, 0.0_dp), kc)
    1095        15254 :          CALL parallel_gemm("N", "N", nao, homo, nao, (2.0_dp, 0.0_dp), smat, cmos, (0.0_dp, 0.0_dp), sc)
    1096        15254 :          CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), sc, kc, (0.0_dp, 0.0_dp), new_errors)
    1097        15254 :          CALL parallel_gemm("N", "C", nao, nao, homo, (1.0_dp, 0.0_dp), kc, sc, (-1.0_dp, 0.0_dp), new_errors)
    1098              :       END IF
    1099              : 
    1100        48832 :       CALL timestop(handle)
    1101              : 
    1102              :    END SUBROUTINE qs_diis_b_calc_err_kp
    1103              : 
    1104              : ! **************************************************************************************************
    1105              : !> \brief Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points
    1106              : !> \param diis_buffer ...
    1107              : !> \param coeffs ...
    1108              : !> \param ib ...
    1109              : !> \param nb ...
    1110              : !> \param delta ...
    1111              : !> \param error_max ...
    1112              : !> \param diis_step ...
    1113              : !> \param eps_diis ...
    1114              : !> \param nspin ...
    1115              : !> \param nkp ...
    1116              : !> \param nkp_local ...
    1117              : !> \param nmixing ...
    1118              : !> \param scf_section ...
    1119              : !> \param para_env_inter_kp communicator connecting the k-point groups
    1120              : ! **************************************************************************************************
    1121        27498 :    SUBROUTINE qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, &
    1122              :                                 nspin, nkp, nkp_local, nmixing, scf_section, para_env_inter_kp)
    1123              : 
    1124              :       TYPE(qs_diis_buffer_type_kp), POINTER              :: diis_buffer
    1125              :       COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT)      :: coeffs
    1126              :       INTEGER, INTENT(IN)                                :: ib, nb
    1127              :       REAL(KIND=dp), INTENT(IN)                          :: delta
    1128              :       REAL(KIND=dp), INTENT(OUT)                         :: error_max
    1129              :       LOGICAL, INTENT(OUT)                               :: diis_step
    1130              :       REAL(KIND=dp), INTENT(IN)                          :: eps_diis
    1131              :       INTEGER, INTENT(IN)                                :: nspin, nkp, nkp_local
    1132              :       INTEGER, INTENT(IN), OPTIONAL                      :: nmixing
    1133              :       TYPE(section_vals_type), POINTER                   :: scf_section
    1134              :       TYPE(mp_para_env_type), POINTER                    :: para_env_inter_kp
    1135              : 
    1136              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_diis_b_step_kp'
    1137              :       REAL(KIND=dp), PARAMETER :: eigenvalue_threshold = 1.0E-12_dp
    1138              : 
    1139              :       CHARACTER(LEN=2*default_string_length)             :: message
    1140              :       COMPLEX(KIND=dp)                                   :: tmp
    1141        27498 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :)     :: a, b
    1142              :       INTEGER                                            :: handle, ikp, ispin, jb, my_nmixing, nb1, &
    1143              :                                                             output_unit
    1144              :       LOGICAL                                            :: eigenvectors_discarded
    1145        27498 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ev
    1146              :       TYPE(cp_cfm_type)                                  :: old_errors
    1147              :       TYPE(cp_cfm_type), POINTER                         :: new_errors
    1148              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
    1149              :       TYPE(cp_fm_type)                                   :: ierr, rerr
    1150              :       TYPE(cp_logger_type), POINTER                      :: logger
    1151              : 
    1152        27498 :       NULLIFY (matrix_struct, new_errors, logger)
    1153              : 
    1154        27498 :       CALL timeset(routineN, handle)
    1155              : 
    1156        27498 :       diis_step = .FALSE.
    1157              : 
    1158        27498 :       my_nmixing = 2
    1159        27498 :       IF (PRESENT(nmixing)) my_nmixing = nmixing
    1160              : 
    1161        27498 :       logger => cp_get_default_logger()
    1162              : 
    1163              :       ! Quick return, if no DIIS is requested
    1164        27498 :       IF (diis_buffer%nbuffer < 1) THEN
    1165            0 :          CALL timestop(handle)
    1166            0 :          RETURN
    1167              :       END IF
    1168              : 
    1169              :       ! Check, if a DIIS step is appropriate
    1170        27498 :       diis_step = ((diis_buffer%ncall >= my_nmixing) .AND. (delta < eps_diis))
    1171              : 
    1172              :       ! Calculate the DIIS buffer, and update it if max_error < eps_diis
    1173        27498 :       CALL cp_cfm_get_info(diis_buffer%error(ib, 1, 1), matrix_struct=matrix_struct)
    1174        27498 :       CALL cp_fm_create(ierr, matrix_struct)
    1175        27498 :       CALL cp_fm_create(rerr, matrix_struct)
    1176        27498 :       CALL cp_cfm_create(old_errors, matrix_struct)
    1177       109992 :       ALLOCATE (b(nb, nb))
    1178        27498 :       b = 0.0_dp
    1179       122794 :       DO jb = 1, nb
    1180       243442 :          DO ikp = 1, nkp_local
    1181       407349 :             DO ispin = 1, nspin
    1182       163907 :                new_errors => diis_buffer%error(ib, ispin, ikp)
    1183       163907 :                CALL cp_cfm_to_fm(diis_buffer%error(jb, ispin, ikp), rerr, ierr)
    1184       163907 :                CALL cp_fm_scale(-1.0_dp, ierr)
    1185       163907 :                CALL cp_fm_to_cfm(rerr, ierr, old_errors)
    1186       163907 :                CALL cp_cfm_trace(old_errors, new_errors, tmp)
    1187       312053 :                b(jb, ib) = b(jb, ib) + 1.0_dp/REAL(nkp, dp)*tmp
    1188              :             END DO
    1189              :          END DO
    1190       122794 :          b(ib, jb) = CONJG(b(jb, ib))
    1191              :       END DO
    1192        27498 :       CALL cp_fm_release(ierr)
    1193        27498 :       CALL cp_fm_release(rerr)
    1194        27498 :       CALL cp_cfm_release(old_errors)
    1195              :       ! cp_cfm_trace has already reduced each trace over the ranks of its
    1196              :       ! k-point group. Reduce the group contributions only once, over the
    1197              :       ! communicator that connects the k-point groups.
    1198        27498 :       CALL para_env_inter_kp%sum(b)
    1199              : 
    1200        27498 :       error_max = SQRT(REAL(b(ib, ib))**2 + AIMAG(b(ib, ib))**2)
    1201              : 
    1202              :       output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIIS_INFO", &
    1203        27498 :                                          extension=".scfLog")
    1204        27498 :       IF (output_unit > 0) THEN
    1205              :          WRITE (UNIT=output_unit, FMT="(/,T9,A,I4,/,(T9,A,ES12.3))") &
    1206            0 :             "DIIS | Current SCF DIIS buffer size:         ", nb, &
    1207            0 :             "DIIS | Maximum SCF DIIS error at last step:  ", error_max, &
    1208            0 :             "DIIS | Current SCF convergence:              ", delta, &
    1209            0 :             "DIIS | Threshold value for a DIIS step:      ", eps_diis
    1210            0 :          IF (error_max < eps_diis) THEN
    1211              :             WRITE (UNIT=output_unit, FMT="(T9,A)") &
    1212            0 :                "DIIS | => The SCF DIIS buffer will be updated"
    1213              :          ELSE
    1214              :             WRITE (UNIT=output_unit, FMT="(T9,A)") &
    1215            0 :                "DIIS | => No update of the SCF DIIS buffer"
    1216              :          END IF
    1217            0 :          IF (diis_step .AND. (error_max < eps_diis)) THEN
    1218              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
    1219            0 :                "DIIS | => A SCF DIIS step will be performed"
    1220              :          ELSE
    1221              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
    1222            0 :                "DIIS | => No SCF DIIS step will be performed"
    1223              :          END IF
    1224              :       END IF
    1225              : 
    1226              :       ! Update the SCF DIIS buffer
    1227        27498 :       IF (error_max < eps_diis) THEN
    1228        82412 :          DO jb = 1, nb
    1229        64354 :             diis_buffer%b_matrix(ib, jb) = b(ib, jb)
    1230        82412 :             diis_buffer%b_matrix(jb, ib) = b(jb, ib)
    1231              :          END DO
    1232              :       ELSE
    1233              : 
    1234         9440 :          diis_step = .FALSE.
    1235              :       END IF
    1236        27498 :       DEALLOCATE (b)
    1237              : 
    1238              :       ! Perform DIIS step
    1239        27498 :       IF (diis_step) THEN
    1240              : 
    1241        10694 :          nb1 = nb + 1
    1242              : 
    1243        42776 :          ALLOCATE (a(nb1, nb1))
    1244        32082 :          ALLOCATE (b(nb1, nb1))
    1245        32082 :          ALLOCATE (ev(nb1))
    1246              : 
    1247              :          ! Set up the linear DIIS equation system
    1248       213814 :          b(1:nb, 1:nb) = diis_buffer%b_matrix(1:nb, 1:nb)
    1249              : 
    1250        52004 :          b(1:nb, nb1) = -1.0_dp
    1251        52004 :          b(nb1, 1:nb) = -1.0_dp
    1252        10694 :          b(nb1, nb1) = 0.0_dp
    1253              : 
    1254              :          ! Solve the linear DIIS equation system
    1255        62698 :          ev(1:nb1) = 0.0_dp !eigenvalues
    1256       317822 :          a(1:nb1, 1:nb1) = 0.0_dp !eigenvectors
    1257        10694 :          CALL diag_complex(b(1:nb1, 1:nb1), a(1:nb1, 1:nb1), ev(1:nb1))
    1258       317822 :          b(1:nb1, 1:nb1) = a(1:nb1, 1:nb1)
    1259              : 
    1260        10694 :          eigenvectors_discarded = .FALSE.
    1261              : 
    1262        62698 :          DO jb = 1, nb1
    1263        62698 :             IF (ABS(ev(jb)) < eigenvalue_threshold) THEN
    1264        18164 :                IF (output_unit > 0) THEN
    1265            0 :                   IF (.NOT. eigenvectors_discarded) THEN
    1266              :                      WRITE (UNIT=output_unit, FMT="(T9,A)") &
    1267            0 :                         "DIIS | Checking eigenvalues of the DIIS error matrix"
    1268              :                   END IF
    1269              :                   WRITE (UNIT=message, FMT="(T9,A,I6,A,ES10.1,A,ES10.1)") &
    1270            0 :                      "DIIS | Eigenvalue ", jb, " = ", ev(jb), " is smaller than "// &
    1271            0 :                      "threshold ", eigenvalue_threshold
    1272            0 :                   CALL compress(message)
    1273            0 :                   WRITE (UNIT=output_unit, FMT="(T9,A)") TRIM(message)
    1274            0 :                   eigenvectors_discarded = .TRUE.
    1275              :                END IF
    1276       108700 :                a(1:nb1, jb) = 0.0_dp
    1277              :             ELSE
    1278       198428 :                a(1:nb1, jb) = a(1:nb1, jb)/ev(jb)
    1279              :             END IF
    1280              :          END DO
    1281              : 
    1282        10694 :          IF ((output_unit > 0) .AND. eigenvectors_discarded) THEN
    1283              :             WRITE (UNIT=output_unit, FMT="(T9,A,/)") &
    1284            0 :                "DIIS | The corresponding eigenvectors were discarded"
    1285              :          END IF
    1286              : 
    1287       328516 :          coeffs(1:nb) = -MATMUL(a(1:nb, 1:nb1), CONJG(b(nb1, 1:nb1)))
    1288              :       ELSE
    1289              : 
    1290        70790 :          coeffs(:) = 0.0_dp
    1291        16804 :          coeffs(ib) = 1.0_dp
    1292              :       END IF
    1293              : 
    1294              :       CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
    1295        27498 :                                         "PRINT%DIIS_INFO")
    1296              : 
    1297        27498 :       CALL timestop(handle)
    1298              : 
    1299        82494 :    END SUBROUTINE qs_diis_b_step_kp
    1300        10694 : END MODULE qs_diis
        

Generated by: LCOV version 2.0-1